Есть ли более простой способ вручную рассчитать оценочные средние значения группы с помощью model.matrix?

Я хочу рассчитать предполагаемые средние баллы группы в гауссовой регрессии 2x2 после получения коэффициентов регрессии. Вот данные игрушки. По 100 наблюдений региона - a и b - и пола - m и f. Я разработал оценки таким образом, чтобы разница между регионами a и b в среднем составляла 5 баллов, но между m и f разницы не было.

set.seed(1234)

d <- data.frame(region = factor(rep(letters[1:2],each=100)),
                sex = factor(rep(c("m", "f"),times=100)),
                score = round(x = c(rnorm(100, mean = 5, sd = 1),
                                    rnorm(100, mean = 10, sd = 1)),
                              digits = 1))

Теперь я буду использовать функцию model.matrix() для получения коэффициентов контраста для каждого наблюдения в зависимости от его принадлежности к группе. Я буду использовать кодирование лечения, то есть [0,1] с указанием региона a и пола m в качестве контрольных уровней для каждого из них.

model.matrix(object = score ~ region*sex,
             data = d,
             contrasts.arg = list(region = contr.treatment(nlevels(d$region)),
                                  sex = contr.treatment(nlevels(d$region)))) -> cmTreat

Теперь мы можем использовать матрицу модели непосредственно в регрессии с помощью функции lm(). Мы указываем 0 + terms, потому что матрица модели уже содержит перехват.

(lm(d$score ~ 0 + cmTreat) -> lmTreat)

# output
# Call:
#   lm(formula = d$score ~ 0 + cmTreat)
# 
# Coefficients:
# cmTreat(Intercept)       cmTreatregion2          cmTreatsex2  cmTreatregion2:sex2  
#              4.814                5.132                0.056                0.140 

Регрессия выявила основные эффекты и взаимодействия. Но что, если мы хотим получить оценочные предельные средние значения, а именно оценочное среднее значение в каждой «ячейке» 2 x 2: регион a — женский, регион a — мужской, регион b — женский, регион b — мужской.

Мы можем сделать это вручную с помощью атрибутов матрицы модели.

treatCoefs <- coef(lmTreat) # assign the vector of coefficients a name

# mean in region a female: intercept[1] + region[0] + sex[0] + region[0]*sex[0]
regionA_f <- treatCoefs[1] + treatCoefs[2]*attr(cmTreat, which = "contrasts")$region[,1][1] + treatCoefs[3]*attr(cmTreat, which = "contrasts")$sex[,1][1] + treatCoefs[4]*attr(cmTreat, which = "contrasts")$region[,1][1]*attr(cmTreat, which = "contrasts")$sex[,1][1]

# mean in region a male: intercept[1] + region[0] + sex[1] + region[0]*sex[1]
regionA_m <- treatCoefs[1] + treatCoefs[2]*attr(cmTreat, which = "contrasts")$region[,1][1] + treatCoefs[3]*attr(cmTreat, which = "contrasts")$sex[,1][2] + treatCoefs[4]*attr(cmTreat, which = "contrasts")$region[,1][1]*attr(cmTreat, which = "contrasts")$sex[,1][2]

# mean in region b female: : intercept[1] + region[1] + sex[0] + region[1]*sex[0]
regionB_f <- treatCoefs[1] + treatCoefs[2]*attr(cmTreat, which = "contrasts")$region[,1][2] + treatCoefs[3]*attr(cmTreat, which = "contrasts")$sex[,1][1] + treatCoefs[4]*attr(cmTreat, which = "contrasts")$region[,1][2]*attr(cmTreat, which = "contrasts")$sex[,1][1]

# mean in group b male: intercept[1] + region[1] + sex[1] + region[1]*sex[1]
regionB_m <-treatCoefs[1] + treatCoefs[2]*attr(cmTreat, which = "contrasts")$region[,1][2] + treatCoefs[3]*attr(cmTreat, which = "contrasts")$sex[,1][2] + treatCoefs[4]*attr(cmTreat, which = "contrasts")$region[,1][2]*attr(cmTreat, which = "contrasts")$sex[,1][2]

Теперь, если мы сравним фактические средние значения группы с расчетными средними (извиняюсь, люди, не занимающиеся этим вопросом)…

(library(dplyr)
d %>%
  group_by(region, sex) %>%
    summarise(actualMean = mean(score)) %>%
      add_column(estMeans = c(regionA_f, regionA_m, regionB_f, regionB_m))

# # A tibble: 4 × 4
# # Groups:   region [2]
# region  sex    actualMean estMeans
# <fct>   <fct>        <dbl>    <dbl>
# 1 a      f           4.81     4.81
# 2 a      m           4.87     4.87
# 3 b      f           9.95     9.95
# 4 b      m           10.1     10.1

Так что это отлично работает. "В чем проблема?" Я слышал, ты спрашиваешь. Итак, вы видели, сколько кода потребовалось для получения оценочных средних значений для каждой группы. И я могу это сделать. Но мне было интересно: «Есть ли более простой способ сделать это вручную?».

Я знаю, что могу использовать превосходный пакет emmeans Расса Лента и использую его часто, но мне хотелось научиться делать это вручную более элегантным способом. Я ничего не знаю о матричной алгебре и мало что знаю о контрастных матрицах. Я просто не могу избавиться от ощущения, что существует лучший способ (тот, метод которого может лучше адаптироваться к различным конструкциям и уровням).

п.с. этот вопрос, возможно, лучше подходил для перекрестной проверки, но я решил сначала попробовать здесь, поскольку он достаточно специфичен, чтобы гарантировать публикацию на SO.