Skip to content

Commit

Permalink
survival analysis truncated F1 to 217 dof
Browse files Browse the repository at this point in the history
  • Loading branch information
@SamGurr
SamGurr committed Aug 29, 2024
1 parent 06a4c79 commit a57e64a
Show file tree
Hide file tree
Showing 10 changed files with 878 additions and 11 deletions.
103 changes: 103 additions & 0 deletions RAnalysis/Data/Physiology/Survival/F2/cumulative_raw/F2_survival_raw_growout.csv
View file
Original file line number Diff line number Diff line change
@@ -0,0 +1,103 @@
Date,Age,Treatment,Rep,Count_corrected,Survival
10/5/2022,50,Low,A,107,1
10/5/2022,50,Low,B,108,1
10/5/2022,50,Low,C,119,1
10/5/2022,50,Low,D,115,1
10/5/2022,50,Moderate,A,114,1
10/5/2022,50,Moderate,B,95,1
10/5/2022,50,Moderate,C,106,1
10/5/2022,50,Moderate,D,113,1
10/5/2022,50,High,A,109,1
10/5/2022,50,High,B,147,1
10/5/2022,50,High,C,118,1
10/5/2022,50,High,D,116,1
10/26/2022,71,Low,A,107,1
10/26/2022,71,Low,B,108,1
10/26/2022,71,Low,C,113,0.949579832
10/26/2022,71,Low,D,114,0.991304348
10/26/2022,71,Moderate,A,114,1
10/26/2022,71,Moderate,B,95,1
10/26/2022,71,Moderate,C,106,1
10/26/2022,71,Moderate,D,110,0.973451327
10/26/2022,71,High,A,109,1
10/26/2022,71,High,B,147,1
10/26/2022,71,High,C,117,0.991525424
10/26/2022,71,High,D,116,1
11/14/2022,90,Low,A,105,0.981308411
11/14/2022,90,Low,B,104,0.962962963
11/14/2022,90,Low,C,111,0.932773109
11/14/2022,90,Low,D,109,0.947826087
11/14/2022,90,Moderate,A,114,1
11/14/2022,90,Moderate,B,95,1
11/14/2022,90,Moderate,C,106,1
11/14/2022,90,Moderate,D,108,0.955752212
11/14/2022,90,High,A,105,0.963302752
11/14/2022,90,High,B,143,0.972789116
11/14/2022,90,High,C,106,0.898305085
11/14/2022,90,High,D,114,0.982758621
12/5/2022,111,Low,A,96,0.897196262
12/5/2022,111,Low,B,84,0.777777778
12/5/2022,111,Low,C,94,0.789915966
12/5/2022,111,Low,D,92,0.8
12/5/2022,111,Moderate,A,101,0.885964912
12/5/2022,111,Moderate,B,86,0.905263158
12/5/2022,111,Moderate,C,102,0.962264151
12/5/2022,111,Moderate,D,99,0.876106195
12/5/2022,111,High,A,99,0.908256881
12/5/2022,111,High,B,131,0.891156463
12/5/2022,111,High,C,93,0.788135593
12/5/2022,111,High,D,104,0.896551724
1/27/2023,164,Low,A,96,0.897196262
1/27/2023,164,Low,B,82,0.759259259
1/27/2023,164,Low,C,94,0.789915966
1/27/2023,164,Low,D,91,0.791304348
1/27/2023,164,Moderate,A,101,0.885964912
1/27/2023,164,Moderate,B,82,0.863157895
1/27/2023,164,Moderate,C,101,0.952830189
1/27/2023,164,Moderate,D,93,0.82300885
1/27/2023,164,High,A,99,0.908256881
1/27/2023,164,High,B,129,0.87755102
1/27/2023,164,High,C,93,0.788135593
1/27/2023,164,High,D,102,0.879310345
2/2/2023,164,Low,A,96,0.897196262
2/2/2023,164,Low,B,82,0.759259259
2/2/2023,164,Low,C,94,0.789915966
2/2/2023,164,Low,D,91,0.791304348
2/2/2023,164,Low,E,,
2/2/2023,164,Low,F,,
2/2/2023,164,Low,G,,
2/2/2023,164,Moderate,A,101,0.885964912
2/2/2023,164,Moderate,B,82,0.863157895
2/2/2023,164,Moderate,C,101,0.952830189
2/2/2023,164,Moderate,D,93,0.82300885
2/2/2023,164,Moderate,E,,
2/2/2023,164,Moderate,F,,
2/2/2023,164,Moderate,G,,
2/2/2023,164,High,A,99,0.908256881
2/2/2023,164,High,B,129,0.87755102
2/2/2023,164,High,C,93,0.788135593
2/2/2023,164,High,D,102,0.879310345
2/2/2023,164,High,E,,
2/2/2023,164,High,F,,
2/2/2023,164,High,G,,
2/27/2023,195,Low,A,90,0.841121495
2/27/2023,195,Low,B,76,0.703703704
2/27/2023,195,Low,C,87,0.731092437
2/27/2023,195,Low,D,87,0.756521739
2/27/2023,195,Low,E,,
2/27/2023,195,Low,F,,
2/27/2023,195,Low,G,,
2/27/2023,195,Moderate,A,98,0.859649123
2/27/2023,195,Moderate,B,72,0.757894737
2/27/2023,195,Moderate,C,98,0.924528302
2/27/2023,195,Moderate,D,93,0.82300885
2/27/2023,195,Moderate,E,,
2/27/2023,195,Moderate,F,,
2/27/2023,195,Moderate,G,,
2/27/2023,195,High,A,99,0.908256881
2/27/2023,195,High,B,129,0.87755102
2/27/2023,195,High,C,90,0.762711864
2/27/2023,195,High,D,99,0.853448276
2/27/2023,195,High,E,,
2/27/2023,195,High,F,,
2/27/2023,195,High,G,,
Binary file added RAnalysis/Output/Survival/F1/F1_Survival_217dpf_ggsurvfit.pdf
View file
Binary file not shown.
Binary file added RAnalysis/Output/Survival/F2/F2_Survival_ggsurvfit.pdf
View file
Binary file not shown.
Binary file added RAnalysis/Output/Survival/F2/F2_Survival_ggsurvfit_LvH.pdf
View file
Binary file not shown.
Binary file added RAnalysis/Output/Survival/F2/F2_Survival_ggsurvfit_LvM.pdf
View file
Binary file not shown.
Binary file added RAnalysis/Output/Survival/F2/F2_Survival_hazardratio_LvH.pdf
View file
Binary file not shown.
Binary file added RAnalysis/Output/Survival/F2/F2_Survival_hazardratio_LvM.pdf
View file
Binary file not shown.
Binary file added RAnalysis/Output/Survival/F2/F2_Survival_meanSE.pdf
View file
Binary file not shown.
172 changes: 161 additions & 11 deletions RAnalysis/Scripts/F1_Phys/F1_Survival.Rmd
View file
Original file line number Diff line number Diff line change
Expand Up @@ -171,13 +171,26 @@ binary_surv_df <- binary_surv_df %>%
Replicate== 'E'~4,
Replicate== 'F'~5,
Replicate== 'G'~6,
Replicate== 'H'~7)) %>%
Replicate== 'H'~7),
TankID = as.numeric(
as.factor(
paste0(
Replicate, Treatment)
)
)
) %>%
na.omit()
unique(binary_surv_df$Replicate) # 1 2 3 4 5 6 7
# convert to numeric to be sure
binary_surv_df$Replicate <- as.numeric(binary_surv_df$Replicate)
# For meaningful comparison to F2, truncate binary_surv_df to 217 dpf (F2 ended at 219 dpf)
binary_surv_df_217dpf <- binary_surv_df %>% filter(Age <= 217)
unique(binary_surv_df_217dpf$Age) # 23 51 64 92 129 147 182 217
```
# Create surv objects
* binary_surv_df - all data with replicate as random factor
Expand All @@ -196,6 +209,29 @@ surv_obj.means <- survfit2(Surv(Age, Count_dead) ~ Treatment, data = binary_surv
# n events median 0.95LCL 0.95UCL
# Treatment=Low 5686 651 318 303 346
# Treatment=Moderate 5739 862 303 273 303
# For meaningful comparison to F2, truncate binary_surv_df to 217 dpf (F2 ended at 219 dpf)
binary_surv_df_217dpf$Age <- as.numeric(binary_surv_df_217dpf$Age)
surv_obj_217dpf.all <- survfit2(Surv(Age, Count_dead) ~ Treatment+(1|Replicate),
data = binary_surv_df_217dpf)
# n events median 0.95LCL 0.95UCL
# Treatment=Low, 1 | Replicate=TRUE 22840 3268 217 217 217
# Treatment=Moderate, 1 | Replicate=TRUE 24192 5307 182 182 182
binary_surv_df_217dpf$Age <- as.numeric(binary_surv_df_217dpf$Age)
surv_obj_217dpf.means <- survfit2(Surv(Age, Count_dead) ~ Treatment,
data = binary_surv_df_217dpf)
# n events median 0.95LCL 0.95UCL
# Treatment=Low 22840 3268 217 217 217
# Treatment=Moderate 24192 5307 182 182 182
```


Expand Down Expand Up @@ -247,6 +283,7 @@ F1_Survplot_MEANS <- ggsurvplot(surv_obj.means,
pdf(paste0("Output/Survival/F1/F1_Survival_plot.pdf"),width=8, height=3)
print(F1_Survplot_MEANS)
dev.off()
```


Expand Down Expand Up @@ -288,6 +325,57 @@ pdf(paste0("Output/Survival/F1/F1_Survival_ggsurvfit.pdf"),width=8, height=4)
print(F1_KaplanMeier.means)
dev.off()
# For meaningful comparison to F2, truncate binary_surv_df to 217 dpf (F2 ended at 219 dpf)
F1_KaplanMeier_217dpf.all <- surv_obj_217dpf.all %>%
ggsurvfit(linetype_aes = TRUE, linewidth = 1) +
scale_color_manual(values = c("forestgreen","orange")) +
scale_fill_manual(values = c("forestgreen","orange")) +
labs(
x = "Days",
y = "Overall survival probability"
) +
add_confidence_interval() +
# add_risktable() +
# add_risktable_strata_symbol() + # (symbol = "\U25CF", size = 10
# scale_ggsurvfit(x_scales = list(breaks = 0:1)) #+
add_pvalue(location = "annotation", x = 8.5)
F1_KaplanMeier_217dpf.all
F1_KaplanMeier_217dpf.means <- surv_obj_217dpf.means %>%
ggsurvfit(linetype_aes = TRUE, linewidth = 1) +
scale_color_manual(values = c("forestgreen","orange")) +
scale_fill_manual(values = c("forestgreen","orange")) +
labs(
x = "Days",
y = "Overall survival probability"
) +
add_confidence_interval() +
# add_risktable() +
# add_risktable_strata_symbol() + # (symbol = "\U25CF", size = 10
# scale_ggsurvfit(x_scales = list(breaks = 0:1)) #+
add_pvalue(location = "annotation", x = 8.5)
F1_KaplanMeier_217dpf.means
# output plot
pdf(paste0("Output/Survival/F1/F1_Survival_217dpf_ggsurvfit.pdf"),width=8, height=4)
print(F1_KaplanMeier_217dpf.means)
dev.off()
```

## Estimating x-time survival
Expand Down Expand Up @@ -379,7 +467,7 @@ F1_logrank_test.means <- survdiff(formula = Surv(Age, Count_dead) ~
# run cox test on all data with replicate as a random factor
F1_coxph_test <- coxph(formula = Surv(Age, Count_dead) ~ Treatment+(1|Replicate), data = binary_surv_df)
F1_coxph_test <- coxph(formula = Surv(Age, Count_dead) ~ Treatment+(1|TankID), data = binary_surv_df)
summary(F1_coxph_test)
# coef exp(coef) se(coef) z Pr(>|z|)
# TreatmentModerate 0.24904 1.28280 0.02429 10.25 <2e-16 ***
Expand All @@ -388,26 +476,23 @@ summary(F1_coxph_test)
ggforest(F1_coxph_test, data = binary_surv_df)
# Likelihood ratio test=106 on 1 df, p=< 2.2e-16
# n= 47946, number of events= 6902
F1_coxph_test_table <- coxph(Surv(Age, Count_dead) ~ Treatment+(1|Replicate), data = binary_surv_df) %>%
F1_coxph_test_table <- coxph(Surv(Age, Count_dead) ~ Treatment+(1|TankID), data = binary_surv_df) %>%
tbl_regression(exp = TRUE)
# COXME IS THE PROPER WAY TO ADDRESS THE TANK REPLICATE AS A RANDOM FACTOR - ALITHOU GGFOREST IS NOT COMPTATIBLE
# THEREFORE, USE THE OUTPUT EXP AND c95 CI FROM THE COXME BUT PLOT THE COXPH HAZARD RATIO - ADJUST PDF MANUALLY OFF R
# run mixed cox model with replicate as a random factor
F1_coxme_mixed_test <- coxme(Surv(Age, Count_dead) ~ # using coxme, notice ther results are the same as coxpH
Treatment + (1 | Replicate), # note we cannot make a hazrd ratio plot using coxme
Treatment + (1 | TankID), # note we cannot make a hazrd ratio plot using coxme
data = binary_surv_df)
# coef exp(coef) se(coef) z p
# TreatmentModerate 0.1211895 1.128839 0.01470208 8.24 2.2e-16
# Random effects
# Group Variable Std Dev Variance
# Replicate Intercept 0.2404987 0.0578396
# coef exp(coef) se(coef) z p
# TreatmentModerate 0.2097335 1.233349 0.1537181 1.36 0.17
F1_coxme_mixed_HR <- coxme(Surv(Age, Count_dead) ~ # using coxme, notice ther results are the same as coxpH
Treatment + (1 | Replicate), # note we cannot make a hazrd ratio plot using coxme
Treatment + (1 | TankID), # note we cannot make a hazrd ratio plot using coxme
data = binary_surv_df) %>%
tbl_regression(exp = TRUE) # thi gives us what we need
# 1.13 1.10, 1.16 <0.001 # THESE ARE THE CORRECT VALUES - OUTPUT THE COXPH HAZARD GGFIREST BUT MANUALLY ADJUST FOR THESE EXPONENTS!!
# 1.23 0.91, 1.67 0.2 # THESE ARE THE CORRECT VALUES - OUTPUT THE COXPH HAZARD GGFIREST BUT MANUALLY ADJUST FOR THESE EXPONENTS!!
Expand All @@ -433,4 +518,69 @@ F1_coxph_test_table.mean <- coxph(Surv(Age, Count_dead) ~ Treatment,
```


```{r 217dpf cox for multiple variables}
# Hazard ratio is equivalent to the exp(coef)
# HR > 1 indicates an increased risk of death, An HR < 1, on the other hand, indicates a decreased risk
# For example, a hazard ratio of 0.25 for treatment groups tells you that individuals who received treatment
# B have a reduced risk of dying compared to patients who received treatment A
# run cox test on all data with replicate as a random factor
F1_coxph_217dpf <- coxph(formula = Surv(Age, Count_dead) ~
Treatment+(1|TankID),
data = binary_surv_df_217dpf)
summary(F1_coxph_217dpf)
# coef exp(coef) se(coef) z Pr(>|z|)
# TreatmentModerate 0.36187 1.43601 0.02225 16.27 <2e-16 ***
# 1 | TankIDTRUE NA NA 0.00000 NA NA
# Hazard ratio = 1.28280
ggforest(F1_coxph_217dpf, data = binary_surv_df_217dpf)
F1_coxph_217dpf_table <- coxph(Surv(Age, Count_dead) ~
Treatment+(1|TankID),
data = binary_surv_df_217dpf) %>%
tbl_regression(exp = TRUE)
# 1.44 1.37, 1.50 <0.001 # tell sus that moderate had signicantly lower survival, however coxpH does NOT address
# the tank factor correctly, need to use coxme below, these data however output out of R (havenot found out how to
# output coxme) therefore output, do not report but use as a template to report the coxme below!!
# COXME IS THE PROPER WAY TO ADDRESS THE TANK REPLICATE AS A RANDOM FACTOR - ALITHOU GGFOREST IS NOT COMPTATIBLE
# THEREFORE, USE THE OUTPUT EXP AND c95 CI FROM THE COXME BUT PLOT THE COXPH HAZARD RATIO - ADJUST PDF MANUALLY OFF R
# run mixed cox model with replicate as a random factor
F1_coxme_mixed_217dpf <- coxme(Surv(Age, Count_dead) ~ # using coxme, notice ther results are the same as coxpH
Treatment + (1 | TankID), # note we cannot make a hazrd ratio plot using coxme
data = binary_surv_df_217dpf)
# coef exp(coef) se(coef) z p
# TreatmentModerate 0.439905 1.55256 0.3328382 1.32 0.19
F1_coxme_mixed_217dpf_HR <- coxme(Surv(Age, Count_dead) ~ # using coxme, notice ther results are the same as coxpH
Treatment + (1 | TankID), # note we cannot make a hazrd ratio plot using coxme
data = binary_surv_df_217dpf) %>%
tbl_regression(exp = TRUE) # thi gives us what we need
# 1.55 0.81, 2.98 0.2 # THESE ARE THE CORRECT VALUES - OUTPUT THE COXPH HAZARD GGFIREST BUT MANUALLY ADJUST FOR THESE EXPONENTS!!
# run cox test on all data as a mean
# * no replicate tank, these are the binary data as a mean
F1_coxph_test.means <- coxph(formula = Surv(Age, Count_dead) ~ Treatment, data = binary_surv_df_MEANS)
summary(F1_coxph_test.means)
# coef exp(coef) se(coef) z Pr(>|z|)
# TreatmentModerate 0.13669 1.14647 0.03066 4.459 8.25e-06 ***
HaxardRatio_plot <- ggforest(F1_coxph_test.means,
data = binary_surv_df_MEANS)
pdf(paste0("Output/Survival/F1/F1_Survival_hazardratio.pdf"),width=6, height=2)
print(HaxardRatio_plot)
dev.off()
# Likelihood ratio test=106 on 1 df, p=< 2.2e-16
# n= 47946, number of events= 6902
F1_coxph_test_table.mean <- coxph(Surv(Age, Count_dead) ~ Treatment,
data = binary_surv_df_MEANS) %>%
tbl_regression(exp = TRUE)
```

Loading

0 comments on commit a57e64a

Please sign in to comment.