Coding interaksi antara nominal dan prediktor berkelanjutan untuk regresi logistik di MATLAB

8

Jadi data kami disusun sebagai berikut:

Kita punya M.peserta, masing-masing peserta dapat dikategorikan menjadi 3 kelompok ( G SEBUAH,B,C), dan untuk setiap peserta yang kami miliki Nsampel variabel kontinu. Dan kami mencoba memprediksi nilai yang 0 atau 1.

Bagaimana kita menggunakan matlab untuk menguji interaksi antara variabel kontinu dan variabel kategori dalam memprediksi nilai-nilai ini?

mpacer
sumber
@Thislstheld Kode apa yang Anda gunakan untuk menyesuaikan regresi logistik Anda?
chl
@ chl - Entah glmfit dengan 'tautan', 'logit' atau mnrfit akan berfungsi, tidak ada preferensi khusus untuk keduanya.
mpacer
Saya tidak berpikir ini adalah pertanyaan statistik, tetapi pertanyaan pemrograman, yang ditempatkan lebih tepat di StackOverflow ...
Manoel Galdino
@Manoel sebagai balasan oleh @chl menunjukkan, pertanyaan ini pada dasarnya bersifat statistik. Jika pertanyaan seperti ini di luar topik, maka akan menjadi sekitar setengah dari pertanyaan di situs ini!
whuber

Jawaban:

9

Cara termudah, IMO, adalah membangun sendiri matriks desain, seperti glmfitmenerima matriks nilai mentah (yang diamati) atau matriks desain. Pengkodean istilah interaksi tidak terlalu sulit setelah Anda menulis model lengkap. Katakanlah kita memiliki dua prediktor,x (kontinu) dan g (kategorikal, dengan tiga tingkat tidak teratur, katakanlah g=1,2,3). Menggunakan notasi Wilkinson, kita akan menulis model ini sebagai y ~ x + g + x:g, mengabaikan sisi kiri (untuk hasil binomial, kita akan menggunakan fungsi tautan logit). Kami hanya membutuhkan dua vektor dummy untuk mengkodekan glevel (sebagai ada / tidak ada untuk pengamatan tertentu), jadi kami akan memiliki 5 koefisien regresi, ditambah istilah intersepsi. Ini dapat diringkas sebagai

β0+β1x+β2sayag=2+β3sayag=3+β4x×sayag=2+β5x×sayag=3,

dimana saya singkatan dari matriks indikator yang mengkode tingkat g.

Di Matlab, menggunakan contoh online, saya akan lakukan sebagai berikut:

x = [2100 2300 2500 2700 2900 3100 3300 3500 3700 3900 4100 4300]';
g = [1 1 1 1 2 2 2 2 3 3 3 3]';
gcat = dummyvar(g);
gcat = gcat(:,2:3); % remove the first column
X = [x gcat x.*gcat(:,1) x.*gcat(:,2)]; 
n = [48 42 31 34 31 21 23 23 21 16 17 21]';
y = [1 2 0 3 8 8 14 17 19 15 17 21]';
[b, dev, stats] = glmfit(X, [y n], 'binomial', 'link', 'probit');

Saya tidak menyertakan kolom yang untuk intersep karena disertakan secara default. Matriks desain terlihat seperti

    2100           0           0           0           0
    2300           0           0           0           0
    2500           0           0           0           0
    2700           0           0           0           0
    2900           1           0        2900           0
    3100           1           0        3100           0
    3300           1           0        3300           0
    3500           1           0        3500           0
    3700           0           1           0        3700
    3900           0           1           0        3900
    4100           0           1           0        4100
    4300           0           1           0        4300

dan Anda dapat melihat bahwa istilah interaksi hanya dikodekan sebagai produk dari xdengan kolom yang sesuai dari g(g = 2 dan g = 3, karena kita tidak memerlukan level pertama).

Hasilnya diberikan di bawah ini, sebagai koefisien, kesalahan standar, statistik dan nilai-p (dari statsstruktur):

   int.   -3.8929    2.0251   -1.9223    0.0546
   x       0.0009    0.0008    1.0663    0.2863
   g2     -3.2125    2.7622   -1.1630    0.2448
   g3     -5.7745    7.5542   -0.7644    0.4446
   x:g2    0.0013    0.0010    1.3122    0.1894
   x:g3    0.0021    0.0021    0.9882    0.3230

Sekarang, menguji interaksi dapat dilakukan dengan menghitung perbedaan penyimpangan dari model penuh di atas dan model yang dikurangi (menghilangkan istilah interaksi, yaitu dua kolom terakhir dari matriks desain). Ini dapat dilakukan secara manual, atau menggunakan lratiotestfungsi yang menyediakan uji hipotesis Likelihood ratio. Penyimpangan untuk model lengkap adalah 4.3122 ( dev), sedangkan untuk model tanpa interaksi adalah 6.4200 (saya menggunakan glmfit(X(:,1:3), [y n], 'binomial', 'link', 'probit');), dan uji LR terkait memiliki dua derajat kebebasan (perbedaan dalam jumlah parameter antara kedua model). Karena penyimpangan skala hanya dua kali log-kemungkinan untuk GLM, kita bisa menggunakan

[H, pValue, Ratio, CriticalValue] = lratiotest(4.3122/2, 6.4200/2, 2)

di mana statistik didistribusikan sebagai a χ2dengan 2 df (nilai kritisnya adalah 5.9915, lihat chi2inv(0.95, 2)). Keluaran menunjukkan hasil yang tidak signifikan: Kami tidak dapat menyimpulkan adanya interaksi antara xdan gdalam sampel yang diamati.

Saya kira Anda dapat menyelesaikan langkah-langkah di atas dalam fungsi pilihan Anda yang nyaman. (Perhatikan bahwa tes LR dapat dilakukan dengan tangan dalam beberapa perintah!)


Saya memeriksa hasil tersebut terhadap keluaran R, yang diberikan selanjutnya.

Berikut adalah kode R:

x <- c(2100,2300,2500,2700,2900,3100,3300,3500,3700,3900,4100,4300)
g <- gl(3, 4)
n <- c(48,42,31,34,31,21,23,23,21,16,17,21)
y <- c(1,2,0,3,8,8,14,17,19,15,17,21)
f <- cbind(y, n-y) ~ x*g
model.matrix(f)  # will be model.frame() for glm()
m1 <- glm(f, family=binomial("probit"))
summary(m1)

Inilah hasilnya, untuk koefisien dalam model lengkap,

Call:
glm(formula = f, family = binomial("probit"))

Deviance Residuals: 
    Min       1Q   Median       3Q      Max  
-1.7124  -0.1192   0.1494   0.3036   0.5585  

Coefficients:
             Estimate Std. Error z value Pr(>|z|)  
(Intercept) -3.892859   2.025096  -1.922   0.0546 .
x            0.000884   0.000829   1.066   0.2863  
g2          -3.212494   2.762155  -1.163   0.2448  
g3          -5.774400   7.553615  -0.764   0.4446  
x:g2         0.001335   0.001017   1.312   0.1894  
x:g3         0.002061   0.002086   0.988   0.3230  

Untuk perbandingan dua model bersarang, saya menggunakan perintah berikut:

m0 <- update(m1, . ~ . -x:g)
anova(m1,m0)

yang menghasilkan "tabel penyimpangan" berikut:

Analysis of Deviance Table

Model 1: cbind(y, n - y) ~ x + g
Model 2: cbind(y, n - y) ~ x * g
  Resid. Df Resid. Dev Df Deviance
1         8     6.4200            
2         6     4.3122  2   2.1078
chl
sumber
Btw, untuk orang lain yang ingin menggunakan ini, untuk referensi Anda, Anda perlu membagi dengan -2 bukan 2, jika tidak maka akan menghasilkan kesalahan. Lihat en.wikipedia.org/wiki/Deviance_%28statistics%29
mpacer
Juga, terima kasih banyak - ini sangat membantu baik dalam pemahaman saya tentang masalah dan masalah umum yang saya hadapi :).
mpacer