GLMMS for Stress-Hardening IPAM Data

Coding

Stats

GLMM

Linear Mixed Models

PAM

Stress-Hardening

glmmTMB

Kevin

Published

January 29, 2023

Based on the PPT from Kevin, I think my Fv/Fm data falls into a Percentage/Proportion data type, where numerical data is bound from 0-1. The assumption for this is that there is overdispersion (the mean is < the variance). The code provided in the PPT as an example was: glmmTMB::glmmTMB(survival ~ treatment + (1|sample), family= list(family=”beta”, link=”logit”)

PAM data during 28-day temperature treatment

#what are the differences in Fv/Fm between treatments? (pre-CBASS)
ipam_tidy_data %>% 
  dplyr::filter(!Treatment_period == "CBASS") %>%
  drop_na(Colony) %>% #this is because there are some corals that never matched up to IPAM values so they're just NA
  dplyr::select(Date, Colony, Puck, Tank, fvfm, Treatment, Species) %>% 
  dplyr::filter(fvfm < 0.75) %>% #filter out the outliers, 0.75 was set in Cunning et al 2021
  dplyr::filter(Date <= "2022-04-20") %>% 
  mutate(Date = as.factor(Date)) %>% 
  dplyr::group_by(Treatment, Date, Species) %>% 
  ggplot(., aes(x=Date, y=fvfm, fill = Treatment)) + 
  geom_boxplot() + 
  facet_wrap(~Species) + 
  scale_x_discrete(labels = c("Mar 16", "Mar 21", "Apr 6", "Apr 20")) + 
  theme_classic() + labs(y = "Fv/Fm") +
  scale_fill_manual(labels=c("Control", "Variable"), values = c( "#60DBDB", "#F54A34"))

#stats for the boxplot of pre-CBASS raw Fv/Fm values
ipam_tidy_data %>% 
  filter(Species == "Pclivosa") %>% 
  drop_na(Tank, Colony, Treatment, fvfm, Date) %>% 
  filter(fvfm < 1) -> Pcli_ipam_data_preCBASS #some fvfm values are >1 which is incorrect....

ipam_tidy_data %>% 
  filter(Species == "Acervicornis") %>% 
  drop_na(Tank, Colony, Treatment) %>% 
  filter(fvfm < 1) -> Acer_ipam_data_preCBASS #some fvfm values are >1 which is incorrect....

summary(glmmTMB::glmmTMB(fvfm ~ Treatment + (1|Colony) + (1|Tank/Date), family=list(family="beta", link="logit"), data=Pcli_ipam_data_preCBASS))

summary(glmmTMB::glmmTMB(fvfm ~ Treatment + (1|Colony) + (1|Tank/Date), family=list(family="beta", link="logit"), data=Acer_ipam_data_preCBASS))

For the glmmTMB function, I get this error: Warning: some components missing from ‘family’: downstream methods may failWarning: specifying ‘family’ as a plain list is deprecated So i think there is something wrong with the family= part of the formula.

For Pclivosa, this was the summary of the model: Family: beta ( logit ) Formula: fvfm ~ Treatment + (1 | Colony) + (1 | Tank/Date) Data: Pcli_ipam_data_preCBASS

 AIC      BIC   logLik deviance df.resid

-2191.1 -2166.7 1101.6 -2203.1 429

Random effects:

Conditional model: Groups Name Variance Std.Dev. Colony (Intercept) 3.924e-03 6.264e-02 Date:Tank (Intercept) 1.915e-02 1.384e-01 Tank (Intercept) 1.688e-09 4.109e-05 Number of obs: 435, groups: Colony, 3; Date:Tank, 12; Tank, 4

Dispersion parameter for beta family (): 787

Conditional model: Estimate Std. Error z value Pr(>|z|) (Intercept) 0.03057 0.06725 0.455 0.649 TreatmentVariable -0.05018 0.08019 -0.626 0.531

For Acervicornis, this was the summary of the model: Family: beta ( logit ) Formula: fvfm ~ Treatment + (1 | Colony) + (1 | Tank/Date) Data: Acer_ipam_data_preCBASS

 AIC      BIC   logLik deviance df.resid

-1524.8 -1501.2 768.4 -1536.8 370

Random effects:

Conditional model: Groups Name Variance Std.Dev. Colony (Intercept) 3.403e-03 5.834e-02 Date:Tank (Intercept) 6.649e-02 2.579e-01 Tank (Intercept) 1.056e-09 3.249e-05 Number of obs: 376, groups: Colony, 3; Date:Tank, 12; Tank, 4

Dispersion parameter for beta family (): 295

Conditional model: Estimate Std. Error z value Pr(>|z|) (Intercept) 0.06647 0.11085 0.600 0.549 TreatmentVariable 0.19180 0.14938 1.284 0.199

To calculate the p-value, there are three options (according to Kevin’s PPT): 1. Likelihood ratio test (anova of model 1 vs. model 2 with and without the fixed effect) 2. Satterthwaite approximation (for normal/gaussian distributions) 3. Wald/chi-square approximation (for binomial, poisson distributions)

I think my data is a binomial distribution, so I will try the Wald/chi-square approximation.

I found this useful PDF as part of the glmmTMB vignette called “model_evaluation”.

https://cran.r-project.org/web/packages/glmmTMB/vignettes/model_evaluation.pdf

It mentions this package called “DHARMa”, which helps create residual diagnostic plots.

Pcli_fvfm_GLMM <- glmmTMB::glmmTMB(fvfm ~ Treatment + (1|Colony) + (1|Tank/Date), family=list(family="beta", link="logit"), data=Pcli_ipam_data_preCBASS)
summary(Pcli_fvfm_GLMM)

Pcli_fvfm_GLMM_simres <- simulateResiduals(Pcli_fvfm_GLMM)
plot(Pcli_fvfm_GLMM_simres)

The variance is significant, so I don’t think you can do a regular anova.

I’m not sure what the assumptions are for the Chi-squared test, but the code for that is:

stats::drop1(Pcli_fvfm_GLMM,test=“Chisq”) #code source: https://cran.r-project.org/web/packages/glmmTMB/vignettes/model_evaluation.pdf

Result: Single term deletions Model: fvfm ~ Treatment + (1 | Colony) + (1 | Tank/Date) Df AIC LRT Pr(>Chi) -2191.1
Treatment 1 -2192.7 0.38531 0.5348

This is all for Pclivosa though so we expect this to not be significant anyways (based on how the boxplot looks).

Ok now let’s try for Acervicornis.

Acer_fvfm_GLMM <- glmmTMB::glmmTMB(fvfm ~ Treatment + (1|Colony) + (1|Tank/Date), family=list(family="beta", link="logit"), data=Acer_ipam_data_preCBASS)

summary(Acer_fvfm_GLMM)

Acer_fvfm_GLMM_simres <- simulateResiduals(Acer_fvfm_GLMM)
plot(Acer_fvfm_GLMM_simres) #things are significant, so I don't think you can do a regular anova

stats::drop1(Acer_fvfm_GLMM,test="Chisq") #not significant

Result: Single term deletions Model: fvfm ~ Treatment + (1 | Colony) + (1 | Tank/Date) Df AIC LRT Pr(>Chi) -1524.8
Treatment 1 -1525.2 1.545 0.2139

My next step I think is to try to normalize the Fv/Fm to the initial value, and then redo these stat tests and see if that helps at all. If not, then none of these are significant. Which is fine.