dat19 = 
  read.csv("2019.biomass.csv") %>%
  mutate(Cage=trimws(Cage))%>%
  group_by(Plot, Quad, Cage, Fire) %>%
  summarise(Biomass = sum(Biomass.g.m.2, na.rm = T)) %>%
  mutate(Plot=as.factor(Plot))%>%
  mutate(Fire=factor(Fire,levels=c("N","L","H")))

1 2019 Community data

We used package glmmADMB, which is built on the open source AD Model Builder nonlinear fitting engine, to fit generalized linear mixed models.

Our model has a gamma response distribution with a log link and includes Plot as a random effect to account for within-plot error.

We modelled 2019 biomass as a function of fire treatment, herbivory treatment, and their interaction

To cite package ‘glmmADMB’ in publications use:

Fournier DA, Skaug HJ, Ancheta J, Ianelli J, Magnusson A, Maunder M, Nielsen A, Sibert J (2012). “AD Model Builder: using automatic differentiation for statistical inference of highly parameterized complex nonlinear models.” Optim. Methods Softw., 27, 233-249.

to cite R use:

R Core Team (2021). R: A language and environment for statistical computing. R Foundation for Statistical Computing, Vienna, Austria. URL https://www.R-project.org/.

#run the model

glmm19=glmmADMB::glmmadmb(Biomass~Cage+Fire+Cage*Fire,  
              random=~1 | Plot, data = dat19, family = "gamma",
              link = "log")

#calculate confidence intervals for the estimates
glmm_terms <-  
  confint(glmm19) %>%
  as.data.frame %>%
  rownames_to_column("term") %>%
  rename(lwr = `2.5 %`, 
         upr = `97.5 %`) %>%
  full_join(
    coef(glmm19) %>%
      as.data.frame() %>%
      rownames_to_column("term") %>%
      rename(estimate = '.') )
              
summary(glmm19)
## 
## Call:
## glmmADMB::glmmadmb(formula = Biomass ~ Cage + Fire + Cage * Fire, 
##     data = dat19, family = "gamma", link = "log", random = ~1 | 
##         Plot)
## 
## AIC: 1824.3 
## 
## Coefficients:
##             Estimate Std. Error z value Pr(>|z|)    
## (Intercept)   4.9926     0.1640   30.45   <2e-16 ***
## CageY         0.7157     0.2319    3.09    0.002 ** 
## FireL         0.0382     0.2319    0.16    0.869    
## FireH        -0.1510     0.2336   -0.65    0.518    
## CageY:FireL  -0.2443     0.3279   -0.74    0.456    
## CageY:FireH   0.3483     0.3285    1.06    0.289    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Number of observations: total=144, Plot=72 
## Random effect variance(s):
## Group=Plot
##             Variance   StdDev
## (Intercept) 1.14e-05 0.003377
## 
## Gamma shape parameter: 1.5502 (std. err.: 0.1688)
## 
## Log-likelihood: -904.147
glmm_terms
##          term        lwr       upr    estimate
## 1 (Intercept)  4.6712242 5.3138964  4.99256032
## 2       CageY  0.2612164 1.1700909  0.71565363
## 3       FireL -0.4162803 0.4925942  0.03815695
## 4       FireH -0.6088861 0.3068483 -0.15101892
## 5 CageY:FireL -0.8869897 0.3984331 -0.24427829
## 6 CageY:FireH -0.2955020 0.9921159  0.34830696
Table 1.1: Results of negative gamma regression for 2019 biomass
  Biomass 2019
Predictors Estimates Conf. Int (95%) p-Value
(Intercept) 147.31 106.83 – 203.14 <0.001
Cage (Yes) 2.05 1.30 – 3.22 0.002
Fire Energy (Low) 1.04 0.66 – 1.64 0.869
Fire Energy (High) 0.86 0.54 – 1.36 0.518
Cage (Yes) * Fire Energy
(Low) 		0.78 0.41 – 1.49 0.456
Cage (Yes) * Fire Energy
(High) 		1.42 0.74 – 2.70 0.289
Observations 144 
ggplot(dat19, aes(y=Biomass, x=Cage)) + 
  geom_boxplot(notch=TRUE, fill=c("goldenrod","darkolivegreen"))+
  geom_jitter(shape=16, position=position_jitter(0.2)) +
  theme_minimal()+
  theme(text = element_text(size = 15))