Memahami pembuatan variabel dummy (manual atau otomatis) di GLM

Jika variabel faktor (misalnya jenis kelamin dengan level M dan F) digunakan dalam rumus glm, variabel dummy dibuat, dan dapat ditemukan dalam ringkasan model glm bersama dengan koefisien yang terkait (misalnya genderM)

Jika, alih-alih mengandalkan R untuk membagi faktor dengan cara ini, faktor tersebut dikodekan dalam serangkaian variabel 0/1 numerik (misalnya genderM (1 untuk M, 0 untuk F), genderF (1 untuk F, 0 untuk M) dan variabel-variabel ini kemudian digunakan sebagai variabel numerik dalam rumus GLM, apakah hasil koefisien akan berbeda?

Pada dasarnya pertanyaannya adalah: apakah R menggunakan perhitungan koefisien yang berbeda ketika bekerja dengan variabel faktor versus variabel numerik?

Pertanyaan tindak lanjut (mungkin dijawab oleh yang di atas): selain hanya efisiensi membiarkan R membuat variabel dummy, apakah ada masalah dengan pengkodean ulang faktor sebagai serangkaian variabel 0,1 numerik dan menggunakan yang ada dalam model sebagai gantinya?

Variabel kategorikal (disebut " faktor " dalam R) perlu diwakili oleh kode numerik dalam model regresi berganda. Ada banyak cara yang memungkinkan untuk membuat kode numerik dengan tepat (lihat daftar hebat ini di situs bantuan statistik UCLA). Secara default, R menggunakan pengkodean level referensi (yang disebut R sebagai "contr.treatment"), dan yang merupakan statistik default. Ini dapat diubah untuk semua kontras untuk seluruh sesi R Anda menggunakan ? Opsi , atau untuk analisis / variabel tertentu menggunakan ? Kontras atau ? C (perhatikan modal). Jika Anda memerlukan informasi lebih lanjut tentang pengkodean tingkat referensi, saya jelaskan di sini: Regresi berdasarkan misalnya pada hari dalam seminggu.

Beberapa orang menganggap pengkodean tingkat referensi membingungkan, dan Anda tidak harus menggunakannya. Jika Anda mau, Anda bisa memiliki dua variabel untuk pria dan wanita; ini disebut level berarti pengkodean. Namun, jika Anda melakukan itu, Anda harus menekan intersep atau matriks model akan tunggal dan regresi tidak dapat sesuai dengan catatan @Affine di atas dan seperti yang saya jelaskan di sini: Pengodean variabel kualitatif mengarah ke singularitas . Untuk menekan intersep, Anda memodifikasi rumus Anda dengan menambahkan -1atau +0seperti itu: y~... -1atau y~... +0.

Menggunakan pengkodean level berarti alih-alih pengkodean level referensi akan mengubah estimasi koefisien dan arti dari tes hipotesis yang dicetak dengan output Anda. Ketika Anda memiliki faktor dua level (misalnya, pria vs wanita) dan Anda menggunakan pengkodean level referensi, Anda akan melihat intersep dipanggil (constant)dan hanya satu variabel yang tercantum dalam output (mungkin sexM). Intersep adalah rata-rata dari kelompok referensi (mungkin perempuan) dan sexMmerupakan perbedaan antara rata-rata laki-laki dan rata-rata perempuan. Nilai p yang terkait dengan intersep adalah uji satu sampel apakah tingkat referensi memiliki rata-rata dan nilai p yang terkait dengan$t$$0$sexMmemberi tahu Anda jika jenis kelamin berbeda pada respons Anda. Tetapi jika Anda menggunakan pengkodean level means, Anda akan memiliki dua variabel terdaftar dan masing-masing nilai p akan sesuai dengan -sampel satu-sampel apakah rata-rata level tersebut adalah . Artinya, tidak ada nilai-p yang akan menjadi ujian apakah jenis kelaminnya berbeda. $t$$0$

set.seed(1)
y    = c(    rnorm(30), rnorm(30, mean=1)         )
sex  = rep(c("Female",  "Male"          ), each=30)
fem  = ifelse(sex=="Female", 1, 0)
male = ifelse(sex=="Male", 1, 0)

ref.level.coding.model   = lm(y~sex)
level.means.coding.model = lm(y~fem+male+0)

summary(ref.level.coding.model)
# ...
# Coefficients:
#             Estimate Std. Error t value Pr(>|t|)    
# (Intercept)  0.08246    0.15740   0.524    0.602    
# sexMale      1.05032    0.22260   4.718 1.54e-05 ***
#   ---
#   Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
# ...
summary(level.means.coding.model)
# ...
# Coefficients:
#      Estimate Std. Error t value Pr(>|t|)    
# fem   0.08246    0.15740   0.524    0.602    
# male  1.13277    0.15740   7.197 1.37e-09 ***
#   ---
#   Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
# ...

Koefisien yang diperkirakan akan menjadi subjek yang sama dengan kondisi yang Anda buat variabel dummy Anda (yaitu yang numerik) konsisten dengan R. Misalnya: mari 'buat data palsu dan muat menggunakan Poisson glm menggunakan faktor. Perhatikan bahwa glfungsi menciptakan variabel faktor.

> counts <- c(18,17,15,20,10,20,25,13,12)
> outcome <- gl(3,1,9)
> outcome
[1] 1 2 3 1 2 3 1 2 3
Levels: 1 2 3
> class(outcome)
[1] "factor"
> glm.1<- glm(counts ~ outcome, family = poisson())
> summary(glm.1)

Call:
glm(formula = counts ~ outcome, family = poisson())

Deviance Residuals: 
    Min       1Q   Median       3Q      Max  
-0.9666  -0.6713  -0.1696   0.8471   1.0494  

Coefficients:
            Estimate Std. Error z value Pr(>|z|)    
(Intercept)   3.0445     0.1260  24.165   <2e-16 ***
outcome2     -0.4543     0.2022  -2.247   0.0246 *  
outcome3     -0.2930     0.1927  -1.520   0.1285    
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

(Dispersion parameter for poisson family taken to be 1)

    Null deviance: 10.5814  on 8  degrees of freedom
Residual deviance:  5.1291  on 6  degrees of freedom
AIC: 52.761

Number of Fisher Scoring iterations: 4

Karena hasil memiliki tiga level, saya membuat dua variabel dummy (dummy.1 = 0 jika hasil = 2 dan dummy.2 = 1 jika hasil = 3) dan mereparasi dengan menggunakan nilai-nilai numerik ini:

> dummy.1=rep(0,9)
> dummy.2=rep(0,9)
> dummy.1[outcome==2]=1
> dummy.2[outcome==3]=1
> glm.2<- glm(counts ~ dummy.1+dummy.2, family = poisson())
> summary(glm.2)

Call:
glm(formula = counts ~ dummy.1 + dummy.2, family = poisson())

Deviance Residuals: 
    Min       1Q   Median       3Q      Max  
-0.9666  -0.6713  -0.1696   0.8471   1.0494  

Coefficients:
            Estimate Std. Error z value Pr(>|z|)    
(Intercept)   3.0445     0.1260  24.165   <2e-16 ***
dummy.1      -0.4543     0.2022  -2.247   0.0246 *  
dummy.2      -0.2930     0.1927  -1.520   0.1285    
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

(Dispersion parameter for poisson family taken to be 1)

    Null deviance: 10.5814  on 8  degrees of freedom
Residual deviance:  5.1291  on 6  degrees of freedom
AIC: 52.761

Number of Fisher Scoring iterations: 4

Seperti yang Anda lihat, koefisien yang diperkirakan sama. Tetapi Anda harus berhati-hati saat membuat variabel dummy Anda jika Anda ingin mendapatkan hasil yang sama. Sebagai contoh jika saya membuat dua variabel dummy sebagai (dummy.1 = 0 jika hasil = 1 dan dummy.2 = 1 jika hasil = 2) maka hasil estimasi berbeda seperti berikut:

> dummy.1=rep(0,9)
> dummy.2=rep(0,9)
> dummy.1[outcome==1]=1
> dummy.2[outcome==2]=1
> glm.3<- glm(counts ~ dummy.1+dummy.2, family = poisson())
> summary(glm.3)

Call:
glm(formula = counts ~ dummy.1 + dummy.2, family = poisson())

Deviance Residuals: 
    Min       1Q   Median       3Q      Max  
-0.9666  -0.6713  -0.1696   0.8471   1.0494  

Coefficients:
            Estimate Std. Error z value Pr(>|z|)    
(Intercept)   2.7515     0.1459   18.86   <2e-16 ***
dummy.1       0.2930     0.1927    1.52    0.128    
dummy.2      -0.1613     0.2151   -0.75    0.453    
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

(Dispersion parameter for poisson family taken to be 1)

    Null deviance: 10.5814  on 8  degrees of freedom
Residual deviance:  5.1291  on 6  degrees of freedom
AIC: 52.761

Number of Fisher Scoring iterations: 4

Ini karena ketika Anda menambahkan outcomevariabel di glm.1, R secara default membuat dua variabel dummy yaitu outcome2dan outcome3dan mendefinisikannya sama dengan dummy.1dan dummy.2di glm.2 yaitu tingkat hasil pertama adalah ketika semua variabel dummy lainnya ( outcome2dan outcome3) ditetapkan menjadi nol.