I have a study where the initial height of seedlings were measured at year 0 and again 11 years after treatments were applied and a net change in height was recorded.

I have a split-plot design as follows:

  • Treatment (whole-plot factor, 4 levels: TMT1, TMT2, TMT3, CONTROL)
  • Unit (replicates of treatment, 20 total: 5 for each treatment, tmts randomly assigned to units)
  • Moisture (split-plot factor, 3 levels: LOW, MOD, HIGH within each Unit)
  • Plot (2 replicate plots within each moisture level)
  • Group (5 levels: OAKS, PINES, MAPLES, ELMS, OTHER, etc.)
  • ID (the tag # of the seedling)
  • HeightChg (response variable, continuous)

I want to know:

  1. Is the effect of treatment type on change in seedling height significant?
  2. Is the effect of moisture on change in seedling height significant?
  3. Is the effect of species group on change in seeding height significant?

It isn't a block design in the traditional sense because all treatments aren't contained within some unit. Rather, each unit contains the same three abstract moisture levels (which will vary some by unit) but treatments are randomly assigned applied to the whole units in the beginning. That is: Treatment < Unit < Moisture Level < Plot < ID

This is my model so far:

library(lmerTest) #load package
mod1 = lmer(HeightChg ~ treatment * moisture * group + (1|unit), mydata) #11yr height difference

However, when I do a Q-Q plot of the model, I have fat tails on both ends:


Here is the plotted distribution of the change in height variable:


My colleague suggested one possibility is to transform the change in height variable since it's the extreme values in the dataset like dieback (abnormally large height loss) and rapid growth (abnormally large height change) that is causing kurtosis in the distribution.

He also said that before transforming the data to get a more linear fit, I should try to properly define my error matrix and random effects in my code. I've scoured the web looking for identical designs that have variable nesting and have found some good information, but haven't found an example that mirrors mine precisely.

Truthfully, the seedling heights were measured multiple times over various intervals but I haven't figured out the proper code for my split-plot design, let alone integrating unbalanced repeated measures of height, so I took the difference between the year 11 measurement and the year 0 measurement to keep it simple.

Any help you can provide how to use lme, lmer, or glmer to properly construct the model as it applies to this specific design (additive, multiplicative, error terms, nesting, etc.) would be greatly appreciated! I plan to run anova(mod1) once it's correct and try to calculate the p-values.


1 Answer 1


It's a little hard to say without being able to look at the data in more detail.

  • you should probably include plot-level variation, i.e. use (1|unit/plot) rather than (1|unit) (assuming there is more than one observation per plot, which I'm guessing there is — presumably several different groups [oaks, elms, pines] occur in each plot?)
  • technically, for the maximal model I think you should include (1+group*moisture|unit) + (1+group|unit:plot), to allow for all possible variation (and covariation) of the effects of predictors that vary within experimental groups. More parsimoniously, (1|unit) + (1|unit:group) + (1|unit:group:moisture) + (1|unit:plot) + (1|unit:plot:group) will allow for all possible variation that can be expressed by a compound symmetric model (i.e., the same correlations between every possible pair of (say) moisture levels within a unit/group combination)
  • it never hurts to try a log transformation; more complex transformations such as Box-Cox aren't always worth the trouble/loss of interpretability
  • the best general advice is to plot your residuals with respect to different predictor variables and see if you can identify a particular group/combination of groups that are particularly odd (although you should be fitting at least the means of these groups well since all of your predictors are categorical and you're fitting the full factorial model). The DHARMa package (and the plotResiduals package with the form argument specified) is good for this, as is broom.mixed::augment() + ggplot2.
  • $\begingroup$ The best analogy is imagine a 6-sided die. The square face is a single unit for which a randomly selected treatment has been applied to it. The 6 dots are seedling plots. The die is oriented parallel to a hillside so 2 plots at the top of the hill, 2 midslope, and 2 at the bottom (3 moisture levels). The log transform I tried improved the fit greatly. [1 / height change] was even better with a lower AIC. Will including group and moisture on the random effects side of the "|" not cause an issue if they are also defined as fixed effect factors in the model? @Ben Bolker $\endgroup$
    – kryogenic7
    Commented Jan 26 at 23:24
  • $\begingroup$ Including group and moisture on the grouping-variable (right) side of | is OK if they are crossed/interacted with (:) or nested within (/) other grouping variables. $\endgroup$
    – Ben Bolker
    Commented Jan 26 at 23:30
  • $\begingroup$ Awesome, I definitely have a learning curve with the coding. You've given me a direction to go with your initial Answer, as well as some packages to check out for the residuals. If the effects are significant, my next step is to use the 'emmeans' package to do some post-hoc testing. Thanks for your help! @Ben Bolker $\endgroup$
    – kryogenic7
    Commented Jan 26 at 23:58

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

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