Not a member of Pastebin yet?
Sign Up,
it unlocks many cool features!
- 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))
- 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")
- 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))
Add Comment
Please, Sign In to add comment