从lmer中提取模型矩阵时出现非相容参数错误

Question

从lmer中提取模型矩阵时出现非相容参数错误

3

我有一些纵向数据，想要在指定时间得到预测平均值。该模型包括2个项、它们的交互和时间变量的样条项。当我试图获取预测平均值时，出现“mm %*% fixef(m4)非一致参数”的错误。

我使用了lmer的睡眠数据集来说明我的问题。首先，我导入数据并创建一个名为“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)

在这一点上，我遇到了错误：Error in mm %*% fixef(m4) : non-conformable arguments

我可以使用预测来得到平均值。

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

这段代码可以正常工作，但我希望能够计算置信区间，因此我需要一个相容的矩阵。

我在某个地方看到过这个问题可能是lmer创建了别名（我找不到那篇文章）。该评论是针对类似任务无法使用effect()函数而发出的。我无法完全理解如何克服这个问题。此外，我记得那篇文章有点旧了，希望别名问题可能已经不再相关。

如果有人有关于我可能做错的建议，我将感激您的反馈。谢谢。

- sianagh

2

在你提供的例子中，你得到了消息 fixed-effect model matrix is rank deficient so dropping 1 column / coefficient。这意味着模型中固定效应的数量少于 mm 中的数量，这就是为什么你会收到错误消息的原因。 - aosmith

是的，aosmith，我收到了那个消息。如果模型仅包括样条，则该消息会消失；即：m4 <- lmer(Reaction ~ ns(Days, df=4) + age + Days:age + (Days | Subject), sleep)。 - sianagh

抱歉，更正一下，在单变量（仅样条）模型中它会消失，但无论如何，这个错误都与使用样条有关。 - sianagh

你看了我的回答吗...？ - Ben Bolker

是的， @BenBolker，我确实看过了——我先需要解决一些错误。我已经直接回复了您的答案。此外，我还发现了这些示例数据的2个奇怪之处，需要新帖子进行讨论，但为了完整起见，在此记录如下：1）我认为使用种子会每次生成相同年龄水平，但这个奇怪的现象不是这样的，因此您有“##下一行只随机抽取了2和3 ...”的评论; 2）当因子变量使用字母顺序与编号顺序不同的标签时，预测和该方法给出的结果不同。 - sianagh

2个回答

0

错误是由于漂移组件引起的，如果您放置

allowdrift=FALSE

将其输入到您的auto.arima预测中，问题就会得到解决。

- Daniel Ghebreial Nieto

网页内容由stack overflow 提供, 点击上面的

可以查看英文原文，
原文链接

- Ben Bolker · Accepted Answer

这里有几个需要注意的地方。

你需要删除一些列，使得你的模型矩阵与实际拟合的固定效应向量相对应（也就是与实际用于拟合的模型矩阵相对应，在删除共线列之后）
更加混乱的是，你只采样了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)))