kimfinale / kmisc Goto Github PK

---
title: "Kim Miscellaneous"
author: "Jong-Hoon Kim"
date: "3/20/2021"
output: html_document
editor_options: 
  chunk_output_type: console
---

```{r}
devtools::document()
```


###  Ordinary differential equation-based SEIR model

![](figs/seir.png){width=60%}

$$
\begin{aligned}
\frac{dS}{dt} &= - \beta S\frac{I}{N}\\ 
\frac{dE}{dt} &= \beta S\frac{I}{N} - \sigma E\\ 
\frac{dI}{dt} &= \sigma E - \gamma I\\ 
\frac{dR}{dt} &= \gamma I
\end{aligned}
$$
## Parameters
- Unit time = day
- $\beta$ = Transmission rate per unit time
- $1/\sigma$ = mean incubation period
- $\gamma$ = transition rate from $E$ to $I$ 

```{r echo=F}
devtools::load_all()
y0 <- c(S = 1e7-10, I = 10, R = 0, CI = 10) # initial values
params <- c(beta = 0.3, gamma = 0.2) # parameter values
times <- seq(0, 50, by = 1) # daily output for 150 days
```

```{r, echo=FALSE, message=FALSE }
library(deSolve)
library(tidyverse)
ode(y = y0, times = times, func = ode_sir, parms = params) %>% 
  as.data.frame() -> out 
```

```{r echo=FALSE}
out_long <- out %>% pivot_longer(- time) 
out_long$name <- factor(out_long$name, levels = c("S", "I", "R", "CI"))
out_long %>% 
  filter(name != "CE") %>%
  ggplot(aes(x = time, y = value, color = name))+
  geom_line(size = 1.2)+
  labs(x = 'Time (day)', y = 'Number of individuals', color = "")+ 
  theme_grey(base_size = 16)

daily_onset <- round(c(0, diff(out$CI)))
dput(daily_onset)
```

```{r echo=F}
y0 <- c(S = 1e3-10, E = 0, I = 10, R = 0, CE = 0, CI = 10) # initial values
params <- c(beta = 2.5/4, sigma = 1/4, gamma = 1/4) # parameter values
times <- seq(0, 150, by = 1) # daily output for 150 days
```

```{r, echo=FALSE, message=FALSE }
library(deSolve)
library(tidyverse)
ode(y = y0, times = times, func = ode_seir, parms = params) %>% 
  as.data.frame() -> out 
```

```{r echo=FALSE}
out_long <- out %>% pivot_longer(- time) 
out_long$name <- factor(out_long$name, levels = c("S", "E", "I", "R", "CE", "CI"))
out_long %>% 
  filter(name != "CE") %>%
  ggplot(aes(x = time, y = value, color = name))+
  geom_line(size = 1.2)+
  labs(x = 'Time (day)', y = 'Number of individuals', color = "")+ 
  theme_grey(base_size = 16)
```


### Stochastic model - tauleap

$$
\begin{align}
\Delta N_{SE} &= \textrm{Binomial}\left( S(t), 1-\textrm{exp}\{-\mu_{SE}\delta \}\right) \\
\Delta N_{EI} &= \textrm{Binomial}\left( E(t), 1-\textrm{exp}\{-\mu_{EI} \delta \}\right) \\
\Delta N_{IR} &= \textrm{Binomial}\left( I(t), 1-\textrm{exp}\{-\mu_{IR} \delta\}\right) \\
\end{align}
$$
$$\begin{align} 
S(t+\delta) &= S(t) - \Delta N_{SE}\\
E(t+\delta) &= E(t) + \Delta N_{SE} - \Delta N_{EI}\\
I(t+\delta) &= I(t) + \Delta N_{EI} - \Delta N_{IR}\\
R(t+\delta) &= R(t) + \Delta N_{IR}\\
\end{align}
$$
#### Run the tauleap model
```{r echo=FALSE, results=FALSE}
nrun <- 100 # number of runs
delta <- 0.1 # 0.1 day
# times <- seq(1, 4)
d <- matrix(NA, nrow = length(times), ncol = nrun)
for (i in 1:nrun) {
  x <- stoch_solve(func = stoch_seir, y = y0, times = times, params = params, delta = delta)
  d[, i]  <- x[,"I"]
}
# d
dmean <- rowMeans(d)
q <- apply(d, 1, quantile, probs = c(0.025, 0.25, 0.5, 0.75, 0.975))
nc <- ncol(q)
q[,(nc-5):nc]
```

#### compare the results between ODE and tauleap model by plotting
```{r echo=FALSE}
out %>%
  gather(state, value, -time) %>%
  filter(state == "I") -> d1
names(d1) <- c("time", "state", "value")
d1$state <- "ODE"

dd <- cbind(time = d1$time, as.data.frame(d) )
dd %>% as.data.frame() %>% 
  gather(state, value, -time) -> d2

dmean <- data.frame(time = d1$time, value = dmean)
dmed <- data.frame(time = d1$time, value = q[3,])

ggplot()+
  geom_line(data=d2, aes(x=time, y=value, group=state), color="grey50", alpha=0.2)+
  geom_line(data=d1, aes(x=time, y=value), color="darkred", size=2, inherit.aes=FALSE)+
  geom_line(data=dmean, aes(x=time, y=value), color="steelblue", size=2, inherit.aes=FALSE)+
  geom_segment(aes(x=rep(95,3), xend=rep(100,3), y=c(150,140,130), yend=c(150,140,130)), 
               color = c("darkred","steelblue","grey"), size=1.2) +
  annotate("text", x=rep(102,3), y=c(150,140,130), 
           label = c("ODE", "Mean of 1000 runs", "Individual runs"), hjust = 0) +
  labs(x='time (day)', y='number of infected individuals (I)') +
  theme_grey(base_size = 16) 
```

## Stochastic model - Gillespie's algorithm
### SEPAIR model 
#### Initial conditions

```{r}
# devtools::install()
devtools::load_all()
R0 <- 1.2
params <- c(beta = R0/5, epsilon = 1/2.5, delta = 1/5, 
            gamma = 1/2.5, rho_p = 1, 
            rho_a = 1, frac_a = 0.3, frac_detect = 1)
y0 <- c(S = 1e7, E = 0, P = 0, A = 0, I = 10, R = 0, CE = 0, CI = 10)
tend <- 175
```

#### Run the ODE version of SEPAIR model
```{r}
times <- seq(0, tend, by = 1) # daily output for 150 days
library(deSolve)
library(tidyverse)
ode(y = y0, times = times, func = ode_sepair, parms = params) %>% 
  as.data.frame() -> out 

out$daily_infected <- c(0, diff(out$CE)) 
out$daily_symptom_onset <- c(0, diff(out$CI)) 
out$daily_confirmed <- c(0, diff(out$R))

tstamp <- format(Sys.time(), "%Y%m%dT%H%M%S")
daily_ode <- list()
daily_ode$daily_infected <- c(0, diff(out$CE)) 
daily_ode$daily_symptom_onset <- c(0, diff(out$CI)) 
daily_ode$daily_confirmed <- c(0, diff(out$R))
# saveRDS(daily_ode, paste0("outputs/daily_ode_", tstamp, ".rds"))

Rt <- rep(NA, length(times))
for (i in seq_along(times)){
  Rt[i] <- get_Rt(times[i])  
}
out <- cbind(out, Rt)

out_long <- out %>% pivot_longer(-time) 
out_long$name <- factor(out_long$name,
                        levels = c("S", "E", "P", "A", "I", "R", 
                                  "CE", "CI", "daily_infected", 
                                  "daily_symptom_onset", "daily_confirmed", "Rt"))
out_long %>% 
  filter(name == "daily_infected" | 
           name == "daily_symptom_onset" | name == "daily_confirmed") %>%
  ggplot(aes(x = time, y = value, color = name)) +
  geom_line(size = 1.2) +
  labs(x='Time (day)', y = 'Number of individuals', color = "") + 
  theme_grey(base_size = 16) 
  # facet_wrap(vars(name), nrow = 2, scales = "free_y")

```

## run model implemented in Gillespie's direct method
```{r}
set.seed(1)
tstart <- Sys.time()
nrun <- 100
res <- gillespie_run(func = gillespie_sepair, 
                     tend = tend, 
                     nrun = nrun, 
                     y0 = y0, 
                     params = params, 
                     report_dt = 1)
Sys.time() - tstart

daily_inf <- data.frame(matrix(0, nrow = tend + 1, ncol = nrun))
daily_ons <- data.frame(matrix(0, nrow = tend + 1, ncol = nrun))
daily_con <- data.frame(matrix(0, nrow = tend + 1, ncol = nrun))

for (i in seq_len(length(res))) {
  inf <- diff(res[[i]]$CE)
  onset <- diff(res[[i]]$CI)
  conf <- diff(res[[i]]$R)
  daily_inf[1:(length(inf)+1), i] <- c(0, inf)
  daily_ons[1:(length(onset)+1), i] <- c(0, onset)
  daily_con[1:(length(conf)+1), i] <- c(0, conf)
}

df <- daily_inf
df$time <- 0:tend
df %>% pivot_longer(cols = -time) -> df

ggplot() +
  geom_line(data = df, aes(time, value, group = name), 
            color = "grey50") +
  geom_line(data = out, mapping = aes(time, daily_infected), color = "darkred",
            size = 2, inherit.aes = FALSE)

daily_stoch <- list()
daily_stoch$confirmed <- daily_con
daily_stoch$onset <- daily_ons
daily_stoch$infected <- daily_inf

# saveRDS(daily_stoch, paste0("outputs/daily_stoch_", tstamp, ".rds"))
```

### Apply a detection rate
```{r}
rr <- nrow(daily_stoch$confirmed)
cc <- ncol(daily_stoch$confirmed)
detected <- data.frame(matrix(0, nrow = rr, ncol = cc))
daily_stoch_partial_obs <- list()
probs <- seq(0.2,0.8,0.1)
daily_stoch_partial_obs$probs <- probs
set.seed(30)
for (k in 1:length(probs)) {
  prob <- probs[k]
  for (i in 1:rr) {
     for (j in 1:cc){
       size <- daily_stoch$confirmed[i, j]
       # message(paste0("i = ", i, ", j = ", j, ", case = ", size, "\n"))
       if (!is.na(size) & size > 0){
         detected[i, j] <- rbinom(1, size = size, prob = prob)
       }
    }
  }
  daily_stoch_partial_obs$detected[[k]] <- detected
}

# saveRDS(daily_stoch_partial_obs, paste0("outputs/daily_stoch_partial_obs_", tstamp, ".rds"))
```
## total dataset
```{r}
Rt_dataset <- list()
Rt_dataset$ode <- daily_ode
Rt_dataset$stoch_perfect_obs <- daily_stoch
Rt_dataset$stoch_partial_obs <- daily_stoch_partial_obs
# saveRDS(Rt_dataset, paste0("outputs/Rt_dataset_", tstamp, ".rds"))

Rt_dataset$stoch_perfect_obs$confirmed$X1[1:20]
Rt_dataset$stoch_partial_obs$detected[[7]]$X1[1:20]

library(ggplot2)

ggplot() +
  geom_line(aes(x=1:101, y = Rt_dataset$ode_perfect_obs$daily_infected)) + 
  geom_line(aes(x=1:101, y = Rt_dataset$stoch_perfect_obs$infected$X1), color = "blue", size = 1, alpha = 0.3)+ 
  geom_line(aes(x=1:101, y = Rt_dataset$stoch_perfect_obs$infected$X3),  color = "blue", size = 1, alpha = 0.6)+
  geom_line(aes(x=1:101, y = Rt_dataset$stoch_partial_obs$detected[[7]]$X1), color = "darkred", size = 1, alpha = 0.3)+
  geom_line(aes(x=1:101, y = Rt_dataset$stoch_partial_obs$detected[[7]]$X3), color = "darkred", size = 1, alpha = 0.6)+
  geom_line(aes(x=1:101, y = Rt_dataset$stoch_partial_obs$detected[[7]]$X1), color = "darkred")+
  geom_line(aes(x=1:101, y = Rt_dataset$stoch_partial_obs$detected[[7]]$X3), color = "darkred")+
  labs(x = "day", y = "daily infection")

plot(Rt_dataset$stoch_perfect_obs$confirmed[,1], type = "l", col = rgb(0.6,0.6,0.6, alpha=0.8), ylim=c(0,100))
for (i in 2:100){
  lines(Rt_dataset$stoch_perfect_obs$confirmed[,i], col = rgb(0.6,0.6,0.6, alpha=0.8))
}
lines(Rt_dataset$ode_perfect_obs$daily_confirmed, lwd=2)
# row_mean <- rowMeans(Rt_dataset$stoch_perfect_obs$confirmed, na.rm=TRUE)
# for (i in 1:nrow(daily_case_list$confirmed)){
#   daily_case_list$confirmed[i, is.na(daily_case_list$confirmed[i, ])] <- 0 
# }
row_mean <- rowMeans(daily_case_list$confirmed, na.rm=TRUE)
lines(row_mean, lwd=2, col=2)
```


```{r eval = F}
# tstamp <- format(Sys.time(), "%Y%m%dT%H%M%S")
# sim_res <- list(daily_infected = daily_inf, 
#                 daily_symptom_onset = daily_ons,
#                 daily_confirmed = daily_con,
#                 ODE = out)
# saveRDS(sim_res, paste0("outputs/sim", tstamp, ".rds"))
# library(tidyverse)
# library(data.table)
# extract_daily_output <- function(df){
#   tend <- max(df$time) + 1
#   rowids <- sapply(1:tend, function(x) length(df$time[(x - df$time) > 0]))
#   return (df[c(1, rowids),]) # first row is the initial condition and included
# }
# df <- extract_daily_output(res[[1]])
# df %>% ggplot(aes(ceiling(time), I)) + geom_point() + geom_line()
```


### Random-walk Metropolis-Hastings MCMC 
```{r}
# data <- rgamma(50, shape = 2.4, rate = 0.4) + rnorm(50, mean = 0, sd = 0.1)
set.seed(1510)
data <- rgamma(10000, shape = 2.4, rate = 0.4)
lb <- c(1e-6, 1e-6)
ub <- c(1e3, 1e3)
# Prior distribution
log_prior <- function (params = NULL, lower = NULL, upper = NULL){
  if (!(length(params) == length(lower) & length(params) == length(lower))){
    stop("The number of elements must be the same for params, lower, and upper")
  }
  shape <- params[1]
  rate <- params[2]
  shape_prior = dunif(shape, lower[1], upper[1], log = T)
  rate_prior = dunif(rate, lower[2], upper[2], log = T)
  return(shape_prior + rate_prior)
}

log_lik <- function(params, data){
  if(length(params) != 2) {
    stop("Two elements are required for the param!")
  }
  return(sum(dgamma(data, shape = params[1], rate = params[2], log = TRUE)))
}

log_posterior <- function (params, data, lower, upper){
    return (log_lik(params, data))
}
# log_prior(params = c(2,3), lower = lb, upper = ub)
# log_posterior(params = c(2,3), data = data, lower = lb, upper = ub)

library(tmvtnorm)
run_mcmc <- function(data, 
                     startvalue,
                     iter, 
                     burnin = round(iter/2), 
                     sigma, 
                     lower,
                     upper,
                     show_progress_bar = TRUE){
  if (show_progress_bar) {
    pb <- txtProgressBar(min = 0, max = iter, style = 3)
  }
  update_step <- max(5, floor(iter/100))
  
  accept <- numeric(iter)
  npar <- length(startvalue)
  chain <- matrix(NA, nrow = iter, ncol = (npar+1)) # unknown parameters + likelihood store
  chain[1, 1:npar] <-  startvalue
  chain[1, (npar+1)] <- log_posterior(chain[1, 1:npar], data, lower, upper)
  
  for (i in 2:iter) {
    # proposal function
    if (show_progress_bar && i%%update_step == 0) {
      setTxtProgressBar(pb, i)
    }
    # random walk
    proposal <- rmvnorm(1, mean = chain[(i-1), 1:npar], sigma = sigma)
    posterior_proposal <- -Inf
    if ((proposal[1] > 0) & (proposal[2] > 0)) {
      posterior_proposal <- log_posterior(proposal, data, lower, upper)
    }
    # cat( proposal, "\n" )
    # # acceptance probability
    alpha <- exp(posterior_proposal - chain[(i-1), (npar+1)])
    # cat( "\niter = ", i, ", log density = ", posterior_proposal, ", log density prev = ", chain[ (i-1), (npar+1) ], ", alpha = ", alpha )
    if (!is.finite(alpha)){ 
      alpha <- 0
    }
    if (runif(1) < alpha){
      chain[i, 1:npar] <- proposal
      chain[i, (npar+1)] <- posterior_proposal
      accept[i] <- 1
      # cat( "\niter = ", i, ", proposal accepted:", chain[ i,],"\n")
    } else { 
      chain[i, 1:npar ] <- chain[(i-1), 1:npar]
      chain[i, (npar+1)] <- chain[(i-1), (npar+1)]
    } 
  } 
  # cat( "\nAcceptance ratio = ", sum(accept[(burnin + 1):iter]) / (iter - burnin), "\n" )
  list(theta = chain[(burnin + 1):iter, (1:npar)],
       log_posterior = chain[(burnin + 1):iter, (npar+1)], 
       acceptance_ratio = sum(accept[(burnin + 1):iter]) / (iter - burnin))
}

run_mcmc2 <- function(data, 
                     startvalue,
                     iter, 
                     burnin = round(iter/2), 
                     sigma, 
                     lower,
                     upper,
                     show_progress_bar = TRUE){
  if (show_progress_bar) {
    pb <- txtProgressBar(min = 0, max = iter, style = 3)
  }
  update_step <- max(5, floor(iter/100))
  
  accept <- numeric(iter)
  npar <- length(startvalue)
  chain <- matrix(NA, nrow = iter, ncol = (npar+1)) # unknown parameters + likelihood store
  chain[1, 1:npar] <-  startvalue
  chain[1, (npar+1)] <- log_posterior(chain[1, 1:npar], data, lower, upper)
  
  for (i in 2:iter) {
    # proposal function
    if (show_progress_bar && i%%update_step == 0) {
      setTxtProgressBar(pb, i)
    }
    # random walk

    proposal <- rtmvnorm(1, mean = chain[(i-1), 1:npar], sigma = sigma,
                         lower = lower, upper = upper)
    # proposal <- rmvnorm(1, mean = chain[(i-1), 1:npar], sigma = sigma)
    posterior_proposal <- -Inf
    if ((proposal[1] > 0) & (proposal[2] > 0)) {
      posterior_proposal <- log_posterior(proposal, data, lower, upper)
    }
    # cat( proposal, "\n" )
    # # acceptance probability
    alpha <- exp(posterior_proposal - chain[(i-1), (npar+1)])
    # acceptbance probability
    old <- chain[(i-1), 1:npar]
    prob_new_given_old <-
      dtmvnorm(proposal, mean = old, sigma = sigma, lower = lower, upper = upper)
    prob_old_given_new <-
      dtmvnorm(old, mean = as.vector(proposal), sigma = sigma, lower = lower, upper = upper)
    alpha <- exp(posterior_proposal - chain[(i-1), (npar+1)]) * prob_old_given_new / prob_new_given_old
    # cat( "\niter = ", i, ", log density = ", posterior_proposal, ", log density prev = ", chain[ (i-1), (npar+1) ], ", alpha = ", alpha )
    if (!is.finite(alpha)){ 
      alpha <- 0
    }
    if (runif(1) < alpha){
      chain[i, 1:npar] <- proposal
      chain[i, (npar+1)] <- posterior_proposal
      accept[i] <- 1
      # cat( "\niter = ", i, ", proposal accepted:", chain[ i,],"\n")
    } else { 
      chain[i, 1:npar ] <- chain[(i-1), 1:npar]
      chain[i, (npar+1)] <- chain[(i-1), (npar+1)]
    } 
  } 
  # cat( "\nAcceptance ratio = ", sum(accept[(burnin + 1):iter]) / (iter - burnin), "\n" )
  list(theta = chain[(burnin + 1):iter, (1:npar)],
       log_posterior = chain[(burnin + 1):iter, (npar+1)], 
       acceptance_ratio = sum(accept[(burnin + 1):iter]) / (iter - burnin))
}


mcmc <- run_mcmc(data = data, 
                     startvalue = c(1, 1),
                     iter =  5e5, 
                     burnin = round(iter/2), 
                     sigma = diag(length(lb)), 
                     lower = lb,
                     upper = ub,
                     show_progress_bar = TRUE)

quantile(mcmc$theta[,1])
quantile(mcmc$theta[,2])
plot(mcmc$theta[,1], type = "l")

```