Guest User

Untitled

a guest
Nov 18th, 2018
98
0
Never
Not a member of Pastebin yet? Sign Up, it unlocks many cool features!
text 1.51 KB | None | 0 0
  1. library(lme4)
  2. library(ggplot2)
  3.  
  4. df <- data.frame(ID = rep(c("a", "b", "c"), each = 20),
  5. Time = c(1:16, -1, -1, -1, -1,
  6. 1:16, -1, -1, -1, -1,
  7. 1:16, -1, -1, -1, -1))
  8. df$y <- 2 + 0.8*df$Time + 1*df$Time^2 + rnorm(30, 0, 3)
  9. df[df$Time < 0,]$y <- rnorm(12, 5, 3)
  10.  
  11. df[df$ID == "b",]$y <- df[df$ID == "b",]$y + 5
  12. df[df$ID == "c",]$y <- df[df$ID == "c",]$y - 5
  13. df$Exposure <- "Before"
  14. df[df$Time > 0,]$Exposure <- "After"
  15. df$Exposure <- factor(df$Exposure, levels = c("Before", "After"))
  16.  
  17. ggplot(df[df$Time > 0,]) +
  18. geom_point(aes(x = Time, y = y, colour = ID)) +
  19. geom_point(data = df[df$Time < 0,], aes(x = -5, y = y, colour = ID))
  20.  
  21. df[df$Time < 0,]$Time <- 0
  22. m <- lmer(y ~ Exposure + poly(Time, 2) + (1|ID), data = df)
  23.  
  24. # output estimates
  25. newdata <- data.frame(Exposure = c("Before", "After", "After", "After", "After", "After"),
  26. Time = c(0, 1, 4, 8, 12, 16))
  27. newdata$Pred <- predict(m, re.form = NA, newdata = newdata)
  28.  
  29. ## plot looks good
  30. ggplot(df[df$Time > 0,]) +
  31. geom_point(aes(x = Time, y = y, colour = ID)) +
  32. geom_point(data = df[df$Time == 0,], aes(x = -5, y = y, colour = ID)) +
  33. geom_line(data = newdata[newdata$Exposure == "After",],
  34. aes(x = Time, y = Pred)) +
  35. geom_point(data = newdata[newdata$Exposure == "Before",],
  36. aes(x = -5, y = Pred), colour = "red")
  37.  
  38. library(contrast)
  39. library(multcomp)
  40.  
  41. cc <- contrast(m,
  42. a = list(Time = 0, Exposure = "Before"),
  43. b = list(Time = c(3, 6, 9), Exposure = "After"))
  44. summary(glht(m, linfct = cc$X))
Add Comment
Please, Sign In to add comment