Apakah ada perbedaan antara lm dan glm untuk keluarga gaussian glm?

45

Secara khusus, saya ingin tahu apakah ada perbedaan antara lm(y ~ x1 + x2)dan glm(y ~ x1 + x2, family=gaussian). Saya pikir kasus glm ini sama dengan lm. Apakah aku salah?

pengguna3682457
sumber
10
Iya dan tidak. Sebagai model statistik, tidak. Sebagai objek yang pas di R, ya; objek kembali berbeda, algoritma berbeda digunakan.
Reinstate Monica - G. Simpson
3
Sepertinya saya ada pertanyaan statistik di sini, serta satu R coding.
Silverfish

Jawaban:

48

Sementara untuk bentuk spesifik dari model yang disebutkan dalam isi pertanyaan (yaitu lm(y ~ x1 + x2)vs glm(y ~ x1 + x2, family=gaussian)), regresi dan GLM adalah model yang sama, pertanyaan judul menanyakan sesuatu yang sedikit lebih umum:

Apakah ada perbedaan antara lm dan glm untuk keluarga gaussian glm?

Yang jawabannya adalah "Ya!".

Alasan mereka bisa berbeda adalah karena Anda juga dapat menentukan fungsi tautan di GLM. Ini memungkinkan Anda untuk mencocokkan bentuk-bentuk tertentu dari hubungan nonlinear antara y (atau lebih tepatnya bersyarat bersyarat) dan variabel- x ; sementara Anda dapat melakukan ini nlsjuga, tidak perlu untuk memulai nilai, terkadang konvergensi lebih baik (juga sintaksnya sedikit lebih mudah).

Bandingkan, misalnya, model-model ini (Anda memiliki R jadi saya menganggap Anda dapat menjalankan ini sendiri):

x1=c(56.1, 26.8, 23.9, 46.8, 34.8, 42.1, 22.9, 55.5, 56.1, 46.9, 26.7, 33.9, 
37.0, 57.6, 27.2, 25.7, 37.0, 44.4, 44.7, 67.2, 48.7, 20.4, 45.2, 22.4, 23.2, 
39.9, 51.3, 24.1, 56.3, 58.9, 62.2, 37.7, 36.0, 63.9, 62.5, 44.1, 46.9, 45.4, 
23.7, 36.5, 56.1, 69.6, 40.3, 26.2, 67.1, 33.8, 29.9, 25.7, 40.0, 27.5)

x2=c(12.29, 11.42, 13.59, 8.64, 12.77, 9.9, 13.2, 7.34, 10.67, 18.8, 9.84, 16.72, 
10.32, 13.67, 7.65, 9.44, 14.52, 8.24, 14.14, 17.2, 16.21, 6.01, 14.23, 15.63, 
10.83, 13.39, 10.5, 10.01, 13.56, 11.26, 4.8, 9.59, 11.87, 11, 12.02, 10.9, 9.5, 
10.63, 19.03, 16.71, 15.11, 7.22, 12.6, 15.35, 8.77, 9.81, 9.49, 15.82, 10.94, 6.53)

y = c(1.54, 0.81, 1.39, 1.09, 1.3, 1.16, 0.95, 1.29, 1.35, 1.86, 1.1, 0.96,
1.03, 1.8, 0.7, 0.88, 1.24, 0.94, 1.41, 2.13, 1.63, 0.78, 1.55, 1.5, 0.96, 
1.21, 1.4, 0.66, 1.55, 1.37, 1.19, 0.88, 0.97, 1.56, 1.51, 1.09, 1.23, 1.2, 
1.62, 1.52, 1.64, 1.77, 0.97, 1.12, 1.48, 0.83, 1.06, 1.1, 1.21, 0.75)

lm(y ~ x1 + x2)
glm(y ~ x1 + x2, family=gaussian) 
glm(y ~ x1 + x2, family=gaussian(link="log")) 
nls(y ~ exp(b0+b1*x1+b2*x2), start=list(b0=-1,b1=0.01,b2=0.1))

Perhatikan bahwa pasangan pertama adalah model yang sama ( ), dan pasangan kedua adalah model yang sama ( dan pasangannya pada dasarnya sama di setiap pasangan.yiN(β0+β1x1i+β2x2i,σ2)yiN(exp(β0+β1x1i+β2x2i),σ2)

Jadi - dalam kaitannya dengan pertanyaan judul - Anda dapat memasukkan variasi model Gaussian yang jauh lebih luas dengan GLM daripada dengan regresi.

Glen_b
sumber
4
+1. Salah satu sisi komputasi dari hal-hal yang saya juga akan berpikir bahwa algoritma GLM akan menggunakan beberapa varian IRWLS (dalam banyak kasus) sementara LM akan menyampaikan pada beberapa varian solusi bentuk-tertutup.
usεr11852 mengatakan Reinstate Monic
@ usεr11852 - Saya akan berpikir itu EM, tetapi mereka mungkin hal yang sama dalam hal ini.
EngrStudent
1
Itu tidak menanggapi melihat "outlier" (kecuali melalui kemungkinan seperti yang dijelaskan di atas); reweighting adalah karena efek dari fungsi varians dan pergeseran dalam pendekatan linier lokal.
Glen_b
1
@ChrisChiasson: +1 ke komentar Glen_b. Seperti yang disebutkan ini tidak terkait dengan kekokohan algoritma dengan adanya outlier. Anda mungkin ingin menjelajahi berbagai keluarga (mis. Distribusi sesuai , atau kehilangan Huber; periksa lebih banyak tentang itu) - maaf baru online setelah beberapa hari libur ..tMASS::rlm
usεr11852 mengatakan Reinstate Monic
1
Anda dapat mencapai jenis ketahanan yang saya pikir Anda inginkan dalam beberapa cara. Namun, dengan glms dan model tipe regresi, Anda harus berhati-hati bukan hanya pada outlier di arah y tetapi juga outlier yang berpengaruh , yang dapat membuat diri mereka tidak terlihat tidak pada tempatnya ..
Glen_b
14

Jawaban singkatnya, mereka persis sama:

# Simulate data:
set.seed(42)
n <- 1000

x1 <- rnorm(n, mean = 150, sd = 3)
x2 <- rnorm(n, mean = 100, sd = 2)
u  <- rnorm(n)
y  <- 5 + 2*x1 + 3*x2 + u

# Estimate with OLS:
reg1 <- lm(y ~ x1 + x2)
# Estimate with GLS
reg2 <- glm(y ~ x1 + x2, family=gaussian)

# Compare:
require(texreg)
screenreg(l = list(reg1, reg2))

=========================================
                Model 1      Model 2     
-----------------------------------------
(Intercept)        6.37 **       6.37 ** 
                  (2.20)        (2.20)   
x1                 1.99 ***      1.99 ***
                  (0.01)        (0.01)   
x2                 3.00 ***      3.00 ***
                  (0.02)        (0.02)   
-----------------------------------------
R^2                0.99                  
Adj. R^2           0.99                  
Num. obs.          1000          1000       
RMSE               1.00                  
AIC                           2837.66    
BIC                           2857.29    
Log Likelihood               -1414.83    
Deviance                       991.82    
=========================================
*** p < 0.001, ** p < 0.01, * p < 0.05

Jawaban yang lebih panjang; Fungsi glm cocok dengan model oleh MLE, namun, karena asumsi yang Anda buat tentang fungsi tautan (dalam hal ini normal), Anda berakhir dengan perkiraan OLS.

Repmat
sumber
+1, salah ketik pada kalimat terakhir. Asumsi normal adalah tentang distribusi kesalahan, bukan tentang fungsi tautan. Dalam contoh Anda, fungsi tautan default adalah "identitas". Bentuk yang lebih lengkap glmadalah glm(y ~ x1 + x2, family = gaussian(link = "identity")).
Paul
14

Dari jawaban @ Repmat, ringkasan model adalah sama, tetapi CI dari koefisien regresi dari confintsedikit berbeda antara lmdan glm.

> confint(reg1, level=0.95)
               2.5 %    97.5 %
(Intercept) 2.474742 11.526174
x1          1.971466  2.014002
x2          2.958422  3.023291
> confint(reg2, level=0.95)
Waiting for profiling to be done...
               2.5 %    97.5 %
(Intercept) 2.480236 11.520680
x1          1.971492  2.013976
x2          2.958461  3.023251

t distribusi digunakan lmsementara distribusi normal digunakan glmketika membangun interval.

> beta <- summary(reg1)$coefficients[, 1]
    > beta_se <- summary(reg1)$coefficients[, 2]
> cbind(`2.5%` = beta - qt(0.975, n - 3) * beta_se, 
        `97.5%` = beta + qt(0.975, n - 3) * beta_se) #t
                2.5%     97.5%
(Intercept) 2.474742 11.526174
x1          1.971466  2.014002
x2          2.958422  3.023291
> cbind(`2.5%` = beta - qnorm(0.975)*beta_se, 
        `97.5%` = beta + qnorm(0.975)*beta_se) #normal
                2.5%     97.5%
(Intercept) 2.480236 11.520680
x1          1.971492  2.013976
x2          2.958461  3.023251
pe-pe-rry
sumber