Dodanie efektu losowego wpływa na oszacowanie współczynnika

Zawsze uczono mnie, że efekty losowe wpływają tylko na wariancję (błąd), a efekty stałe wpływają tylko na średnią. Ale znalazłem przykład, w którym efekty losowe wpływają również na średnią - oszacowanie współczynnika:

require(nlme)
set.seed(128)
n <- 100
k <- 5
cat <- as.factor(rep(1:k, each = n))
cat_i <- 1:k # intercept per kategorie
x <- rep(1:n, k)
sigma <- 0.2
alpha <- 0.001
y <- cat_i[cat] + alpha * x + rnorm(n*k, 0, sigma)
plot(x, y)

# simulate missing data
y[c(1:(n/2), (n*k-n/2):(n*k))] <- NA

m1 <- lm(y ~ x)
summary(m1)

m2 <- lm(y ~ cat + x)
summary(m2)

m3 <- lme(y ~ x, random = ~ 1|cat, na.action = na.omit)
summary(m3)

Widać, że szacowany współczynnik dla xmodelu m1wynosi -0,013780, podczas gdy dla modelu m3wynosi 0,0011713 - oba znacznie różnią się od zera.

Zauważ, że kiedy usuwam linię symulującą brakujące dane, wyniki są takie same (jest to pełna matryca).

Dlaczego?

PS: proszę zauważyć, że nie jestem zawodowym statystykiem, więc jeśli masz zamiar odpowiedzieć z dużą ilością matematyki, proszę również zrobić kilka prostych podsumowań dla manekinów :-)

„Zawsze uczono mnie, że efekty losowe wpływają tylko na wariancję (błąd), a efekty stałe wpływają tylko na średnią”.

Jak odkryłeś, dotyczy to tylko zbilansowanych, kompletnych (tj. Brakujących danych) zestawów danych bez ciągłych predyktorów. Innymi słowy, dla rodzajów danych / modeli omawianych w klasycznych tekstach ANOVA. W tych idealnych okolicznościach ustalone efekty i efekty losowe można oszacować niezależnie od siebie.

Gdy te warunki się nie utrzymują (ponieważ bardzo często nie występują w „świecie rzeczywistym”), ustalone i losowe efekty nie są niezależne. Co ciekawe, „nowoczesne” modele mieszane są szacowane przy użyciu iteracyjnych metod optymalizacji, a nie dokładnie rozwiązywane za pomocą odrobiny algebry macierzowej, jak w przypadku klasycznej mieszanej analizy ANOVA: aby oszacować ustalone efekty, musimy znamy losowe efekty, ale aby oszacować losowe efekty, musimy znać ustalone efekty! Bardziej adekwatne do obecnego pytania, oznacza to również, że gdy dane są niezrównoważone / niekompletne i / lub w modelu występują ciągłe predyktory, wówczas dostosowanie struktury efektów losowych modelu mieszanego może zmienić oszacowania stałej części modelu , i wzajemnie.

Edytuj 2016-07-05. Z komentarzy: „ Czy mógłbyś rozwinąć lub przytoczyć powód, dla którego ciągłe predyktory wpłyną na oszacowania stałej części modelu? ”

Oszacowania dla stałej części modelu będą zależeć od oszacowań losowej części modelu - to znaczy oszacowanych składników wariancji - jeśli (ale nie tylko jeśli) wariancja predyktorów różni się między klastrami. Co prawie na pewno będzie prawdą, jeśli któryś z predyktorów jest ciągły (przynajmniej w danych z „prawdziwego świata” - teoretycznie byłoby to niemożliwe, np. W skonstruowanym zbiorze danych).

Na pierwszym poziomie myślę, że ignorujecie kurczenie się w kierunku wartości populacji; „ nachylenia i przechwyty dla poszczególnych podmiotów z modelu efektów mieszanych są bliższe oszacowaniom populacji niż szacunki najmniejszych kwadratów w obrębie podmiotu. ” [ref. 1]. Poniższy link prawdopodobnie również będzie pomocny ( jakie są odpowiednie opisy, aby szukać moich modeli mieszanych? ), Patrz odpowiedź Mike Lawrence).

Co więcej, myślę, że masz mało szczęścia w przykładzie z zabawką, ponieważ masz idealnie zrównoważony projekt, który powoduje, że masz dokładnie takie same oszacowanie w przypadku braku brakujących wartości.

Wypróbuj następujący kod, który ma ten sam proces bez brakującej wartości:

 cat <- as.factor(sample(1:5, n*k, replace=T) ) #This should be a bit unbalanced.
 cat_i <- 1:k # intercept per kategorie
 x <- rep(1:n, k)
 sigma <- 0.2
 alpha <- 0.001
 y <- cat_i[cat] + alpha * x + rnorm(n*k, 0, sigma) 

 m1 <- lm(y ~ x)  
 m3 <- lme(y ~ x, random = ~ 1|cat, na.action = na.omit) 

 round(digits= 7,fixef(m3)) ==  round(digits=7, coef(m1)) #Not this time lad.
 #(Intercept)           x 
 #      FALSE       FALSE 

Gdzie teraz, ponieważ twój projekt nie jest idealnie zrównoważony, nie masz takich samych oszacowań współczynników.

W rzeczywistości, jeśli grasz w zgubny wzór wartości w głupi sposób (na przykład:), y[ c(1:10, 100 + 1:10, 200 + 1:10, 300 + 1:10, 400 +1:10)] <- NAwięc twój projekt jest nadal idealnie zrównoważony, otrzymasz znowu te same współczynniki.

 require(nlme)
 set.seed(128)
 n <- 100
 k <- 5
 cat <- as.factor(rep(1:k, each = n))
 cat_i <- 1:k # intercept per kategorie
 x <- rep(1:n, k)
 sigma <- 0.2
 alpha <- 0.001
 y <- cat_i[cat] + alpha * x + rnorm(n*k, 0, sigma)
 plot(x, y)

 # simulate missing data in a perfectly balanced way
 y[ c(1:10, 100 + 1:10, 200 + 1:10, 300 + 1:10, 400 +1:10)] <- NA

 m1 <- lm(y ~ x)  
 m3 <- lme(y ~ x, random = ~ 1|cat, na.action = na.omit) 

 round(digits=7,fixef(m3)) ==  round(digits=7, coef(m1)) #Look what happend now...
 #(Intercept)           x 
 #       TRUE        TRUE 

Nieznacznie wprawia Cię w błąd idealny projekt oryginalnego eksperymentu. Po włożeniu NA w sposób niezrównoważony zmieniłeś wzór „siły”, jaką poszczególne osoby mogły pożyczyć od siebie.

Krótko mówiąc, różnice, które widzisz, są spowodowane efektami skurczu, a konkretnie dlatego, że zniekształciłeś swój pierwotnie idealnie zrównoważony projekt z niedokładnie zrównoważonymi brakującymi wartościami.

Ref 1: Douglas Bates lme4: Modelowanie efektów mieszanych z R , strony 71-72