Kesulitan menguji linearitas dalam regresi

Dalam Pemodelan Statistik: Dua Budaya Leo Breiman menulis

Praktik yang diterapkan saat ini adalah untuk memeriksa model data yang sesuai dengan menggunakan uji goodness-of-fit dan analisis residual. Pada satu titik, beberapa tahun yang lalu, saya membuat masalah regresi disimulasikan dalam tujuh dimensi dengan jumlah nonlinier yang terkontrol. Tes standar good-of-fit tidak menolak linearitas sampai nonlinier ekstrem.

Breiman tidak memberikan detail simulasi itu. Dia merujuk makalah yang menurutnya memberikan pembenaran teoretis untuk pengamatannya, tetapi makalah itu tidak diterbitkan.

Adakah yang melihat hasil simulasi yang dipublikasikan atau makalah teoretis untuk mendukung klaim Brieman?

Saya membuat simulasi yang akan menjawab deskripsi Breiman dan hanya menemukan yang jelas: hasilnya tergantung pada konteks dan pada apa yang dimaksud dengan "ekstrim."

Banyak hal yang bisa dikatakan, tetapi izinkan saya membatasinya hanya pada satu contoh yang dilakukan dengan menggunakan Rkode yang mudah dimodifikasi untuk digunakan pembaca yang berminat dalam penyelidikan mereka sendiri. Kode ini dimulai dengan menyiapkan matriks desain yang terdiri dari nilai-nilai independen yang terdistribusi secara merata yang kira-kira ortogonal (sehingga kita tidak masuk ke masalah multikolinieritas). Ini menghitung interaksi kuadratik tunggal (yaitu, nonlinier) antara dua variabel pertama: ini hanya salah satu dari banyak jenis "nonlinier" yang dapat dipelajari, tetapi setidaknya itu adalah yang umum, yang dipahami dengan baik. Kemudian ia menstandarkan semuanya sehingga koefisien akan sebanding:

set.seed(41)
p <- 7                                            # Dimensions
n <- 2^p                                          # Observations
x <- as.matrix(do.call(expand.grid, lapply(as.list(1:p), function(i) c(-1,1))))
x <- x + runif(n*p, min=-1, max=1)
x <- cbind(x, x.12 = x[,1]*x[,2])                 # The nonlinear part
x <- apply(x, 2, function(y) (y - mean(y))/sd(y)) # Standardization

Untuk model OLS dasar (tanpa nonlinier) kita harus menentukan beberapa koefisien dan standar deviasi dari kesalahan residual. Berikut adalah sekumpulan koefisien satuan dan SD yang sebanding:

beta <- rep(c(1,-1), p)[1:p]
sd <- 1

$1/4$$1$$-1$

gamma = 1/4          # The standardized interaction term
df <- data.frame(x)
df$y <- x %*% c(beta, gamma) + rnorm(n, sd=sd)
summary(df)
cor(df)*100
plot(df, lower.panel=function(x,y) lines(lowess(y~x)), 
     upper.panel=function(x,y) points(x,y, pch=".", cex=4))
summary(lm(df$y ~ x))

Daripada mengarungi semua output di sini, mari kita lihat data ini menggunakan output dari plotperintah:

Jejak lowess pada segitiga bawah pada dasarnya tidak menunjukkan hubungan linear antara interaksi ( x.12) dan variabel dependen ( y) dan hubungan linear sederhana antara variabel lain dan y. Hasil OLS mengkonfirmasi bahwa; interaksi ini hampir tidak signifikan:

Coefficients:
            Estimate Std. Error t value Pr(>|t|)    
(Intercept)   0.0263     0.0828    0.32    0.751    
xVar1         0.9947     0.0833   11.94   <2e-16 ***
xVar2        -0.8713     0.0842  -10.35   <2e-16 ***
xVar3         1.0709     0.0836   12.81   <2e-16 ***
xVar4        -1.0007     0.0840  -11.92   <2e-16 ***
xVar5         1.0233     0.0836   12.24   <2e-16 ***
xVar6        -0.9514     0.0835  -11.40   <2e-16 ***
xVar7         1.0482     0.0835   12.56   <2e-16 ***
xx.12         0.1902     0.0836    2.27    0.025 *  

Saya akan mengambil nilai p dari istilah interaksi sebagai tes nonlinier: ketika nilai p ini cukup rendah (Anda dapat memilih seberapa rendah), kami akan mendeteksi nonlinieritas.

(Ada kehalusan di sini tentang apa yang sebenarnya kita cari. Dalam praktiknya kita mungkin perlu memeriksa semua 7 * 6/2 = 21 kemungkinan interaksi kuadratik tersebut, serta mungkin 7 istilah kuadratik yang lain, daripada berfokus pada satu istilah tunggal seperti yang dilakukan di sini. Kami ingin membuat koreksi untuk 28 tes yang saling terkait ini. Saya tidak secara eksplisit membuat koreksi ini di sini, karena sebaliknya saya menampilkan distribusi simulasi nilai-p. Anda dapat membaca tingkat deteksi langsung dari yang histogram pada akhir berdasarkan Anda ambang batas signifikansi.)

Tetapi jangan hanya melakukan analisis ini sekali saja; mari kita lakukan berkali-kali, menghasilkan nilai-nilai baru ydi setiap iterasi sesuai dengan model yang sama dan matriks desain yang sama. Untuk mencapai ini, kami menggunakan fungsi untuk melakukan satu iterasi dan mengembalikan nilai p dari istilah interaksi:

test <- function(gamma, sd=1) {
  y <- x %*% c(beta, gamma) + rnorm(n, sd=sd)
  fit <- summary(lm(y ~ x))
  m <- coef(fit)
  n <- dim(m)[1]
  m[n, 4]
}

Saya memilih untuk menyajikan hasil simulasi sebagai histogram nilai-p, memvariasikan koefisien standar gammadari istilah interaksi. Pertama, histogram:

h <- function(g, n.trials=1000) {
  hist(replicate(n.trials, test(g, sd)), xlim=c(0,1), 
       main=toString(g), xlab="x1:x2 p-value")
}
par(mfrow=c(2,2)) # Draw a 2 by 2 panel of results

Sekarang untuk melakukan pekerjaan. Dibutuhkan beberapa detik untuk 1000 percobaan per simulasi (dan empat simulasi independen, dimulai dengan nilai yang diberikan dari istilah interaksi dan secara bertahap membagi dua setiap kali):

temp <- sapply(2^(-3:0) * gamma, h)

Hasil:

xsdbeta$1/4$$1/8$$1/16$gamma$1/2$

$1/32$$1/4$xsdbetasd

Singkatnya, simulasi seperti ini dapat membuktikan apa pun yang Anda suka jika Anda hanya mengaturnya dan menafsirkannya dengan cara yang benar. Itu menunjukkan ahli statistik individu harus melakukan eksplorasi mereka sendiri, sesuai dengan masalah khusus yang mereka hadapi, untuk sampai pada pemahaman pribadi dan mendalam tentang kemampuan dan kelemahan prosedur yang mereka gunakan.

Tidak yakin itu memberikan jawaban akhir untuk pertanyaan itu, tetapi saya akan memberikan ini . Terutama poin 2. Lihat juga diskusi dalam lampiran A2 dari makalah ini .