Coral Survival Analysis
Daisy Ponce
6/30/2023
To investigate the effect of light spectra and species (fixed factors) on survival over time, a survival analysis, specifically a Cox proportional hazard model, was used. The Kaplan-Meier estimator was utilized to generate survival curves and visualize mortality among treatment groups.
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:
STATUS (0=alive or 1=dead)
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(survival)Import data:
surv <- read_csv("/Users/Day/Downloads/SurvivalFinalData.csv")## New names:
## Rows: 480 Columns: 22
## ── Column specification
## ──────────────────────────────────────────────────────── Delimiter: "," chr
## (6): ID, TANK, SPECIES, TREATMENT, TILE, ST dbl (2): AGE, STATUS lgl (14):
## ...9, ...10, ...11, ...12, ...13, ...14, ...15, ...16, ...17, ...1...
## ℹ 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.
## • `` -> `...9`
## • `` -> `...10`
## • `` -> `...11`
## • `` -> `...12`
## • `` -> `...13`
## • `` -> `...14`
## • `` -> `...15`
## • `` -> `...16`
## • `` -> `...17`
## • `` -> `...18`
## • `` -> `...19`
## • `` -> `...20`
## • `` -> `...21`
## • `` -> `...22`#Define factors
View(surv)
#predictors
surv$TREATMENT<- factor(surv$TREATMENT)
surv$SPECIES<- factor(surv$SPECIES)
#view object type
class(surv$TREATMENT)## [1] "factor"class(surv$SPECIES)## [1] "factor"class(surv$AGE)## [1] "numeric"class(surv$STATUS)## [1] "numeric"attach (surv)COX MODEL
#creates a matrix X that combines the TREATMENT and SPECIES variables into a
#single matrix for use in a Cox proportional hazards model
X<-cbind(TREATMENT,SPECIES)
# right censored:
cox<- coxph (Surv(AGE, STATUS)~ X, method='breslow')Model Summary
summary(cox)## Call:
## coxph(formula = Surv(AGE, STATUS) ~ X, method = "breslow")
##
## n= 480, number of events= 73
##
## coef exp(coef) se(coef) z Pr(>|z|)
## XTREATMENT -0.2463 0.7817 0.1433 -1.719 0.08566 .
## XSPECIES -0.6622 0.5157 0.2407 -2.751 0.00594 **
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## exp(coef) exp(-coef) lower .95 upper .95
## XTREATMENT 0.7817 1.279 0.5902 1.0352
## XSPECIES 0.5157 1.939 0.3217 0.8266
##
## Concordance= 0.615 (se = 0.03 )
## Likelihood ratio test= 10.87 on 2 df, p=0.004
## Wald test = 10.58 on 2 df, p=0.005
## Score (logrank) test = 10.9 on 2 df, p=0.004#Mortality Counts per Week
event_counts <- table(cox$y)
print(event_counts)##
## 6 7 8 9 10 11 12 13 14 14+
## 10 17 11 9 4 6 3 6 7 407At 95% confidence level, survival among light spectra treatments could not be considered statistically significant, but since the same exact ranking of light spectra was observed in both species and the p-value for this test (p-value = 0.086) was very close to rejection threshold (0.05) and would be considered significant if a 90% confidence level had been applied.
KAPLAN-MEIER SURVIVAL CURVE (all species)
Obtaining Percent Survival For Each Curve
# Get the survival probabilities at specific time points (14 weeks/end of experiment)
time_points <- c(14) # Example time points, modify as needed
survival_summary <- summary(km, times = time_points)
# Print the survival probabilities for each group at the specified time points
print(survival_summary)## Call: survfit(formula = Surv(AGE, STATUS) ~ TREATMENT + SPECIES)
##
## TREATMENT=BLUE, SPECIES=PCLI
## time n.risk n.event survival std.err lower 95% CI
## 14.0000 62.0000 12.0000 0.8378 0.0428 0.7579
## upper 95% CI
## 0.9262
##
## TREATMENT=BLUE, SPECIES=PSTR
## time n.risk n.event survival std.err lower 95% CI
## 14.0000 73.0000 13.0000 0.8488 0.0386 0.7764
## upper 95% CI
## 0.9280
##
## TREATMENT=NATURAL, SPECIES=PCLI
## time n.risk n.event survival std.err lower 95% CI
## 14.0000 54.0000 21.0000 0.7162 0.0524 0.6205
## upper 95% CI
## 0.8267
##
## TREATMENT=NATURAL, SPECIES=PSTR
## time n.risk n.event survival std.err lower 95% CI
## 14.0000 74.0000 13.0000 0.8488 0.0386 0.7764
## upper 95% CI
## 0.9280
##
## TREATMENT=SUNLIGHT, SPECIES=PCLI
## time n.risk n.event survival std.err lower 95% CI
## 14.0000 67.0000 12.0000 0.8378 0.0428 0.7579
## upper 95% CI
## 0.9262
##
## TREATMENT=SUNLIGHT, SPECIES=PSTR
## time n.risk n.event survival std.err lower 95% CI
## 14.0000 84.0000 2.0000 0.9767 0.0163 0.9454
## upper 95% CI
## 1.0000Plot all survival curves
plot(km, lwd = 2.5, xlab = 'Age (weeks post-settlement)', ylab = 'Probability of survival',
conf.int = FALSE, lty = c(2, 1, 2, 1, 2, 1),
col = c('#1811ac', '#1811ac', '#bdd07d', '#bdd07d', '#f7b58b', '#f7b58b'),
main= "",
ylim = c(0.5, 1),
xlim = c(4, max(km$time)))
legend('bottomleft',
c('Blue PCLI', 'Blue PSTR', 'Natural PCLI', 'Natural PSTR', 'Sunlight PCLI', 'Sunlight PSTR'),
lty = c(2, 1, 2, 1, 2, 1),
col =c('#1811ac', '#1811ac', '#bdd07d', '#bdd07d', '#f7b58b', '#f7b58b'),
bty = 'n',
ncol = 2, cex = 0.5)
detach(surv)Read me:
The models and figures presented below were created to generate separate survival curves in order to enhance the clarity and visual appeal of the figures. However, it is important to note that variables demonstrating significance may not necessarily indicate influence due to the reduction in sample size.
P.strigosa Cox Model and Survival Curves
#CSV contains Pstr data only
#library(readxl)
surv2 <- read_csv("/Users/Day/Downloads/SurvPstr.csv") ## New names:
## Rows: 258 Columns: 22
## ── Column specification
## ──────────────────────────────────────────────────────── Delimiter: "," chr
## (6): ID, TANK, SPECIES, TREATMENT, TILE, ST dbl (2): AGE, STATUS lgl (14):
## ...9, ...10, ...11, ...12, ...13, ...14, ...15, ...16, ...17, ...1...
## ℹ 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.
## • `` -> `...9`
## • `` -> `...10`
## • `` -> `...11`
## • `` -> `...12`
## • `` -> `...13`
## • `` -> `...14`
## • `` -> `...15`
## • `` -> `...16`
## • `` -> `...17`
## • `` -> `...18`
## • `` -> `...19`
## • `` -> `...20`
## • `` -> `...21`
## • `` -> `...22`#Define factors
View(surv2)
#predictors
surv2$TREATMENT<- factor(surv2$TREATMENT)
surv2$SPECIES<- factor(surv2$SPECIES)
surv2$ST<- factor(surv2$ST)
#view object type
class(surv2$TREATMENT)## [1] "factor"class(surv2$SPECIES)## [1] "factor"class(surv2$AGE)## [1] "numeric"class(surv2$ST)## [1] "factor"class(surv2$STATUS)## [1] "numeric"attach (surv2)Pstr Cox Model
# right censored:
cox2<- coxph (Surv(AGE, STATUS)~ TREATMENT, method='breslow')Model Summary
summary(cox2)## Call:
## coxph(formula = Surv(AGE, STATUS) ~ TREATMENT, method = "breslow")
##
## n= 258, number of events= 28
##
## coef exp(coef) se(coef) z Pr(>|z|)
## TREATMENTNATURAL 0.0164 1.0165 0.3923 0.042 0.9666
## TREATMENTSUNLIGHT -1.9166 0.1471 0.7596 -2.523 0.0116 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## exp(coef) exp(-coef) lower .95 upper .95
## TREATMENTNATURAL 1.0165 0.9837 0.47119 2.193
## TREATMENTSUNLIGHT 0.1471 6.7980 0.03319 0.652
##
## Concordance= 0.642 (se = 0.041 )
## Likelihood ratio test= 11.86 on 2 df, p=0.003
## Wald test = 6.88 on 2 df, p=0.03
## Score (logrank) test = 9.28 on 2 df, p=0.01KAPLAN-MEIER SURVIVAL CURVE (PSTR)
km2 <-survfit(Surv(AGE,STATUS)~ TREATMENT)Plot Curves
plot(km2, lwd = 2.5, xlab = 'Time (weeks post-settlement)', ylab = 'Probability of survival',
conf.int = FALSE, lty = c(1, 1, 1),
col = c('#1811ac', '#bdd07d', '#f7b58b'),
main= " ", ylim = c(0.5, 1),xlim = c(5, 14))
legend('bottomleft',
c( 'LED Blue-shifted spectrum', 'LED Reef-depth spectrum', 'Sunlight at surface depth '),
lty = c(1, 1, 1),
col =c('#1811ac', '#bdd07d', '#f7b58b'),
bty = 'n')
detach (surv2)P.clivosa Cox Model and Survival Curves
#CSV contains Pcli data only
#Import data
surv3 <- read_csv("/Users/Day/Downloads/SurvPcli.csv")## New names:
## Rows: 222 Columns: 22
## ── Column specification
## ──────────────────────────────────────────────────────── Delimiter: "," chr
## (6): ID, TANK, SPECIES, TREATMENT, TILE, ST dbl (2): AGE, STATUS lgl (14):
## ...9, ...10, ...11, ...12, ...13, ...14, ...15, ...16, ...17, ...1...
## ℹ 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.
## • `` -> `...9`
## • `` -> `...10`
## • `` -> `...11`
## • `` -> `...12`
## • `` -> `...13`
## • `` -> `...14`
## • `` -> `...15`
## • `` -> `...16`
## • `` -> `...17`
## • `` -> `...18`
## • `` -> `...19`
## • `` -> `...20`
## • `` -> `...21`
## • `` -> `...22`#Define factors
View(surv3)
#predictors
surv3$TREATMENT<- factor(surv3$TREATMENT)
surv3$SPECIES<- factor(surv3$SPECIES)
#view object type
class(surv3$TREATMENT)## [1] "factor"class(surv3$SPECIES)## [1] "factor"class(surv3$AGE)## [1] "numeric"class(surv3$STATUS)## [1] "numeric"attach (surv3)PCLI Cox Model
# right censored:
cox3<- coxph (Surv(AGE, STATUS)~ TREATMENT, method='breslow')Model Summary
summary(cox3)## Call:
## coxph(formula = Surv(AGE, STATUS) ~ TREATMENT, method = "breslow")
##
## n= 222, number of events= 45
##
## coef exp(coef) se(coef) z Pr(>|z|)
## TREATMENTNATURAL 0.64800 1.91171 0.36195 1.790 0.0734 .
## TREATMENTSUNLIGHT -0.04344 0.95749 0.40827 -0.106 0.9153
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## exp(coef) exp(-coef) lower .95 upper .95
## TREATMENTNATURAL 1.9117 0.5231 0.9404 3.886
## TREATMENTSUNLIGHT 0.9575 1.0444 0.4301 2.131
##
## Concordance= 0.591 (se = 0.04 )
## Likelihood ratio test= 4.87 on 2 df, p=0.09
## Wald test = 5.03 on 2 df, p=0.08
## Score (logrank) test = 5.22 on 2 df, p=0.07KAPLAN-MEIER SURVIVAL CURVE (PCLI)
km3 <-survfit(Surv(AGE,STATUS)~ TREATMENT)Plot Curves
plot(km3, lwd = 2.5, xlab = 'Time (weeks post-settlement)', ylab = 'Probability of survival',
conf.int = FALSE, lty = c(1, 1, 1),
col = c('#1811ac', '#bdd07d', '#f7b58b'),
main= "", ylim = c(0.5, 1),xlim = c(5, 14))
legend('bottomleft',
c( 'LED Blue-shifted spectrum', 'LED Reef-depth spectrum', 'Sunlight at surface depth '),
lty = c(1, 1, 1),
col =c('#1811ac', '#bdd07d', '#f7b58b'),
bty = 'n')
detach(surv3)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] survival_3.5-5 dplyr_1.1.2 readr_2.1.3 readxl_1.4.1
##
## loaded via a namespace (and not attached):
## [1] highr_0.10 cellranger_1.1.0 bslib_0.4.2 compiler_4.2.2
## [5] pillar_1.9.0 jquerylib_0.1.4 tools_4.2.2 bit_4.0.5
## [9] digest_0.6.31 lattice_0.20-45 jsonlite_1.8.4 evaluate_0.20
## [13] lifecycle_1.0.3 tibble_3.2.1 pkgconfig_2.0.3 rlang_1.1.1
## [17] Matrix_1.5-3 cli_3.6.0 rstudioapi_0.14 parallel_4.2.2
## [21] yaml_2.3.6 xfun_0.36 fastmap_1.1.0 stringr_1.5.0
## [25] knitr_1.41 generics_0.1.3 vctrs_0.6.2 sass_0.4.4
## [29] hms_1.1.2 bit64_4.0.5 grid_4.2.2 tidyselect_1.2.0
## [33] glue_1.6.2 R6_2.5.1 fansi_1.0.3 vroom_1.6.0
## [37] rmarkdown_2.20 tzdb_0.3.0 magrittr_2.0.3 htmltools_0.5.4
## [41] ellipsis_0.3.2 splines_4.2.2 utf8_1.2.2 stringi_1.7.12
## [45] cachem_1.0.6 crayon_1.5.2