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

```{r, echo = F, warning = F, include = T, fig.cap = “Fv/Fm values of corals within the control or variable treatment during the one-month treatment period (March 23-April 20).”} #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”))

```{r}
#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)