pathways_paper_fitting.Rmd

---
title: "Patient pathways paper - code for LoS fitting"
author: "Quentin J. Leclerc, Gwenan M. Knight"
date: "`r format(Sys.time(), '%A %d %B %Y')`"
output:
  bookdown::pdf_document2:
    fig_width: 7.5
    fig_height: 4
---


# Setup 

We load packages, pathways data, and read in the movement data from Excel.

```{r setup, include = F}
knitr::opts_chunk$set(echo = FALSE, message = F, warning = F)

library(openxlsx)
library(dplyr)
library(linelist)
library(knitr)
library(reshape2)
library(rlist)
library(here)
library(Hmisc)

source(here::here("code", "core_model_functions.R"))
source(here::here("code", "paths_table_to_blocks.R"))

set.seed(3210)

```

```{r load_hospital_data}

#latest version of hosp data can be generated by sourcing the corresponding code:
#(generates a ton of objects, recommend clearing global envir after!)
#source(here::here("code", "import_hosp_data.R"))

hospital_data = read.csv(here::here("data", "uk_hospital_data.csv")) %>%
  mutate(date = as.Date(date))

regions_hospital_data = hospital_data %>%
  filter(!(geography %in% c("England", "Scotland", "Wales", "Northern Ireland")))

cov_curve = hospital_data %>%
  filter(geography == "England") %>%
  filter(value_type == "hospital_inc")


```

```{r calc_region_props}
#calculate proportion of cumulative hospitalisations in each region
#used later to calculate fitted weighted average LoS for England

hosp_props = hospital_data %>%
  filter(value_type == "hospital_inc") %>%
  filter(geography != "England") %>%
  group_by(geography) %>%
  summarise(n = sum(value)) %>%
  ungroup %>%
  mutate(n = n/(sum(n)))

```

## Fitting LoS by region

```{r fit_by_region}
all_regions_results = data.frame()
all_regions_true_prev = data.frame()
all_regions_LoS = data.frame()


optim_func_ward = function(par, total_LoS_all, cov_curve, true_prev){
  
  message("Try Ward LoS ", par)
  
  total_LoS_all$Ward$LoS$mean = par
  
  #giving it a spin with total LoS assumptions
  res_total_LoS_all = multi_pathways_model(10, cov_curve, total_LoS_all, length(cov_curve))
  
  sq_diff = sum((res_total_LoS_all$Ward - true_prev$Ward)^2, na.rm = T)
  
  sq_diff
  
}

optim_func_CC = function(par, total_LoS_all, cov_curve, true_prev){
  
  message("Try CC LoS ", par)
  
  total_LoS_all$CC$LoS$mean = par
  
  #giving it a spin with total LoS assumptions
  res_total_LoS_all = multi_pathways_model(10, cov_curve, total_LoS_all, length(cov_curve))
  
  sq_diff = sum((res_total_LoS_all$CC - true_prev$CC)^2, na.rm = T)
  
  sq_diff
  
}


for (i in unique(regions_hospital_data$geography)) {
  
  message("Currently on region ", i)
  
  data_reg = regions_hospital_data %>%
    filter(geography == i)
  
  cov_curve = data_reg %>%
    filter(value_type == "hospital_inc") %>%
    select(value) %>%
    pull
  
  true_reg = data.frame(
    time = c(1:length(unique(data_reg$date))),
    date = unique(data_reg$date),
    Ward = data_reg %>% filter(value_type == "hospital_prev") %>% select(value) %>% pull,
    CC = data_reg %>% filter(value_type == "icu_prev") %>% select(value) %>% pull,
    region = i
  )
  
  opt_res_ward = optimize(optim_func_ward, lower = 1, upper = 10,
                          tol = 0.01, total_LoS_all = total_LoS_all,
                          cov_curve = cov_curve, true_prev = true_reg)
  opt_res_CC = optimize(optim_func_CC, lower = 8, upper = 20,
                        tol = 0.01, total_LoS_all = total_LoS_all,
                        cov_curve = cov_curve, true_prev = true_reg)
  
  total_LoS_all$Ward$LoS$mean = opt_res_ward$minimum
  total_LoS_all$CC$LoS$mean = opt_res_CC$minimum
  
  #giving it a spin with total LoS assumptions
  res_total_LoS_all = multi_pathways_model(100, cov_curve, total_LoS_all, length(cov_curve))
  res_total_LoS_all$region = i
  res_total_LoS_all$date = true_reg$date
  
  all_regions_results = rbind(all_regions_results, res_total_LoS_all)
  all_regions_true_prev = rbind(all_regions_true_prev, true_reg)
  
  
  res_region_LoS = data.frame(
    Ward = opt_res_ward$minimum,
    CC = opt_res_CC$minimum,
    region = i
  )
  
  all_regions_LoS = rbind(all_regions_LoS, res_region_LoS)
  
}

eng_fitted = c(Hmisc::wtd.mean(all_regions_LoS$Ward, hosp_props$n),
               Hmisc::wtd.mean(all_regions_LoS$CC, hosp_props$n),
               "England")

saveRDS(all_regions_results, here::here("outputs", "results", "fit_regions_eng.rds"))

all_regions_LoS = rbind(all_regions_LoS, eng_fitted)

write.csv(all_regions_LoS, here::here("outputs", "tables", "best_fit_LoS.csv"), row.names = F)


```




## Fitting prop to CC by region

```{r fit_by_region_proba}
all_regions_results2 = data.frame()
all_regions_true_prev2 = data.frame()
all_regions_proba2 = data.frame()


optim_func2 = function(par, total_LoS_all, cov_curve, true_prev){
  
  message("Try CC proba ", par)
  
  total_LoS_all$CC$proba = par
  total_LoS_all$Ward$proba = 1-par
  
  
  #giving it a spin with total LoS assumptions
  res_total_LoS_all = multi_pathways_model(10, cov_curve, total_LoS_all, length(cov_curve))
  
  sq_diff = sum((res_total_LoS_all$Ward - true_prev$Ward)^2, na.rm = T)
  sq_diff2 = sum((res_total_LoS_all$CC - true_prev$CC)^2, na.rm = T)
  
  sq_diff + sq_diff2
  
}


for (i in unique(regions_hospital_data$geography)) {
  
  message("Currently on region ", i)
  
  data_reg = regions_hospital_data %>%
    filter(geography == i)
  
  cov_curve = data_reg %>%
    filter(value_type == "hospital_inc") %>%
    select(value) %>%
    pull
  
  true_reg = data.frame(
    time = c(1:length(unique(data_reg$date))),
    date = unique(data_reg$date),
    Ward = data_reg %>% filter(value_type == "hospital_prev") %>% select(value) %>% pull,
    CC = data_reg %>% filter(value_type == "icu_prev") %>% select(value) %>% pull,
    region = i
  )
  
  opt_res = optimize(optim_func2, lower = 0.01, upper = 0.5,
                     tol = 0.01, total_LoS_all = total_LoS_all,
                     cov_curve = cov_curve, true_prev = true_reg)
  
  total_LoS_all$CC$proba = opt_res$minimum
  total_LoS_all$Ward$proba = 1-opt_res$minimum
  
  #giving it a spin with total LoS assumptions
  res_total_LoS_all = multi_pathways_model(100, cov_curve, total_LoS_all, length(cov_curve))
  res_total_LoS_all$region = i
  res_total_LoS_all$date = true_reg$date
  
  all_regions_results2 = rbind(all_regions_results2, res_total_LoS_all)
  all_regions_true_prev2 = rbind(all_regions_true_prev2, true_reg)
  
  
  res_region_proba = data.frame(
    Ward = 1-opt_res$minimum,
    CC = opt_res$minimum,
    region = i
  )
  
  all_regions_proba2 = rbind(all_regions_proba2, res_region_proba)
  
}

eng_fitted2 = c(Hmisc::wtd.mean(all_regions_proba2$Ward, hosp_props$n),
                Hmisc::wtd.mean(all_regions_proba2$CC, hosp_props$n),
                "England")

saveRDS(all_regions_results2, here::here("outputs", "results", "fit_regions_eng_proba.rds"))

all_regions_proba2 = rbind(all_regions_proba2, eng_fitted2)



write.csv(all_regions_proba2, here::here("outputs", "tables", "best_fit_proba.csv"), row.names = F)

```