Existe alguma diferença entre lm e glm para a família gaussiana de glm?

45

Especificamente, quero saber se há uma diferença entre lm(y ~ x1 + x2)e glm(y ~ x1 + x2, family=gaussian). Penso que este caso particular de glm é igual a lm. Estou errado?

user3682457
fonte
10
Sim e não. Como modelo estatístico, não. Como um objeto ajustado em R, sim; objetos retornados diferentes, algoritmo diferente usado.
Reintegrar Monica - G. Simpson
3
Parece-me que há uma questão estatística aqui, assim como uma questão de codificação R.
Silverfish

Respostas:

48

Enquanto para a forma específica de modelo mencionada no corpo da pergunta (ou seja, lm(y ~ x1 + x2)vs glm(y ~ x1 + x2, family=gaussian)), regressão e GLMs são o mesmo modelo, a pergunta do título faz algo um pouco mais geral:

Existe alguma diferença entre lm e glm para a família gaussiana de glm?

Para o qual a resposta é "Sim!".

A razão pela qual eles podem ser diferentes é porque você também pode especificar uma função de link no GLM. Isso permite que você ajuste formas particulares de relacionamento não linear entre (ou melhor, sua média condicional) e as variáveis ; Embora você possa fazer isso também, não há necessidade de iniciar valores, às vezes a convergência é melhor (também a sintaxe é um pouco mais fácil).yxnls

Compare, por exemplo, esses modelos (você tem R, então eu presumo que você possa executá-los):

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))

Observe que o primeiro par é o mesmo modelo ( ) e o segundo par é o mesmo modelo ( e os ajustes são essencialmente os mesmos em cada par.yiN(β0+β1x1i+β2x2i,σ2)yiN(exp(β0+β1x1i+β2x2i),σ2)

Portanto - em relação à questão do título - você pode ajustar uma variedade substancialmente mais ampla de modelos gaussianos com um GLM do que com regressão.

Glen_b
fonte
4
+1. No lado computacional, eu também pensaria que um algoritmo GLM usaria alguma variante IRWLS (na maioria dos casos), enquanto um LM retransmitiria alguma variante de solução de forma fechada.
usεr11852 diz Reinstate Monic
@ usεr11852 - Eu pensaria que era EM, mas eles podem ser a mesma coisa neste caso.
EngrStudent - Restabelece Monica
1
Não responde a ver "outliers" (exceto pela probabilidade como descrito acima); a reponderação se deve ao efeito da função de variação e à mudança na aproximação linear local.
#
1
@ChrisChiasson: +1 ao comentário de Glen_b. Como mencionado, isso não está relacionado à robustez do algoritmo na presença de outliers. Você pode querer explorar famílias diferentes (por exemplo, adequadamente dimensionados. ; -distributions, ou uma perda Huber cheque - desculpe ficou on-line depois de alguns dias de folga .. para saber mais sobre isso)tMASS::rlm
usεr11852 diz Reintegrar Monic
1
Você poderia obter o tipo de robustez que acho que pretende de várias maneiras. No entanto, com MLG e modelos tipo de regressão, você tem que tomar cuidado não apenas de outliers em direção y mas de influentes de outliers, que podem tornar-se não olhar para fora do lugar ..
Glen_b
14

Resposta curta, eles são exatamente os mesmos:

# 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

Resposta mais longa; A função glm se ajusta ao modelo pelo MLE, no entanto, devido à suposição que você fez sobre a função de link (neste caso, normal), você acaba com as estimativas do OLS.

Repmat
fonte
+1, um erro de digitação na última frase. A suposição normal é sobre a distribuição de erros, não sobre a função de link. No seu exemplo, a função de link padrão é "identidade". Um formulário mais completo para glmé glm(y ~ x1 + x2, family = gaussian(link = "identity")).
Paul
14

Na resposta de @ Repmat, o resumo do modelo é o mesmo, mas os ICs dos coeficientes de regressão de confintsão ligeiramente diferentes entre lme 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 distribuição é usada lmenquanto a distribuição normal é usada na glmconstrução dos intervalos.

> 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
fonte