2015-03-25 17 views
5

Próbuję dopasować model efektów mieszanych, a następnie użyć tego modelu do wygenerowania szacunków dla nowego zestawu danych, który może mieć różne poziomy. Spodziewałem się, że szacunki dotyczące nowego zestawu danych wykorzystają średnią wartość oszacowanych parametrów, ale wydaje się, że tak nie jest. Oto przykład minimalna pracy:Prognozy z lme4 na nowych poziomach

library(lme4) 
d = data.frame(x = rep(1:10, times = 3), 
       y = NA, 
       grp = rep(1:3, each = 10)) 
d$y[d$grp == 1] = 1:10 + rnorm(10) 
d$y[d$grp == 2] = 1:10 * 1.5 + rnorm(10) 
d$y[d$grp == 3] = 1:10 * 0.5 + rnorm(10) 
fit = lmer(y ~ (1+x)|grp, data = d) 
newdata = data.frame(x = 1:10, grp = 4) 
predict(fit, newdata = newdata, allow.new.levels = TRUE) 

W tym przykładzie, ja zasadniczo zdefiniowania trzech grup z różnymi równania regresji (zboczach 1, 1,5 i 0,5). Jednak gdy próbuję przewidzieć nowy zestaw danych z niewidocznym poziomem, otrzymuję stałą wartość szacunkową. Spodziewałbym się, że oczekiwana wartość nachylenia i przechwycenia zostanie wykorzystana do wygenerowania prognoz dla tych nowych danych. Czy oczekuję niewłaściwej rzeczy? Albo, co robię źle z moim kodem?

+2

Wierzę, że 'predict.merMod' po prostu używa współczynników ze stałych części efektów dla nowych poziomów. 'y ~ x + (x | grp)' jest bardziej sensowną specyfikacją modelu. – Roland

+0

Ach, to ma sens! Jeśli dodasz to jako odpowiedź, zaakceptuję to. –

Odpowiedz

8

Generalnie nie uwzględniłbym losowego nachylenia bez uwzględnienia stałego nachylenia. Wygląda na to, że zgadza się ze mną predict.merMod, ponieważ wydaje się, że używa tylko stałych efektów do przewidywania nowych poziomów. Dokumentacja mówi, że "prognozy będą wykorzystywać bezwarunkowe (wartości na poziomie populacji) dane z wcześniej nieobserwowanymi poziomami", ale te wartości nie wydają się być oszacowane przy użyciu specyfikacji modelu.

Zatem proponuję tego modelu:

fit = lmer(y ~ x + (x|grp), data = d) 
newdata = data.frame(x = 1:10, grp = 4) 
predict(fit, newdata = newdata, allow.new.levels = TRUE) 
#  1   2   3   4   5   6   7   8   9  10 
#1.210219 2.200685 3.191150 4.181616 5.172082 6.162547 7.153013 8.143479 9.133945 10.124410 

To jest taka sama, jak tylko za pomocą efektów stałych części modelu

t(cbind(1, newdata$x) %*% fixef(fit)) 
#   [,1]  [,2] [,3]  [,4]  [,5]  [,6]  [,7]  [,8]  [,9] [,10] 
#[1,] 1.210219 2.200685 3.19115 4.181616 5.172082 6.162547 7.153013 8.143479 9.133945 10.12441 
5

Może to nie jest wystarczająco jasne, ale myślę, że dokumentacja dla ?predict.merMod stwierdza (rozsądnie) wyraźnie, co się dzieje, gdy allow.new.levels=TRUE. Chyba dwuznaczność może być w tym, co „bezwarunkowe (populacja poziomu) wartości” oznacza ...

allow.new.levels: logiczny, czy nowe poziomy (lub wartości NA) „NewData” są dozwolone. Jeśli FALSE (wartość domyślna), takie nowe wartości w 'newdata' spowodują błąd; jeśli PRAWDA, to predykcja użyje wartości bezwarunkowych (na poziomie populacji) dla danych z wcześniej niezaobserwowanymi poziomami (lub NA) z wcześniejszymi .

Powiązane problemy