Warum erhalte ich in OL die gleichen Ergebnisse für OLS und GLS?

8

Wenn ich diesen Code ausführe:

require(nlme)

a <- matrix(c(1,3,5,7,4,5,6,4,7,8,9))

b <- matrix(c(3,5,6,2,4,6,7,8,7,8,9))

res <- lm(a ~ b)

print(summary(res))

res_gls <- gls(a ~ b)

print(summary(res_gls))

Ich bekomme die gleichen Koeffizienten und die gleiche statistische Signifikanz für die Koeffizienten:

Loading required package: nlme

Call:
lm(formula = a ~ b)

Residuals:
    Min      1Q  Median      3Q     Max 
-2.7361 -1.1348 -0.2955  1.2463  3.8234 

Coefficients:
            Estimate Std. Error t value Pr(>|t|)  
(Intercept)   2.0576     1.8732   1.098   0.3005  
b             0.5595     0.2986   1.874   0.0937 .
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1   1 

Residual standard error: 2.088 on 9 degrees of freedom
Multiple R-squared: 0.2807, Adjusted R-squared: 0.2007 
F-statistic: 3.512 on 1 and 9 DF,  p-value: 0.09371 

Generalized least squares fit by REML
  Model: a ~ b 
  Data: NULL 
      AIC      BIC    logLik
  51.0801 51.67177 -22.54005

Coefficients:
                Value Std.Error  t-value p-value
(Intercept) 2.0576208 1.8731573 1.098477  0.3005
b           0.5594796 0.2985566 1.873948  0.0937

 Correlation: 
  (Intr)
b -0.942

Standardized residuals:
       Min         Q1        Med         Q3        Max 
-1.3104006 -0.5434780 -0.1415446  0.5968911  1.8311781 

Residual standard error: 2.087956 
Degrees of freedom: 11 total; 9 residual

Warum passiert dies? In welchen Fällen stimmen die OLS-Schätzungen mit den GLS-Schätzungen überein?

Akavall
quelle
5
Ein GLS-Modell ermöglicht es, dass die Fehler korreliert werden und / oder ungleiche Varianzen aufweisen. Wenn Sie nicht eine solche Korrelation oder Differenz von Residualvarianz mit den Optionen angeben correlationoder weightsin der glsFunktion sind die Ergebnisse von GLS gleich jenen aus lm.
COOLSerdash
2
OK, danke das macht Sinn. Im Grunde hatte ich die gleichen Ergebnisse, weil ich sagte gls, ich solle mich so verhalten lm. Eine andere Frage ist, was ich für correlationund setzen sollte weights.
Akavall

Antworten:

13

Sie haben die gleichen Ergebnisse erhalten, weil Sie in der glsFunktion keine spezielle Varianz oder Korrelationsstruktur angegeben haben . Ohne solche Optionen verhält sich ein GLS wie ein OLS. Der Vorteil eines GLS-Modells gegenüber einer normalen Regression besteht in der Möglichkeit, eine Korrelationsstruktur anzugeben (Option correlation) oder zuzulassen, dass sich die Restvarianz unterscheidet (Option weights). Lassen Sie mich dies anhand eines Beispiels zeigen.

library(nlme)

set.seed(1500)

x <- rnorm(10000,100,12) # generate x with arbitrary values

y1 <- 10 + 15*x + rnorm(10000,0,5) # the first half of the dataset

y2 <-  -2 - 5*x + rnorm(10000,0,15) # the 2nd half of the data set with 3 times larger residual SD (15 vs. 5)

y <- c(y1, y2)
x.new <- c(x, x)

dummy.var <- c(rep(0, length(y1)), rep(1, length(y2))) # dummy variable to distinguish the first half of the dataset (y1) from the second (y2)

# Calculate a normal regression model   

lm.mod <- lm(y~x.new*dummy.var)

summary(lm.mod)

Coefficients:
                 Estimate Std. Error   t value Pr(>|t|)    
(Intercept)      10.27215    0.94237    10.900   <2e-16 ***
x.new            14.99691    0.00935  1603.886   <2e-16 ***
dummy.var       -12.07076    1.33272    -9.057   <2e-16 ***
x.new:dummy.var -19.99891    0.01322 -1512.387   <2e-16 ***

# Calculate a GLS without any options

gls.mod.1 <- gls(y~x.new*dummy.var)

summary(gls.mod.1)

Coefficients:
                    Value Std.Error    t-value p-value
(Intercept)      10.27215 0.9423749    10.9003       0
x.new            14.99691 0.0093504  1603.8857       0
dummy.var       -12.07076 1.3327194    -9.0572       0
x.new:dummy.var -19.99891 0.0132234 -1512.3868       0

# GLS again, but allowing different residual variance for y1 and y2

gls.mod.2 <- gls(y~x.new*dummy.var, weights=varIdent(form=~1|dummy.var))

summary(gls.mod.2)

 Parameter estimates:
       0        1 
1.000000 2.962565 

Coefficients:
                    Value Std.Error   t-value p-value
(Intercept)      10.27215 0.4262268    24.100       0
x.new            14.99691 0.0042291  3546.144       0
dummy.var       -12.07076 1.3327202    -9.057       0
x.new:dummy.var -19.99891 0.0132234 -1512.386       0

# Perform a likelihood ratio test

anova(gls.mod.1, gls.mod.2)

          Model df      AIC      BIC    logLik   Test  L.Ratio p-value
gls.mod.1     1  5 153319.4 153358.9 -76654.69                        
gls.mod.2     2  6 143307.2 143354.6 -71647.61 1 vs 2 10014.15  <.0001

Das erste GLS-Modell ( gls.mod.1) und das normale lineare Regressionsmodell ( lm.mod) liefern genau die gleichen Ergebnisse. Das GLS-Modell, das unterschiedliche Standardabweichungen ( gls.mod.2) berücksichtigt, schätzt, dass die Rest-SD y2etwa dreimal so groß ist wie die Rest-SD, y1die genau das ist, was wir bei der Generierung der Daten angegeben haben. Die Regressionskoeffizienten sind praktisch gleich, aber die Standardfehler haben sich geändert. Der Likelihood-Ratio-Test (und AIC) legen nahe, dass das GLS-Modell mit den verschiedenen Restvarianzen ( gls.mod.2) besser zu den Daten passt als das normale Modell ( lm.mododer gls.mod.1).


Varianz- und Korrelationsstrukturen in gls

Sie können in der glsFunktion und der Option mehrere Varianzstrukturen angeben weights. Siehe hier für eine Liste. Eine Liste der Korrelationsstrukturen für die Option finden correlationSie hier .

COOLSerdash
quelle
Was bestimmt die zu wählende Varianzstruktur?
Rafael
@ Rafael In diesem Fall habe ich die Daten simuliert und wusste, welche Varianzstruktur zu nehmen ist. In der Praxis würde ich verschiedene Varianzstrukturen ausprobieren, die auf Fachwissen und explorativen Grafiken basieren. Die verschiedenen Modelle mit unterschiedlichen Varianzstrukturen können dann unter Verwendung von Likelihood-Ratio-Tests verglichen werden. Ich weiß nicht, ob es ein empfohlenes "Goldstandard" -Verfahren zur Auswahl der Varianzstruktur gibt.
COOLSerdash
Hallo COOLSerdash, danke für deine Antwort. Ich werde verschiedene Strukturen und Modellvergleiche mit dem LR-Test ausprobieren.
Rafael
1

und um es klar zu machen, können Sie im Falle einer seriellen Korrelation der Residuen einfach die OLS-Schätzung verwenden, z. B. gls(..., cor=corAR1(0.6))hier die 0,6 sowie die Reihenfolge von OLS, können Sie sie mit der arFunktion für die Residuen berechnen von OLS

Wiktor Olszowy
quelle