¿Hay alguna diferencia entre lm y glm para la familia gaussiana de glm?

45

Específicamente, quiero saber si hay una diferencia entre lm(y ~ x1 + x2)y glm(y ~ x1 + x2, family=gaussian). Creo que este caso particular de glm es igual a lm. ¿Me equivoco?

usuario3682457
fuente
10
Si y no. Como modelo estadístico, no. Como un objeto ajustado en R, sí; diferentes objetos devueltos, diferentes algoritmos utilizados.
Restablece a Monica - G. Simpson el
3
Me parece que hay una pregunta estadística aquí, así como una pregunta de codificación R.
Silverfish

Respuestas:

48

Mientras que para la forma específica de modelo mencionada en el cuerpo de la pregunta (es decir, lm(y ~ x1 + x2)vs glm(y ~ x1 + x2, family=gaussian)), la regresión y los GLM son el mismo modelo, la pregunta del título hace algo un poco más general:

¿Hay alguna diferencia entre lm y glm para la familia gaussiana de glm?

A lo que la respuesta es "¡Sí!".

La razón por la que pueden ser diferentes es porque también puede especificar una función de enlace en el GLM. Esto le permite ajustar formas particulares de relación no lineal entre (o más bien su media condicional) y las variables ; Si bien también puede hacer esto , no hay necesidad de valores iniciales, a veces la convergencia es mejor (también la sintaxis es un poco más fácil).yxnls

Compare, por ejemplo, estos modelos (tiene R, así que supongo que puede ejecutarlos usted mismo):

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

Tenga en cuenta que el primer par es el mismo modelo ( ), y el segundo par es el mismo modelo ( y los ajustes son esencialmente los mismos dentro de cada par.yiN(β0+β1x1i+β2x2i,σ2)yiN(exp(β0+β1x1i+β2x2i),σ2)

Entonces, en relación con la pregunta del título, puede ajustar una variedad sustancialmente más amplia de modelos gaussianos con un GLM que con una regresión.

Glen_b
fuente
44
+1. En el lado computacional de las cosas, también pensaría que un algoritmo GLM usaría alguna variante IRWLS (en la mayoría de los casos) mientras que un LM transmitiría una variante de solución de forma cerrada.
usεr11852 dice Reinstate Monic
@ usεr11852 - Pensé que era EM, pero podrían ser lo mismo en este caso.
EngrStudent - Restablece a Monica
1
No responde a ver "valores atípicos" (excepto a través de la probabilidad descrita anteriormente); la reponderación se debe al efecto de la función de varianza y al cambio en la aproximación lineal local.
Glen_b
1
@ChrisChiasson: +1 al comentario de Glen_b. Como se mencionó, esto no está relacionado con la robustez del algoritmo en presencia de valores atípicos. Es posible que desee explorar las diferentes familias (por ejemplo, escala adecuada. ; -distributions, o una pérdida Huber cheque - lo siento acabo de línea después de un par de días de descanso .. Para más información sobre esto)tMASS::rlm
usεr11852 dice Restablecer Monic
1
Podría lograr el tipo de robustez que creo que tiene la intención de varias maneras. Sin embargo, con los modelos de tipo glms y de regresión, debe tener cuidado no solo de los valores atípicos en la dirección y sino de los valores atípicos influyentes , que pueden hacer que no se vean fuera de lugar ...
Glen_b
14

Respuesta corta, son exactamente lo mismo:

# 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

Respuesta más larga; La función glm se ajusta al modelo de MLE, sin embargo, debido a la suposición que hizo sobre la función de enlace (en este caso normal), termina con las estimaciones de OLS.

Repmat
fuente
+1, un error tipográfico en la última oración. La suposición normal es sobre la distribución de errores, no sobre la función de enlace. En su ejemplo, la función de enlace predeterminada es "identidad". Una forma más completa para glmes glm(y ~ x1 + x2, family = gaussian(link = "identity")).
Paul
14

De la respuesta de @ Repmat, el resumen del modelo es el mismo, pero los IC de los coeficientes de regresión de confintson ligeramente diferentes entre lmy 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 distribución se usa lmmientras que la distribución normal se usa glmal construir los 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
fuente