Dlaczego wprowadzenie losowego efektu nachylenia zwiększyło SE nachylenia?

Próbuję przeanalizować wpływ roku na zmienną logInd dla konkretnej grupy osób (mam 3 grupy). Najprostszy model:

> fix1 = lm(logInd ~ 0 + Group + Year:Group, data = mydata)
> summary(fix1)

Call:
lm(formula = logInd ~ 0 + Group + Year:Group, data = mydata)

Residuals:
    Min      1Q  Median      3Q     Max 
-5.5835 -0.3543 -0.0024  0.3944  4.7294 

Coefficients:
              Estimate Std. Error t value Pr(>|t|)    
Group1       4.6395740  0.0466217  99.515  < 2e-16 ***
Group2       4.8094268  0.0534118  90.044  < 2e-16 ***
Group3       4.5607287  0.0561066  81.287  < 2e-16 ***
Group1:Year -0.0084165  0.0027144  -3.101  0.00195 ** 
Group2:Year  0.0032369  0.0031098   1.041  0.29802    
Group3:Year  0.0006081  0.0032666   0.186  0.85235    
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1 

Residual standard error: 0.7926 on 2981 degrees of freedom
Multiple R-squared: 0.9717,     Adjusted R-squared: 0.9716 
F-statistic: 1.705e+04 on 6 and 2981 DF,  p-value: < 2.2e-16

Widzimy, że Grupa 1 znacznie spada, Grupy 2 i 3 rosną, ale nie tak znacząco.

Oczywiście jednostka powinna być efektem losowym, dlatego wprowadzam losowy efekt przechwytywania dla każdej osoby:

> mix1a = lmer(logInd ~ 0 + Group + Year:Group + (1|Individual), data = mydata)
> summary(mix1a)
Linear mixed model fit by REML 
Formula: logInd ~ 0 + Group + Year:Group + (1 | Individual) 
   Data: mydata 
  AIC  BIC logLik deviance REMLdev
 4727 4775  -2356     4671    4711
Random effects:
 Groups     Name        Variance Std.Dev.
 Individual (Intercept) 0.39357  0.62735 
 Residual               0.24532  0.49530 
Number of obs: 2987, groups: Individual, 103

Fixed effects:
              Estimate Std. Error t value
Group1       4.6395740  0.1010868   45.90
Group2       4.8094268  0.1158095   41.53
Group3       4.5607287  0.1216522   37.49
Group1:Year -0.0084165  0.0016963   -4.96
Group2:Year  0.0032369  0.0019433    1.67
Group3:Year  0.0006081  0.0020414    0.30

Correlation of Fixed Effects:
            Group1 Group2 Group3 Grp1:Y Grp2:Y
Group2       0.000                            
Group3       0.000  0.000                     
Group1:Year -0.252  0.000  0.000              
Group2:Year  0.000 -0.252  0.000  0.000       
Group3:Year  0.000  0.000 -0.252  0.000  0.000

Miał oczekiwany efekt - SE nachyleń (współczynniki Grupa 1-3: Rok) są teraz niższe, a resztkowe SE jest również niższe.

Poszczególne osoby różnią się także nachyleniem, dlatego wprowadziłem również losowy efekt nachylenia:

> mix1c = lmer(logInd ~ 0 + Group + Year:Group + (1 + Year|Individual), data = mydata)
> summary(mix1c)
Linear mixed model fit by REML 
Formula: logInd ~ 0 + Group + Year:Group + (1 + Year | Individual) 
   Data: mydata 
  AIC  BIC logLik deviance REMLdev
 2941 3001  -1461     2885    2921
Random effects:
 Groups     Name        Variance  Std.Dev. Corr   
 Individual (Intercept) 0.1054790 0.324775        
            Year        0.0017447 0.041769 -0.246 
 Residual               0.1223920 0.349846        
Number of obs: 2987, groups: Individual, 103

Fixed effects:
              Estimate Std. Error t value
Group1       4.6395740  0.0541746   85.64
Group2       4.8094268  0.0620648   77.49
Group3       4.5607287  0.0651960   69.95
Group1:Year -0.0084165  0.0065557   -1.28
Group2:Year  0.0032369  0.0075105    0.43
Group3:Year  0.0006081  0.0078894    0.08

Correlation of Fixed Effects:
            Group1 Group2 Group3 Grp1:Y Grp2:Y
Group2       0.000                            
Group3       0.000  0.000                     
Group1:Year -0.285  0.000  0.000              
Group2:Year  0.000 -0.285  0.000  0.000       
Group3:Year  0.000  0.000 -0.285  0.000  0.000

Ale teraz, wbrew oczekiwaniom, SE nachyleń (współczynniki Grupa 1-3: Rok) są teraz znacznie wyższe, nawet wyższe niż bez żadnego losowego efektu!

Jak to jest możliwe? Spodziewałbym się, że losowy efekt „zje” niewyjaśnioną zmienność i zwiększy „pewność” oszacowania!

Resztkowa SE zachowuje się jednak zgodnie z oczekiwaniami - jest niższa niż w modelu przechwytywania losowego.

Oto dane w razie potrzeby.

Edytować

Teraz zrozumiałem zdumiewający fakt. Jeśli wykonam regresję liniową dla każdej osoby osobno, a następnie uruchomię ANOVA na powstałych zboczach, otrzymam dokładnie taki sam wynik jak losowy model zbocza! Czy wiesz dlaczego?

indivSlope = c()
for (indiv in 1:103) {
    mod1 = lm(logInd ~ Year, data = mydata[mydata$Individual == indiv,])
    indivSlope[indiv] = coef(mod1)['Year']
}

indivGroup = unique(mydata[,c("Individual", "Group")])[,"Group"]


anova1 = lm(indivSlope ~ 0 + indivGroup)
summary(anova1)

Call:
lm(formula = indivSlope ~ 0 + indivGroup)

Residuals:
      Min        1Q    Median        3Q       Max 
-0.176288 -0.016502  0.004692  0.020316  0.153086 

Coefficients:
              Estimate Std. Error t value Pr(>|t|)
indivGroup1 -0.0084165  0.0065555  -1.284    0.202
indivGroup2  0.0032369  0.0075103   0.431    0.667
indivGroup3  0.0006081  0.0078892   0.077    0.939

Residual standard error: 0.04248 on 100 degrees of freedom
Multiple R-squared: 0.01807,    Adjusted R-squared: -0.01139 
F-statistic: 0.6133 on 3 and 100 DF,  p-value: 0.6079

Oto dane w razie potrzeby.

r mixed-model lme4-nlme random-effects-model

— Ciekawy
źródło

Potrzebujesz efektu ustalonego na rok, jeśli masz mieć efekt ustalony na rok: interakcja w grupie. Zasadniczo nie można dołączyć terminu interakcji bez uwzględnienia głównych efektów. Czy naprawdę uważasz, że nie ma stałego elementu efektu roku? A jeśli tak, to jak mógłby istnieć ustalony rok: interakcja w grupie?

— John

A dlaczego nie ma stałego przechwytywania? Możesz mieć zarówno stałe, jak i losowe.

— John

@John, ten model jest całkowicie ważny. Jest to tylko kwestia pożądanego kodowania zmiennej jakościowej. W ten sposób jest punktem przecięcia grupy , a jest nachyleniem w grupie . Jeśli uwzględniony zostanie główny efekt Roku i przecięcia, wówczas szacunki będą różnicami między przecięciem dla Grupy i Grupy 1, i podobnie ze spadkami. Group

i

$i$

i

$i$ Group

i

$i$ :Year

i

$i$

i

$i$

— Aniko,

@John, to nie jest temat na moje pytanie, niemniej: uwierz mi, to jest OK, przeprowadziłem z tym wiele eksperymentów. Mój pierwszy model lm jest całkowicie równoważny logInd ~ Year*Group, tylko współczynniki mają inny kształt, nic więcej. Wszystko zależy od Twojego gustu i kształtu współczynników, które lubisz. Nie ma wykluczenia „Głównego efektu roku” w moim pierwszym modelu, gdy piszesz ... logInd ~ Year*Grouprobi dokładnie to samo, wówczas Yearwspółczynnik nie jest głównym efektem, ale Group1: Year.

— Ciekawy

OK, schludnie, nie uważałem, że zarówno punkt przecięcia 0, jak i grupa są kategoryczne.

— John

Myślę, że problem dotyczy twoich oczekiwań :) Zauważ, że kiedy dodałeś losowy przechwycenie dla każdej osoby, standardowy błąd przechwyceń wzrósł. Ponieważ każda osoba może mieć własne przechwytywanie, średnia grupy jest mniej pewna. To samo stało się z losowym nachyleniem: nie szacuje się już jednego wspólnego (wewnątrz grupy) nachylenia, ale średnią różnych nachyleń.

EDYCJA: Dlaczego lepszy model nie daje dokładniejszej oceny?

Zastanówmy się nad tym na odwrót: dlaczego początkowy model nie docenia standardowego błędu? Zakłada niezależność obserwacji, które nie są niezależne. Drugi model rozluźnia to założenie (w sposób, który wpływa na przechwyty), a trzeci rozluźnia je dalej.

EDYCJA 2: związek z wieloma modelami specyficznymi dla pacjenta

Twoja obserwacja jest znaną właściwością (a gdybyś miał tylko dwa lata, wówczas model efektów losowych byłby równoważny sparowanemu testowi t). Nie sądzę, że dam sobie radę z prawdziwym dowodem, ale być może spisanie dwóch modeli sprawi, że relacja będzie wyraźniejsza. Zignorujmy zmienną grupującą, ponieważ skomplikowałoby to notację. Będę używać liter greckich do efektów losowych i liter łacińskich do efektów stałych.

Model efektów losowych to ( - podmiot, - replikacja w obrębie podmiotu): gdzie i . $i$ $j$

Y_{i j} = a + α_{i} + (b + β_{i}) x_{i j} + ϵ_{i j},

$Y_{ij} = a + \alpha_i + (b+\beta_i)x_{ij} + \epsilon_{ij},$

(α_{i}, β_{i})^{'} \sim N (0, Σ)

$(\alpha_i,\beta_i)'\sim N(0,\Sigma)$

ϵ_{i j} \sim N (0, σ^{2})

$\epsilon_{ij}\sim N(0,\sigma^2)$

Jeśli dopasujesz osobne modele dla każdego przedmiotu, to gdzie .

Y_{i j} = a_{i} + b_{i} x_{i j} + ϵ_{i j},

$Y_{ij} = a_i + b_i x_{ij}+ \epsilon_{ij},$

ϵ_{i j} \sim N (0, σ_{i}^{2})

$\epsilon_{ij}\sim N(0,\sigma_i^2)$

[Uwaga: poniższe czynności to po prostu ręczne tworzenie fal:]

Możesz zobaczyć wiele podobieństw między tymi dwoma modelami, gdzie odpowiada i do . Średnia odpowiada , ponieważ efekty losowe wynoszą średnio 0. Nieograniczona korelacja losowego przecięcia i nachylenia prowadzi do tego, że modele można po prostu dopasować osobno. Nie jestem pewien, w jaki sposób założenie single zazębia się ze specyficznym dla tematu , ale zakładam, że odbiera różnicę. $a_i$ $a+\alpha_i$ $b_i$ $b+\beta_i$ $b_i$ $b$ $\sigma$ $\sigma_i$ $\alpha_i$

— Aniko
źródło

Dzięki Aniko. Masz rację, moje obliczenia to potwierdzają, ale chciałbym zobaczyć, dlaczego ... Wydaje się to sprzeczne z intuicją. Ulepszyłem modele - wprowadzając efekty losowe lepiej opisałem strukturę błędów. Błąd resztkowy to potwierdza - jest coraz niższy. Więc z tymi lepszymi, bardziej precyzyjnymi modelami spodziewałbym się bardziej precyzyjnego nachylenia ... Wiem, że gdzieś się mylę, pomóżcie mi to zobaczyć.

— Ciekawy

Dzięki Aniko, to ciekawy punkt widzenia! Interesują mnie tylko stoki (Grupa *: Rok), nie przechwytuj tutaj ... więc mój pierwszy krok wprowadzenia losowego efektu itcept złagodził to założenie niezależności i doprowadził do niższego SE .. (stoku ..), a następnie następny krok było prawdopodobnie za dużo (??) i zrobiło się wręcz przeciwnie (jeszcze gorzej SE ..) .. może muszę o tym pomyśleć, dzięki.

— Ciekawy

Teraz jestem również zaskoczony bardzo interesującym faktem - zobacz moją edycję. Czy wiesz, dlaczego tak jest?

— Ciekawy

Nie sądzę, by założenie niepodległości było zbyt rozluźnione! Na początku było źle.

— Aniko,

Tomas, „precyzyjny” model nie oznacza, że szacunki będą bardziej precyzyjne. Jako skrajny przykład weź dowolny model wolny od danych, na przykład taki, który przewiduje, że wszystkie odpowiedzi są zerowe. Ten model jest absolutnie pewny w swojej ocenie zerowej. Dlatego jest tak precyzyjny, jak to tylko możliwe - ale prawdopodobnie jest również tak błędny, jak to możliwe. Zwiększenie zakresu modelu w celu dopasowania parametrów oznacza zwykle, że parametry te są dopasowywane z mniejszą dokładnością, a nie z większą dokładnością. Lepszy model, ponieważ może oszacować niepewność nie uchwyconą przez gorszy model, często ma większe błędy standardowe.

— whuber