Gibt es einen Unterschied zwischen lm und glm für die Gaußsche Familie von glm?

45

Insbesondere möchte ich wissen, ob es einen Unterschied zwischen lm(y ~ x1 + x2)und gibt glm(y ~ x1 + x2, family=gaussian). Ich denke, dass dieser spezielle Fall von glm gleich lm ist. Liege ich falsch?

user3682457
quelle
10
Ja und nein. Als statistisches Modell, nein. Als angepasstes Objekt in R, ja; Verschiedene zurückgegebene Objekte, verschiedene verwendete Algorithmen.
Setzen Sie Monica - G. Simpson
3
Mir scheint, es gibt hier eine statistische und eine R-codierende Frage.
Silverfish

Antworten:

48

Während für die spezifische Form des Modells, die im Hauptteil der Frage (dh lm(y ~ x1 + x2)vs glm(y ~ x1 + x2, family=gaussian)) erwähnt wird, Regression und GLMs dasselbe Modell sind, werden in der Titelfrage etwas allgemeinere Fragen gestellt:

Gibt es einen Unterschied zwischen lm und glm für die Gaußsche Familie von glm?

Darauf lautet die Antwort "Ja!".

Der Grund dafür, dass sie unterschiedlich sein können, besteht darin, dass Sie im GLM auch eine Verknüpfungsfunktion angeben können . Auf diese Weise können Sie bestimmte Formen der nichtlinearen Beziehung zwischen (oder vielmehr dem bedingten Mittelwert) und den Variablen anpassen. Auch wenn Sie dies tun können, brauchen Sie keine Startwerte, manchmal ist die Konvergenz besser (auch die Syntax ist etwas einfacher).yxnls

Vergleichen Sie zum Beispiel diese Modelle (Sie haben R, ich gehe also davon aus, dass Sie diese selbst ausführen können):

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

Beachten Sie, dass das erste Paar dasselbe Modell ist ( ) und das zweite Paar dasselbe Modell ist ( und die Anpassungen sind im Wesentlichen innerhalb jedes Paares gleich.yiN(β0+β1x1i+β2x2i,σ2)yiN(exp(β0+β1x1i+β2x2i),σ2)

In Bezug auf die Titelfrage können Sie also eine wesentlich größere Vielfalt von Gaußschen Modellen mit einer GLM ausrüsten als mit einer Regression.

Glen_b
quelle
4
+1. Eine rechnerische Seite von Dingen ist, dass ein GLM-Algorithmus (in den meisten Fällen) eine IRWLS-Variante verwendet, während ein LM eine geschlossene Lösungsvariante verwendet.
usεr11852 sagt Reinstate Monic
@ usεr11852 - Ich hätte gedacht, dass es EM ist, aber in diesem Fall könnten sie dasselbe sein.
EngrStudent - Wiedereinsetzung von Monica
1
Es reagiert nicht auf das Erkennen von "Ausreißern" (außer über die oben beschriebene Wahrscheinlichkeit). Die Neugewichtung beruht auf dem Effekt der Varianzfunktion und der Verschiebung der lokalen linearen Approximation.
Glen_b
1
@ ChrisChiasson: +1 zu Glen_bs Kommentar. Wie bereits erwähnt, hängt dies nicht mit der Robustheit des Algorithmus bei Ausreißern zusammen. Sie könnten verschiedene Familien erkunden wollen (zB in geeigneter Weise skaliert. ; -distributions oder einen Huber Verlust Check - sorry nur bekomme Online nach ein paar freien Tagen .. für mehr dazu)tMASS::rlm
usεr11852 sagen wieder einzusetzen Monic
1
Sie könnten die Robustheit erreichen, die Sie in vielerlei Hinsicht beabsichtigen. Bei Glms- und Regressionstyp-Modellen müssen Sie jedoch nicht nur auf Ausreißer in y-Richtung achten, sondern auch auf einflussreiche Ausreißer, die sich nicht
verstellen können
14

Kurze Antwort, sie sind genau das gleiche:

# 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

Längere Antwort; Die glm-Funktion passt zum Modell von MLE. Aufgrund der Annahme, die Sie bezüglich der Verknüpfungsfunktion getroffen haben (in diesem Fall normal), erhalten Sie jedoch die OLS-Schätzungen.

Repmat
quelle
+1, ein Tippfehler im letzten Satz. Die normale Annahme betrifft die Fehlerverteilung, nicht die Verknüpfungsfunktion. In Ihrem Beispiel lautet die Standardverbindungsfunktion "Identität". Ein vollständigeres Formular für glmist glm(y ~ x1 + x2, family = gaussian(link = "identity")).
Paul
14

Ausgehend von der Antwort von @ Repmat ist die Modellzusammenfassung dieselbe, aber die CIs der Regressionskoeffizienten von unterscheiden confintsich geringfügig zwischen lmund 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 -Verteilung wird verwendet , lmwährend eine normale Verteilung wird verwendet , in , glmwenn die Intervalle zu konstruieren.

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