pathways_paper_fitting_nhstrust.Rmd

---
title: "Patient pathways paper - code for LoS fitting by NHS Trust"
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_nhstrust2.csv")) %>%
  mutate(date = as.Date(date))


```


## 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(hospital_data$geography)) {
  
  cat("\nCurrently on region", i)
  
  data_reg = 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 = 30,
                          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 = 1, upper = 30,
                        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)
  
}

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

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


```