Из-за отсутствия лучшей формулировки: мне нужно создать новый фрейм данных из существующего фрейма данных, чтобы запустить модель lm/lmer, чтобы получить прогнозы модели.
Предположим, у меня есть:
x <- as.numeric(rep(1:6,5))
y <- as.numeric(rep(1:5,6))
int1 <- factor(rep(c("a","b"),15))
int2 <- factor(rep(c("11","12","13"),10))
g <- rep(c("f","m"),15)
df <- data.frame(x,y,int1,int2,g)
head(df)
# x y int1 int2 g
#1 1 1 a 11 f
#2 2 2 b 12 m
#3 3 3 a 13 f
#4 4 4 b 11 m
#5 5 5 a 12 f
#6 6 1 b 13 m
library(lme4)
mod <- lmer(y ~ x + int1*int2 + (1|g), data=df) #model I need
#OR:
#mod <- lm(y ~ x + int1*int2, data=df) #might be enough, don't need
#random effects in new df.
Теперь мне нужно создать новый df для запуска подобранной модели, чтобы получить прогнозы. Меня волнуют только термины взаимодействия (взаимодействие 2x3), но, видимо, новый df должен иметь все фиксированные эффекты, иначе выдает ошибку. Случайные эффекты исключены. Новый df имеет стандартные ошибки и границы доверительного интервала.
Это должно выглядеть примерно так:
> foo
int1/int2fit se lwr upr
11 68.86 2.03 64.91 72.86 #main effect
12 43.44 5.78 32.50 55.10 #main effect
13 38.77 4.14 31.12 47.19 #main effect
a 36.81 5.87 26.05 48.72 #main effect
b 34.58 3.59 27.55 41.71 #main effect
11a 28.04 4.40 19.87 37.31 #interaction term
11b 32.69 3.92 25.28 40.48 #interaction term
11c more numbers … … #interaction term
12a … … … … #interaction term
12b … … … … #interaction term
12c … … … … #interaction term
13a … … … … #interaction term
13b … … … … #interaction term
13c … … … … #interaction term
Код, который я использовал, приведен ниже. Не работает, ошибки в первой строке.
newdata <- data.frame(int1 = levels(df$int1), int2 = levels(df$int2),
x = range(df$x)) #wrong. How change it?
fitmod = fitted(mod, newdata = newdata, re_formula = NA, summary =
T)*100 #convert to %
colnames(fitmod) = c('fit', 'se', 'lwr', 'upr')
foo = cbind(newdata, fitmod)
Создание дополнительного столбца, объединяющего int1 и int2, также не работает. Если бы модель имела только один предиктор, который был бы фактором, это было бы:
newdata <- data.frame(int1 = levels(df$int1))
Как правильно настроить этот df, чтобы получить все правильные коэффициенты? Большое спасибо
Самый простой в использовании modelr
:
newdat <- modelr::data_grid(df, .model = mod)
newdat$pr <- predict(mod, newdat)
Но похоже, что вместо этого вы используете emmeans
?
emmeans::emmeans(mod, ~int1:int2)
Я попробовал ваше предложение с modelr
, но newdat
это data frame with 0 columns and 0 rows
Затем вам нужно вручную указать нужные вам столбцы в data_grid
Но emmeans, казалось, работал отлично! Спасибо еще раз! Однако один вопрос: учитывая, что это байесовский анализ, и мне может понадобиться извлечь цепочки MCMC аналогичным образом, нет способа сделать это с помощью emmeans, не так ли? В моем коде выше все, что нужно изменить, это summary = F
(из (summary = T)*100
), по-видимому, и цепочки извлекаются вместо сводки.
Только что увидела ваш ответ, отлично, тоже попробую.
Спасибо большое, попробую! На самом деле я использую пакет
brms
для байесовской регрессии, но я думаю, что сlme4
все точно так же.