1

这个问题是由此处发布的另一个问题引起的：non-conformable arguments error from lmer when trying to extract information from the model matrix

当尝试从包含因子变量的 lmer 模型中获取预测均值时，输出会根据因子变量的指定方式而有所不同。

我有一个可变年龄组，可以使用“15 岁以下儿童”、“15-49 岁成人”、“50 岁以上老人”或“0-15 岁”、“15-49 岁”、“50 岁以上”组来指定你”。我的选择很重要，因为对于前者，标签的字母顺序不同于级别的数字顺序。为了说明这一点，我再次使用了睡眠数据。 

library(lme4)
sleep <- as.data.frame(sleepstudy)   #import the sleep data

我必须为年龄创建一个变量。

set.seed(13)  #set a seed for creating a new variable, age
sleep$age <- sample(1:3,length(sleep),rep=TRUE) #create a new variable, age
sleep$agegroup1 <- factor(sleep$age, levels = c(1,2,3), 
        labels = c("Children <15 years", "Adults 15-49 years", "Elderly 50+ years"))
table(sleep$agegroup)  #should have 3 age groups

运行模型

m1 <- lmer(Reaction ~ Days + agegroup1 + Days:agegroup1 + (Days | Subject), sleep) 
summary(m1)

# New data frame for predicted means
d <- seq(0,9,1)  # make a vector of days = 0 to 9
newdat1 <- data.frame(Days=d,      
                          agegroup1=factor(rep(levels(sleep$agegroup1),length(d))))
newdat1 <- newdat1[order(newdat1$Days,newdat1$agegroup1),]   #order by Days 
mm <- model.matrix(formula(m1,fixed.only=TRUE)[-2], newdat1)  #create the matrix

现在，我尝试使用模型矩阵和预测函数输出预测均值：

newdat1$mm <- mm%*%fixef(m1)    
newdat1$predict <- predict(m1, newdata=newdat1, re.form=NA)
head(newdat1)

这里，来自模型矩阵的预测均值和预测函数是不同的；成人和儿童年龄组是倒置的。 

   Days          agegroup1       mm  predict
11    0 Adults 15-49 years 252.2658 252.8241
1     0 Children <15 years 252.8241 252.2658
21    0  Elderly 50+ years 249.1254 249.1254
2     1 Adults 15-49 years 262.3326 263.2674
22    1 Children <15 years 263.2674 262.3326
12    1  Elderly 50+ years 260.0171 260.0171

如果我使用字母顺序与级别的数字顺序相同的因子标签再次运行此脚本，我会得到不同的结果：

#set new labels for agegroup
sleep$agegroup2 <- factor(sleep$age, levels = c(1,2,3), 
                        labels = c("0-15y", "15-49y", "50+y"))
m2 <- lmer(Reaction ~ Days + agegroup2 + Days:agegroup2 + (Days | Subject), sleep) 
summary(m2)

# New data frame for predicted means
d <- seq(0,9,1)  # make a vector of days = 0 to 9
newdat2 <- data.frame(Days=d,
                    agegroup2=factor(rep(levels(sleep$agegroup2),length(d))))
newdat2 <- newdat2[order(newdat2$Days,newdat2$agegroup2),]   #order by Days
mm <- model.matrix(formula(m2,fixed.only=TRUE)[-2], newdat2)
newdat2$mm <- mm%*%fixef(m2)   
newdat2$predict <- predict(m2, newdata=newdat2, re.form=NA)
head(newdat2)

这里，来自模型矩阵的预测均值和预测函数是相同的。

   Days agegroup2       mm  predict
1     0     0-15y 252.2658 252.2658
11    0    15-49y 252.8241 252.8241
21    0      50+y 249.1254 249.1254
22    1     0-15y 262.3326 262.3326
2     1    15-49y 263.2674 263.2674
12    1      50+y 260.0171 260.0171

Predict 似乎忽略了标签并专注于级别，而直接访问模型矩阵则正确地专注于标签。那么，我的问题是，在尝试使用模型矩阵时是否总是需要确保因子水平和标签具有相同的顺序？还是有其他方法可以克服这个问题？

4

1 回答 1

2

模型矩阵的列顺序和模型的固定效应的顺序必须匹配，以便正确地进行矩阵乘法以“手动”计算预测值。这意味着，是的，新数据集中因子水平的顺序必须与要使用的原始数据集中的顺序相同，model.matrix并且fixef与您所做的一样。

您可以通过在新数据集中设置因子水平的顺序来实现这一点。这是最容易做到的，只需使用原始数据集中的因子水平即可。例如，newdat1你可以这样做：

factor(rep(levels(sleep$agegroup1), length(d)), levels = levels(sleep$agegroup1)))

于 2015-12-18T16:35:17.990 回答