12
$\begingroup$

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)..."

Specifically,
*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!

$\endgroup$
  • $\begingroup$ can you explain more regarding what you mean by 'centered' on the last day of the study? why not just pump in day as 0 through number of days on the last day of the study? Also- in your summary table- where are the interaction effects? $\endgroup$ – colin Jun 24 '14 at 20:46
1
$\begingroup$

This may not help answer your question, but I noticed that you have a repeated measure (Day) in your experiment, but you did not indicate that this was a repeated measure in your model. I would have thought the random term in your model to be as such:

mymodel <- lme(dv ~ Treatment*Day, random = ~1|Subject/Day, 
               data = mydf, na.action = na.omit,
               correlation = corAR1(form = ~1 |Subject/Day), method = "REML")

As for reporting the results, did you intend to report on the day on which you start seeing significant differences between the treatments? If so, then I think you'll need to look at/report on the contrasts on the interaction term as well. I'm a stats novice myself and basically have the same question as you do :-)

Andy Field's "Discovering Statistics Using R" explains how to report results from a linear mixed effects model in Ch14. I don't have the book at hand but can edit this post once I get my hands on it again.

$\endgroup$
  • $\begingroup$ that depends, if the Day is the same for all Subjects, I guess it is a crossed rather than a nested random factor, right? Then it's $random=~Day|Subject$ because day is continuous it is placed in front of the vertical line. Comments appreciated! $\endgroup$ – Aura Borealis Feb 25 '16 at 16:08

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.