In dchodge/hero-study: HERO study statistical analysis

setwd("~/documents/research/covid/hero/serobayes")

save_all_datafit()
data(datafit) # loads all data-sets for fitting into environment
load_fit_mcmc_all()   # loads all mcmc fits into environment

save_all_dataplt()
data(dataplt) # loads all data-sets for plotting figures

Check MCMC

# m1
# says bad for n eff for a few generated quantities but these are fixed values so their sd is 0.

check_hmc_diagnostics(m1_sero[[1]][[1]]) # fine
check_hmc_diagnostics(m1_sero[[1]][[2]]) # fine
check_hmc_diagnostics(m1_sero[[2]][[1]]) # fine
check_hmc_diagnostics(m1_sero[[2]][[2]]) # fine
check_hmc_diagnostics(m1_sero[[3]][[1]]) # fine
check_hmc_diagnostics(m1_sero[[3]][[2]]) # fine

check_hmc_diagnostics(m1_asymp[[1]][[1]]) # fine
check_hmc_diagnostics(m1_asymp[[1]][[2]]) # fine
check_hmc_diagnostics(m1_asymp[[2]][[1]]) # fine
check_hmc_diagnostics(m1_asymp[[2]][[2]]) # fine
check_hmc_diagnostics(m1_asymp[[3]][[1]]) # fine
check_hmc_diagnostics(m1_asymp[[3]][[2]]) # fine

check_hmc_diagnostics(m2_start[[1]]) # fine
check_hmc_diagnostics(m2_start[[2]]) # fine
check_hmc_diagnostics(m2_start[[3]]) # fine
check_hmc_diagnostics(m2_start[[4]]) # fine

check_hmc_diagnostics(m2_change[[1]]) # fine
check_hmc_diagnostics(m2_change[[2]]) # fine
check_hmc_diagnostics(m2_change[[3]]) # fine
check_hmc_diagnostics(m2_change[[4]]) # fine


data_full_more <- data_full$full %>% rename(IgG_S_unit_log2 = IgG_S_unit,
                          IgG_N_unit_log2 = IgG_N_unit,
                          IgG_S_unit_hero_log2 = IgG_S_unit_hero,
                          IgG_N_unit_hero_log2 = IgG_N_unit_hero,
                          IgA_Serum_S_unit_log2 = IgA_Serum_S_unit,
                          IgA_Serum_N_unit_log2 = IgA_Serum_N_unit,
                          IgA_Oral_S_unit_log2 = IgA_Oral_S_unit,
                          IgA_Oral_N_unit_log2 = IgA_Oral_N_unit) %>% 
  mutate(IgG_S_unit = case_when(IgG_S_unit_log2>0~(2^IgG_S_unit_log2), IgG_S_unit_log2 == 0~0),
         IgG_N_unit = case_when(IgG_N_unit_log2>0~(2^IgG_N_unit_log2), IgG_N_unit_log2 == 0~0),
         IgG_S_unit_hero = case_when(IgG_S_unit_hero_log2>0~(2^IgG_S_unit_hero_log2), IgG_S_unit_hero_log2 == 0~0), 
         IgG_N_unit_hero = case_when(IgG_N_unit_hero_log2>0~(2^IgG_N_unit_hero_log2), IgG_N_unit_hero_log2 == 0~0), 
         IgA_Serum_S_unit = case_when(IgA_Serum_S_unit_log2>0~(2^IgA_Serum_S_unit_log2), IgA_Serum_S_unit_log2 == 0~0), 
         IgA_Serum_N_unit = case_when(IgA_Serum_N_unit_log2>0~(2^IgA_Serum_N_unit_log2), IgA_Serum_N_unit_log2 == 0~0),
         IgA_Oral_S_unit = case_when(IgA_Oral_S_unit_log2>0~(2^IgA_Oral_S_unit_log2)*2000, IgA_Oral_S_unit_log2 == 0~0), 
         IgA_Oral_N_unit = case_when(IgA_Oral_N_unit_log2>0~(2^IgA_Oral_N_unit_log2)*2000, IgA_Oral_N_unit_log2 == 0~0))


data_full_more$IgG_S_unit %>% max

raw_data <- read.csv(system.file('extdata', 'HERO_data_clean_21Dec2020.csv', package='serobayes'))
raw_data$IgG_V1_spike_antibody_units %>% max(na.rm = TRUE)
raw_data$IgG_V2_spike_antibody_units %>% max(na.rm = TRUE)

write.table(data_full_more, here::here("inst", "extdata", "clean_hero.csv"), sep = ",", row.names = FALSE)

Some interesting stats

data_all <- data_full$full %>% filter(sample_no == 1)
data_complete_r <- data_full$full[!(data_full$full$IgG_S_OD %>% is.na),]
data_complete <- data_complete_r %>% filter(sample_no == 2)
data_pos <- data_full$sero_pos %>% filter(sample_no == 1)
data_symp <- data_pos %>% filter(symp_pos == 0)

all <- data_all %>% filter(sample_no == 1)
allpos <- data_all %>% filter(sample_no == 1, sero_pos == 1)
allasymp <- data_all %>% filter(sample_no == 1, sero_pos == 1, symp_pos == 0)

(allasymp %>% nrow)
(allpos %>% nrow)

# loc_coded, job_grouped, c19zones, gender, ethnic, age_group
values <- all %>% pull(loc_coded) %>% table
valuespos <- allpos %>% pull(gender) %>% table
valuesasymp <- allasymp %>% pull(gender) %>% table
data.frame(number = valuesasymp, prec = valuesasymp / valuespos)

allneg <- data_full$full %>% filter(sample_no == 1, sero_pos == 0)
allposV1 <- data_full$full %>% filter(sample_no == 2, sero_pos == 1)
allrevert <- allposV1 %>% filter(Sample_ID %in% allneg$Sample_ID)

valuesneg <- allneg %>% pull(loc_coded) %>% table
valuesrevert <- allrevert %>% pull(loc_coded) %>% table 
data.frame(number = valuesrevert, prec = valuesrevert / valuesneg)


allsinci %>% nrow
(allneg %>% nrow)


vals <- seq_len(length(valuespos)) %>% map(~allasymp %>% filter(age_group == valuespos[.x]) %>% nrow %>% setNames(valuespos[.x])) %>% unlist
c(24, 18, 24, 25, 10) / c(51, 60, 61, 67, 25)

(data_symp$gender %>% table)/(data_pos$gender %>% table) * 100

(data_pos$c19zone %>% table)/(data_complete$c19zone %>% table)
(data_pos$gender %>% table)/(data_complete$gender %>% table)
data_pos$gender %>% length

data_complete$ethnic %>% table

# number of samples recruited
length(data_full$full$HERO_ID %>% unique)

# number of samples recruited
length(data_full$full$IgG_S_OD %>% na.omit)
sum(data_full$full$IgG_S_OD %>% is.na)

# number seropositive on first bleed
data_full$full %>%
  pipeline_exclude_strange_dates %>%
  filter(sero_pos == 1, sample_no == 1) %>%
  nrow

data_full$full %>%
  pipeline_exclude_strange_dates %>%
  filter(sero_pos == 1, sample_no == 1, c19_diagnosis_category == "No") %>%
  nrow

# number who seroconverted
length(data_full$sero_con$HERO_ID %>% unique)

# number who seroreverted
length(data_full$sero_revert$HERO_ID %>% unique)

# get dates
  df_dates <- data.frame(
    "Symptom Onset" = as.Date(data_full$symp$inf_date, origin = "2020-01-01")[seq(1, 326, 2)],
    "First bleed" = as.Date(data_full$symp$date_of_sample , origin = "2020-01-01")[seq(1, 326, 2)],
    "Second bleed" = as.Date(data_full$symp$date_of_sample , origin = "2020-01-01")[seq(2, 326, 2)]
  ) %>% tidyr::gather(event, value)
df_dates %>% filter(event == "First.bleed") %>% pull(value) %>% sort
df_dates %>% filter(event == "Second.bleed") %>% pull(value) %>% sort

# Seroprevalence overall
post_sero <- extract(m1_sero[[1]][[1]])
props <- (datafit_prev$c19zones %>% table %>% as.numeric)/nrow(datafit_prev)
(post_sero$marginal[,,1]*props) %>% apply(1, sum) %>% mean
(post_sero$marginal[,,1]*props) %>% apply(1, sum) %>% quantile(c(0.025, 0.975))

(post_sero$marginal[,6,1]) %>% mean
(post_sero$marginal[,6,1]) %>% quantile(c(0.025, 0.975))


post_sero <- extract(m1_sero[[2]][[2]])
# Mean and ci for OT
post_sero$marginal[,9,1] %>% mean
post_sero$marginal[,9,1] %>% quantile(c(0.025, 0.975))
# Mean and ci for Allied
post_sero$marginal[,2,1] %>% mean
post_sero$marginal[,2,1] %>% quantile(c(0.025, 0.975))

post_sero <- extract(m1_sero[[3]][[2]])
# Mean and ci for OT
post_sero$marginal[,2,1] %>% mean
post_sero$marginal[,2,1] %>% quantile(c(0.025, 0.975))


# Seroprevalence overall
post_sero <- extract(m1_asymp[[1]][[1]])
props <- (datafit_prev$c19zones %>% table %>% as.numeric)/nrow(datafit_prev)
(post_sero$marginal[,,1]*props) %>% apply(1, sum) %>% mean
(post_sero$marginal[,,1]*props) %>% apply(1, sum) %>% quantile(c(0.025, 0.975))

post1 <- extract(m2_change[[2]])
get_mean_ci(post1$marginal_s)

inter_start <- c("marginal_a_i", "marginal_e_i", "marginal_g_i", "marginal_s_i")
slope_change <- c("marginal_a_s", "marginal_e_s", "marginal_g_s", "marginal_s_s")

start <- init_model_ab_kin(name = "start",
                           stan_file = "start",
                           datafit = datafit_start,
                           marginal_vars = inter_start,
                           y_axis = "log(AU) at first bleed")
change <- init_model_ab_kin(name = "change",
                            stan_file = "change",
                            datafit = datafit_change,
                            marginal_vars = slope_change,
                            y_axis = "Change in log(AU) after 28 days")

start_sample <- clean_samples(start)
change_sample <- clean_samples(change)

all_start_S <- 1:4 %>% map(~rnorm(1000, start_sample[.x, ]$`50%`, (start_sample[.x, ]$`50%` - start_sample[.x, ]$`2.5%`) / 2) ) %>% unlist
get_mean_ci(all_start_S)

all_start_N <- 14:18 %>% map(~rnorm(1000, start_sample[.x, ]$`50%`, (start_sample[.x, ]$`50%` - start_sample[.x, ]$`2.5%`) / 2) ) %>% unlist
get_mean_ci(all_start_N)

all_change_S <- 1:4 %>% map(~rnorm(1000, change_sample[.x, ]$`50%`, (change_sample[.x, ]$`50%` - change_sample[.x, ]$`2.5%`) / 2) ) %>% unlist
get_mean_ci(all_change_S)

all_change_N <- 14:18 %>% map(~rnorm(1000, change_sample[.x, ]$`50%`, (change_sample[.x, ]$`50%` - change_sample[.x, ]$`2.5%`) / 2) ) %>% unlist
get_mean_ci(all_change_N)

bind_rows(
  get_mean_ci(all_start_S),
  get_mean_ci(all_start_N),
  get_mean_ci(all_change_S),
  get_mean_ci(all_change_N)
)

data(fit_mcmc_start)

x <- str2lang(paste0(start$marginal[1], "[j]"))
sd <- medianvalues$sigma %>% unique %>% `*`(10)
full <- m2_start[[2]] %>%
        as_draws_df %>%
        tidybayes::spread_draws(sigma, !!x)
medianvalues <- m2_start[[2]] %>%
        as_draws_df %>%
        tidybayes::spread_draws(sigma, !!x) %>%
        ggdist::median_qi() %>%
        as.data.frame

fullsamples <- list()
for (i in 1:5) {
  age <- full %>% filter(j == i) %>% pull(marginal_a_i)
  agefull <- 1:4000 %>% map(~rnorm(1, age[.x], sd)) %>% unlist

  fullsamples[[i]] <- data.frame(
    age = i,
    titre_start = agefull
  )
}

fullsamples %>% bind_rows %>% group_by(age) %>%
  dplyr::summarise(mean = mean(titre_start), sd = sd(titre_start))





medianvalues %>%