Linear mixed-effects models aren't commonly used in my corner of biology, and I need to report the statistical test I used in a paper I'm trying to write. I know that awareness of multilevel modeling is starting to appear in some areas of the biosciences (A solution to dependency: using multilevel analysis to accommodate nested data), but I'm still trying to learn how to report my results!
My experimental design, in brief:
*Subjects were assigned to one of four treatment groups
*Measurements of the dependent variable were taken on various days after the start of treatment
*The design is unbalanced (unequal numbers of subjects in the treatment groups, and missing measurements for some subjects on some days)
*Treatment A is the reference category
*I centered the data at the final day of treatment
I want to know whether Treatment A (the reference category) yields significantly better results than the other treatments (at the end of treatment).
I did my analysis in R, using nlme:
mymodel <- lme(dv ~ Treatment*Day, random = ~1|Subject, data = mydf, na.action = na.omit, + correlation = corAR1(form = ~1 |Subject), method = "REML")
And the output (in part; truncated for brevity) is:
>anova(mymodel) numDF denDF F-value p-value (Intercept) 1 222 36173.09 <.0001 Treat 3 35 16.61 <.0001 Day 7 222 3.43 0.0016 Treat:Day 21 222 3.62 <.0001 >summary(mymodel) Fixed effects: dv ~ Treatment * Day Value Std.Error DF t-value p-value (Intercept) 7.038028 0.1245901 222 56.48945 0.0000 TreatmentB 0.440560 0.1608452 35 2.73903 0.0096 TreatmentC 0.510214 0.1761970 35 2.89570 0.0065 TreatmentD 0.106202 0.1637436 35 0.64859 0.5208
So, I know that the effect of Day differs by Treatment, and that, on the final day of treatment (where the data are centered), dv is significantly different in Treatment A than Treatments B or C.
What I want to say is: "As predicted, we found that Dependent Variable was significantly lower in subjects receiving Treatment A (mean +/- SE) than in subjects receiving Treatment B (mean +/- SE, p=0.0096) or Treatment C (mean +/- SE, p=0.0065), as measured on the final day of treatment."
But, I have to indicate what statistical test was done. Would this be an acceptable way to describe the analysis? "[Measurement Method] was performed on days indicated and Dependent Variable (units) was determined; we analyzed the log-transformed data using a linear mixed-effects model centered at [final day of treatment]. Symbols represent mean dv; error bars are standard error. On the final day of treatment, dv was significantly lower in Treatment A (mean +/- SE) than Treatment B (mean +/- SE, p=0.0096)..."
*Does that say enough about the statistical test used? (Readers are used to seeing something more like "mean +/- SE, p=0.0096, Student's t-test," but it seems weird to write "p=0.0096, coefficient for Treatment B vs. Treatment A from linear mixed-effects model at [final day of treatment].")
*Is there a better way to put this?
(The methods section will include more information about the stats: "[Measurement Method] data were analyzed using R and the R packages... We analyzed the log-transformed Dependent Variable data by using linear mixed effects models using Subjects as random effects and an autocorrelation structure of order 1 (AR1). As fixed effects, we included Treatment and Day, and the interaction of Treatment and Day. We checked for normality and homogeneity by visual inspections of plots of residuals against fitted values. To assess the validity of the mixed effects analyses, we performed likelihood ratio tests comparing the models with fixed effects to the null models with only the random effects.")
Any advice on how to report results of a linear mixed effects model for an often stats-averse audience (and written by a relative stats novice) would be greatly appreciated!