DICE_basic function.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## This Jupyter Notebook contains basic functions for the implementation of algorithms in the paper [_The Development and Deployment of a Model for Hospital-level COVID-19 Associated Patient Demand Intervals from Consistent Estimators (DICE)_](https://arxiv.org/abs/2011.09377)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Libraries\n",
    "library(ggplot2)\n",
    "library(dplyr)\n",
    "library(hrbrthemes)\n",
    "library(changepoint)\n",
    "Sys.setenv(LANGUAGE = \"en\")\n",
    "library(reshape2)\n",
    "library(forecast)\n",
    "library(gridExtra)\n",
    "library(invgamma)\n",
    "rm(list=ls())\n",
    "#work_dir='...'\n",
    "#setwd(work_dir)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Actual hospital-level hospitalizations\n",
    "shc.data=read.csv('...')\n",
    "\n",
    "# Actual county-level hospitalizations\n",
    "scc.data=read.csv('...')\n",
    "\n",
    "# county-level hospitalization forecasts\n",
    "county.forecast=read.csv('...')\n",
    "\n",
    "\n",
    "N_list = scc.data$Hospitalizations\n",
    "A_list = shc.data$covid_AAU\n",
    "B_list = shc.data$covid_ICU\n",
    "\n",
    "hat_p = sum(A_list)/sum(N_list)\n",
    "hat_q = sum(B_list)/sum(N_list)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Perfect Forcasts: Bootstrapping intervals"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "bootstrap_perfect_forecast<-function(n,delta,alpha,b0,hat_p,hat_q,F_list,F_r){\n",
    "    N_lambda=F_list\n",
    "    A_lower_plugin_list=c()\n",
    "    A_upper_plugin_list=c()\n",
    "        \n",
    "    B_lower_plugin_list=c()\n",
    "    B_upper_plugin_list=c()  \n",
    "    for(b in 1:b0){\n",
    "        A_star_s=c()\n",
    "        B_star_s=c()\n",
    "        N_star_s=c()\n",
    "\n",
    "    \n",
    "         N_star_s=rpois(n=length(N_lambda),N_lambda)\n",
    "        for(j in 1:n){\n",
    "            multi_result = sapply(N_star_s[j], rmultinom,n=1, prob=c(hat_p,hat_q,(1-hat_p-hat_q)))\n",
    "            A_j = multi_result[1,]\n",
    "            B_j = multi_result[2,]\n",
    "            A_star_s=c(A_star_s,A_j)\n",
    "            B_star_s=c(B_star_s,B_j)           \n",
    "        }\n",
    "        p_star = sum(A_star_s)/sum(N_star_s)\n",
    "        q_star = sum(B_star_s)/sum(N_star_s)\n",
    "        \n",
    "        A_lower_plugin_b=qpois(delta/2,p_star*F_r)\n",
    "        A_upper_plugin_b=qpois(1-delta/2,p_star*F_r)\n",
    "        \n",
    "        B_lower_plugin_b=qpois(delta/2,q_star*F_r)\n",
    "        B_upper_plugin_b=qpois(1-delta/2,q_star*F_r)\n",
    "        \n",
    "        A_lower_plugin_list=c(A_lower_plugin_list,A_lower_plugin_b)\n",
    "        A_upper_plugin_list=c(A_upper_plugin_list,A_upper_plugin_b)\n",
    "        \n",
    "        B_lower_plugin_list=c(B_lower_plugin_list,B_lower_plugin_b)\n",
    "        B_upper_plugin_list=c(B_upper_plugin_list,B_upper_plugin_b)        \n",
    "    }\n",
    "    A_lower_plugin=qpois(delta/2,hat_p*F_r)\n",
    "    A_upper_plugin=qpois(1-delta/2,hat_p*F_r)\n",
    "    B_lower_plugin=qpois(delta/2,hat_q*F_r)\n",
    "    B_upper_plugin=qpois(1-delta/2,hat_q*F_r)\n",
    "    \n",
    "    z_al=ceiling(quantile(A_lower_plugin_list-A_lower_plugin,1-alpha))\n",
    "    z_bl=ceiling(quantile(B_lower_plugin_list-B_lower_plugin,1-alpha))\n",
    "    \n",
    "    z_ar=floor(quantile(A_upper_plugin_list-A_upper_plugin,alpha))\n",
    "    z_br=floor(quantile(B_upper_plugin_list-B_upper_plugin,alpha))\n",
    "    \n",
    "    A_lower_bootstrap=A_lower_plugin-z_al\n",
    "    B_lower_bootstrap=B_lower_plugin-z_bl\n",
    "    A_upper_bootstrap=A_upper_plugin-z_ar\n",
    "    B_upper_bootstrap=B_upper_plugin-z_br\n",
    "    \n",
    "    return(list(c(A_lower_plugin,A_upper_plugin,B_lower_plugin,B_upper_plugin),\n",
    "               c(A_lower_bootstrap,A_upper_bootstrap,B_lower_bootstrap,B_upper_bootstrap),\n",
    "               A_lower_plugin_list,A_upper_plugin_list,B_lower_plugin_list,B_upper_plugin_list))\n",
    "}"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Unbiased Forecasts with Lognormal Errors"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "bootstrap_unbiased_forecast<-function(n,delta,alpha,b0,m,hat_p,hat_q,hat_mu,hat_sigma2,hat_rho,F_list,F_r){\n",
    "    A_lower_plugin_list=c()\n",
    "    A_upper_plugin_list=c()\n",
    "        \n",
    "    B_lower_plugin_list=c()\n",
    "    B_upper_plugin_list=c()  \n",
    "    for(b in 1:b0){\n",
    "        A_star_s=c()\n",
    "        B_star_s=c()\n",
    "        N_star_s=c()\n",
    "        Y=rep(0,n)\n",
    "        Y[1]=rnorm(n=1,mean=hat_mu/(1-hat_rho),sd=sqrt(hat_sigma2/(1-hat_rho^2)))\n",
    "        for(j in 2:n){\n",
    "            Y[j]=hat_rho*Y[j-1]+rnorm(n=1,mean=hat_mu,sd=sqrt(hat_sigma2))      \n",
    "        }\n",
    "        N_star_s=rpois(n=n,F_list*exp(Y))\n",
    "        fraction1 = mean((N_star_s^2-N_star_s)/F_list^2)\n",
    "        fraction2 = mean(N_star_s[1:(length(N_star_s)-1)]*N_star_s[2:length(N_star_s)]/(F_list[1:(length(F_list)-1)]*F_list[2:length(F_list)]))\n",
    "\n",
    "        fn <- function(param){\n",
    "            rho=param[1]\n",
    "            sigma2=param[2]\n",
    "            mu=-sigma2/(2+2*rho)\n",
    "            obj = 0\n",
    "            obj = obj+(fraction1-exp(2*mu/(1-rho)+2*sigma2/(1-rho^2)))^2\n",
    "            obj = obj+(fraction2-exp(2*mu/(1-rho)+sigma2/(1-rho)))^2\n",
    "            obj = obj+10000*max(0,(rho^2-1))+10000*max(0,(-sigma2))\n",
    "            return(obj)\n",
    "        }\n",
    "        result = optim(rep(0,2),lower=c(-0.99,0),upper=c(0.99,1),fn,method=\"L-BFGS-B\")\n",
    "        hat_rho_star=result$par[1]\n",
    "        hat_sigma2_star=result$par[2]\n",
    "        hat_mu_star=-hat_sigma2_star/(2+2*hat_rho_star^2)\n",
    "        \n",
    "        for(j in 1:n){\n",
    "            multi_result = sapply(N_star_s[j], rmultinom,n=1, prob=c(hat_p,hat_q,(1-hat_p-hat_q)))\n",
    "            A_j = multi_result[1,]\n",
    "            B_j = multi_result[2,]\n",
    "            A_star_s=c(A_star_s,A_j)\n",
    "            B_star_s=c(B_star_s,B_j)           \n",
    "        }\n",
    "        p_star = sum(A_star_s)/sum(N_star_s)\n",
    "        q_star = sum(B_star_s)/sum(N_star_s)\n",
    "        \n",
    "        A_r_star=rep(0,m)\n",
    "        B_r_star=rep(0,m)\n",
    "        for(i in 1:m){\n",
    "            Yr=rnorm(n=1,mean=hat_mu_star/(1-hat_rho_star),\n",
    "                     sd=sqrt(hat_sigma2_star/(1-hat_rho_star^2)))\n",
    "            A_r_star[i]=rpois(n=1,lambda=p_star*F_r*exp(Yr))\n",
    "            B_r_star[i]=rpois(n=1,lambda=q_star*F_r*exp(Yr))\n",
    "        }\n",
    "               \n",
    "        A_lower_plugin_b=ceiling(quantile(A_r_star,delta/2))\n",
    "        A_upper_plugin_b=floor(quantile(A_r_star,1-delta/2))\n",
    "        \n",
    "        B_lower_plugin_b=ceiling(quantile(B_r_star,delta/2))\n",
    "        B_upper_plugin_b=floor(quantile(B_r_star,1-delta/2))\n",
    "        \n",
    "        A_lower_plugin_list=c(A_lower_plugin_list,A_lower_plugin_b)\n",
    "        A_upper_plugin_list=c(A_upper_plugin_list,A_upper_plugin_b)\n",
    "        \n",
    "        B_lower_plugin_list=c(B_lower_plugin_list,B_lower_plugin_b)\n",
    "        B_upper_plugin_list=c(B_upper_plugin_list,B_upper_plugin_b)        \n",
    "    }\n",
    "    A_r=rep(0,m)\n",
    "    B_r=rep(0,m)\n",
    "    for(i in 1:m){\n",
    "        Yr=rnorm(n=1,mean=hat_mu/(1-hat_rho),\n",
    "                     sd=sqrt(hat_sigma2/(1-hat_rho^2)))\n",
    "        A_r[i]=rpois(n=1,lambda=hat_p*F_r*exp(Yr))\n",
    "        B_r[i]=rpois(n=1,lambda=hat_q*F_r*exp(Yr))\n",
    "    }\n",
    "               \n",
    "    \n",
    "    A_lower_plugin=ceiling(quantile(A_r,delta/2))\n",
    "    A_upper_plugin=floor(quantile(A_r,1-delta/2))\n",
    "    B_lower_plugin=ceiling(quantile(B_r,delta/2))\n",
    "    B_upper_plugin=floor(quantile(B_r,1-delta/2))\n",
    "\n",
    "    \n",
    "    z_al=ceiling(quantile(A_lower_plugin_list-A_lower_plugin,1-alpha))\n",
    "    z_bl=ceiling(quantile(B_lower_plugin_list-B_lower_plugin,1-alpha))\n",
    "    \n",
    "    z_ar=floor(quantile(A_upper_plugin_list-A_upper_plugin,alpha))\n",
    "    z_br=floor(quantile(B_upper_plugin_list-B_upper_plugin,alpha))\n",
    "    \n",
    "    A_lower_bootstrap=A_lower_plugin-z_al\n",
    "    B_lower_bootstrap=B_lower_plugin-z_bl\n",
    "    A_upper_bootstrap=A_upper_plugin-z_ar\n",
    "    B_upper_bootstrap=B_upper_plugin-z_br\n",
    "    \n",
    "    return(list(c(A_lower_plugin,A_upper_plugin,B_lower_plugin,B_upper_plugin),\n",
    "               c(A_lower_bootstrap,A_upper_bootstrap,B_lower_bootstrap,B_upper_bootstrap),\n",
    "               A_lower_plugin_list,A_upper_plugin_list,B_lower_plugin_list,B_upper_plugin_list))\n",
    "}"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Biased Forecasts with Lognormal Errors"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "bootstrap_biased_forecast <- function(n,delta,alpha,b0,m,hat_p,hat_q,hat_mu,hat_sigma2,hat_rho,F_list,F_r){\n",
    "    A_lower_plugin_list=c()\n",
    "    A_upper_plugin_list=c()\n",
    "        \n",
    "    B_lower_plugin_list=c()\n",
    "    B_upper_plugin_list=c()  \n",
    "    for(b in 1:b0){\n",
    "        A_star_s=c()\n",
    "        B_star_s=c()\n",
    "        N_star_s=c()\n",
    "        Y=rep(0,n)\n",
    "        Y[1]=rnorm(n=1,mean=hat_mu/(1-hat_rho),sd=sqrt(hat_sigma2/(1-hat_rho^2)))\n",
    "        for(j in 2:n){\n",
    "            Y[j]=hat_rho*Y[j-1]+rnorm(n=1,mean=hat_mu,sd=sqrt(hat_sigma2))      \n",
    "        }\n",
    "        N_star_s=rpois(n=n,F_list*exp(Y))\n",
    "        \n",
    "        fraction1 = mean(N_star_s/F_list)\n",
    "        fraction2 = mean((N_star_s^2-N_star_s)/F_list^2)\n",
    "        fraction3 = mean(N_star_s[1:(length(N_star_s)-1)]*N_star_s[2:length(N_star_s)]/(F_list[1:(length(F_list)-1)]*F_list[2:length(F_list)]))\n",
    "\n",
    "        fn <- function(param){\n",
    "            rho=param[1]\n",
    "            sigma2=param[2]\n",
    "            mu=param[3]\n",
    "            obj = 0\n",
    "            obj = obj+(fraction1-exp(mu/(1-rho)+0.5*sigma2/(1-rho^2)))^2\n",
    "            obj = obj+(fraction2-exp(2*mu/(1-rho)+2*sigma2/(1-rho^2)))^2\n",
    "            obj = obj+(fraction3-exp(2*mu/(1-rho)+sigma2/(1-rho)))^2\n",
    "            obj = obj+10000*max(0,(rho^2-1))+10000*max(0,(-sigma2))\n",
    "            return(obj)\n",
    "        }\n",
    "        result = optim(rep(0,3),lower=c(-0.99,0.0001,-4),upper=c(0.99,4,4),fn,method=\"L-BFGS-B\")\n",
    "\n",
    "        hat_rho_star=result$par[1]\n",
    "        hat_sigma2_star=result$par[2]\n",
    "        hat_mu_star=result$par[3]\n",
    "        \n",
    "        for(j in 1:n){\n",
    "            multi_result = sapply(N_star_s[j], rmultinom,n=1, prob=c(hat_p,hat_q,(1-hat_p-hat_q)))\n",
    "            A_j = multi_result[1,]\n",
    "            B_j = multi_result[2,]\n",
    "            A_star_s=c(A_star_s,A_j)\n",
    "            B_star_s=c(B_star_s,B_j)           \n",
    "        }\n",
    "        p_star = sum(A_star_s)/sum(N_star_s)\n",
    "        q_star = sum(B_star_s)/sum(N_star_s)\n",
    "        \n",
    "        A_r_star=rep(0,m)\n",
    "        B_r_star=rep(0,m)\n",
    "        for(i in 1:m){\n",
    "            Yr=rnorm(n=1,mean=hat_mu_star/(1-hat_rho_star),\n",
    "                     sd=sqrt(hat_sigma2_star/(1-hat_rho_star^2)))\n",
    "            A_r_star[i]=rpois(n=1,lambda=p_star*F_r*exp(Yr))\n",
    "            B_r_star[i]=rpois(n=1,lambda=q_star*F_r*exp(Yr))\n",
    "        }\n",
    "               \n",
    "        A_lower_plugin_b=ceiling(quantile(A_r_star,delta/2))\n",
    "        A_upper_plugin_b=floor(quantile(A_r_star,1-delta/2))\n",
    "        \n",
    "        B_lower_plugin_b=ceiling(quantile(B_r_star,delta/2))\n",
    "        B_upper_plugin_b=floor(quantile(B_r_star,1-delta/2))\n",
    "        \n",
    "        A_lower_plugin_list=c(A_lower_plugin_list,A_lower_plugin_b)\n",
    "        A_upper_plugin_list=c(A_upper_plugin_list,A_upper_plugin_b)\n",
    "        \n",
    "        B_lower_plugin_list=c(B_lower_plugin_list,B_lower_plugin_b)\n",
    "        B_upper_plugin_list=c(B_upper_plugin_list,B_upper_plugin_b)        \n",
    "    }\n",
    "    A_r=rep(0,m)\n",
    "    B_r=rep(0,m)\n",
    "    for(i in 1:m){\n",
    "        Yr=rnorm(n=1,mean=hat_mu/(1-hat_rho),\n",
    "                     sd=sqrt(hat_sigma2/(1-hat_rho^2)))\n",
    "        A_r[i]=rpois(n=1,lambda=hat_p*F_r*exp(Yr))\n",
    "        B_r[i]=rpois(n=1,lambda=hat_q*F_r*exp(Yr))\n",
    "    }\n",
    "               \n",
    "    \n",
    "    A_lower_plugin=ceiling(quantile(A_r,delta/2))\n",
    "    A_upper_plugin=floor(quantile(A_r,1-delta/2))\n",
    "    B_lower_plugin=ceiling(quantile(B_r,delta/2))\n",
    "    B_upper_plugin=floor(quantile(B_r,1-delta/2))\n",
    "\n",
    "    \n",
    "    z_al=ceiling(quantile(A_lower_plugin_list-A_lower_plugin,1-alpha))\n",
    "    z_bl=ceiling(quantile(B_lower_plugin_list-B_lower_plugin,1-alpha))\n",
    "    \n",
    "    z_ar=floor(quantile(A_upper_plugin_list-A_upper_plugin,alpha))\n",
    "    z_br=floor(quantile(B_upper_plugin_list-B_upper_plugin,alpha))\n",
    "    \n",
    "    A_lower_bootstrap=A_lower_plugin-z_al\n",
    "    B_lower_bootstrap=B_lower_plugin-z_bl\n",
    "    A_upper_bootstrap=A_upper_plugin-z_ar\n",
    "    B_upper_bootstrap=B_upper_plugin-z_br\n",
    "    \n",
    "    return(list(c(A_lower_plugin,A_upper_plugin,B_lower_plugin,B_upper_plugin),\n",
    "               c(A_lower_bootstrap,A_upper_bootstrap,B_lower_bootstrap,B_upper_bootstrap),\n",
    "               A_lower_plugin_list,A_upper_plugin_list,B_lower_plugin_list,B_upper_plugin_list))\n",
    "}"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "R",
   "language": "R",
   "name": "ir"
  },
  "language_info": {
   "codemirror_mode": "r",
   "file_extension": ".r",
   "mimetype": "text/x-r-source",
   "name": "R",
   "pygments_lexer": "r",
   "version": "3.6.1"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}