我有一些纵向数据，希望在指定的时间获得预测的均值。该模型包括2个项，它们的相互作用和时间变量的样条项。当我尝试获取预测的均值时，出现“ mm％*％fixef（m4）中的错误：不合格的参数”
我使用来自lmer的sleep数据集来说明我的问题。首先，我导入数据并为交互创建一个变量“ age”
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）中的错误：参数不一致
我可以使用预测来获取手段
newdat$pred <- predict(m4, newdata=newdat, re.form=NA)
效果很好，但是我希望能够计算出置信区间，因此我需要一个合适的矩阵。
我在某处读到问题可能是lmer创建了别名（我找不到该帖子）。该评论是关于不能将effect（）用于类似任务的。我不太明白如何克服这个问题。此外，我记得帖子有点旧，希望别名问题不再重要。
如果有人对我可能做错的事情有任何建议，我将非常感谢您的反馈。谢谢。
这里有几件事。
我已经整理了一下代码...
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)))
本文收集自互联网，转载请注明来源。
如有侵权，请联系[email protected] 删除。
我来说两句