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

У меня есть некоторые продольные данные, из которых я хотел бы получить прогнозируемые средние значения в определенное время. Модель включает 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 создает псевдонимы (не могу найти этот пост). Этот комментарий был сделан по поводу невозможности использования эффекта() для аналогичной задачи. Не совсем понял как победить эту проблему. Более того, я помню, что этот пост был немного старым, и я надеялся, что проблема псевдонима больше не актуальна.

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


person sianagh    schedule 11.12.2015    source источник
comment
В приведенном вами примере вы получаете сообщение fixed-effect model matrix is rank deficient so dropping 1 column / coefficient. Это означает, что количество фиксированных эффектов в модели меньше числа в mm, поэтому вы получаете сообщение об ошибке.   -  person aosmith    schedule 12.12.2015
comment
Да, аосмит, я понял это сообщение. Он исчезает, если модель включает только сплайн; то есть: m4 ‹- lmer(Реакция ~ ns(Days, df=4) + age + Days:age + (Days | Subject), sleep)   -  person sianagh    schedule 14.12.2015
comment
Извините, поправочка, в одномерной (только сплайновой) модели она исчезает, но все равно ошибка связана с использованием сплайнов.   -  person sianagh    schedule 14.12.2015
comment
ты читал мой ответ...?   -  person Ben Bolker    schedule 14.12.2015
comment
Да, @BenBolker, я читал это - у меня были некоторые ошибки, которые мне нужно было сначала исправить. Я ответил прямо на ваш ответ. Кроме того, я выявил 2 странности этих примерных данных, которые нуждаются в новых сообщениях, но я отмечаю их здесь для полноты: 1) Я думал, что использование семени будет генерировать одинаковые уровни возраста каждый раз, но это странно не так, следовательно ваш комментарий «## следующая строка встречается только с образцами 2 и 3…»; 2) прогнозировать, и этот метод дает разные результаты, когда факторная переменная использует метки, которые имеют алфавитный порядок, отличный от порядка нумерации.   -  person sianagh    schedule 14.12.2015


Ответы (2)


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

  • вам нужно удалить столбцы, чтобы ваша модельная матрица соответствовала вектору фиксированного эффекта, который был фактически подобран (т.е. соизмерим с модельной матрицей, которая фактически использовалась для подбора, после удаления коллинеарных столбцов)
  • для дополнительной путаницы у вас оказались только выборки возрастов 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)))
person Ben Bolker    schedule 12.12.2015
comment
Большое спасибо за ваш ответ. Действительно, это решение работает. Я планировал использовать bootMer для генерации интервалов прогнозирования, см. (r-bloggers.com/confidence-intervals-for-prediction-in-glmms). Однако эти данные, похоже, не подтверждают это. Если я запускаю следующее: predFun<-function(.) mm%*%fixef(.) bb<-bootMer(m4,FUN=predFun,nsim=200), я получаю 2 предупреждения: unable to evaluate scaled gradient и Model failed to converge: degenerate Hessian with 1 negative eigenvalues. Я думаю, проблема в том, что данные не могут поддерживать эту функцию. Ты согласен? - person sianagh; 14.12.2015
comment
На основе этого примера я исправил ошибку в lme4. Насколько сложно будет вам установить разрабатываемую версию с Github, чтобы протестировать ее (например, library("devtools"); install_github("lme4/lme4") ? - person Ben Bolker; 15.12.2015
comment
При использовании разрабатываемой версии при первом запуске появляется та же ошибка, но по-прежнему выдаются оценки. При последующих запусках нет предупреждения об ошибке. Оценки каждый раз разные (как и должно быть). Мне непонятно, почему ошибка исчезает после первого запуска. - person sianagh; 17.12.2015
comment
Для тех, кому интересно, я столкнулся с другой проблемой этого анализа, которая указывала на необходимость тщательного построения матрицы модели (если она используется). Пост здесь. stackoverflow.com/questions/34346755/ - person sianagh; 23.12.2015

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

allowdrift=FALSE

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

person Daniel Ghebreial Nieto    schedule 18.03.2019