I need to compare total concentrations of a chemical across different tissues of a plant species. I sampled multiple tissues from each individual, and sampled at multiple sites. I do not know how to handle the individual and site covariates. There is good reason to expect a site effect on chemical concentrations. I need to control for the fact that concentrations will be correlated across tissues within individuals, but am only interested in whether concentrations differ consistently across tissues. I.e. are leaf concentrations always the lowest, etc.
I have two questions –
Question 1 what it the most appropriate analysis to use?
I do not think I have enough data for a mixed effects model, and would like to avoid using one if possible.
So far, I've used: (1) a three way ANOVA with tissue, individual, and site as factors.
3wy_mod <- (lm(total_conc ~ tissue+indv+site, data=data))
However, I feel like that isn't correct. I think that the model structure should have tissue nested within individual nested within site.
(2) to address the nesting issue, I've also used this model:
nest <- anova(lm(total_conc ~ tissue/indv/site, data=data))
The nested model seems correct. However, I'm not sure how to validate this model. This brings me to
Question 2
How do I assess model fit, check assumptions, etc with whichever model structure turns out to be appropriate?
Here is an attempt at a reproducible example:
data <- structure(list(indv = c(1L, 1L, 1L, 1L, 2L, 2L, 2L, 2L, 3L, 3L,
3L, 3L, 4L, 4L, 4L, 4L, 5L, 5L, 5L, 5L, 6L, 6L, 6L, 6L, 7L, 7L,
7L, 8L, 8L, 8L, 8L, 9L, 9L, 9L, 9L, 10L, 10L, 10L, 10L, 11L,
11L, 11L, 11L, 12L, 12L, 12L, 12L, 13L, 13L, 13L, 13L, 14L, 14L,
14L, 14L, 15L, 15L, 15L, 15L, 16L, 16L, 16L, 16L, 17L, 17L, 17L,
17L, 18L, 18L, 18L, 18L, 19L, 19L, 19L, 19L, 20L, 20L, 20L, 20L,
21L, 21L, 21L, 21L, 22L, 22L, 22L, 22L, 23L, 23L, 23L, 23L, 24L,
24L, 24L, 24L, 25L, 25L, 25L, 25L, 25L, 26L, 26L, 26L, 26L, 26L,
27L, 27L, 27L, 27L, 28L, 28L, 28L, 28L), tissue = structure(c(3L,
5L, 1L, 4L, 3L, 5L, 1L, 4L, 3L, 5L, 1L, 4L, 3L, 5L, 1L, 4L, 3L,
5L, 1L, 4L, 3L, 5L, 1L, 4L, 3L, 5L, 1L, 3L, 1L, 5L, 4L, 3L, 1L,
5L, 4L, 3L, 1L, 5L, 4L, 3L, 1L, 5L, 4L, 3L, 1L, 5L, 4L, 3L, 1L,
5L, 4L, 3L, 1L, 5L, 4L, 3L, 5L, 1L, 4L, 3L, 5L, 1L, 4L, 3L, 5L,
1L, 4L, 3L, 5L, 1L, 4L, 3L, 5L, 1L, 4L, 3L, 5L, 1L, 4L, 3L, 5L,
1L, 4L, 3L, 1L, 5L, 2L, 3L, 1L, 5L, 4L, 3L, 1L, 5L, 4L, 4L, 3L,
1L, 5L, 4L, 3L, 1L, 5L, 2L, 4L, 3L, 1L, 5L, 4L, 3L, 1L, 5L, 4L
), .Label = c("flower", "fruit", "leaf", "pollen", "stem"), class = "factor"),
site = structure(c(1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L,
1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L,
1L, 1L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L,
2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L,
3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L,
3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 4L, 4L,
4L, 4L, 4L, 4L, 4L, 4L, 4L, 4L, 4L, 4L, 4L, 4L, 4L, 4L, 4L,
4L, 4L, 4L, 4L, 4L, 4L, 4L, 4L, 4L, 4L, 4L, 4L, 4L), .Label = c("a",
"b", "c", "d"), class = "factor"), total_conc = c(8251.26,
3845.59, 6681.23, 504.47, 21494.77, 13514.46, 14244.08, 1091,
4156.25, 8301.19, 12968.21, 5536.14, 9034.05, 11401.89, 12010.11,
2013.21, 14997.58, 12266.25, 23068.28, 4784.39, 9326.1, 5522.93,
9814.8, 3467.77, 7107.22, 6129.18, 7330.3, 7242.1, 8668.38,
5214.82, 1891.02, 5918.55, 11268.32, 3267.63, 3648.9, 13406.13,
15847.31, 6775.97, 1781.21, 10943.45, 14477.16, 7890.34,
1738.36, 14467.25, 11777.99, 5551.34, 1136, 5537.86, 18114.18,
2895.41, 999.91, 3633.19, 7243.96, 4856.73, 8317.76, 15885.53,
6811.57, 15818.11, 1316.81, 15347.71, 8068.19, 14688.62,
2247.32, 17762.4, 6628.74, 16494.26, 11255.5, 19459.68, 5960.92,
20607.73, 3226.06, 11572.8, 7924.07, 12091.07, 2934.25, 12734.17,
5158.48, 16520.8, 1431.31, 14007.17, 4689.99, 23774.68, 5799.91,
12903.34, 18982.81, 9125.87, 14869.07, 20022.23, 28254.47,
9459.41, 1977.36, 22562.52, 22195.51, 12376.81, 2284.3, 1358.63,
24866.6, 15789.27, 9544.07, 275.32, 17669.95, 18935.79, 9357.06,
18818.68, 1676.3, 23224.24, 25432.66, 12328.59, 1367.51,
24176.56, 26710.84, 21184.79, 1642.89)), .Names = c("indv",
"tissue", "site", "total_conc"), class = "data.frame", row.names = c(NA,
-113L))
3wy_mod <- anova(lm(total_conc ~ tissue+indv+site, data=data))
nest <- anova(lm(total_conc ~ tissue/indv/site, data=data))
Best Answer
The crossed, non-nested model is incorrect, as you suspect, and seems to be over-specified for your posted data. It's important to get the nesting correct, which can be confusing. This page is a good guide. With multiple individuals per site and each individual restricted to one site, the individuals are nested within site. If you were to use your
anova(lm())
approach you should be writing:which is expanded to:
So this model examines the main effect of
tissue
plus its interactions withsite
and withsite:indiv
, without main effects forsite
orindiv
. The result of youranova(lm())
on your posted data* would then be:You only have 1 residual degree of freedom because most cases have only one measurement on each tissue from an individual. So the three-way interaction is fairly meaningless. Although results from the Type I sums of squares used by
anova()
, when applied to unbalanced data, can depend on the order of entry of variables into the model formula, in this case there really is only one order of entry that makes sense, so that's not such an issue.These data are highly unbalanced, with observations on
fruit
restricted to 2indiv
at one particularsite
. If that's characteristic of your real data, then you probably should consider removing such tissues with low numbers of observations from your analyses. Would you (or a skeptical reviewer) really trust analyses based on 2 observations from a singlesite
? Removing thefruit
observations from these data, at least, would leave a reasonably balanced design that would probably be good enough for standard nested ANOVA.Note that both the nested ANOVA approach and mixed models with intercept-only terms like
(1|indiv)
are making an implicit assumption that all effects of individuals and sites ontotal_conc
are additive without regard to thetissue
being evaluated. You must use your knowledge of the subject matter as to the validity of that assumption. With the raw means oftotal_conc
ranging from 2804 (pollen
) to 16844 (fruit
) I would be a bit worried about that assumption. In principle the formalism of a mixed model could incorporate random effects that differ amongtissue
values, but you don't seem to have enough data for that. The suggestion from @MarkWhite to try a Bayesian multilevel model might help, but I have no experience with such models.If there is some irreducible imbalance in the data then I think that many would argue for a mixed-model approach. My guess is that the errors that you discover with mixed models due to small numbers of observations would also be seen in the
anova(lm())
approach if you dug more deeply. For example, that approach on the crossed model applied to your data provides what seems to be a perfectly reasonable ANOVA table.Nevertheless, the underlying crossed model can't define 3 regression coefficients because of singularities:
*Note that the
indiv
values in the posted data need to be converted to factors first.