R – Representing Repeated Measurements Within Sample Plots in Linear Mixed Effects Models

lme4-nlmerrepeated measures

I have collected data on gas fluxes from plots of soil subjected to 5 different treatments ("D2", "K2", "M", "N", and "O2"), which also possessed variable clay contents. The experiment was laid out in a randomized complete block design, with 4 replications. Within each plot, two separate measurements of flux were performed. The resulting data.frame resembles the following:

block   treatment   subsample       flux        clay
1           D2          1           112068.6003 14.8
1           D2          2           129223.1641 14.8
1           K2          1           256712.4712 15.5
1           K2          2           113343.9756 15.5
1           M2          1           85794.47834 16.4
1           M2          2           -33620.6990 16.4
1           N           1           70283.98133 18.2
1           N           2           49569.84621 18.2
1           O2          1           100553.1116 13.4
1           O2          2           38885.99674 13.4
2           D2          1           96968.58451 15.8

I want to build a linear mixed effect model that takes account of this subsampling, and have come up with:

flux.lme <- lme(flux ~ treatment + block, random = ~1|subsample)

Which produces an ANOVA table:

> anova(flux.lme)
               numDF denDF   F-value p-value
(Intercept)        1    26 158.15781  <.0001
treatment          4    26   8.88691  0.0001
clay               1    26   1.72640  0.2003
block              3    26   1.59188  0.2153
treatment:clay     4    26   1.73011  0.1736

This output seems a bit strange to me, as the denominator degrees of freedom should not be 26, which is taking each repeated sampling as independent experimental unit. It should instead be based on the number of “main plots”, which is 20. In my case the denominator d.f. should be 20-1-3-1-4-4=7. Is is possible to instruct the lme() function to use this value?

Best Answer

I would recommend getting a few resources if you're just getting started fitting mixed models. I think Mixed Effects Models and Extensions in Ecology in R (Zuur et al. 2009) is a nice (approachable) place to start. I often peruse http://glmm.wikidot.com/faq, as well. It is R specific but I always learn a lot about mixed models in general.

I think it will help you to write out your study design explicitly to make sure you are accounting for all the variation in your model.
You have 4 blocks (account for variation among blocks).
You have 5 plots in each block for 20 plots total (account for variation within blocks).
You have 2 subsamples in each plot for 40 observations total (account for within plot variation).

I had to make a fake dataset, next time please make any coding questions reproducible by providing data (possibly with dput()).

block = factor(rep(1:4, each = 10))
treatment = factor(rep(rep(1:5, each = 2), 4))
subsample = factor(rep(1:2, 20))
flux = rnorm(40, 10, 7)
clay = rep(rnorm(20, 10, 2), each = 2)

dat1 = data.frame(block, treatment, subsample, flux, clay)

Here is your first model. You account for block to block variation by using block as a fixed effect. You use the two-level variable "subsample" as a random effect.

require(nlme)

fluxlme = lme(flux ~ treatment*clay + block, random = ~1|subsample, data = dat1)
summary(fluxlme)
anova(fluxlme)

At the bottom of the summary you will see this:

Number of Observations: 40
Number of Groups: 2

This is where I check if my model reflects my design. Do you have two groups with 20 observations per group? Nope, you have 20 groups (plots) with 2 observations per group. You failed to account for one of the levels of variation in your design in your model.

Make a variable to represent your plots nested in blocks. Because treatment represents each plot within each block, you can combine the block variable and the treatment variable to make a new variable with a unique identifier for each plot. This is called explicit nesting. You can read more about implicit vs explicit nesting online, starting with the website I listed above. I've found that using explicit nesting avoids a lot of confusion and mistakes on my part.

dat1$plot = with(dat1, interaction(block, treatment) )

fluxlme2 = lme(flux ~ treatment*clay + block, random = ~1|plot, data = dat1)
summary(fluxlme2)
anova(fluxlme2)

The subsamples are the observation-level measurement, which is represented by the residual error term in linear mixed models.

Related Solutions

Solved – Repeated-measures linear mixed effect model

okay should work out okay then-- so

Yes or you can use the lmer() and lme4. There is another one but I don't remember off the top of my head I think it is just lme?

2.You have a nested structure so yes you need (1|sample/participant)

Did you plot rating score vs stim.level to see evidence of a quadratic relationship? If not try plotting to see-- i you do see a quadratic pattern then yes you should add stim.level quadratic effect by

model <- lmer (rating.score ~ stim.level + I(stim.level^2) factor + stim.level*factor +(1|sample/participant) , mydata)

To reply to the comment

so if you are fitting a parabola and not a line you are fitting the generic

 y= a + b*x + c*x^2

so you need the linear and quadratic term so stim.level is the linear term and I(stim.level^2) would be the quadratic term. Try writing your model out on paper in equation form like

rating score = b_0 + b_1 * stim.level + b_2 *stim.level^2 + b.3*(stim.level*factor)

and if wanted you could have an interaction with the squared term but that might not be directly interpretable. 2) Remember you are fitting a parabola NOT a line http://biopt.ub.edu/_/rsrc/1257372769717/force-detection/equipartition-theorem/ex4%20X%20potential.jpg?height=315&width=420 something that looks like that. Did you plot stim.level against rating.score and see a quadratic relationship?

Solved – Maximal model for linear mixed-effects model for repeated mesaures design

The maximal structure would need to include also a random effect for the interaction between color and shape, that is:

Y ~ color * shape + (color + shape + color:shape | subject)

This will result in all your predictors (color, shape and their interaction) having a fixed effect (constant for all subjects), and a random effect (individual fluctuations around the estimated fixed effect). In this sense the model is the maximal one. Note that it might not be fully equivalent to a repeated-measures ANOVA as it doesn't make equally strict assumptions on the correlational structure (see Tom's answer).

If you don't include the interaction in the random effect part of the formula, individual variation in the interaction effect will not be considered as "random", and the model will not be equivalent to a repeated-measures ANOVA. Of course, the variance of the random deviates for the interaction (or any other random effect) might be so small that including it in the model do not improve much the fit. You can check this not only with the AIC, but with a likelihood ratio test, as model with vs without one random effect are nested one another. In principle if the likelihood ratio test is not significant, it means that you can safely remove that random effect. Simplifying the random effect structures by removing negligible components would be an example of what in the article you linked is called data-driven approach.

You can simplify the model in this way, and it would still be equivalent to a repeated-measures ANOVA:

Y ~ color*shape + (1|subject) + (0+color|subject) + (0+shape|subject) + (0+color:shape|subject)

This syntax tells lmer to not estimate the correlations of random deviates across subjects. The drawback here is that, for example, you won't be able to tell whether subjects that have a large effect of color tend to have also a larger effect of shape (or smaller effect, in case of negative correlation).

You can easily include a between-subjects predictor, the only difference is that you can't add a random effect for it. "gender" for example cannot have a random effect grouped according to subject, but it can interact with the other fixed effects, e.g.:

Y ~ color * shape * gender + (color + shape + color:shape | subject)

Best Answer

Related Solutions

Solved – Repeated-measures linear mixed effect model

Solved – Maximal model for linear mixed-effects model for repeated mesaures design

Related Question