У меня есть модель, которая комбинирует фиктивную и непрерывную переменную для описания результата после нарушения. Так что, если было нарушение, у меня есть измерения времени в моменты 1:16 после нарушения. Если в недавнем прошлом не было никаких нарушений, результат закодирован до поддельного значения времени -1. Вот представление набора данных:
library(lme4)
library(ggplot2)
df <- data.frame(ID = rep(c("a", "b", "c"), each = 20),
Time = c(1:16, -1, -1, -1, -1,
1:16, -1, -1, -1, -1,
1:16, -1, -1, -1, -1))
df$y <- 2 + 0.8*df$Time + 1*df$Time^2 + rnorm(30, 0, 3)
df[df$Time < 0,]$y <- rnorm(12, 5, 3)
df[df$ID == "b",]$y <- df[df$ID == "b",]$y + 5
df[df$ID == "c",]$y <- df[df$ID == "c",]$y - 5
df$Exposure <- "Before"
df[df$Time > 0,]$Exposure <- "After"
df$Exposure <- factor(df$Exposure, levels = c("Before", "After"))
ggplot(df[df$Time > 0,]) +
geom_point(aes(x = Time, y = y, colour = ID)) +
geom_point(data = df[df$Time < 0,], aes(x = -5, y = y, colour = ID))
То, что я ищу, сравнивает оценку «без помех» с различными периодами времени после нарушения, чтобы увидеть, когда разница станет существенной.
Перед моделированием назначьте данные "без помех" на время 0.
df[df$Time < 0,]$Time <- 0
m <- lmer(y ~ Exposure + poly(Time, 2) + (1|ID), data = df)
# output estimates
newdata <- data.frame(Exposure = c("Before", "After", "After", "After", "After", "After"),
Time = c(0, 1, 4, 8, 12, 16))
newdata$Pred <- predict(m, re.form = NA, newdata = newdata)
## plot looks good
ggplot(df[df$Time > 0,]) +
geom_point(aes(x = Time, y = y, colour = ID)) +
geom_point(data = df[df$Time == 0,], aes(x = -5, y = y, colour = ID)) +
geom_line(data = newdata[newdata$Exposure == "After",],
aes(x = Time, y = Pred)) +
geom_point(data = newdata[newdata$Exposure == "Before",],
aes(x = -5, y = Pred), colour = "red")
Как бы я сравнил, скажем, До оценки с оценками После, например, Time==3, Time == 6 и Time == 9? Примерно так было бы замечательно, но я не могу понять, как устранить ошибку, которую я получаю.
library(contrast)
library(multcomp)
cc <- contrast(m,
a = list(Time = 0, Exposure = "Before"),
b = list(Time = c(3, 6, 9), Exposure = "After"))
summary(glht(m, linfct = cc$X))
### ОБНОВЛЕНИЕ
После отличных изменений в rvl я провел пробную проверку своих фактических данных и столкнулся с новой проблемой. Моя фактическая переменная времени не является целым числом, но я хочу делать прогнозы в целочисленном масштабе. Когда я обновляю игрушечный пример, вложение, кажется, нарушается:
df$Time <- df$Time + rnorm(60, 0, 0.5)
df[df$Exposure == "Before",]$Time <- -1.12
m <- lmer(y ~ Exposure + poly(Time, 2) + (1|ID), data = df)
# freshly installed emmeans from github
emm = emmeans(m, "Time", at = list(Time = c(0,3,6,9)))
emm ## no longer get the nesting info, and the preds aren't nested
По моим собственным данным (и используя спецификацию at, я на самом деле получаю только одну строку для Time == 0 и Exposure == Before, и все - ничего больше в выводе ... какие-либо предложения ??
## ОБНОВЛЕНИЕ2
По какой-то причине решение работает с игрушечным примером, но не с моими собственными данными ... Вот небольшое подмножество моего набора данных. Подгонка модели единственная, но проблемы, которые я получаю для emmeans, такие же, как и для всего моего набора данных ... help?
df <- structure(list(ID = structure(c(2L, 1L, 2L, 2L, 1L, 2L, 2L, 1L,
2L, 1L, 2L, 1L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 1L, 1L, 1L, 1L, 2L,
1L, 1L, 2L, 1L, 1L, 2L, 1L, 2L, 1L, 2L, 2L, 2L, 2L, 2L, 2L, 2L,
2L, 2L, 1L, 2L, 2L), .Label = c("B", "A"), class = "factor"),
Exposure = structure(c(1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L,
1L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L,
2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L,
2L, 2L, 2L, 2L, 2L, 2L), .Label = c("No exposure", "Exposure"
), class = "factor"), Time = c(0, 0, 0, 0, 0, 0, 0, 0, 0,
0, 4.78757545912946, 9.63531173739354, 5.47889766247861,
7.17017886302881, 1.43155423003375, 3.72391354120779, 2.56353688399906,
8.29779117320654, 9.52304006615339, 9.48174174807695, 0.859601950498583,
4.63141168677387, 7.92347302279951, 7.92067346608815, 5.23250024053785,
5.57671787587839, 1.85126003367584, 3.1097216702916, 7.72389534567839,
9.36144591805227, 2.70213603445334, 1.84811002303022, 6.82448971585652,
7.88336338096561, 3.84031339520175, 5.62874085650497, 4.0972590990481,
2.09535527965164, 2.22160757456982, 7.35862943664427, 7.41826702411403,
8.24309337727667, 4.7943847267765, 5.8840472004994, 7.02963322046381
), Response = c(-7.16922413711838, 143.482571506177, 16.45347120693,
25.022565770909, -55.8024015971315, -124.925019624537, -16.4000310854958,
40.9499232825204, 2.46651714407957, -34.3558611547229, -80.1711009500979,
-58.5220697399603, 17.6390452197579, -11.2077688506688, 87.0618648836916,
113.611468732, -27.1400972587652, -30.0256851366867, -111.149731873181,
-24.2689502403869, -16.2737794106996, -125.618994529607,
95.9640135688539, 46.4163972081548, 6.72470222784859, -0.148508667228167,
-118.897875455802, 28.6093848128793, -57.5632050845714, 31.390260468939,
27.6826377837027, -40.7112943346364, -53.5934755706868, 27.0754421268185,
165.146183257597, 39.6762439690417, -9.74912218853661, 18.3454700992841,
33.8006770750647, -18.6013173700368, 12.7360264627221, 178.646948999019,
93.5496871933183, -8.68468960982507, 2.86668462850576)), row.names = c(1L,
3L, 5L, 7L, 9L, 11L, 13L, 15L, 17L, 19L, 21L, 23L, 25L, 27L,
29L, 31L, 33L, 35L, 37L, 39L, 41L, 43L, 45L, 47L, 49L, 51L, 53L,
55L, 57L, 59L, 61L, 63L, 65L, 67L, 69L, 71L, 73L, 75L, 77L, 79L,
81L, 83L, 85L, 87L, 89L), class = c("tbl_df", "tbl", "data.frame"))
Запуск модели и emmeans:
m <- lmer(Response ~ Exposure + poly(Time, 2) + (1|ID), data = df)
## this only gives one row instead of 8?
emmeans(m, c("Time", "Exposure"), at = list(Time = c(0,3,6,9)))
## when I specify the nesting myself, I get a "multiple actual arguments" error...
emmeans(m, c("Time", "Exposure"), at = list(Time = c(0,3,6,9)),
nesting = "Time %in% Exposure")