Duża różnica zdań w oszacowaniu nachylenia, gdy grupy są traktowane jako losowe vs. ustalone w modelu mieszanym

Rozumiem, że używamy modeli efektów losowych (lub efektów mieszanych), gdy uważamy, że niektóre parametry modelu zmieniają się losowo w zależności od czynnika grupującego. Chcę dopasować model, w którym odpowiedź została znormalizowana i wyśrodkowana (nie idealnie, ale całkiem blisko) w obrębie czynnika grupującego, ale zmienna niezależna xnie została w żaden sposób dostosowana. Doprowadziło mnie to do następującego testu (przy użyciu sfabrykowanych danych), aby upewnić się, że znajdę efekt, którego szukałem, jeśli rzeczywiście tam jest. Uruchomiłem jeden model efektów mieszanych z przypadkowym przechwytywaniem (pomiędzy grupami zdefiniowanymi przez f) i drugi model efektu stałego z czynnikiem f jako predyktorem efektu stałego. Użyłem pakietu R lmerdla modelu efektu mieszanego i funkcji podstawowejlm()dla modelu z efektem stałym. Poniżej znajdują się dane i wyniki.

Zauważ, że yniezależnie od grupy, waha się w okolicach 0. I xróżni się w zależności od ygrupy, ale różni się znacznie w zależności od grupy niży

> data
      y   x f
1  -0.5   2 1
2   0.0   3 1
3   0.5   4 1
4  -0.6  -4 2
5   0.0  -3 2
6   0.6  -2 2
7  -0.2  13 3
8   0.1  14 3
9   0.4  15 3
10 -0.5 -15 4
11 -0.1 -14 4
12  0.4 -13 4

Jeśli jesteś zainteresowany pracą z danymi, oto dput()wynik:

data<-structure(list(y = c(-0.5, 0, 0.5, -0.6, 0, 0.6, -0.2, 0.1, 0.4, 
-0.5, -0.1, 0.4), x = c(2, 3, 4, -4, -3, -2, 13, 14, 15, -15, 
-14, -13), f = structure(c(1L, 1L, 1L, 2L, 2L, 2L, 3L, 3L, 3L, 
4L, 4L, 4L), .Label = c("1", "2", "3", "4"), class = "factor")), 
.Names = c("y","x","f"), row.names = c(NA, -12L), class = "data.frame")

Dopasowanie modelu efektów mieszanych:

> summary(lmer(y~ x + (1|f),data=data))
Linear mixed model fit by REML 
Formula: y ~ x + (1 | f) 
   Data: data 
   AIC   BIC logLik deviance REMLdev
 28.59 30.53  -10.3       11   20.59
Random effects:
 Groups   Name        Variance Std.Dev.
 f        (Intercept) 0.00000  0.00000 
 Residual             0.17567  0.41913 
Number of obs: 12, groups: f, 4

Fixed effects:
            Estimate Std. Error t value
(Intercept) 0.008333   0.120992   0.069
x           0.008643   0.011912   0.726

Correlation of Fixed Effects:
  (Intr)
x 0.000 

Zwracam uwagę, że komponent wariancji przechwytywania jest szacowany na 0 i, co dla mnie ważne, xnie jest znaczącym predyktorem y.

Następnie dopasowuję model z efektem stałym fjako predyktor zamiast czynnika grupującego dla przypadkowego przechwytywania:

> summary(lm(y~ x + f,data=data))

Call:
lm(formula = y ~ x + f, data = data)

Residuals:
     Min       1Q   Median       3Q      Max 
-0.16250 -0.03438  0.00000  0.03125  0.16250 

Coefficients:
            Estimate Std. Error t value Pr(>|t|)    
(Intercept) -1.38750    0.14099  -9.841 2.38e-05 ***
x            0.46250    0.04128  11.205 1.01e-05 ***
f2           2.77500    0.26538  10.457 1.59e-05 ***
f3          -4.98750    0.46396 -10.750 1.33e-05 ***
f4           7.79583    0.70817  11.008 1.13e-05 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1 

Residual standard error: 0.1168 on 7 degrees of freedom
Multiple R-squared: 0.9484, Adjusted R-squared: 0.9189 
F-statistic: 32.16 on 4 and 7 DF,  p-value: 0.0001348 

Teraz zauważam, że zgodnie z oczekiwaniami xjest znaczącym predyktorem y.

To, czego szukam, to intuicja dotycząca tej różnicy. W jaki sposób błędne jest moje myślenie? Dlaczego niepoprawnie oczekuję znalezienia znaczącego parametru dla xobu tych modeli, ale faktycznie widzę go tylko w modelu ze stałym efektem?

Tutaj dzieje się kilka rzeczy. To interesujące kwestie, ale wyjaśnienie tego zajmie sporo czasu / miejsca.

Po pierwsze, wszystko to staje się o wiele łatwiejsze do zrozumienia, jeśli wykreślimy dane . Oto wykres punktowy, w którym punkty danych są pokolorowane według grup. Dodatkowo dla każdej grupy mamy osobną linię regresji specyficzną dla grupy, a także prostą linię regresji (ignorowanie grup) pogrubioną czcionką:

plot(y ~ x, data=dat, col=f, pch=19)
abline(coef(lm(y ~ x, data=dat)), lwd=3, lty=2)
by(dat, dat$f, function(i) abline(coef(lm(y ~ x, data=i)), col=i$f))

Model z efektem stałym

$x$$x$$x$$x$$x$$x$$x$$y$$t$

$x$$x$$x$lm()

Model mieszany

$x$$x$$x$$x$

$x$

Oto współczynniki dla prostego modelu regresji (przerywana pogrubiona linia na wykresie):

> lm(y ~ x, data=dat)

Call:
lm(formula = y ~ x, data = dat)

Coefficients:
(Intercept)            x  
   0.008333     0.008643  

Jak widać, współczynniki tutaj są identyczne z uzyskanymi w modelu mieszanym. Dokładnie tego się spodziewaliśmy, ponieważ, jak już zauważyłeś, mamy oszacowanie 0 wariancji dla losowych przechwyceń, co czyni wspomnianą wcześniej korelację stosunku / wewnątrz klasy 0. Więc szacunki modelu mieszanego w tym przypadku są tylko proste oszacowania regresji liniowej, a jak widać na wykresie, nachylenie tutaj jest znacznie mniej wyraźne niż nachylenia wewnątrz gromady.

To prowadzi nas do jednego końcowego zagadnienia koncepcyjnego ...

Dlaczego wariancja losowych przechwyceń jest szacowana na 0?

Odpowiedź na to pytanie może stać się trochę techniczna i trudna, ale postaram się, aby była tak prosta i nietechniczna, jak to tylko możliwe (ze względu na nas samych!). Ale może nadal będzie trochę rozwlekły.

$y$(lub bardziej poprawnie błędy modelu) wywołane przez strukturę klastrowania. Korelacja wewnątrz klasy mówi nam, jak przeciętnie podobne są dwa błędy wyciągnięte z tego samego klastra, w stosunku do średniego podobieństwa dwóch błędów wyciągniętych z dowolnego miejsca w zbiorze danych (tj. Może, ale nie musi być w tym samym klastrze). Dodatnia korelacja wewnątrz klasy mówi nam, że błędy z tego samego klastra wydają się być względnie bardziej do siebie podobne; jeśli narysuję jeden błąd z klastra i ma on wysoką wartość, to mogę spodziewać się ponad szansą, że następny błąd, który narysuję z tego samego klastra, również będzie miał wysoką wartość. Chociaż nieco mniej powszechne, korelacje wewnątrz klasy mogą być również ujemne; dwa błędy wyciągnięte z tego samego klastra są mniej podobne (tj. bardziej oddalone od siebie pod względem wartości) niż zwykle można oczekiwać w całym zbiorze danych.

Rozważany przez nas model mieszany nie korzysta z wewnątrzklasowej metody korelacji reprezentującej zależność w danych. Zamiast tego opisuje zależność w kategoriach składników wariancji . Wszystko jest w porządku, dopóki korelacja wewnątrz klasy jest dodatnia. W takich przypadkach korelację wewnątrz klasy można łatwo napisać w kategoriach składników wariancji, szczególnie jako wcześniej wspomniany stosunek losowej wariancji przechwytywania do całkowitej wariancji. (Zobacz stronę wiki dotyczącą korelacji wewnątrzklasowejwięcej informacji na ten temat.) Niestety modele składowe wariancji mają trudności z radzeniem sobie z sytuacjami, w których mamy ujemną korelację wewnątrz klasy. W końcu pisanie korelacji wewnątrzklasowej pod względem składników wariancji wymaga zapisania jej jako proporcji wariancji, a proporcje nie mogą być ujemne.

$y$$y$$y$, podczas gdy błędy zaczerpnięte z różnych klastrów będą miały bardziej umiarkowaną różnicę.) Tak więc twój model mieszany robi to, co w praktyce często robią modele mieszane w tym przypadku: daje oszacowania, które są równie spójne z ujemną korelacją wewnątrz klasy ponieważ może zbierać, ale zatrzymuje się na dolnej granicy 0 (to ograniczenie jest zwykle programowane w algorytmie dopasowania modelu). W rezultacie otrzymujemy szacunkową wariancję losowego przechwytywania wynoszącą 0, co nadal nie jest bardzo dobrym oszacowaniem, ale jest tak bliskie, jak to możliwe z tym typem modelu wariancji-składników.

Więc co możemy zrobić?

$x$

$x$

Aby to zrobić, bierzemy nasze $x$$x_b$$x$$x_w$$x$

> dat <- within(dat, x_b <- tapply(x, f, mean)[paste(f)])
> dat <- within(dat, x_w <- x - x_b)
> dat
      y   x f x_b x_w
1  -0.5   2 1   3  -1
2   0.0   3 1   3   0
3   0.5   4 1   3   1
4  -0.6  -4 2  -3  -1
5   0.0  -3 2  -3   0
6   0.6  -2 2  -3   1
7  -0.2  13 3  14  -1
8   0.1  14 3  14   0
9   0.4  15 3  14   1
10 -0.5 -15 4 -14  -1
11 -0.1 -14 4 -14   0
12  0.4 -13 4 -14   1
> 
> mod <- lmer(y ~ x_b + x_w + (1|f), data=dat)
> mod
Linear mixed model fit by REML 
Formula: y ~ x_b + x_w + (1 | f) 
   Data: dat 
   AIC   BIC logLik deviance REMLdev
 6.547 8.972  1.726   -23.63  -3.453
Random effects:
 Groups   Name        Variance Std.Dev.
 f        (Intercept) 0.000000 0.00000 
 Residual             0.010898 0.10439 
Number of obs: 12, groups: f, 4

Fixed effects:
            Estimate Std. Error t value
(Intercept) 0.008333   0.030135   0.277
x_b         0.005691   0.002977   1.912
x_w         0.462500   0.036908  12.531

Correlation of Fixed Effects:
    (Intr) x_b  
x_b 0.000       
x_w 0.000  0.000

$x_w$$x_b$$y$$x$$x$$x_b$$t$-statystyka jest większa. Nie jest to również zaskakujące, ponieważ wariancja rezydualna jest znacznie mniejsza w tym modelu mieszanym ze względu na przypadkowe efekty grupowe pochłaniające dużą wariancję, z którą musiał sobie poradzić prosty model regresji.

Wreszcie, nadal mamy oszacowanie 0 dla wariancji przypadkowych przechwyceń, z powodów, które rozwinąłem w poprzedniej sekcji. Nie jestem do końca pewien, co możemy z tym zrobić, przynajmniej bez przełączania się na oprogramowanie inne niż lmer(), i nie jestem też pewien, w jakim stopniu będzie to miało negatywny wpływ na nasze szacunki w tym ostatecznym mieszanym modelu. Być może inny użytkownik może wtrącić się z przemyśleniami na temat tego problemu.

Bibliografia

• Bell, A. i Jones, K. (2014). Objaśnienie efektów stałych: Modelowanie efektów losowych danych przekroju i panelu szeregów czasowych. Badania i metody nauk politycznych. PDF
• Bafumi, J., i Gelman, AE (2006). Dopasowanie modeli wielopoziomowych, gdy korelują predyktory i efekty grupowe. PDF

Po dłuższej kontemplacji wydaje mi się, że odkryłem własną odpowiedź. Wierzę, że ekonometryczny zdefiniowałby moją zmienną niezależną jako endogenną, a zatem byłby skorelowany zarówno ze zmiennymi niezależnymi, jak i zależnymi. W takim przypadku zmienne te są pomijane lub nie są obserwowane . Obserwuję jednak grupy, między którymi zmienna pominięta powinna się zmieniać.

Wierzę, że ekonometryczny sugerowałby model z efektem stałym . Oznacza to, że model zawiera atrapę dla każdego poziomu grupowania (lub równoważną specyfikację, która warunkuje model tak, że wiele manekinów grupujących nie jest wymaganych) w tym przypadku. Dzięki modelowi o stałym efekcie można mieć nadzieję, że wszystkie nieobserwowane i niezmienne w czasie zmienne mogą być kontrolowane przez warunkowanie w obrębie grupy (lub pomiędzy poszczególnymi) odmianami. Rzeczywiście, drugi model w moim pytaniu jest modelem o ustalonym efekcie i jako taki daje oszacowanie, którego się spodziewam.

Z zadowoleniem przyjmuję komentarze, które dodatkowo wyjaśnią tę okoliczność.