0

I have repeated measurements on a continuous outcome that looks as follows:

library(magrittr)
library(ggplot2)
library(nlme)
mydata <- structure(list(ID = structure(c(1L, 1L, 1L, 1L, 1L, 1L, 2L, 2L, 
3L, 3L, 3L, 3L, 3L, 4L, 4L, 5L, 5L, 5L, 6L, 6L, 6L, 6L, 6L, 7L, 
7L, 7L, 7L, 7L, 8L, 8L, 9L, 9L, 10L, 10L, 11L, 11L, 11L, 11L, 
12L, 12L, 13L, 13L, 14L, 14L, 14L, 14L, 14L, 15L, 15L, 16L, 16L, 
17L, 17L, 17L, 17L, 17L, 18L, 18L, 19L, 19L, 20L, 20L, 21L, 21L, 
22L, 22L, 22L, 22L, 22L, 23L, 23L, 24L, 24L, 24L, 24L), .Label = c("2", 
"3", "4", "7", "8", "13", "14", "20", "21", "22", "24", "25", 
"27", "29", "30", "31", "34", "36", "37", "38", "39", "40", "48", 
"49", "50", "51", "52", "54", "58", "60", "61", "65", "74", "75", 
"76", "77", "80", "81", "82", "83", "84", "86", "87", "88", "92", 
"94", "95", "96", "103", "104", "105", "114", "115", "116", "117", 
"119", "125", "126", "127", "132", "134", "135", "137", "138", 
"141", "142", "145", "152", "153", "154", "157", "159", "160", 
"162", "164", "165", "171", "172", "179", "180", "184", "185", 
"189", "194", "195", "197", "198", "202", "203", "205", "209", 
"213", "221", "253", "255", "258", "262", "271", "273", "277", 
"279", "310", "315", "320"), class = "factor"), date_ct = structure(c(15923, 
16122, 16715, 16902, 17086, 18003, 16150, 16841, 16421, 16764, 
16951, 17135, 18011, 16622, 18247, 16582, 16752, 18045, 16729, 
16862, 17042, 17226, 18102, 16568, 16736, 16916, 17100, 18040, 
16743, 16841, 16589, 16729, 16526, 16729, 16619, 16862, 17042, 
17226, 16407, 18437, 16512, 16953, 16457, 16946, 17112, 17310, 
17989, 16573, 16841, 15923, 16752, 16505, 16729, 16909, 17107, 
18038, 16540, 16743, 15951, 16122, 16624, 18202, 16623, 18221, 
16694, 16715, 16902, 17086, 18037, 16451, 16743, 16421, 16736, 
16909, 17100), class = "Date"), age = c(56.6, 57.1, 58.8, 59.3, 
59.8, 62.3, 43.2, 45.1, 52, 52.9, 53.4, 53.9, 56.3, 58.5, 63, 
57.4, 57.9, 61.4, 57.8, 58.2, 58.7, 59.2, 61.6, 52.4, 52.8, 53.3, 
53.8, 56.4, 70.8, 71.1, 61.4, 61.8, 59.2, 59.8, 61.5, 62.2, 62.7, 
63.2, 48.9, 54.5, 54.2, 55.4, 50.1, 51.4, 51.8, 52.4, 54.3, 55.4, 
56.1, 48.6, 50.9, 64.2, 64.8, 65.3, 65.8, 68.4, 68.3, 68.8, 66.7, 
67.1, 60.5, 64.8, 56.5, 60.9, 62.7, 62.8, 63.3, 63.8, 66.4, 49, 
49.8, 61, 61.8, 62.3, 62.8), continuous_outcome = c(1636.4, 544.1, 
1408, 1594.7, 1719.4, 2345.9, 115.3, 226, 2678.2, 3451.6, 3702.7, 
3632.7, 5805, 155.2, 1095, 992.2, 296.6, 2020.4, 3708.6, 2710.7, 
2934.2, 3080.4, 4489.7, 3459.4, 4965.3, 5553.1, 5037.8, 7315.7, 
29980.8, 35407.5, 2263.2, 2060.6, 3220.7, 4467.1, 5902.3, 6407.2, 
5947.1, 6271.6, 306, 689.3, 1430.6, 1672.1, 9.9, 58.7, 69.9, 
125.3, 39.5, 3842.5, 5136.3, 216.6, 332.4, 5719.3, 5386, 5490.7, 
5268.2, 6166.7, 12520.6, 12981.8, 2896.1, 2976.8, 5495.6, 6470.6, 
4235.5, 7603.5, 3887, 3344.5, 2885.7, 3324.1, 6401, 1942.2, 2000.9, 
2401.7, 2231.5, 2749.7, 2741.7)), row.names = c(NA, -75L), class = c("tbl_df", 
"tbl", "data.frame"))

Through exploratory analyses I've found that a spline with 2 knots run the following linear mixed model to determine if age is a predictor of my outcome:

model1 <- 
   lme(data=mydata, 
      fixed=continuous_outcome ~ age, 
      random=~1|ID)

I want to find a way to plot the predicted values of this model over the observed values for my outcome. Any help?

Ben Bolker
  • 211,554
  • 25
  • 370
  • 453
tcvdb1992
  • 413
  • 3
  • 12
  • 1
    You should be able to use `predict`. For example, to get the plot: `plot(mydata$continuous_outcome, predict(model1))` – dario Sep 28 '21 at 13:45

1 Answers1

4

predict() is the basic answer. There are more sophisticated approaches using the ggeffects package (but it doesn't seem to work with lme models?) or sjPlot (but it's not super-easy to overlay the data?)

model1 <-
   lme(data=mydata,
      fixed=continuous_outcome ~ splines::ns(age, 2),
      random=~1|ID)

Construct prediction data frame: this step isn't always necessary, but it's useful if your data points are non-unique, or out of order, or if you want to plot a nonlinear response, or confidence intervals [which are generally curved], or if you have a very large data set (in which case predictions for the whole data set are probably unnecessary); it's also usually necessary if you have a more complex model with non-focal predictors (the effects/ggeffects packages are designed to automate this step).

pframe <- data.frame(age=seq(min(mydata$age), max(mydata$age), length.out=101))

Use newdata to specify the values to make predictions for, and level=0 to predict at the population level (this also means we don't need to give ID values in newlevels)

pframe$continuous_outcome <- predict(model1, newdata = pframe, level = 0)
ggplot(mydata, aes(age, continuous_outcome)) + 
  geom_point() +
  geom_line(data=pframe)

enter image description here

Ben Bolker
  • 211,554
  • 25
  • 370
  • 453
  • Hi @ Ben Bolker. Thanks for your answer. I am confused about two lines of your code: 1. What exactly does `pframe <- data.frame(age=seq(min(mydata$age), max(mydata$age), length.out=101))` do? Make a new dataframe based upon only the minimum and maximum values of age? Can't I just add the predicted values to my original dataset `mydata`? 2. What do the newdata and level arguments do in `pframe$continuous_outcome <- predict(model1, newdata = pframe, level = 0)`? Thanks again! – tcvdb1992 Sep 28 '21 at 14:51
  • I think I've answered your questions in my edits. – Ben Bolker Sep 28 '21 at 16:17
  • Thanks so much for that explanation. A few things are still not completely clear to me: so what I do is I make a new dataframe with random numbers between the minimum and maximum age of my original dataset (with `seq`) and then I add the predicted values to that data. What if I also have for example the variable sex. Do I add a random sequence of the same sex factor in my model to the prediction dataset `pframe`, for example by adding `sex=sample(mydata$sex, size=101`?. A second question is, what does the `level=0` in the `predict` function do? – tcvdb1992 Oct 02 '21 at 15:17
  • (1) `seq()` doesn't use random numbers (although you could say "arbitrary"); (2) If you have non-focal variables, everything gets more complicated. You could draw two separate lines; you could average the values for male and female; etc.. Addressing these complexities is what packages like `emmeans`, `sjPlot`, `ggeffects`, etc. are for: you might not want to reinvent it all. (3) I already answered this - it's to predict at the population level ("average" group) rather than at the level of individual groups – Ben Bolker Oct 02 '21 at 19:08
  • Hi @BenBolker , I was wondering if you could give me some ideas to construct the new data frame? I tried the hints here but not successfully hehe Could you have a look, please ? https://stackoverflow.com/questions/71210932/broom-mixed-exp-model-predictions – Rosa Maria Feb 21 '22 at 18:33