severeC19_3_model_over60.R

rm(list = ls())
invisible(gc(full = T))

library(data.table)
library(dplyr)
library(lubridate)
library(stringr)
library(ggplot2)
library(scales)
library(survival)
library(broom)
library(tidyr)
library(rlang)
library(fst)
library(future.apply)

setwd("N:/output/scripts/acorbetta/")
na_to_0 = function(x){ifelse(is.na(x), 0,x)}
tresh22 = ymd("2022-06-01")


df = read_fst("N:/output/data/COV-CVD/working_severeC19_static_QCed.fst")
setDT(df)

df = df[age_2019 + 2 >= 60]

#### create tv datasets ####

setDTthreads(percent = 100)        # let data.table use all CPU threads
options(future.globals.maxSize = 12 * 1024^3)  # 6 GiB
plan(multisession, workers = 2)  # choose number of workers

# --- source data as data.table ---
# df is your original data.frame

# keep only columns we need (if some are absent, intersect handles it)
need <- c("COD_SOGGETTO","AGE_CLA","SEX","MCS_CLA","URBAN_DOM","CONT", 'PRESYMP',
          paste0("OUTCOME", 1:4),
          paste0("DATA_END", 1:4), paste0("DATE_END_", 1:4),
          paste0("DATA_VACCINO", 1:4), "begin")
DT <- df[, intersect(names(df), need), with = FALSE]

rm(df)
gc(full=T)
# If you have a constant study start (e.g. tresh22), set it here once.
DT[, begin := tresh22]   # uncomment if you use a constant start

# --- helper: map end columns (supports DATA_END* or DATE_END*) ---
map_end <- function(i) {
  # prefer DATA_ENDi; else DATE_ENDi; else DATE_END_i (with underscore) if present
  cands <- c(sprintf("DATA_END_%d", i), sprintf("DATE_END_%d", i), sprintf("DATE_END_%d", i))
  cands[cands %in% names(DT)][1]
}

END_COLS <- vapply(1:4, map_end, character(1))
if (any(is.na(END_COLS))) stop("Couldn't find all 4 end-date columns.")

# --- vaccination backfill / monotone cleanup (vectorised) ---
# rule: if dose k exists, ensure dose k-1 exists; if missing or later than k, copy from k
for (k in 2:4) {
  v  <- sprintf("DATA_VACCINO%d", k)
  vp <- sprintf("DATA_VACCINO%d", k - 1)
  if (!v  %in% names(DT)) next
  if (!vp %in% names(DT)) DT[, (vp) := as.IDate(NA)]
  DT[!is.na(get(v)) & (is.na(get(vp)) | get(vp) > get(v)), (vp) := get(v)]
}
# ensure non-decreasing sequence across doses
for (k in 2:4) {
  v  <- sprintf("DATA_VACCINO%d", k)
  vp <- sprintf("DATA_VACCINO%d", k - 1)
  if (all(c(v, vp) %in% names(DT)))
    DT[!is.na(get(v)) & !is.na(get(vp)) & get(v) < get(vp), (v) := get(vp)]
}


expand_one_outcome <- function(z, out_col, end_col, dose_cols, baseline_cols) {
  # z is a one-row data.table
  
  start <- z$begin
  end   <- z[[end_col]]
  stat  <- z[[out_col]]
  
  # dose dates in canonical order (DATA_VACCINO1..4)
  all_doses <- if (length(dose_cols)) as.Date(unlist(z[, ..dose_cols])) else as.Date(character())
  
  n_pre <- sum(!is.na(all_doses) & all_doses <= start)
  onstudy_doses <- sort(all_doses[!is.na(all_doses) & all_doses > start & all_doses <= end])
  
  changes <- sort(unique(c(start, onstudy_doses, end)))
  if (length(changes) < 2L) return(NULL)
  
  tstart_date <- head(changes, -1L)
  tstop_date  <- tail(changes, -1L)
  
  n_on_by_tstart <- findInterval(tstart_date, onstudy_doses)
  exposure <- pmax(4L - (n_pre + n_on_by_tstart), 0L)
  
  ans <- data.table(
    COD_SOGGETTO = z$COD_SOGGETTO,
    tstart_date  = tstart_date,
    tstop_date   = tstop_date,
    tstart       = as.numeric(tstart_date - start),
    tstop        = as.numeric(tstop_date  - start),
    exposure     = as.integer(exposure),
    status       = as.integer(tstop_date == end & isTRUE(stat == 1)),
    begin        = start,
    end          = end
  )
  
  if (length(baseline_cols))
    ans[, (baseline_cols) := z[, ..baseline_cols][rep(1L, .N)]]
  
  ans
}

#  one-time setup (do this once, not inside build_one)
setDTthreads()                               # let data.table use all threads
setkeyv(DT, "COD_SOGGETTO")                  # fast grouping

dose_cols     <- intersect(names(DT), paste0("DATA_VACCINO", 1:4))
baseline_cols <- intersect(names(DT), c("AGE_CLA","SEX","MCS_CLA","URBAN_DOM","CONT","PRESYMP"))

#  faster build_one
build_one <- function(i) {
  out_col <- sprintf("OUTCOME%d", i)
  end_col <- END_COLS[i]
  
  sdcols <- c("COD_SOGGETTO", "begin", out_col, end_col, dose_cols, baseline_cols)
  
  res <- DT[
    , expand_one_outcome(.SD, out_col, end_col, dose_cols, baseline_cols),
    by = COD_SOGGETTO,
    .SDcols = sdcols
  ]
  
  res[, outcome_no := i]
  setkey(res, NULL)
  res[, COD_SOGGETTO := NULL]
  setcolorder(res, c("COD_SOGGETTO","outcome_no"))
  
  # Fast write (optional) — roughly 10–20× faster than RDS for big tables:
  write_fst(res, sprintf("N:/output/data/COV-CVD/OVER60_tv_dataset_outcome%d.fst", i), compress = 50)
  
  # If you must stick to .rds, at least turn off compression to save CPU:
  #saveRDS(res, sprintf("N:/output/data/COV-CVD/tmp_tv_dataset_outcome%d.rds", i), compress = FALSE)
  
  res
}


gc(full=T)
# Parallel build (four independent jobs)
future_lapply(3:4, build_one)


#### clean up ####
rm(DT)
gc(full=T)
#### RUN parallel model ####

library(future.apply)
options(future.globals.maxSize = 8 * 1024^3)  # 7 GiB
plan(multisession, workers = 4)  # choose number of workers

fit_tv_model_basic = function(df_adj){
  tv_model <- coxph(Surv(time = tstart, time2 = tstop, event = status) ~ 
                      VACCINE_DEFICIT + SEX + AGE_CLA + CONT + MCS_CLA + URBAN_DOM,
                    data = df_adj)
  return(tv_model)
}

files = list(c("N:/output/data/COV-CVD/OVER60_tv_dataset_outcome1.fst"),
             c("N:/output/data/COV-CVD/OVER60_tv_dataset_outcome2.fst"),
             c("N:/output/data/COV-CVD/OVER60_tv_dataset_outcome3.fst"),
             c("N:/output/data/COV-CVD/OVER60_tv_dataset_outcome4.fst"))


models <- future_lapply(files,function(out){
  df_adj = read_fst(out) %>% setDT()
  if (grepl("3|4", out)) {
    df_adj = df_adj[PRESYMP == 0]
  }
  df_adj[, `:=`(
    VACCINE_DEFICIT = as.factor(exposure),
    AGE_CLA = relevel(factor(AGE_CLA), ref = "[60,65)"))
    ]
  mod = fit_tv_model_basic(df_adj)
  return(mod)
})

names(models) = paste0("OUTCOME",1:4)
plan(sequential)
gc(full=T)

saveRDS(models,"N:/output/data/COV-CVD/OVER_60_tv_models_outcome.rds", compress = F)



#### parallel more adj ####
df = read_fst("N:/output/data/COV-CVD/working_severeC19_static_QCed.fst")
setDT(df)

df = df[age_2019 + 2 >= 60]

gc(full=T)

library(future.apply)
options(future.globals.maxSize = 8* 1024^3)  # 7 GiB
plan(multisession, workers = 4)  # choose number of workers

fit_tv_model = function(df_adj){
  tv_model <- coxph(Surv(time = tstart, time2 = tstop, event = status) ~ 
                      VACCINE_DEFICIT + SEX + AGE_CLA + CONT + MCS_CLA + URBAN_DOM + 
                      HOSP + CHOSP + N_PCR_POS + N_PCR + LAST_PCR + ATS,
                    data = df_adj)
  return(tv_model)
}

files_adj = list(c("N:/output/data/COV-CVD/OVER60_tv_dataset_outcome1.fst"),
             c("N:/output/data/COV-CVD/OVER60_tv_dataset_outcome2.fst"),
             c("N:/output/data/COV-CVD/OVER60_tv_dataset_outcome3.fst"),
             c("N:/output/data/COV-CVD/OVER60_tv_dataset_outcome4.fst"))

add_covs = c("COD_SOGGETTO","HOSP", "CHOSP", "N_PCR_POS", "N_PCR", "LAST_PCR", "ATS")
df_cov = df[,..add_covs]
(miss = apply(df_cov, 2, function(x){sum(is.na(x))}))
# 19480 missing CONT
# 40923 missing domicile
tmp = na.omit(df_cov,invert = T)
nrow(tmp) #40923
df_cov = df_cov[!COD_SOGGETTO %in% tmp$COD_SOGGETTO]
rm(tmp,df)
gc(full=T)

models <- future_lapply(files_adj,function(out){
  df_adj = read_fst(out) %>% setDT()
  if (grepl("3|4", out)) {
    df_adj = df_adj[PRESYMP == 0]
  }
  df_adj[, `:=`(
    VACCINE_DEFICIT = as.factor(exposure),
    AGE_CLA = relevel(factor(AGE_CLA), ref = "[60,65)"))
  ]
  df_adj = merge(df_adj, df_cov, by = "COD_SOGGETTO")
  mod = fit_tv_model(df_adj)
  return(mod)
})

names(models) = paste0("OUTCOME",1:4)
plan(sequential)

gc(full=T)
saveRDS(models, file = "N:/output/data/COV-CVD/OVER60_tv_models_outcome_add.rds",
        compress =F)


#### EXTRA ####

#fit_tv = coxph(Surv(time = tstart, time2 = tstop, event = OUTCOME2) ~ 
#                 VACCINE_DEFICIT + SEX + AGE_CLA + CONT + MCS_CLA + URBAN_DOM + 
#                 HOSP + CHOSP + N_PCR_POS + N_PCR + LAST_PCR + ATS,
#               data = df_adj ,weights = weight, x =T)

#new = df_adj %>%
#  filter(COD_SOGGETTO %in% SAMPLE_ALL)

#adjsurv <- adjustedsurv(data=new,
#                        variable="VACCINE_DEFICIT",
#                        ev_time="tstop",
#                        event="OUTCOME2",
#                        method="direct",
#                        outcome_model=fit_tv,
#                        conf_int=TRUE)
#
#p_adj = plot(adjsurv,conf_int=TRUE)
#ggsave("N:/output/plots/SEVERC19_ADJ_CURVES.png", plot = p_adj, width = 10, height = 10, dpi = 300)


test_ph_tv <- cox.zph(models$OUTCOME1)
test_ph_tv$table
plot(test_ph_tv, size = 0.1)