Ошибка несогласованных аргументов от lmer при попытке извлечь информацию из модельной матрицы

Question

У меня есть некоторые продольные данные, из которых я хотел бы получить прогнозируемые средние значения в указанное время. Модель включает в себя 2 члена, их взаимодействие и сплайн-член для переменной времени. Когда я пытаюсь получить предсказанные средние значения, я получаю "Ошибка в мм%*% fixef(m4): несоответствующие аргументы"

Я использовал данные сна из lmer, чтобы проиллюстрировать свою проблему. Сначала я импортирую данные и создаю переменную "возраст" для моего взаимодействия

sleep <- as.data.frame(sleepstudy)  #get the sleep data
# create fake variable for age with 3 levels
set.seed(1234567)
sleep$age <- as.factor(sample(1:3,length(sleep),rep=TRUE))

Затем я запускаю свою модель Lmer

library(lme4)
library(splines)
m4 <- lmer(Reaction ~ Days + ns(Days, df=4) + age + Days:age + (Days | Subject), sleep) 

Наконец, я создаю данные и матрицу, необходимые для получения предсказанных средств

#new data frame for predicted means
d <- c(0:9)  # make a vector of days = 0 to 9 to obtain predictions for each day
newdat <- as.data.frame(cbind(Days=d, age=rep(c(1:3),length(d))))
newdat$Days <- as.numeric(as.character(newdat$Days))
newdat$age <- as.factor(newdat$age)

# create a matrix 
mm<-model.matrix(~Days + ns(Days, df=4) + age + Days:age, newdat)  
newdat$pred<-mm%*%fixef(m4) 

В этот момент я получаю сообщение об ошибке: Ошибка в мм%*% fixef(m4): несоответствующие аргументы

Я могу использовать прогнозирование, чтобы получить средства

newdat$pred <- predict(m4, newdata=newdat, re.form=NA)

это прекрасно работает, но я хочу иметь возможность рассчитать доверительный интервал, поэтому мне нужна соответствующая матрица.

Я где-то читал, что проблема может быть в том, что lmer создает псевдонимы (я не могу найти этот пост). Этот комментарий был сделан относительно невозможности использовать effect() для аналогичной задачи. Я не мог понять, как преодолеть эту проблему. Более того, я вспоминаю, что пост был немного старым и надеялся, что проблема с псевдонимами, возможно, больше не актуальна.

Если у кого-то есть предложения по поводу того, что я могу делать неправильно, я буду признателен за отзыв. Благодарю.

Answer 1

Здесь есть пара вещей.

• вам нужно отбросить столбцы, чтобы сделать вашу матрицу модели соизмеримой с вектором фактического эффекта, который был фактически подогнан (т.е. соизмерим с матрицей модели, которая фактически использовалась для подгонки, после сброса коллинеарных столбцов)
• для дополнительной путаницы, вы случайно выбрали только возраст 2 и 3 (из возможных {1,2,3})

Я немного почистил код...

library("lme4")
library("splines")
sleep <- sleepstudy  #get the sleep data
set.seed(1234567)
## next line happens to sample only 2 and 3 ...
sleep$age <- as.factor(sample(1:3,length(sleep),rep=TRUE))
length(levels(sleep$age))  ## 2

Подходящая модель:

m4 <- lmer(Reaction ~ Days + ns(Days, df=4) +
    age + Days:age + (Days | Subject), sleep)
## message; fixed-effect model matrix is 
##    rank deficient so dropping 1 column / coefficient

Проверьте фиксированные эффекты:

f1 <- fixef(m4)
length(f1)  ## 7
f2 <- fixef(m4,add.dropped=TRUE)
length(f2)  ## 8

Мы могли бы использовать эту расширенную версию фиксированных эффектов (которая имеет NA ценность в этом), но это только испортило бы нас, распространяя NA значения через вычисления...

Проверьте матрицу модели:

X <- getME(m4,"X")
ncol(X)  ## 7
(which.dropped <- attr(getME(m4,"X"),"col.dropped"))
## ns(Days, df = 4)4 
##             6

Новый фрейм данных для прогнозируемых средних

d <- 0:9  
## best to use data.frame() directly, avoid cbind()
##   generate age based on *actual* levels in data
newdat <- data.frame(Days=d,
   age=factor(rep(levels(sleep$age),length(d))))

Создать матрицу:

mm <- model.matrix(formula(m4,fixed.only=TRUE)[-2], newdat)
mm <- mm[,-which.dropped]   ## drop redundant columns
## newdat$pred <- mm%*%fixef(m4)    ## works now

Добавлено sianagh: Код для получения доверительных интервалов и построения данных:

predFun <- function(x) predict(x,newdata=newdat,re.form=NA)
newdat$pred <- predFun(m4)
bb <- bootMer(m4,
   FUN=predFun,
    nsim=200)  
## nb. this produces an error message on its first run, 
## but not on subsequent runs (using the development version of lme4)
bb_ci <- as.data.frame(t(apply(bb$t,2,quantile,c(0.025,0.975))))
names(bb_ci) <- c("lwr","upr")
newdat <- cbind(newdat,bb_ci)

Участок:

plot(Reaction~Days,sleep)
with(newdat,
    matlines(Days,cbind(pred,lwr,upr),
            col=c("red","green","green"),
            lty=2,
            lwd=c(3,2,2)))

Answer 2

Ошибка вызвана из-за дрейфа компонента, если поставить

allowdrift=FALSE

в вашем прогнозе auto.arima это будет исправлено.