# How to state my model in R in words (linear mixed effect model, general linear hypothesis, Tukey, Benjamini Hochberg)

I have a dataset that has four treatment groups (AX, AY, BX, BY). I am trying to determine if there is a significant difference in a response variable (Length) between any of these four treatment groups. I am also considering a random variable called Experiment. All this information is in the R object d below:

set.seed(1)
d = data.frame(Length = c(rnorm(9,0), rnorm(9,-0.5)), Drug = c(rep("A",9), rep("B",9)), Location = rep(c(rep("X",3), rep("Y",3)),3), Treatment = paste0(c(rep("A",9), rep("B",9)), rep(c(rep("X",3), rep("Y",3)),3)), Experiment = sample(c(1,2), 18, replace = TRUE))

# Linear mixed effect model
lmeOut = lme(Length ~ Treatment, data=d, random = ~1|Experiment)
# General linear hypothesis and multiple comparison


I am unsure how to describe my code above. There is a linear mixed effect model (lme), general linear hypothesis and multiple comparison (glht), and both Tukey and Benjamini-Hochburg used (BH).

Would it be fair for me to say:

"A linear mixed effect model was used with the Experiment
treated as a random variable. A pairwise combination was
performed using a general linear model with
Benjamini-Hochburg correction".


I am planning to have letters over a bar chart for groups AX, AY, BX, BY. As can be seen by running my code above, none of the six pairwise combinations reach significance (adjusted p-value < 0.05). So, I would just have letter "a" over each bar. Should I call this Tukey HSD letters? If not, what should I call these letters?

Very confused about the wording here and any advice would be much appreciated.

• How many experiments do you have? If just 2, it seems a bit over-kill to treat Experiment as a random grouping factor. Also, what is the rationale for adjusting the Tukey p-values (which are already adjusted for multiplicity) further by BH? – Isabella Ghement Sep 15 '18 at 23:51
• Thanks, @IsabellaGhement. I am mostly trying to use "lme" function to do a pairwise combination between the 4 treatment groups (AX, AY, BX, BY), which should result in 6 p-values. Then, I want to use BH correct between these six p-values using p.adjust(pvals, "BH"). However, I am unable to get the six p-values out of the "lme" function, so I am using glht to do so. Is there a way to do this more directly without involving Tukey? – luckButtered Sep 16 '18 at 16:42
• Thank you! Please see my response below. The issue of possibly having too few experiments for fitting an lme model remains. Also, Length tends to be a skewed variable in many cases so make sure that's not an issue in your case, or else you may need to log-transform the data perhaps. – Isabella Ghement Sep 16 '18 at 17:41

One roundabout way to get the unadjusted p-values for all the pairwise comparisons would be to fit the model several times. Each time you fit the model, you set the reference level of Treatment to something else.
Something like this:

d$Treatment <- factor(d$Treatment)


First model fit:

d$Treatment <- relevel(d$Treatment, ref = "AX")
lme1 <- lme(Length ~ Treatment, data =d, random = ~ 1|Experiment)
summary(lme1)


The summary output for the first model fit will give you the unadjusted p-values for the treatment comparisons AY vs AX, BX vs AX, BY vs AX.

Second model fit:

d$Treatment <- relevel(d$Treatment, ref = "AY")
lme2 <- lme(Length ~ Treatment, data =d, random = ~ 1|Experiment)
summary(lme2)


The summary output for the second model fit will give you the unadjusted p-values for the treatment comparisons BX vs AY and BY vs AY.

Third model fit:

d$Treatment <- relevel(d$Treatment, ref = "BX")
lme3 <- lme(Length ~ Treatment, data =d, random = ~ 1|Experiment)
summary(lme3)


The summary output for the third model fit will give you the unadjusted p-values for the treatment comparison BY vs BX.

Once you have the 6 unadjusted p-values produced by the above, you can store them into a vector p:

p <- c(0.002, 0.124, 0.317, 0.276, 0.439, 0.721)


(Note that I made up the above p-values for illustration purposes - you would have to insert the actual p-values derived from your data in the above code line.)

p.adjust(p, "BH")


• A linear mixed effects model was used to relate the continuous outcome variable Length to the categorical explanatory variable Treatment, with the latter variable having the categories AX, AY, BX and BY;

• The model included a random (intercept) effect for the grouping variable Experiment;

• The model was fitted to the data by REML (or ML);

• Post-hoc pairwise comparisons of the 4 levels of the explanatory Treatment variable with respect to mean Length were performed (for the typical experiment) after fitting the model to the data and Benjamini-Hochberg adjusted p-values were reported for the 6 resulting comparisons;

• The post-hoc pairwise comparisons found no statistically significant differences between treatments with respect to mean Length for the typical experiment; their results are visualized in Figure 1.

Something specific along the lines suggested here is more informative for the reader than a vague paragraph which asks the reader to do all the heavy lifting and fill in the missing information. You can add something on the fact that you use BH to control for the false-discovery rate or something to that effect.

As for the letter system to be used in your Figure 1, as long as you explain in the figure legend what that system is, readers should be able to make sense of it.

• Thanks for the correction on relevel() - that function should indeed be applied to d\$Treatment. 😊 – Isabella Ghement Sep 16 '18 at 19:51