5
$\begingroup$

I am trying to assess the effects of an experimental treatment on the insect fauna of artificial ponds. The treatment is applied to the entire pond.

I was able to sample each pond four times, each sample collected from a different spatial location within the pond, but at the same moment in time. So I have replicated similar samples from each pond.

I have used a mixed effects model approach using each sample as my lowest-level experimental unit (rows in a matrix), treatment as a fixed effect, and pond identification as a random effect.

However, recently I was advised to pool all samples from the same pond together and use the pond as my lowest-level experimental unit and treatment as a fixed effect. I feel that this approach is the simplest one, but that I would be losing statistical power just for the sake of simplicity.

Which approach should I use?

$\endgroup$
1
  • 1
    $\begingroup$ Simplicity isn't very convincing as an argument if there is no problem to run and interpret a more complex analysis that makes better use of the available information. (I can't tell you what you should run without seeing the data though; both approaches rely on assumptions and I would want to have an impression regarding whether these are reasonable.) $\endgroup$ Commented May 9, 2023 at 13:26

2 Answers 2

7
$\begingroup$

Murtaugh (2007) makes the point that if you have a nested (all samples within a group share the same predictor variables), balanced (all groups are sampled at the same intensity) design and a linear (responses treated as Gaussian) model, then you will get exactly the same inferences at the population level (i.e. about differences between treatment) if you average the samples each group as if you fit a multi-level model (see example below). In this case he recommends averaging because you're less likely to make a mistake in the analysis.

However, it is best not to pool if your situation deviates from this case.

  • If your responses are pooled but not balanced you may be able to use weights to fit the aggregated data
  • If your design is not pooled (e.g. a randomized complete block design) but you only have two levels of the treatment you may be able to use a paired $t$-test to analyze within-block differences (this test is a special case of a LMM estimated via REML)
  • If your responses are non-Gaussian (and you want to handle that by using a GLM rather than by transforming) you should almost certainly be using the multi-level/mixed GLMM
  • If you're interested in quantifying the variability at different levels (between-pond vs. within-pond) then you need the multilevel model

For example (the simulation here fits the assumptions of the Gaussian LMM perfectly because that was easiest to do, but Murtaugh's point about the equivalence of pooled vs multilevel linear models holds regardless of the conditional distribution of the data as long as the design is nested and balanced):

library(lmerTest)

## simulate data
dd <- expand.grid(pond = factor(1:10),
                  sample = factor(1:4))
dd$treat <- ifelse(as.numeric(dd$pond) <= 5, "C", "T")
dd$y <- simulate(~ treat + (1|pond),
                 seed = 101,
                 newdata = dd,
                 newparams = list(beta = c(1, 2),
                                  theta = 1,
                                  sigma = 1))[[1]]

## aggregate data to pond level
dd2 <- aggregate(dd, y ~ pond, FUN = mean)
dd2 <- merge(dd2, unique(dd[c("pond", "treat")]))

## fit both models (could also use nlme::lme())
mod1 <- lmer(y ~ treat + (1|pond), dd)
mod2 <- lm(y ~ treat, dd2)

## results

printCoefmat(coef(summary(mod1)), digits = 3)

            Estimate Std. Error    df t value Pr(>|t|)   
(Intercept)    0.899      0.330 8.000    2.72   0.0262 * 
treatT         2.260      0.467 8.000    4.84   0.0013 **
---

printCoefmat(coef(summary(mod2)), digits = 3)
            Estimate Std. Error t value Pr(>|t|)   
(Intercept)    0.899      0.330    2.72   0.0262 * 
treatT         2.260      0.467    4.84   0.0013 **

Murtaugh, Paul A. 2007. “Simplicity and Complexity in Ecological Data Analysis.” Ecology 88 (1): 56–62.

$\endgroup$
0
0
$\begingroup$

Thank you very much for this complete answer @BenBolker. I guess things do get messy with my data because my responses are non-Gaussian. Your example has quite the design I have (just simpler).

Because my responses are abundance counts I am using a negative binomial distribution to model abundances (negative binomial was a better fit than poison because of overdispersion). I made a few tests and indeed, if I use the Gaussian distribution I get the exact same Wald statistics and p-values for the mixed model (but not always) and the one adding up abundances from multiple samples of the same pond. But for the negative binomial distribution, it is not the same. In this case, p-values are smaller for the mixed model approach.

See this example with simulated data:

> dd <- expand.grid(pond = factor(1:10),
+                   sample = factor(1:4))
> dd$treat <- ifelse(as.numeric(dd$pond) <= 5, "C", "T")
> 
> set.seed(2)
> s1_C <- rnbinom(n = 5, mu = 140, size = 50)
> s1_T <- rnbinom(n = 5, mu = 130, size = 50)
> s2_C <- rnbinom(n = 5, mu = 140, size = 50)
> s2_T <- rnbinom(n = 5, mu = 130, size = 50)
> s3_C <- rnbinom(n = 5, mu = 140, size = 50)
> s3_T <- rnbinom(n = 5, mu = 130, size = 50)
> s4_C <- rnbinom(n = 5, mu = 140, size = 50)
> s4_T <- rnbinom(n = 5, mu = 130, size = 50)
> 
> dd$y <- c(s1_C, s1_T, s2_C, s2_T, s3_C, s3_T, s4_C, s4_T)
> 
> dd_summaryzed <- aggregate(dd, y ~ pond, FUN = sum)
> dd_summaryzed <- merge(dd_summaryzed, unique(dd[c("pond", "treat")]))
> 
> library(glmmTMB)
> library(car)
> 
> mod_gaussian <- glmmTMB(y ~ treat + (1|pond), data = dd, family = "gaussian")
> Anova(mod_gaussian)
Analysis of Deviance Table (Type II Wald chisquare tests)

Response: y
       Chisq Df Pr(>Chisq)
treat 0.5421  1     0.4616
> 
> mod_gaussian_sum <- glmmTMB(y ~ treat, data = dd_summaryzed, family = "gaussian")
> Anova(mod_gaussian_sum)
Analysis of Deviance Table (Type II Wald chisquare tests)

Response: y
       Chisq Df Pr(>Chisq)
treat 0.5421  1     0.4616
> 
> mod_negbin <- glmmTMB(y ~ treat + (1|pond), data = dd, family = "nbinom1")
> Anova(mod_negbin)
Analysis of Deviance Table (Type II Wald chisquare tests)

Response: y
      Chisq Df Pr(>Chisq)
treat  0.48  1     0.4884
> 
> mod_negbin_sum <- glmmTMB(y ~ treat, data = dd_summaryzed, family = "nbinom1")
> Anova(mod_negbin_sum)
Analysis of Deviance Table (Type II Wald chisquare tests)

Response: y
       Chisq Df Pr(>Chisq)
treat 0.6612  1     0.4161 

Now an example with real data:

> Abundances <- rowSums(com_SS4[isolation_SS2_3_4 == "120",])
> Abundances_summarized <- tapply(Abundances, INDEX = as.character(ID_SS2_3_4[isolation_SS2_3_4 == "120"]), FUN = sum)
> 
> data <- data.frame(ID = ID_SS2_3_4[isolation_SS2_3_4 == "120"], treatments = treatments_SS2_3_4[isolation_SS2_3_4 == "120"], ab = Abundances)
> data_summaryzed <- data.frame(ID = ID_SS1[isolation_SS1 == "120"], treatments = treatments_SS1[isolation_SS1 == "120"], ab = Abundances_summarized)
> 
> library(glmmTMB)
> library(car)
> 
> mod1 <- glmmTMB(ab ~ treatments + (1|ID), data = data, family = "gaussian")
> Anova(mod1)
Analysis of Deviance Table (Type II Wald chisquare tests)

Response: ab
            Chisq Df Pr(>Chisq)  
treatments 4.9448  2    0.08438 .
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
> 
> mod2 <- glmmTMB(ab ~ treatments, data = data_summaryzed, family = "gaussian")
> Anova(mod2)
Analysis of Deviance Table (Type II Wald chisquare tests)

Response: ab
            Chisq Df Pr(>Chisq)  
treatments 4.9448  2    0.08438 .
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
> 
> mod1 <- glmmTMB(ab ~ treatments + (1|ID), data = data, family = "nbinom1")
> Anova(mod1)
Analysis of Deviance Table (Type II Wald chisquare tests)

Response: ab
            Chisq Df Pr(>Chisq)  
treatments 5.4888  2    0.06429 .
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
> 
> mod2 <- glmmTMB(ab ~ treatments, data = data_summaryzed, family = "nbinom1")
> Anova(mod2)
Analysis of Deviance Table (Type II Wald chisquare tests)

Response: ab
            Chisq Df Pr(>Chisq)
treatments 4.3773  2     0.1121

Note that in this case, p-values are 'almost' significant only for the mixed model. Indeed, if run all of my analyses without the mixed-model approach I lose a lot of significant results. Maybe I do have some variation among samples of the same pond. For instance, some taxa may be more prone to be sampled in the first samples while others in the last ones.

$\endgroup$

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.