Existe-t-il une différence entre lm et glm pour la famille gaussienne de glm?

45

Plus précisément, je veux savoir s’il existe une différence entre lm(y ~ x1 + x2)et glm(y ~ x1 + x2, family=gaussian). Je pense que ce cas particulier de glm est égal à lm. Ai-je tort?

utilisateur3682457
la source
10
Oui et non. En tant que modèle statistique, non. En tant qu’objet ajusté dans R, oui; différents objets retournés, différents algorithmes utilisés.
Réintégrer Monica - G. Simpson
3
Il me semble qu’il s’agit là d’une question statistique, ainsi que d’une question portant le code R.
Silverfish

Réponses:

48

Alors que pour la forme spécifique de modèle mentionnée dans le corps de la question (c.-à-d. lm(y ~ x1 + x2)Vs glm(y ~ x1 + x2, family=gaussian)), régression et GLM sont le même modèle, la question du titre pose une question légèrement plus générale:

Existe-t-il une différence entre lm et glm pour la famille gaussienne de glm?

A quoi la réponse est "Oui!".

Ils peuvent être différents parce que vous pouvez également spécifier une fonction de lien dans le GLM. Cela vous permet d’adapter des formes particulières de relation non linéaire entre y (ou plutôt sa moyenne conditionnelle) et les variables x ; bien que vous puissiez le faire nlségalement, vous n'avez pas besoin de valeurs de départ, parfois la convergence est meilleure (la syntaxe aussi est un peu plus simple).

Comparez, par exemple, ces modèles (vous avez R, donc je suppose que vous pouvez les exécuter vous-même):

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

yiN(β0+β1x1i+β2x2i,σ2)yiN(exp(β0+β1x1i+β2x2i),σ2) et les ajustements sont essentiellement les mêmes dans chaque paire.

Donc, en ce qui concerne la question du titre, vous pouvez adapter un modèle GLM à une variété beaucoup plus large de modèles gaussiens qu’une régression.

Glen_b
la source
4
+1 Pour ce qui est du calcul, je pense aussi qu'un algorithme GLM utiliserait une variante IRWLS (dans la plupart des cas), tandis qu'un LM relèverait d'une variante de solution fermée.
usεr11852 dit Rétablir Monic
@ usεr11852 - J'aurais pensé que c'était un pays émergent, mais ce pourrait être la même chose dans ce cas.
EngrStudent
1
Il ne répond pas à voir des "valeurs aberrantes" (sauf via la probabilité décrite ci-dessus); La repondération est due à l'effet de la fonction de variance et au décalage de l'approximation linéaire locale.
Glen_b
1
tMASS::rlm
1
Vous pouvez obtenir le type de robustesse que vous souhaitez, à mon avis, à bien des égards. Cependant, avec les modèles de type glms et de régression, vous devez vous méfier non seulement des valeurs aberrantes dans la direction y, mais également des valeurs aberrantes influentes , qui risquent de ne pas sembler
déplacées
14

Réponse courte, ils sont exactement les mêmes:

# 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

Réponse plus longue; La fonction glm correspond au modèle de MLE, cependant, en raison de l'hypothèse que vous avez faite à propos de la fonction de liaison (dans ce cas, normale), vous vous retrouvez avec les estimations MLS.

Repmat
la source
+1, une faute de frappe dans la dernière phrase. L'hypothèse normale concerne la distribution des erreurs, pas la fonction de lien. Dans votre exemple, la fonction de lien par défaut est "identité". Un formulaire plus complet pour glmis glm(y ~ x1 + x2, family = gaussian(link = "identity")).
Paul
14

D'après la réponse de @ Repmat, le résumé du modèle est identique, mais les CI des coefficients de régression de confintsont légèrement différents entre lmet 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

tlmglm

> 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
la source