I would like to ask for some help with plotting prediction values from my model as well as the equations generated by the estimation of the lmer().
So, the data that I have is the mass volume of different rats across different days. Each rat has different time points where they took the measurement of that volume.
So, then the model that I use is :
m1 <- lmer(lVolume ~ Country*Day + (1|Rat))
I do this because I am interested in exp(fitted) values and then obtaining an exponential approach for this model instead of using a nonlinear mixed effect model (for the moment)
To plot the predictions from this model, my attempt was:
m1%>%
augment() %>%
clean_names() %>%
ggplot(data = .,
mapping = aes(x = day,
y = exp(l_volume),
group = rat)) +
geom_point(alpha = 0.5) +
geom_line(alpha = 0.5) +
geom_point(aes(y = exp(fitted)),
color = "red") +
geom_line(aes(y = exp(fitted)),
color = "red") +
expand_limits(x = 0 , y = 0)
Here I plotted more rats but, as you can see the (0,0) is too far from the predictions of the lmer. I was wondering how I plot the prediction that my model is generating to see points from (0,200). I tried a hint from here by creating a new data frame and then plot using again predict(m1, newdata = new_df), but I am clueless how to create this data frame since I have 20 rats and I do not know how to expand this to the predict().
My attempt:
pframe <- data.frame(Day=seq(0, 200, length.out=101))
pframe$continuous_outcome <- predict(m1, newdata = pframe, level = 0)
ggplot(data, aes(Day,lVolume)) +
geom_point() +
geom_line(data=pframe)
but I got an error:
Error in eval(predvars, data, env) : object 'Rat' not found
And, also is there a way to plot also the equations that you generate from each estimation, i.e, from each rat you have a set of estimators fixed and random, how can I plot the equation (red curves) that the lmer is generating for each of the rats?
It turned out to be easier to use predict than broom.mixed::augment.
(all combinations of Rat/Country/Days 0-150 (taking day up to 200 led to some extreme predictions that blew the vertical scale)
library(tidyverse)
dc <- distinct(dplyr::select(dat1, Rat, Country))
pframe <- (with(dat1,
expand_grid(Rat = unique(Rat),
Day = 0:150))
%>% full_join(dc, by = "Rat")
%>% mutate(lVolume = predict(m1, newdata = .))
)
Combine data and predictions into a single data frame (you don't have to do this but it makes the legend easy)
comb <- dplyr::bind_rows(list(data = dat1, model = pframe),
.id = "type")
Plot:
ggplot(comb, aes(Day, exp(lVolume), colour = type)) +
geom_point(alpha = 0.2) +
geom_line(aes(group = interaction(type, Rat))) +
scale_colour_manual(values = c("black", "red"))
Reconstructing data:
dat0 <- list(
list("rat1", vol=c(78,304,352,690,952,1250), days = c(89,110,117,124,131,138), country = "Chile"),
list("rat2", vol=c(202,440,520,870,1380), days = c(75,89,96,103,110), country = "Chile"),
list("rat3", vol=c(186,370,620,850,1150), days = c(75,89,96,103,110), country = "Chile"),
list("rat4", vol=c(92,250,430,450,510,850,1000,1200), days = c(47,61,75,82,89,97,103,110), country = "England"),
list("rat5", vol=c(110,510,710,1200), days = c(47,61,75,82), country = "England"),
list("rat6", vol=c(115,380,480,540,560,850,1150,1350), days = c(47,61,75,82,89,97,103,110), country = "England"))
dat1 <- purrr::map_dfr(dat0,
~ data.frame(Rat = .[[1]],
lVolume = log(.$vol), Day = .$days,
Country = .$country))
m1 <- lmer(lVolume ~ Country*Day + (1|Rat), data = dat1)