I am doing my first data analysis and I have a hard time translating the experiment design to the model I want to fit. I have a couple of basic questions about the overall coding of the model, and a more complicated one that has been giving me headaches.

We have 6 amounts of the same fertiliser added to pots containing the same soil. Each quantity is added to 3 pots. Inside each pot, there are 2 plants of the same kind. From each plant, we measure photosynthesis, soil temperature, and soil moisture once a month for 6 months. The pots are in-situ, and the effect of time is known for photosynthesis.

how data is labelled

(fertil. quantity) 6 x (pots) 3 x (plants) 2 x (time) 6 = 216 observations.

> str(photo)
'data.frame':   216 obs. of  9 variables:
 $ quantity: int  169 169 169 169 169 169 76 76 76 76 ...
 $ pot     : Factor w/ 18 levels "a","b","c","d",..: 1 1 2 2 3 3 4 4 5 5 ...
 $ plant   : Factor w/ 36 levels "10e","11f","12f",..: 11 22 30 31 32 33 34 
35 36 1 ...
 $ month   : int  5 5 5 5 5 5 5 5 5 5 ...
 $ co2     : num  0.101 0.0669 0.1075 0.0893 0.0846 ...
 $ tsoil   : num  9 8.75 11.05 9.4 10.65 ...
 $ msoil   : num  16.4 18.8 14.4 7.8 15.3 ...

Can I consider that plant is nested in pot and pot is nested in fertilizer quantity, with model as follows?

full<-lmer(co2~quantity+factor(month)+tsoil+msoil+(quantity|pot/plant), data=photo, REML=FALSE)

If it’s the case, to avoid pseudoreplication, should I consider plant as subsamples and average their values per pot?

full2<-lmer(co2~quantity+factor(month)+tsoil+msoil+(1|quantity/pot), data=photo, REML=FALSE)

If I am not interested in time, is it ok to consider it as a random effect? A random slope comes from the relationship between Y and X interacting with months.

full3<-lmer(co2~quantity+tsoil+msoil+(1|quantity/pot)+(1+quantity|month), data=photo, REML=FALSE)

Now a little bit more complicated part. There are other variables of interest, that I hope will explain the response, that were measured from the soil after the experiment (after 6 months). So, these variables contain 18 observations (1 per pot) as time and plant aren’t taken into account. Can I even dream of adding them to the model? One idea at the moment is add the variables to the model based on their relationship to the intercepts of pots (n=18) from one of those models above. I haven't found any information about it, so I'd appreciate comments and suggestions.


First question:

You can look at the ICC to determine if the extra nested levels help your analysis. For example, if the variance of error at level 1 is .2 and the variance of error at level 2 is .3, then the ICC is .3/(.2+.3)=.6, meaning that level 2 matters a lot. This is in my field and you can determine a meaningful ICC based on the literature in your field. You can do an ICC for other levels, too, to determine if those levels help explain more. A four level model, plants in pots in fertilizer in time may not strike a balance between parsimony and information.

Second question: you could consider them subsamples, but that may complicate things such that there is a lack of balance between parsimony and information.

Third question: if you are not interested in time then you wouldn't include it at all in the model. Including it as a random effect means it has an effect.

Fourth question: you might include them as fixed effects at the pot level.

Final note 1: if you make your model longitudinal, make sure for autocorrelation by either including an autoregressive predictor or constraining the error matrix to follow an AR(1) process.

Final note 2: As an alternative to a classic longitudinal model, consider a dynamic multilevel model (Hamaker and Wichers, 2017), such that parameter estimates at time 1 are predictors of a future outcome variable. For example, perhaps a higher tsoil parameter estimate at time 1 is correlated with higher co2 at time 2.


A popular source in my field in Raudenbush and Bryk (2002). Others include OConnel and McCoach (2008), Hox (2010), and Snijders and Bosker (2012). In the Raudenbush and Bryk text, the procedure I refer to in #4 is in chapter 8 in the 2nd edition.

I would say that doing ICC by hand for all the levels will be easier because extracting the values from the functions is not always straightforward in R. Doing ICCs for additional levels is also in chapter 8 of the Raudenbush and Bryk text.

Deciding on whether to include another level can be done with ICC, but also with proportional variance explain (Raudenbush and Bryk, 2002, p. 75 and 79), chi-squared differences test, likelihood ratio test, AIC, BIC, and others (McCoach nad Black, 2008). This is really the beauty of HLM because you can set up many alternative models and see which one best fits the data. Of course, there are an infinite number of alternative models, so make sure to test ones that are based on theory!

  • $\begingroup$ Thank you for the nice suggestions! I'd very much like to make this model as simple as possible! I calculated the ICC "by hand" the way you mentioned and then confirmed it with icc function from sjstats library for a model with (1|pot/plant) and it was quite low. For (quantity|pot/plant) I can only get a icc for the whole structure icc(model, adjusted = TRUE), and it's quite high. Can I use the icc for pot and plant calculated from (1|pot/plant) to decide? On another note, could you kindly direct me to a source to do what you recommended on the 4th question? $\endgroup$ – Tomás Rimoas Oct 5 '18 at 8:26
  • $\begingroup$ @TomásRimoas I made edits to my answer in response to your question. Let me know if I can help further. If this answers your questions, please mark it as answered. Thanks! $\endgroup$ – Jay Schyler Raadt Oct 5 '18 at 17:20

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.