I am trying to graph the results of a week by condition interaction in a mixed-effects model with multiple covariates. Data is nested within-subject. I am predicting steps by week (factor with 5 levels), condition (2 levels), age, sex, education, health, days using the app (which was the intervention). The model, fitted using lme4 package, is pasted below.
m1<-lmer(Steps~ as.factor(Week) + Condition + Age + Gender + Edu + Health + DaysUse + Condition*Week + (1|SubjectID), Dataset1, REML=FALSE)
summary(m1)
Originally to graph this, I had aggregated steps by week and condition, to get the mean and standard error of steps at each week, by condition. Then I graphed this in ggplot using the following syntax.
ggplot(data=stepgraph, aes(x=Week,y=Step.Mean, fill=Condition, group=Condition)) +
geom_point(stat = "identity", aes(shape = Condition), size = 3)+
geom_errorbar(aes(ymin=Step.Mean-Step.SE, ymax=Step.Mean+Step.SE, linetype = Condition), width=.2)+
geom_line(aes(linetype = Condition))+
theme_bw()+
coord_cartesian(ylim = c(1000, 10000))+
labs(title ="Average weekly steps by Condition", x = "Study Week", y = "Weekly Step Average") +
theme(plot.title = element_text(hjust = 0.5))+
theme(
legend.position = c(.95, .05),
legend.justification = c("right", "bottom"),
legend.box.just = "right",
legend.margin = margin(6, 6, 6, 6)
)
However, my advisor noted that this would not technically be a graph of my model results. Not only does this not account for the covariates in my model, I don't think it takes into account the fact that the data is nested within-person.
I also thought about adjusting the y axis using 'position_nudge', creating the nudge by multiplying the mean of each of my covariates by the beta of these covariates, e.g.,
AgeMB1<-m1@beta[6]*mean(Dataset1$Age, na.rm=T)
and so on for all my other covariates... then I summed these values and nudged everything on the y axis accordingly.
However, this just moves all the points up or down by a set amount. I am not sure that this is the correct way to do it either.
Finally, based on another suggestion I found online, I tried constraining the intercept to 0, to output model estimates of mean steps at each week. I did this using the following model... however, to do this separated by condition, I needed to subset the data by condition and run the model twice.
m1C<-lmer(Steps~ 0 + as.factor(Week) + Age + Gender + Edu + Health + DaysUse + (1|SubjectID), Dataset1.1, REML=FALSE)
summary(m1C)
When I take the betas and standard errors of the estimates here, the standard errors are HUGE. I am wondering if this is because the standard error is estimated between-person, not within.
I guess my main question, is whether there is a best-practice for graphing an interaction like this, with data that's nested within person, AND accounting for the covariates in my model. Or, is there a way to extract the predicted step means and standard errors directly from my model, at each of the different weeks by condition?
Thank you in advance for any assistance.