2

我有一些纵向数据，我想从中获得指定时间的预测均值。该模型包括 2 项，它们的交互作用和时间变量的样条项。当我尝试获得预测的平均值时，我得到“错误在 mm %*% fixef(m4) : non-conformable arguments”

我使用了 lmer 的 sleep 数据集来说明我的问题。首先，我导入数据并为我的交互创建一个变量“年龄”

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)

正是在这一点上，我得到了错误： mm %*% fixef(m4) 中的错误：不合格的参数

我可以使用 predict 来获取手段

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

效果很好，但我希望能够计算置信区间，所以我需要一个一致的矩阵。

我在某处读到问题可能是 lmer 创建了别名（我找不到那个帖子）。此评论是关于无法将 effect() 用于类似任务的问题。我不太明白如何克服这个问题。此外，我记得那个帖子有点旧，希望别名问题可能不再相关。

如果有人对我可能做错的事情提出建议，我将不胜感激。谢谢。

4

2 回答 2

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)))
于 2015-12-12T02:59:22.217 回答
0

错误是由于漂移分量引起的，如果你把

allowdrift=FALSE

进入您的 auto.arima 预测，它将被修复。

于 2019-03-18T16:59:43.573 回答