Предикат и model.matrix дают разные прогнозируемые средства в пределах уровней факторной переменной

Question

Этот вопрос возник в результате другого вопроса, размещенного здесь: ошибка несогласованных аргументов от 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 игнорирует метки и фокусируется на уровнях, а прямой доступ к матрице модели правильно фокусируется на метках. Мой вопрос, поэтому, всегда ли необходимо гарантировать, что уровни факторов и метки имеют одинаковый порядок при попытке использовать матрицу модели? Или есть какой-то другой способ преодолеть эту проблему?

Answer 1

Порядок столбцов матрицы модели и фиксированных эффектов модели должен совпадать, чтобы правильно выполнить умножение матрицы для расчета прогнозируемых значений "вручную". Это означает, что да, порядок уровней фактора в новом наборе данных должен быть таким же, как и в исходном наборе данных для использования model.matrix а также fixef как ты.

Вы можете достичь этого, установив порядок уровней факторов в вашем новом наборе данных. Это проще всего сделать, просто используя уровни фактора из исходного набора данных. Например, в newdat1 ты можешь сделать:

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