Coral Size Prediction Model
Daisy Ponce
7/1/2023
Data
Variables
Predictors:
- TREATMENT represents light spectrum
- LED Blue-shifted (BLUE)
- LED reef-depth (NATURAL)
- Sunlight at surface depth (SUNLIGHT)
- SPECIES -Pseudodiploria strigosa (PSTR) -Pseudodiploria clivosa (PCLI)
- AGE represents time (weeks post-settlement)
Response:
GROWTH represents surface area (mm^2)
Random effect variables:
- TANK = The tank in which the tile were placed. (E or F)
- TILE = Tile with 3-8 recruits of either pstr or pcli
- ID = Individual coral recruit on tile
Libraries
library(readxl)
library(readr)
library(dplyr)
library(lme4)
library(car)
library(glmmTMB)
library(DHARMa)
library(emmeans)
library(performance)Import Data
FGS <- read_csv("/Users/Day/Downloads/GrowthFinalData.csv")## Rows: 4365 Columns: 7
## ── Column specification ────────────────────────────────────────────────────────
## Delimiter: ","
## chr (5): ID, TANK, SPECIES, TREATMENT, TILE
## dbl (2): AGE, GROWTH
##
## ℹ Use `spec()` to retrieve the full column specification for this data.
## ℹ Specify the column types or set `show_col_types = FALSE` to quiet this message.#Specifying class of variables
FGS$TREATMENT<- factor(FGS$TREATMENT)
FGS$TANK<- factor(FGS$TANK)
FGS$TILE<- factor(FGS$TILE)
FGS$ID<- factor(FGS$ID)
FGS$SPECIES<- factor(FGS$SPECIES)
#Checking class
class(FGS$TREATMENT)## [1] "factor"class(FGS$TANK)## [1] "factor"class(FGS$TILE)## [1] "factor"class(FGS$ID)## [1] "factor"class(FGS$SPECIES)## [1] "factor"class (FGS$AGE)## [1] "numeric"class (FGS$GROWTH)## [1] "numeric"Descriptive stats
Mean,standard deviation & standard error for all corals at 5 weeks post-settlement
#Mean for all corals at 5 weeks post-settlement
FGS_age_5_14 <- FGS %>% filter(AGE %in% c(5, 14))
mean5 <- FGS_age_5_14 %>% filter( AGE == 5 )
mean(mean5$GROWTH) #0.5279123## [1] 0.5279123sd(mean5$GROWTH) #0.2481307## [1] 0.2481307n <- length(mean5$GROWTH)
standard_error <- 0.2481307/ sqrt(n)
standard_error #0.01133738## [1] 0.01133738Mean, standard deviation & standard error for each species at 5 weeks
P.strigosa
FGS_age_5_14 <- FGS %>% filter(AGE %in% c(5, 14))
meanpstr5 <- FGS_age_5_14 %>% filter( AGE == 5 & SPECIES=="PSTR" )
mean(meanpstr5$GROWTH) #0.5692636## [1] 0.5692636sd(meanpstr5$GROWTH) #0.2237613## [1] 0.2237613n <- length(meanpstr5$GROWTH)
standard_error <- 0.2237613/ sqrt(n)
standard_error #0.01393077## [1] 0.01393077P.clivosa
FGS_age_5_14 <- FGS %>% filter(AGE %in% c(5, 14))
meanpcli5 <- FGS_age_5_14 %>% filter( AGE == 5 & SPECIES=="PCLI" )
mean(meanpcli5$GROWTH) #0.479638## [1] 0.479638sd(meanpcli5$GROWTH) #0.2663386## [1] 0.2663386n <- length(meanpcli5$GROWTH)
standard_error <- 0.2663386/ sqrt(n)
standard_error #0.01791587## [1] 0.01791587Mean, standard deviation & standard error by treatment at 5 weeks
##### by treatment at 5 weeks
blue <- FGS_age_5_14 %>% filter(TREATMENT == "BLUE" & AGE == 5 )
mean(blue$GROWTH)## [1] 0.4538365sd(blue$GROWTH)## [1] 0.210466# > mean(blue$GROWTH)
# [1] 0.4538365
# > sd(blue$GROWTH)
# [1] 0.210466
nat <- FGS_age_5_14 %>% filter(TREATMENT == "NATURAL" & AGE == 5 )
mean(nat$GROWTH)## [1] 0.5906875sd(nat$GROWTH)## [1] 0.2576603mean(nat$GROWTH)## [1] 0.5906875# [1] 0.5906875
# > sd(nat$GROWTH)
# [1] 0.2576603
sun <- FGS_age_5_14 %>% filter(TREATMENT == "SUNLIGHT" & AGE == 5 )
mean(sun$GROWTH)## [1] 0.53875sd(sun$GROWTH)## [1] 0.2551797# > mean(sun$GROWTH)
# [1] 0.53875
# > sd(sun$GROWTH)
# [1] 0.2551797Mean, standard deviation & standard error by treatment at 14 weeks post-settlement
blue14 <- FGS_age_5_14 %>% filter(TREATMENT == "BLUE" & AGE == 14 )
mean(blue14$GROWTH)## [1] 1.129743#1.129743
sd(blue14$GROWTH)## [1] 0.7974191#0.7974191
nat14 <- FGS_age_5_14 %>% filter(TREATMENT == "NATURAL" & AGE == 14 )
mean(nat14$GROWTH)## [1] 1.155159sd(nat14$GROWTH)## [1] 0.8159495# > mean(nat14$GROWTH)
# [1] 1.155159
# > sd(nat14$GROWTH)
# [1] 0.8159495
sun14 <- FGS_age_5_14 %>% filter(TREATMENT == "SUNLIGHT" & AGE == 14 )
mean(sun14$GROWTH)## [1] 0.7660274sd(sun14$GROWTH)## [1] 0.507882# > mean(sun14$GROWTH)
# [1] 0.7660274
# > sd(sun14$GROWTH)
# [1] 0.507882Mean, standard deviation & standard error by treatment and species at 14 weeks post-settlement
blue_df_14c <- FGS_age_5_14 %>% filter(TREATMENT == "BLUE" & AGE == 14 & SPECIES=="PCLI")
sd(blue_df_14c$GROWTH)#0.6831951## [1] 0.6831951mean(blue_df_14c$GROWTH)#0.9024603## [1] 0.9024603summary(blue_df_14c)## ID TANK SPECIES TREATMENT AGE GROWTH
## 3EB2-12_1: 1 E:29 PCLI:63 BLUE :63 Min. :14 Min. :0.1000
## 3EB2-12_2: 1 F:34 PSTR: 0 NATURAL : 0 1st Qu.:14 1st Qu.:0.4050
## 3EB2-12_3: 1 SUNLIGHT: 0 Median :14 Median :0.7300
## 3EB2-5_1 : 1 Mean :14 Mean :0.9025
## 3EB2-5_2 : 1 3rd Qu.:14 3rd Qu.:1.2000
## 3EB2-5_3 : 1 Max. :14 Max. :3.0900
## (Other) :57
## TILE
## 7EB2-2 : 7
## 7FB2-4 : 7
## 6FB2-1 : 6
## 5EB2-9 : 5
## 4FB2-1 : 4
## 4FB2-6 : 4
## (Other):30blue_df_14s <- FGS_age_5_14 %>% filter(TREATMENT == "BLUE" & AGE == 14 & SPECIES=="PSTR")
sd(blue_df_14s$GROWTH) #0.8403171## [1] 0.8403171mean(blue_df_14s$GROWTH)#1.32589## [1] 1.32589summary(blue_df_14s)## ID TANK SPECIES TREATMENT AGE GROWTH
## 3EB1-11_1: 1 E:33 PCLI: 0 BLUE :73 Min. :14 Min. :0.270
## 3EB1-11_2: 1 F:40 PSTR:73 NATURAL : 0 1st Qu.:14 1st Qu.:0.680
## 3EB1-11_3: 1 SUNLIGHT: 0 Median :14 Median :1.120
## 3EB1-17_1: 1 Mean :14 Mean :1.326
## 3EB1-17_2: 1 3rd Qu.:14 3rd Qu.:1.670
## 3EB1-17_3: 1 Max. :14 Max. :3.670
## (Other) :67
## TILE
## 8FB1-6 : 8
## 7FB1-4 : 6
## 5FB1-3 : 5
## 6EB1-5 : 5
## 6FB1-4 : 5
## 7EB1-6 : 5
## (Other):39natural_df_14c <- FGS_age_5_14 %>% filter(TREATMENT == "NATURAL" & AGE == 14 & SPECIES=="PCLI")
sd(natural_df_14c$GROWTH) #0.8303951## [1] 0.8303951mean(natural_df_14c$GROWTH) #1.063019## [1] 1.063019summary(natural_df_14c)## ID TANK SPECIES TREATMENT AGE GROWTH
## 3EN2-2_1 : 1 E:26 PCLI:53 BLUE : 0 Min. :14 Min. :0.120
## 3EN2-2_2 : 1 F:27 PSTR: 0 NATURAL :53 1st Qu.:14 1st Qu.:0.570
## 3EN2-6_1 : 1 SUNLIGHT: 0 Median :14 Median :0.860
## 3EN2-6_2 : 1 Mean :14 Mean :1.063
## 3EN2-6_3 : 1 3rd Qu.:14 3rd Qu.:1.310
## 3FN2-10_1: 1 Max. :14 Max. :3.960
## (Other) :47
## TILE
## 5EN2-7 : 5
## 5FN2-5 : 5
## 4EN2-10: 4
## 4FN2-7 : 4
## 5EN2-4 : 4
## 5FN2-10: 4
## (Other):27natural_df_14s <- FGS_age_5_14 %>% filter(TREATMENT == "NATURAL" & AGE == 14 & SPECIES=="PSTR")
sd(natural_df_14s$GROWTH) # 0.8043989## [1] 0.8043989mean(natural_df_14s$GROWTH) #1.222055## [1] 1.222055summary(natural_df_14s)## ID TANK SPECIES TREATMENT AGE GROWTH
## 3EN1-1_1 : 1 E:40 PCLI: 0 BLUE : 0 Min. :14 Min. :0.200
## 3EN1-1_2 : 1 F:33 PSTR:73 NATURAL :73 1st Qu.:14 1st Qu.:0.680
## 3EN1-1_3 : 1 SUNLIGHT: 0 Median :14 Median :0.960
## 3EN1-18_1: 1 Mean :14 Mean :1.222
## 3EN1-18_2: 1 3rd Qu.:14 3rd Qu.:1.670
## 3EN1-18_3: 1 Max. :14 Max. :4.850
## (Other) :67
## TILE
## 7EN1-3 : 7
## 6EN1-1 : 6
## 6FN1-3 : 6
## 7FN1-1 : 6
## 8EN1-3 : 6
## 5EN1-5 : 5
## (Other):37sunlight_df_14c <- FGS_age_5_14 %>% filter(TREATMENT == "SUNLIGHT" & AGE == 14 & SPECIES=="PCLI")
sd(sunlight_df_14c$GROWTH) #0.3915653## [1] 0.3915653mean(sunlight_df_14c$GROWTH)#0.5643548## [1] 0.5643548summary(sunlight_df_14c)## ID TANK SPECIES TREATMENT AGE GROWTH
## 3ES2-1_1: 1 E:31 PCLI:62 BLUE : 0 Min. :14 Min. :0.1300
## 3ES2-1_2: 1 F:31 PSTR: 0 NATURAL : 0 1st Qu.:14 1st Qu.:0.2650
## 3ES2-1_3: 1 SUNLIGHT:62 Median :14 Median :0.4550
## 3ES2-3_1: 1 Mean :14 Mean :0.5644
## 3ES2-3_2: 1 3rd Qu.:14 3rd Qu.:0.7700
## 3ES2-3_3: 1 Max. :14 Max. :1.9500
## (Other) :56
## TILE
## 7FS2-3 : 7
## 6ES2-4 : 6
## 6FS2-3 : 6
## 5FS2-3 : 5
## 4ES2-2 : 4
## 4ES2-3 : 4
## (Other):30sunlight_df_14s <- FGS_age_5_14 %>% filter(TREATMENT == "SUNLIGHT" & AGE == 14 & SPECIES=="PSTR")
sd(sunlight_df_14s$GROWTH) #0.5339814## [1] 0.5339814mean(sunlight_df_14s$GROWTH)#0.914881## [1] 0.914881summary(sunlight_df_14s)## ID TANK SPECIES TREATMENT AGE GROWTH
## 3ES1-15_1: 1 E:41 PCLI: 0 BLUE : 0 Min. :14 Min. :0.1400
## 3ES1-15_2: 1 F:43 PSTR:84 NATURAL : 0 1st Qu.:14 1st Qu.:0.5250
## 3ES1-15_3: 1 SUNLIGHT:84 Median :14 Median :0.8000
## 3ES1-16_1: 1 Mean :14 Mean :0.9149
## 3ES1-16_2: 1 3rd Qu.:14 3rd Qu.:1.1775
## 3ES1-16_3: 1 Max. :14 Max. :2.9300
## (Other) :78
## TILE
## 8FS1-5 : 8
## 7ES1-5 : 7
## 7FS1-2 : 7
## 8ES1-1 : 7
## 6FS1-2 : 6
## 5ES1-1 : 5
## (Other):44Model Selection
A generalized linear mixed effect model (GLMM) with a Gamma distribution and a log link function was used to compare the growth of juvenile corals (surface area over nine weeks) among different light spectra treatments. A general linear mixed effect model was chosen to account for the effects of random factors, including the random intercept and slope of tile and individual coral, which may cause variation over time. The most complex model was built, including all fixed predictors (light spectrum, species, and time) and their interactions. A backward stepwise method was used to remove non-significant and multicollinear variables and interactions from the model.
# MOST COMPLEX MODEL: Interaction between fixed predictors and all random effects included
mod0 <- glmmTMB (GROWTH~TREATMENT * SPECIES * AGE + (1|TANK)+(1|TILE)+(1|ID), family = Gamma(log),
data = FGS)## Warning in Matrix::sparseMatrix(dims = c(0, 0), i = integer(0), j =
## integer(0), : 'giveCsparse' is deprecated; setting repr="T" for you
## Warning in Matrix::sparseMatrix(dims = c(0, 0), i = integer(0), j =
## integer(0), : 'giveCsparse' is deprecated; setting repr="T" for you
## Warning in Matrix::sparseMatrix(dims = c(0, 0), i = integer(0), j =
## integer(0), : 'giveCsparse' is deprecated; setting repr="T" for yousummary(mod0)## Family: Gamma ( log )
## Formula:
## GROWTH ~ TREATMENT * SPECIES * AGE + (1 | TANK) + (1 | TILE) + (1 | ID)
## Data: FGS
##
## AIC BIC logLik deviance df.resid
## -594.5 -492.4 313.3 -626.5 4349
##
## Random effects:
##
## Conditional model:
## Groups Name Variance Std.Dev.
## TANK (Intercept) 9.589e-11 9.792e-06
## TILE (Intercept) 6.598e-02 2.569e-01
## ID (Intercept) 1.858e-01 4.310e-01
## Number of obs: 4365, groups: TANK, 2; TILE, 102; ID, 481
##
## Dispersion estimate for Gamma family (sigma^2): 0.0945
##
## Conditional model:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -1.044844 0.091723 -11.391 < 2e-16 ***
## TREATMENTNATURAL 0.024990 0.130114 0.192 0.847692
## TREATMENTSUNLIGHT 0.450939 0.129904 3.471 0.000518 ***
## SPECIESPSTR -0.114103 0.125655 -0.908 0.363842
## AGE 0.032639 0.004163 7.841 4.47e-15 ***
## TREATMENTNATURAL:SPECIESPSTR 0.536080 0.178182 3.009 0.002625 **
## TREATMENTSUNLIGHT:SPECIESPSTR 0.209942 0.178083 1.179 0.238440
## TREATMENTNATURAL:AGE 0.008111 0.006111 1.327 0.184459
## TREATMENTSUNLIGHT:AGE -0.049228 0.005933 -8.297 < 2e-16 ***
## SPECIESPSTR:AGE 0.047246 0.005673 8.328 < 2e-16 ***
## TREATMENTNATURAL:SPECIESPSTR:AGE -0.047284 0.008234 -5.742 9.34e-09 ***
## TREATMENTSUNLIGHT:SPECIESPSTR:AGE 0.001765 0.008059 0.219 0.826603
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1glmmTMB:::Anova.glmmTMB(mod0)## Analysis of Deviance Table (Type II Wald chisquare tests)
##
## Response: GROWTH
## Chisq Df Pr(>Chisq)
## TREATMENT 3.8559 2 0.1454
## SPECIES 42.4525 1 7.242e-11 ***
## AGE 476.7247 1 < 2.2e-16 ***
## TREATMENT:SPECIES 2.0223 2 0.3638
## TREATMENT:AGE 146.3545 2 < 2.2e-16 ***
## SPECIES:AGE 97.8936 1 < 2.2e-16 ***
## TREATMENT:SPECIES:AGE 44.7653 2 1.903e-10 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1r2(mod0)## Warning: Can't compute random effect variances. Some variance components equal
## zero. Your model may suffer from singularity (see `?lme4::isSingular`
## and `?performance::check_singularity`).
## Solution: Respecify random structure! You may also decrease the
## `tolerance` level to enforce the calculation of random effect variances.## Random effect variances not available. Returned R2 does not account for random effects.## # R2 for Mixed Models
##
## Conditional R2: NA
## Marginal R2: 0.433Multicollinearity was observed between the variables “TANK” and “ID.” To address this, the random effect variable with the lowest variance should be removed. Therefore, a model was constructed by including “TANK” as a fixed predictor to evaluate its significance and determine if there is a tank effect or if it can be eliminated from the model.
#Drop tank (singularity issue with id) check for tank effect by adding it as fixed
mod0.1 <-glmmTMB (GROWTH~TREATMENT * SPECIES * AGE +TANK+(1|TILE)+(1|ID), family = Gamma(log),
data = FGS)## Warning in Matrix::sparseMatrix(dims = c(0, 0), i = integer(0), j =
## integer(0), : 'giveCsparse' is deprecated; setting repr="T" for you
## Warning in Matrix::sparseMatrix(dims = c(0, 0), i = integer(0), j =
## integer(0), : 'giveCsparse' is deprecated; setting repr="T" for you
## Warning in Matrix::sparseMatrix(dims = c(0, 0), i = integer(0), j =
## integer(0), : 'giveCsparse' is deprecated; setting repr="T" for yousummary(mod0.1)## Family: Gamma ( log )
## Formula:
## GROWTH ~ TREATMENT * SPECIES * AGE + TANK + (1 | TILE) + (1 | ID)
## Data: FGS
##
## AIC BIC logLik deviance df.resid
## -594.5 -492.4 313.3 -626.5 4349
##
## Random effects:
##
## Conditional model:
## Groups Name Variance Std.Dev.
## TILE (Intercept) 0.06605 0.257
## ID (Intercept) 0.18574 0.431
## Number of obs: 4365, groups: TILE, 102; ID, 481
##
## Dispersion estimate for Gamma family (sigma^2): 0.0945
##
## Conditional model:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -1.041467 0.097537 -10.678 < 2e-16 ***
## TREATMENTNATURAL 0.025012 0.130141 0.192 0.847591
## TREATMENTSUNLIGHT 0.450949 0.129931 3.471 0.000519 ***
## SPECIESPSTR -0.114113 0.125681 -0.908 0.363900
## AGE 0.032638 0.004163 7.841 4.47e-15 ***
## TANKF -0.006734 0.066071 -0.102 0.918815
## TREATMENTNATURAL:SPECIESPSTR 0.536080 0.178219 3.008 0.002630 **
## TREATMENTSUNLIGHT:SPECIESPSTR 0.209930 0.178120 1.179 0.238563
## TREATMENTNATURAL:AGE 0.008110 0.006111 1.327 0.184530
## TREATMENTSUNLIGHT:AGE -0.049229 0.005933 -8.297 < 2e-16 ***
## SPECIESPSTR:AGE 0.047248 0.005673 8.328 < 2e-16 ***
## TREATMENTNATURAL:SPECIESPSTR:AGE -0.047289 0.008235 -5.743 9.31e-09 ***
## TREATMENTSUNLIGHT:SPECIESPSTR:AGE 0.001765 0.008059 0.219 0.826672
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1glmmTMB:::Anova.glmmTMB(mod0.1)## Analysis of Deviance Table (Type II Wald chisquare tests)
##
## Response: GROWTH
## Chisq Df Pr(>Chisq)
## TREATMENT 3.8536 2 0.1456
## SPECIES 42.4302 1 7.325e-11 ***
## AGE 476.7045 1 < 2.2e-16 ***
## TANK 0.0104 1 0.9188
## TREATMENT:SPECIES 2.0211 2 0.3640
## TREATMENT:AGE 146.3558 2 < 2.2e-16 ***
## SPECIES:AGE 97.8942 1 < 2.2e-16 ***
## TREATMENT:SPECIES:AGE 44.7722 2 1.896e-10 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1#drop tank (even as fixed pred not sig)###Fitting Random effects as intercept, slopes and both.
#Random intercept
mod1 <- glmmTMB (GROWTH~TREATMENT * SPECIES * AGE + (1|TILE/ID), family = Gamma(log),
data = FGS)## Warning in Matrix::sparseMatrix(dims = c(0, 0), i = integer(0), j =
## integer(0), : 'giveCsparse' is deprecated; setting repr="T" for you
## Warning in Matrix::sparseMatrix(dims = c(0, 0), i = integer(0), j =
## integer(0), : 'giveCsparse' is deprecated; setting repr="T" for you
## Warning in Matrix::sparseMatrix(dims = c(0, 0), i = integer(0), j =
## integer(0), : 'giveCsparse' is deprecated; setting repr="T" for you#Random Slope
mod2 <- glmmTMB (GROWTH~TREATMENT * SPECIES * AGE + (AGE|TILE/ID), family = Gamma(log),
data = FGS)## Warning in Matrix::sparseMatrix(dims = c(0, 0), i = integer(0), j =
## integer(0), : 'giveCsparse' is deprecated; setting repr="T" for you
## Warning in Matrix::sparseMatrix(dims = c(0, 0), i = integer(0), j =
## integer(0), : 'giveCsparse' is deprecated; setting repr="T" for you
## Warning in Matrix::sparseMatrix(dims = c(0, 0), i = integer(0), j =
## integer(0), : 'giveCsparse' is deprecated; setting repr="T" for you#Random slope AND intercept
mod3 <- glmmTMB(GROWTH ~ TREATMENT * SPECIES * AGE + (1 + AGE|TILE/ID), family = Gamma(log), data = FGS)## Warning in Matrix::sparseMatrix(dims = c(0, 0), i = integer(0), j =
## integer(0), : 'giveCsparse' is deprecated; setting repr="T" for you
## Warning in Matrix::sparseMatrix(dims = c(0, 0), i = integer(0), j =
## integer(0), : 'giveCsparse' is deprecated; setting repr="T" for you
## Warning in Matrix::sparseMatrix(dims = c(0, 0), i = integer(0), j =
## integer(0), : 'giveCsparse' is deprecated; setting repr="T" for youanova(mod1,mod2,mod3)## Data: FGS
## Models:
## mod1: GROWTH ~ TREATMENT * SPECIES * AGE + (1 | TILE/ID), zi=~0, disp=~1
## mod2: GROWTH ~ TREATMENT * SPECIES * AGE + (AGE | TILE/ID), zi=~0, disp=~1
## mod3: GROWTH ~ TREATMENT * SPECIES * AGE + (1 + AGE | TILE/ID), zi=~0, disp=~1
## Df AIC BIC logLik deviance Chisq Chi Df Pr(>Chisq)
## mod1 15 -596.52 -500.8 313.26 -626.52
## mod2 19 -1256.65 -1135.4 647.32 -1294.65 668.13 4 <2e-16 ***
## mod3 19 -1256.65 -1135.4 647.32 -1294.65 0.00 0 1
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1#Mod2 in better than mod1 but mod2 and mod3 are the same... mod3 is selected because it make the most sense based on the experimental set up (individual across treatments start at different sizes and do not have the same slope(surface area/week).
glmmTMB:::Anova.glmmTMB(mod3) #drop 3 way interaction## Analysis of Deviance Table (Type II Wald chisquare tests)
##
## Response: GROWTH
## Chisq Df Pr(>Chisq)
## TREATMENT 15.6946 2 0.0003908 ***
## SPECIES 29.0037 1 7.224e-08 ***
## AGE 60.5388 1 7.214e-15 ***
## TREATMENT:SPECIES 3.9541 2 0.1384741
## TREATMENT:AGE 17.2565 2 0.0001790 ***
## SPECIES:AGE 18.6106 1 1.603e-05 ***
## TREATMENT:SPECIES:AGE 5.1500 2 0.0761554 .
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1r2(mod3) #Conditional R2: 0.834 , Marginal R2: 0.168## # R2 for Mixed Models
##
## Conditional R2: 0.834
## Marginal R2: 0.168#Dropped 3way interaction
mod4 <- glmmTMB(GROWTH ~ TREATMENT + SPECIES + AGE+TREATMENT:SPECIES+
TREATMENT:AGE + SPECIES:AGE + (1 + AGE|TILE/ID), family = Gamma(log), data = FGS)## Warning in Matrix::sparseMatrix(dims = c(0, 0), i = integer(0), j =
## integer(0), : 'giveCsparse' is deprecated; setting repr="T" for you
## Warning in Matrix::sparseMatrix(dims = c(0, 0), i = integer(0), j =
## integer(0), : 'giveCsparse' is deprecated; setting repr="T" for you
## Warning in Matrix::sparseMatrix(dims = c(0, 0), i = integer(0), j =
## integer(0), : 'giveCsparse' is deprecated; setting repr="T" for youglmmTMB:::Anova.glmmTMB(mod4) #drop TREATMENT:SPECIES## Analysis of Deviance Table (Type II Wald chisquare tests)
##
## Response: GROWTH
## Chisq Df Pr(>Chisq)
## TREATMENT 15.7436 2 0.0003814 ***
## SPECIES 28.9235 1 7.530e-08 ***
## AGE 57.0699 1 4.206e-14 ***
## TREATMENT:SPECIES 3.8801 2 0.1436986
## TREATMENT:AGE 16.2303 2 0.0002990 ***
## SPECIES:AGE 17.5386 1 2.815e-05 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1r2(mod3) #Conditional R2: 0.835,Marginal R2: 0.168## # R2 for Mixed Models
##
## Conditional R2: 0.834
## Marginal R2: 0.168#Dropped TREATMENT:SPECIES
mod5 <- glmmTMB(GROWTH ~ TREATMENT + SPECIES + AGE+
TREATMENT:AGE + SPECIES:AGE + (1 + AGE|TILE/ID), family = Gamma(log), data = FGS)## Warning in Matrix::sparseMatrix(dims = c(0, 0), i = integer(0), j =
## integer(0), : 'giveCsparse' is deprecated; setting repr="T" for you
## Warning in Matrix::sparseMatrix(dims = c(0, 0), i = integer(0), j =
## integer(0), : 'giveCsparse' is deprecated; setting repr="T" for you
## Warning in Matrix::sparseMatrix(dims = c(0, 0), i = integer(0), j =
## integer(0), : 'giveCsparse' is deprecated; setting repr="T" for yousummary(mod5)## Family: Gamma ( log )
## Formula:
## GROWTH ~ TREATMENT + SPECIES + AGE + TREATMENT:AGE + SPECIES:AGE +
## (1 + AGE | TILE/ID)
## Data: FGS
##
## AIC BIC logLik deviance df.resid
## -1255.9 -1160.1 642.9 -1285.9 4350
##
## Random effects:
##
## Conditional model:
## Groups Name Variance Std.Dev. Corr
## ID:TILE (Intercept) 0.109695 0.33120
## AGE 0.001639 0.04049 -0.21
## TILE (Intercept) 0.106284 0.32601
## AGE 0.001188 0.03447 -0.67
## Number of obs: 4365, groups: ID:TILE, 481; TILE, 102
##
## Dispersion estimate for Gamma family (sigma^2): 0.0721
##
## Conditional model:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -1.138298 0.079092 -14.392 < 2e-16 ***
## TREATMENTNATURAL 0.306580 0.095560 3.208 0.001336 **
## TREATMENTSUNLIGHT 0.508687 0.095085 5.350 8.80e-08 ***
## SPECIESPSTR 0.121999 0.078049 1.563 0.118030
## AGE 0.033332 0.008650 3.853 0.000116 ***
## TREATMENTNATURAL:AGE -0.018491 0.010507 -1.760 0.078438 .
## TREATMENTSUNLIGHT:AGE -0.041615 0.010341 -4.024 5.71e-05 ***
## SPECIESPSTR:AGE 0.035809 0.008547 4.190 2.79e-05 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1glmmTMB:::Anova.glmmTMB(mod5) #all predictors & interactions are significant!!!## Analysis of Deviance Table (Type II Wald chisquare tests)
##
## Response: GROWTH
## Chisq Df Pr(>Chisq)
## TREATMENT 14.230 2 0.0008128 ***
## SPECIES 29.052 1 7.045e-08 ***
## AGE 57.277 1 3.786e-14 ***
## TREATMENT:AGE 16.271 2 0.0002930 ***
## SPECIES:AGE 17.554 1 2.793e-05 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1r2(mod4) #Conditional R2: 0.834, Marginal R2: 0.161## # R2 for Mixed Models
##
## Conditional R2: 0.835
## Marginal R2: 0.168 #The conditional R2 of 0.834 indicates that 83.4% of the total variance in the dependent variable is explained by the predictors and random effects, while the marginal R2 of 0.161 suggests that the fixed predictors alone account for 16.1% of the varianceConfirming selected model is the best
#Lowest AIC value is selected
AIC(mod1,mod2,mod3,mod4,mod5)## df AIC
## mod1 15 -596.5215
## mod2 19 -1256.6486
## mod3 19 -1256.6486
## mod4 17 -1255.6592
## mod5 15 -1255.8702Residuals
library(DHARMa)
simulationOutput <- simulateResiduals(mod5, n = 1000)## Warning in Matrix::sparseMatrix(dims = c(0, 0), i = integer(0), j =
## integer(0), : 'giveCsparse' is deprecated; setting repr="T" for you
## Warning in Matrix::sparseMatrix(dims = c(0, 0), i = integer(0), j =
## integer(0), : 'giveCsparse' is deprecated; setting repr="T" for you
## Warning in Matrix::sparseMatrix(dims = c(0, 0), i = integer(0), j =
## integer(0), : 'giveCsparse' is deprecated; setting repr="T" for you
## Warning in Matrix::sparseMatrix(dims = c(0, 0), i = integer(0), j =
## integer(0), : 'giveCsparse' is deprecated; setting repr="T" for you
## Warning in Matrix::sparseMatrix(dims = c(0, 0), i = integer(0), j =
## integer(0), : 'giveCsparse' is deprecated; setting repr="T" for you
## Warning in Matrix::sparseMatrix(dims = c(0, 0), i = integer(0), j =
## integer(0), : 'giveCsparse' is deprecated; setting repr="T" for youplot(simulationOutput)
plotResiduals(simulationOutput)
plotResiduals(simulationOutput, FGS$TREATMENT)
#Leven test for homogeniety of variance suggests blue and natural have sig
#within-group deviantions from uniformity.
#If the deviation from uniformity is small and the sample sizes are large,
#it may not have a significant impact on the results.Pairwise Comparison
emtrends(mod5, pairwise ~ TREATMENT:SPECIES, var = "AGE")## $emtrends
## TREATMENT SPECIES AGE.trend SE df lower.CL upper.CL
## BLUE PCLI 0.03333 0.00865 4350 0.01637 0.05029
## NATURAL PCLI 0.01484 0.00881 4350 -0.00242 0.03211
## SUNLIGHT PCLI -0.00828 0.00859 4350 -0.02513 0.00856
## BLUE PSTR 0.06914 0.00837 4350 0.05274 0.08554
## NATURAL PSTR 0.05065 0.00849 4350 0.03400 0.06730
## SUNLIGHT PSTR 0.02753 0.00829 4350 0.01127 0.04378
##
## Confidence level used: 0.95
##
## $contrasts
## contrast estimate SE df t.ratio p.value
## BLUE PCLI - NATURAL PCLI 0.01849 0.01051 4350 1.760 0.4921
## BLUE PCLI - SUNLIGHT PCLI 0.04162 0.01034 4350 4.024 0.0008
## BLUE PCLI - BLUE PSTR -0.03581 0.00855 4350 -4.190 0.0004
## BLUE PCLI - NATURAL PSTR -0.01732 0.01352 4350 -1.281 0.7958
## BLUE PCLI - SUNLIGHT PSTR 0.00581 0.01341 4350 0.433 0.9981
## NATURAL PCLI - SUNLIGHT PCLI 0.02312 0.01046 4350 2.211 0.2328
## NATURAL PCLI - BLUE PSTR -0.05430 0.01357 4350 -4.003 0.0009
## NATURAL PCLI - NATURAL PSTR -0.03581 0.00855 4350 -4.190 0.0004
## NATURAL PCLI - SUNLIGHT PSTR -0.01268 0.01352 4350 -0.938 0.9366
## SUNLIGHT PCLI - BLUE PSTR -0.07742 0.01343 4350 -5.767 <.0001
## SUNLIGHT PCLI - NATURAL PSTR -0.05893 0.01350 4350 -4.367 0.0002
## SUNLIGHT PCLI - SUNLIGHT PSTR -0.03581 0.00855 4350 -4.190 0.0004
## BLUE PSTR - NATURAL PSTR 0.01849 0.01051 4350 1.760 0.4921
## BLUE PSTR - SUNLIGHT PSTR 0.04162 0.01034 4350 4.024 0.0008
## NATURAL PSTR - SUNLIGHT PSTR 0.02312 0.01046 4350 2.211 0.2328
##
## P value adjustment: tukey method for comparing a family of 6 estimates#Plotting the estimated marginal means
emmip(mod5, TREATMENT:SPECIES ~ AGE, cov.reduce = range)
This output is the result of running a pairwise comparison using emmeans R package of the effects of three different treatments (LED blue-shifted spectrum, LED reef-depth spectrum, and sunlight at surface spectrum) on two different species (P. clivosa and P. strigosa). Contrast that shows significant differences between the treatment and species pairing.
Predictive Size Model
predictions<-(predict(mod5,newdata=FGS,type='response',re.form=NA)) #not including re
all.data <- cbind(FGS, predictions)
par(mfrow=c(1,1))
#color scheme
col_scheme <- c("BLUE" = "#1811ac", "NATURAL" = "#bdd07d", "SUNLIGHT" = "#f7b58b")
# Set the x-axis limits and tick marks
xlimits <- c(5, 14)
xtickmarks <- seq(xlimits[1], xlimits[2], by = 1)
#size
par(mar = c(5, 5, 4, 2))
# Plot the data with custom tick marks
plot(GROWTH ~ AGE, data = FGS, pch = 20,lwd = 4, xlab = 'Time (weeks post-settlement)', ylab = expression(paste("Surface Area (mm"^2*")")), col = factor(SPECIES), type = "n", ylim = c(0, 1), xlim = xlimits)
axis(side = 1, at = xtickmarks, labels = xtickmarks)
# Add the lines and legends
sub1 <- subset(all.data, all.data$SPECIES == "PCLI" & all.data$TREATMENT == "BLUE")
sub2 <- subset(all.data, all.data$SPECIES == "PCLI" & all.data$TREATMENT == "NATURAL")
sub3 <- subset(all.data, all.data$SPECIES == "PCLI" & all.data$TREATMENT == "SUNLIGHT")
sub4 <- subset(all.data, all.data$SPECIES == "PSTR" & all.data$TREATMENT == "BLUE")
sub5 <- subset(all.data, all.data$SPECIES == "PSTR" & all.data$TREATMENT == "NATURAL")
sub6 <- subset(all.data, all.data$SPECIES == "PSTR" & all.data$TREATMENT == "SUNLIGHT")
sub1 <- sub1[order(sub1$predictions), ]
sub2 <- sub2[order(sub2$predictions), ]
sub3 <- sub3[order(sub3$predictions), ]
sub4 <- sub4[order(sub4$predictions), ]
sub5 <- sub5[order(sub5$predictions), ]
sub6 <- sub6[order(sub6$predictions), ]
points(sub1$predictions ~ sub1$AGE, type = 'l', col = col_scheme["BLUE"], lty = 2)
points(sub2$predictions ~ sub2$AGE, type = 'l', col = col_scheme["NATURAL"], lty = 2)
points(sub3$predictions ~ sub3$AGE, type = 'l', col = col_scheme["SUNLIGHT"], lty = 2)
points(sub4$predictions ~ sub4$AGE, type = 'l', col = col_scheme["BLUE"], lty = 1)
points(sub5$predictions ~ sub5$AGE, type = 'l', col = col_scheme["NATURAL"], lty = 1)
points(sub6$predictions ~ sub6$AGE, type = 'l', col = col_scheme["SUNLIGHT"], lty = 1)
legend("bottomleft", legend = c('Pseudodiploria strigosa ', 'Pseudodiploria clivosa'), lty = c(1, 2), col = c("black", "black"))
legend("bottomright", legend = c('LED Blue-shifted spectrum', 'LED Reef-depth spectrum', 'Sunlight at surface depth'), lty = 1, col = col_scheme)
###Growth rate Growth rate, as measured by the slope of surface area over time, is a key indicator of coral growth, with a steeper slope indicating faster growth. Code below returns the slope of each line:
# Calculate the slopes for each line
slope_sub1 <- lm(predictions ~ AGE, data = sub1)$coefficients[2]
slope_sub2 <- lm(predictions ~ AGE, data = sub2)$coefficients[2]
slope_sub3 <- lm(predictions ~ AGE, data = sub3)$coefficients[2]
slope_sub4 <- lm(predictions ~ AGE, data = sub4)$coefficients[2]
slope_sub5 <- lm(predictions ~ AGE, data = sub5)$coefficients[2]
slope_sub6 <- lm(predictions ~ AGE, data = sub6)$coefficients[2]
# Print the slopes
print(slope_sub1)# 0.01467065 BLUE PCLI## AGE
## 0.01467059print(slope_sub2)# 0.00742943 NATURAL PCLI## AGE
## 0.007429386print(slope_sub3)# -0.004082206 SUNLIGHT PCLI## AGE
## -0.004082271print(slope_sub4)# 0.04863315 BLUE PSTR## AGE
## 0.04863306print(slope_sub5)# 0.04043051 NATURAL PSTR## AGE
## 0.04043045print(slope_sub6)# 0.02155518 SUNLIGHT PSTR## AGE
## 0.02155508Checking predictions with raw data
Dots represent the observations and the line represents the best fit model.
par(mfrow = c(2, 3)) # Set the layout to 2 rows and 3 columns
par(mar = c(2, 2, 0.9, 1)) # Adjust the margin to make the graphs taller
# Function to create plots with consistent y scale
create_plot <- function(sub_data, col_scheme, line_type) {
plot(GROWTH ~ AGE, data = sub_data, pch = 16, lwd = 1,
xlab = 'Time (weeks post-settlement)',
ylab = expression(paste("Surface Area (mm"^2*")")),
col = col_scheme, main = "", # Empty main title
ylim = c(0, 5)) # Set the y scale to 0-5
lines(predictions ~ AGE, data = sub_data, col = "black", lty = line_type)
}
# Create plots for each subset with the modified y scale
create_plot(sub3, col_scheme["SUNLIGHT"], 2) # Dashed line (lty = 2)
create_plot(sub2, col_scheme["NATURAL"], 2) # Dashed line (lty = 2)
create_plot(sub1, col_scheme["BLUE"], 2) # Dashed line (lty = 2)
create_plot(sub6, col_scheme["SUNLIGHT"], 1) # Solid line (lty = 1)
create_plot(sub5, col_scheme["NATURAL"], 1) # Solid line (lty = 1)
create_plot(sub4, col_scheme["BLUE"], 1) # Solid line (lty = 1)
Axis label created separately. x axis= Surface area(mm^2) ; y axis= Time (weeks post-settlement)
Session Info
sessionInfo()## R version 4.2.2 (2022-10-31)
## Platform: x86_64-apple-darwin17.0 (64-bit)
## Running under: macOS Mojave 10.14.1
##
## Matrix products: default
## BLAS: /Library/Frameworks/R.framework/Versions/4.2/Resources/lib/libRblas.0.dylib
## LAPACK: /Library/Frameworks/R.framework/Versions/4.2/Resources/lib/libRlapack.dylib
##
## locale:
## [1] en_US.UTF-8/en_US.UTF-8/en_US.UTF-8/C/en_US.UTF-8/en_US.UTF-8
##
## attached base packages:
## [1] stats graphics grDevices utils datasets methods base
##
## other attached packages:
## [1] performance_0.10.2 emmeans_1.8.4-1 DHARMa_0.4.6 glmmTMB_1.0.2.9000
## [5] car_3.1-1 carData_3.0-5 lme4_1.1-31 Matrix_1.5-3
## [9] dplyr_1.1.2 readr_2.1.3 readxl_1.4.1
##
## loaded via a namespace (and not attached):
## [1] Rcpp_1.0.9 mvtnorm_1.1-3 lattice_0.20-45 zoo_1.8-11
## [5] digest_0.6.31 utf8_1.2.2 R6_2.5.1 cellranger_1.1.0
## [9] gap.datasets_0.0.5 evaluate_0.20 ggplot2_3.4.0 highr_0.10
## [13] pillar_1.9.0 rlang_1.1.1 multcomp_1.4-23 rstudioapi_0.14
## [17] minqa_1.2.5 nloptr_2.0.3 jquerylib_0.1.4 gap_1.5
## [21] rmarkdown_2.20 labeling_0.4.2 splines_4.2.2 stringr_1.5.0
## [25] TMB_1.9.1 munsell_0.5.0 bit_4.0.5 compiler_4.2.2
## [29] xfun_0.36 pkgconfig_2.0.3 htmltools_0.5.4 insight_0.18.8
## [33] tidyselect_1.2.0 tibble_3.2.1 codetools_0.2-18 fansi_1.0.3
## [37] withr_2.5.0 crayon_1.5.2 tzdb_0.3.0 MASS_7.3-58.2
## [41] grid_4.2.2 gtable_0.3.1 nlme_3.1-162 jsonlite_1.8.4
## [45] xtable_1.8-4 lifecycle_1.0.3 magrittr_2.0.3 scales_1.2.1
## [49] estimability_1.4.1 cli_3.6.0 stringi_1.7.12 vroom_1.6.0
## [53] cachem_1.0.6 farver_2.1.1 bslib_0.4.2 ellipsis_0.3.2
## [57] generics_0.1.3 vctrs_0.6.2 boot_1.3-28 sandwich_3.0-2
## [61] TH.data_1.1-1 tools_4.2.2 bit64_4.0.5 glue_1.6.2
## [65] hms_1.1.2 parallel_4.2.2 abind_1.4-5 fastmap_1.1.0
## [69] survival_3.5-5 yaml_2.3.6 colorspace_2.0-3 knitr_1.41
## [73] sass_0.4.4