Incubation Phase

Parmetric survival model for incubation phase durataion

Clinical Onset

The find_clinical_onset() function identifies first occurrence of score > 0 then creates new “Event” column with 1 at this date, 0’s before this date, and a value 3 after that date.

Hide code

clin_start_df <- as.data.frame(
  find_clinical_onset(antem_df)
)

Prepare Data

Quick stats to check data

Hide code

unique(clin_start_df$Event)

[1] 0 1 3

Hide code

check_stats <- clin_start_df %>%
  filter(Event == 1 & censor_status == 1)

dim(check_stats)

[1] 13 16

Hide code

range(check_stats$dpe)

[1] 2 5

Hide code

mean(check_stats$dpe)

[1] 3.076923

Hide code

median(check_stats$dpe)

[1] 3

Remove Group 1 Group 1 was not infected and has no expectation of developing disease. Animals that were infected, verified by nasal, serum, or fever but never developed lesions are censored.

Hide code

clin_start_df <- clin_start_df %>%
 filter(group != "Group 1",
         Event == 1 | censor_status == 0 & censor_k == 1)

check data

Hide code

clin_start_df %>%
  filter(Event == 1) %>%
  summarise(mean_dpe = mean(dpe),
            median_dpe = median(dpe))

mean_dpe	median_dpe
3.076923	3

Scale Time
Add an arbitrarily small value to eliminate zeros. Really not need with his specific data set, by an important step.

Hide code

clin_start_df$scaled_duration <- clin_start_df$dpe + 0.0001

Weighted Contacts
Contact with groups by donors was sequential at 24hr intervals, shedding rates (nasal swabs) varied over this period.

Hide code

shed_rates <- antem_df %>%
  filter(group == "donor") %>%
  group_by(dpi) %>%
  summarise(tot_shed = sum(ifelse(nasal == 45, 0, nasal), na.rm=T)) %>%
  mutate(group = paste("Group", dpi, sep=" "))

# match cumulative shed from donors based on time of contact
clin_start_df$donor_shed <- with(shed_rates,
                                 tot_shed[match(
                                   clin_start_df$group,
                                                 group)])

# an integer index is needed for modeling purposes, group contacts are 1:1 correlated with time
clin_start_df$shed_time <- as.integer(as.factor(clin_start_df$group))

Models

Creating a survival model object to ensure time and censored animals are correctly indicated.

Hide code

surv_obj <- inla.surv(clin_start_df$scaled_duration, clin_start_df$Event)

Parametric Suvival
Accelerated Failure Time (AFT)

Parametric model for study-wide average duration.

Hide code

return_quants <- c(0.025, 0.05, 0.25, 0.5, 0.75, 0.95, 0.975)

pc_prec_iid <- list(theta = list(prior="pc.prec", 
                                 param=c(0.5, 0.001)))

incubation_dur <- inla(surv_obj ~ 1 + 
                         f(shed_time, donor_shed,
                           model = "rw1",
                           constr=TRUE,
                           scale.model = TRUE,
                           hyper=pc_prec_iid),
                        data = clin_start_df, 
                        verbose=FALSE,
                        quantiles = return_quants,
                        family = "exponential.surv",
                        control.fixed = list(prec = 1, prec.intercept  = 0.0001),
                        control.compute=list(dic = TRUE, cpo = FALSE, waic = TRUE))

Treatment is specific to group therefore using group below to identify important differences.

Hide code

return_quants <- c(0.025, 0.05, 0.25, 0.5, 0.75, 0.95, 0.975)

pc_prec_iid <- list(theta = list(prior="pc.prec", 
                                 param=c(0.5, 0.001)))

incubation_aft <- inla(surv_obj ~ 1 +
                       f(group,
                           model = "iid",
                           constr=FALSE,
                           hyper=pc_prec_iid),
                        data = clin_start_df, 
                        verbose=FALSE,
                        quantiles = return_quants,
                        family = "lognormal.surv",
                        control.fixed = list(prec = 1, prec.intercept  = 0.0001),
                        control.compute=list(dic = TRUE, cpo = FALSE, waic = TRUE))

Warning in sqrt(ee): NaNs produced

Sample Marginals

Performing sampling on the parametric model results

Hide code

incubation_samples <- compute_survival_marginals(incubation_dur, 14)

Check estimates at the 0.5 probability (median)

Hide code

median_incubation <- find_closest_quant(incubation_samples, 0.5)
median_incubation

$quant0.025
[1] 1.501

$quant0.25
[1] 2.301

$quant0.5
[1] 2.901

$quant0.75
[1] 3.701

$quant0.975
[1] 6.101

Survival Curve

Exceedance survival curve fro parametric model.

Hide code

plot_survival_marginals(incubation_samples, x_max = 14, xlabel = "Incubation Phase Duration")

AFT Treament Effects

Get treatment-level estimated duration

Hide code

aft_incubation <- incubation_aft$summary.random$group[,c(1,4,6,7,8, 10)]
names(aft_incubation) <- c("Group","Q_0.025","Q_0.25", "Q_0.5", "Q_0.75", "Q_0.975")
mean_aft <- incubation_aft$summary.fixed$mean

aft_incubation[,2:6] <- exp(aft_incubation[,2:6] + mean_aft)
study_wide <- exp(incubation_aft$summary.fixed[,c(3,5:7,9)])
rownames(study_wide) <- NULL
names(study_wide) <- c("Q_0.025","Q_0.25", "Q_0.5", "Q_0.75", "Q_0.975")
study_wide$Group <- "study"
aft_incubation <- rbind(study_wide, aft_incubation)
aft_incubation

Q_0.025	Q_0.25	Q_0.5	Q_0.75	Q_0.975	Group
2.289518	2.858287	3.170443	3.518581	4.405664	study
2.653615	3.380300	3.808423	4.301140	5.558763	donor
3.506350	4.422156	4.962956	5.582167	7.161666	Group 2
1.464660	1.880941	2.116881	2.376889	2.997234	Group 3
1.765721	2.255789	2.530923	2.832588	3.563349	Group 4

Compare treatment groups

Hide code

plot_aft_linerange(incubation_aft, ylabel = "Duration (days)")

Warning: Using `size` aesthetic for lines was deprecated in ggplot2 3.4.0.
ℹ Please use `linewidth` instead.

save incubation samples

Hide code

saveRDS(incubation_samples, here("assets/incubation_samples.rds"))
saveRDS(median_incubation, here("assets/incubation_median.rds"))
saveRDS(aft_incubation, here("assets/incubation_aft_median.rds"))