Hello to the CrossValidated community. I could use some help.

I am struggling with the analysis of an experiment, specifically with deciding and defining in code the structure of a mixed model for analysing it.

The experimental design is as follows:

We wanted to see if there are any differences on the growth of lesions caused by different combinations of pathogens (Treatment in the data).

To investigate this we used detached leaves and inoculated them with said Treatments. Each leaf was inoculated in two spots which were averaged to give us the Area of each lesion. Measurements were done per day for all leaves in the experiment, for days 4-9 after the inoculation (Day in the data below).

Two difference cultivars of potato plants were used in the experiment (Cultivar in the data).

Leaves were arranged in boxes (Box in the data) for a complete randomized block design. Each box contained 7 leaves of each cultivar, each leaf inoculated with one of the Treatment, so a total of 14 leaves in each box.

LeafID is a unique identifier for each of leaves in the experiment.

My understanding is that Box, Cultivar, Treatment and Day are crossed factors. LeafID is crossed with Day and Box and nested within Treatment and Cultivar. I am mentioning this to check my understanding of the experimental design and furthermore to understand if it's important for defining the model. I seem to remember reading that lmer accounts for nested or crossed designs internally but can't quote a source unfortunately.

Here is a pick at my data structure in R, Day and Area are coded as double, all other variables are coded as factors

Box Cultivar Treatment  Day  LeafID Area
1   Arsenal  Control    4    121    0.0000000
1   Arsenal  Control    5    121    0.0000000
1   Arsenal  Control    6    121    0.0000000
1   Arsenal  Control    7    121    0.0000000
1   Arsenal  Control    8    121    0.0000000
1   Arsenal  Control    9    121    0.0000000
1   Arsenal  A.solani   4    21     0.1745344
1   Arsenal  A.solani   5    21     0.3062256
1   Arsenal  A.solani   6    21     0.3721086
1   Arsenal  A.solani   7    21     0.4217482

My understanding is that this can be analysed as either a repeated measures or a time-series design, and the choice really depends on the research question. I chose to look at this as a repeated measures design, mostly because it's easier for me to grasp the concepts and statistics involved. I am more than happy to hear differing opinions of course!

I tried to fit my model the same way Crawley does in the R Book (p,695, "Mixed-effects models with temporal pseudoreplication")

Initially we came up with the following model:

m0.full <- lme(data=DLA2, Area ~ Treatment * Cultivar * Day  ,
random = ~Day | LeafID, na.action = na.exclude ,method = "ML")

The difference from what Crawley does is that our model has extra fixed effects and their interactions Treatment * Cultivar * Day and we are using Day and its interaction as we care about the progress over time and how this affects our growth.

What I am struggling with is deciding if Cultivar should be a fixed or random effect and if including Day in the fixed effects also gives the wrong interpretation.

I've spent countless hours going through mixed model tutorials and I know that going for a very complex model is not the way to go. I guess the most confusing part is deciding what to have as fixed and random. I've seen some people suggesting you can use a factor in both fixed and random parts of the model but i have trouble grasping the concept.

I would be happy to share the entire dataset if that helps, any suggestions of the best way to do this would be welcome.

I am looking forward for some suggestions on how to handle my analysis and what could be improved.

I would be happy using either lme or lmer.

Thank you all for your time!

EDIT: I added a plot of my data for each different Treatment*Cultivar combination. Each panel shows the full replicates corresponding to the combination

Daily growth in lesion radius of each Treatment*Cultivar combination


I came up with these two models, mod.rep takes into account the repeated measures of the Day variable, while mod.time uses a time-series approach as used by Crawley(The R Book). I also accounted for heteroscedasticity using the varIdent function.

Are those models correct ?

mod.rep <- lme(data=DLA2, Radius ~ Treatment * Cultivar, random = ~ 
Day|LeafID, na.action = na.omit ,method = "ML", weights = 

mod.time <- lme(data=DLA2, Radius ~  Day * Treatment * Cultivar, random = 
~1|LeafID, na.action = na.exclude ,method = "ML", weights = 

I still have some trouble translating the random effects formulas into words...

  • $\begingroup$ Does this Q&A help you stats.stackexchange.com/questions/4700/… $\endgroup$ – mdewey Aug 18 '17 at 12:12
  • $\begingroup$ Thank you for pointing it out. I have read it in the past and it does not help much with my particular questions regarding my dataset (like how to define the models and deal with the Day variable) $\endgroup$ – Ioannis Baltzakis Aug 18 '17 at 12:22
  • 1
    $\begingroup$ Thanks for the answer. Regarding the cultivar as a fixed effect, I would accept this without question if my main questions was how Cultivar affects the lesion growth together with the treatment. $\endgroup$ – Ioannis Baltzakis Aug 18 '17 at 12:44
  • 2
    $\begingroup$ In my case, I use 2 different Cultivars to get information on the effect of the Cultivar might have on the Treatment but not directly caring about it. And since the 2 cultivars are drawn from a very diverse choice of many different cultivars, let me quote Crawley : Random effects have factor levels that are drawn from a large (potentially very large) population in which the individuals differ in many ways, but we do not know exactly how or why they differ. $\endgroup$ – Ioannis Baltzakis Aug 18 '17 at 12:50
  • 2
    $\begingroup$ It is philosophically reasonable but not practical to treat a two-level factor as random (when using modern mixed models): see last paragraph of stats.stackexchange.com/questions/4700/… $\endgroup$ – Ben Bolker Aug 18 '17 at 16:16

In section 1.5 of Pinheiro and Bates (2000) Mixed Effects Models in S and S-Plus, you can find the reference for analyzing nested factors with the nlme package, which is related to lmer. The syntax for defining the nested factors can also be found in the same passage. The book is written about S, but these functions mostly work in R without problems.

For example to introduce a free intercept for a nested factor you could write random = ~1|Box/LeafID (This is just an example, I don't say that this is relevant to the specific experiment.) You can then gradually introduce slopes in the random effects, like update(your_model, random = ~Treatment|Box/LeafID), so that to get random effects for the slope of Treatment as well, and compare model fit. (You can similarly build up your way to the triple interaction term in the fixed effects.)

| cite | improve this answer | |
  • $\begingroup$ So how does the use of random=~Day|LeafID account for the repeated measure pseudoreplication in my first model ? I don't quite get the syntax (this is what Crawley does...) $\endgroup$ – Ioannis Baltzakis Aug 29 '17 at 14:51
  • $\begingroup$ It is interpreted as a different slope ($b_{day,i}$ for each leaf. (Having in mind that Y = ... random = ~1|Box means 'assume a different mean for each Box'. github.com/giswqs/Learning-R/blob/master/… -- see lines 105-6 for using an lme as repeated measures. I believe it is something in the lines of random = ~1|Box/Day -- different intercept for each Box and each Day, because Box is the entity and Day the repeated measure. $\endgroup$ – nijk Aug 29 '17 at 15:12
  • $\begingroup$ I am not taking Box into account at all. LeafID is just the identity of each leaf. Day is where pseudoreplication occurs, unless I try and analyse as a time series. So I am assuming you mean random = ~1|LeafID/Day ? $\endgroup$ – Ioannis Baltzakis Aug 29 '17 at 15:14
  • $\begingroup$ If your theoretical entity is the Leaf then yes. $\endgroup$ – nijk Aug 29 '17 at 15:24
  • $\begingroup$ So is random = ~Day|LeafID the same as random = ~1|LeafID/Day ? $\endgroup$ – Ioannis Baltzakis Aug 29 '17 at 15:55

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.