I have an unbalanced repeated measures data set to analyse, and I've read that the way most statistical packages handle this with ANOVA (i.e. type III sum of squares) is wrong. Therefore, I would like to use a mixed effects model to analyse these data. I have read a lot about mixed models in
R, but I am still very new to
R and mixed effect models and not very confident I am doing things right. Note that I can't yet entirely divorce myself of "traditional" methods, and still need $p$-values and post hoc tests.
I would like to know if the following approach makes sense, or if I am doing something horribly wrong. Here's my code:
# load packages library(lme4) library(languageR) library(LMERConvenienceFunctions) library(coda) library(pbkrtest) # import data my.data <- read.csv("data.csv") # create separate data frames for each DV & remove NAs region.data <- na.omit(data.frame(time=my.data$time, subject=my.data$subject, dv=my.data$dv1)) # output summary of data data.summary <- summary(region.data) # fit model # "time" is a factor with three levels ("t1", "t2", "t3") region.lmer <- lmer(dv ~ time + (1|subject), data=region.data) # check model assumptions mcp.fnc(region.lmer) # remove outliers (over 2.5 standard deviations) rm.outliers <- romr.fnc(region.lmer, region.data, trim=2.5) region.data <- rm.outliers$data region.lmer <- update(region.lmer) # re-check model assumptions mcp.fnc(region.lmer) # compare model to null model region.lmer.null <- lmer(dv ~ 1 + (1|subject), data=region.data) region.krtest <- KRmodcomp(region.lmer, region.lmer.null) # output lmer summary region.lmer.summary <- summary(region.lmer) # run post hoc tests t1.pvals <- pvals.fnc(region.lmer, ndigits=10, withMCMC=TRUE) region.lmer <- lmer(dv ~ relevel(time,ref="t2") + (1|subject), data=region.data) t2.pvals <- pvals.fnc(region.lmer, ndigits=10, withMCMC=TRUE) region.lmer <- lmer(dv ~ relevel(time,ref="t3") + (1|subject), data=region.data) t3.pvals <- pvals.fnc(region.lmer, ndigits=10, withMCMC=TRUE) # Get mcmc mean and 50/95% HPD confidence intervals for graphs # repeated three times and stored in a matrix (not shown here for brevity) as.numeric(t1.pvals$fixed$MCMCmean) as.numeric(t1.pvals$fixed$HPD95lower) as.numeric(t1.pvals$fixed$HPD95upper) HPDinterval(as.mcmc(t1.pvals$mcmc),prob=0.5) HPDinterval(as.mcmc(t1.pvals$mcmc),prob=0.5)
Some specific questions I have:
- Is this a valid way of analysing mixed effects models? If not, what should I be doing instead.
- Are the criticism plots output by mcp.fnc good enough for verifying model assumptions, or should I be taking additional steps.
- I get that for mixed models to be valid, the data need respect assumptions of normality and homoscedasticity. How to I judge what is "approximately normal" and what is not by looking at the criticism plots generated by mcp.fnc? Do I just need to get a feel for this, or is their a prescribed way of doing things? How robust are mixed models in respect to these assumptions?
- I need to assess differences between the three time points for ~20 characteristics (biomarkers) of the subjects in my sample. Is fitting and testing separate models for each acceptable so long as I report all undertaken tests (significant or not), or do I need any form of correction for multiple comparisons.
To be a little more precise in regards to the experiment, here are some more details. We followed a number of participants longitudinally as they underwent a treatment. We measured a number of biomarkers before the start of the treatment and at two time points after. What I'd like to see is if there are difference in these biomarkers between the three time points.
I am basing most of what I am doing here on this tutorial, but made some changes based on my needs and things I read. The changes I made are:
- relevel the "time" factor to obtain t1-t2, t2-t3, and t1-t3 comparisons with pvals.fnc (from the languageR package)
- compare my mixed model to the null model using an approximate F-test based on a Kenward-Roger's approach (using the pbkrtest package) rather than a likelihood ratio test (because I read, that Kenward-Roger's is better regarded right now)
- Use the LMERConvenienceFunctions package to check assumptions and remove outliers (because I read that mixed models are very sensitive to outliers)