Czy jest jakaś różnica między lm a glm dla gaussowskiej rodziny glm?

45

W szczególności chcę wiedzieć, czy istnieje różnica między lm(y ~ x1 + x2)i glm(y ~ x1 + x2, family=gaussian). Myślę, że ten konkretny przypadek glm jest równy lm. Czy się mylę?

użytkownik3682457
źródło
10
Tak i nie. Jako model statystyczny nie. Jako obiekt dopasowany w R tak; różne zwracane obiekty, inny algorytm.
Przywróć Monikę - G. Simpson
3
Wydaje mi się, że jest tu pytanie statystyczne, a także kodujące R.
Silverfish,

Odpowiedzi:

48

Podczas gdy dla konkretnej formy modelu wymienionej w treści pytania (tj. lm(y ~ x1 + x2)Vs glm(y ~ x1 + x2, family=gaussian)) regresja i GLM są tym samym modelem, pytanie tytułowe wymaga czegoś nieco bardziej ogólnego:

Czy jest jakaś różnica między lm a glm dla gaussowskiej rodziny glm?

Na które odpowiedź brzmi „Tak!”.

Powodem, dla którego mogą się różnić, jest to, że można również określić funkcję łącza w GLM. Pozwala to dopasować określone formy nieliniowej zależności między y (a raczej jego średnią warunkową) a zmiennymi x ; chociaż można to również zrobić nls, nie ma potrzeby uruchamiania wartości początkowych, czasami zbieżność jest lepsza (również składnia jest nieco łatwiejsza).

Porównaj na przykład te modele (masz R, więc zakładam, że możesz sam je uruchomić):

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

yjaN.(β0+β1x1ja+β2)x2)ja,σ2))yjaN.(exp(β0+β1x1ja+β2)x2)ja),σ2)) i dopasowania są zasadniczo takie same w obrębie każdej pary.

Tak więc - w odniesieniu do pytania tytułowego - można dopasować znacznie szerszą gamę modeli Gaussa z GLM niż z regresją.

Glen_b
źródło
4
+1. Po jednej stronie obliczeniowej rzeczy, pomyślałbym również, że algorytm GLM użyłby jakiegoś wariantu IRWLS (w większości przypadków), podczas gdy LM polegałby na pewnym wariancie rozwiązania w formie zamkniętej.
usεr11852 mówi: Przywróć Monic
@ usεr11852 - Myślałem, że to EM, ale w tym przypadku mogą być takie same.
EngrStudent - Przywróć Monikę
1
Nie reaguje na widzenie „wartości odstających” (z wyjątkiem prawdopodobieństwa opisanego powyżej); ponowne ważenie wynika z efektu funkcji wariancji i przesunięcia w lokalnym przybliżeniu liniowym.
Glen_b,
1
tMASS::rlm
1
Możesz osiągnąć taki poziom solidności, który moim zdaniem zamierzasz na wiele sposobów. Jednak w przypadku modeli typu glms i regresji musisz wystrzegać się nie tylko wartości odstających w kierunku Y, ale także wpływowych wartości odstających, które mogą sprawić, że nie będą wyglądać nie na miejscu.
Glen_b 16.09.17
14

Krótka odpowiedź, są dokładnie takie same:

# 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

Dłuższa odpowiedź; Funkcja glm pasuje do modelu według MLE, jednak ze względu na założenie dotyczące funkcji łącza (w tym przypadku normalne), otrzymujesz oszacowania OLS.

Repmat
źródło
+1, literówka w ostatnim zdaniu. Normalne założenie dotyczy rozkładu błędów, a nie funkcji łącza. W twoim przykładzie domyślną funkcją łącza jest „tożsamość”. Bardziej kompletny formularz glmto glm(y ~ x1 + x2, family = gaussian(link = "identity")).
Paul
14

Z odpowiedzi @ Repmat, podsumowanie modelu jest takie samo, ale współczynniki CI współczynników regresji z confintsą nieco różne pomiędzy lmi 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
źródło