---
title: "Week 6-7 Model"
author: "Bill Foote"
date: "`r Sys.Date()`"
output: html_document
---

```{r , include=FALSE}
knitr::opts_chunk$set( echo = TRUE, warning = FALSE, message = TRUE )
library(deSolve)
library(tidyverse)
library(bbmle)
library(mvtnorm)
library(plotly)
options( digits=5, scipen=9999)
```

## Operational Forced Outages and Restoration

The reliability of many processes are preyed upon by unexpected forced outages. Restoration of regular operations is a time and resource allocation problem. We use a version of the predator-prey model to model to a first approximation the timing of restoration of regular operations.

Our stock variables are hours of operation $H$, hours of forced outage $F$, and hour of restored operations $R$, all across time $t$, here in weeks.  

$$
\begin{align}
\frac{dH}{dt} & = -\beta H F \\
\frac{dF}{dt} &  = \beta H F - \nu F \\
\frac{dR}{dt} &  = \nu F
\end{align}
$$

Here is the [Vensim model](models/hours-forced-restored.mdl) of this formulation.

![**Fig. 1. Hours of operation, forced outage, and restored hours.** ](images/hours-forced-restored.png)
We can export the `xmile` file to recreate this model using **readsdr**`::read_xmile()`. From the file we pull the function rendered in R.

```{r}
options( digits=5, scipen=9999999)
#
hfr <- function(t,x,parms){
  H <- x[1]
  F <- x[2]
  R <- x[3]
  with(as.list(parms),
       {
         dH <- -beta*H*F
         dF <- beta*H*F - nu*F
         dR <- nu*F
         out <- c(dH,dF,dR)
         list(out)
       })
}
#
H0 <- 1e4
parms <- c(H0=H0,beta=0.0001, nu = 1/7)
times <- seq(0,30,0.1)
y0 <- c(H0,1,0)
out_HFR <- ode(y=y0, times, hfr, parms)
colnames(out_HFR) <- c("time","H","F","R")
n_times <- length(times)
t <- rep( 1:n_times, 3)
states <- c( rep( "Hours Operated", n_times), rep( "Forced Outage", n_times), rep( "Restored", n_times))
hours <- c( out_HFR[,2],
            out_HFR[,3],
            out_HFR[,4]
            )
d_plot_out <- tibble(
  weeks=out_HFR[,1],
  H=out_HFR[,2],
  F=out_HFR[,3],
  R=out_HFR[,4]
)
p <- d_plot_out |>
  ggplot( aes( x=weeks, y=H ) ) +
  geom_line( color = "blue") +
  geom_line( aes( x=weeks, y=F ), color="red" ) +
  geom_line( aes( x=weeks, y=R ), color="black" ) +
  labs( title="Hours of operation, forced outage, restored hours",
        x="weeks",
        y="hours"
  )
ggplotly( p )
#
```

We recently measured the number of hours per week restored.

```{r}
set.seed(42)
plant_0 <- c(0, 4, 10, 15, 18, 21, 31, 51, 53, 97, 125, 183, 292, 390, 448,
            641, 771, 701, 696, 867, 925, 801, 580, 409, 351, 210, 113, 65, 
            52, 51, 39, 33)
plant <- abs(plant_0*(1 + rlogis(length(plant_0))))
cat(plant)
plant <- c( 0,1,11,5,13,0,12,67,129,353,327,170,472,703,1079,2211,705,9, 1939,185,2161,1286,403,513,1636,236,124,20,100,215,42,18)
cumplant <- cumsum(plant)
weeks <- 0:31
d_plot_obs <- tibble(
  weeks=weeks,
  cumplant=cumplant
)
d_plot_obs |>
  ggplot( aes( x=weeks, y=cumplant) ) +
  geom_point( color = "darkred") +
  labs( title="Cumulative Restored Hours Forced Out", 
        x="Weeks", 
        y="Cumulative Restoration (hours)")
#
```

Now comes the hard work of estimating the parameters $\beta \,\text{and}\, \nu$. Our strategy here will be to keep the estimation entirely in \R for the time being. To that end we will use Laplace's quadratic approximation routine. Effectively this means we will maximize the log posterior likelihood of observing restoration data given that $\beta \,\text{and}\, \nu$ range fairly widely as priors.

The estimator will work if we transform the parameters, as well as input initial values to the ODE, into log scores. When we pull the parameters back we will take advantage of obvious logistic structure of the ODE and of the data. Thus we will use the `plogis()` instead of `exp()`. We will make good use of the `mle2()` function from the **bbmle** package.

```{r}
# Form likelihood function Pr( y | beta, nu)Pr(beta)Pr(nu)
hfrLL <- function(lbeta, lnu, logH0, logF0) {
  parms <- c(beta=plogis(lbeta), nu=plogis(lnu))
  y0 <- c(H=exp(logH0), F=exp(logF0), R=0)
  out <- ode(y=y0, weeks, hfr, parms)
  SD <- sqrt(sum( (cumplant-out[,4])^2)/length(weeks) )
  -sum(dnorm(cumplant, mean=out[,4], sd=SD, log=TRUE))
}
#
lbeta <- qlogis(1e-5); lnu <- qlogis(0.2) 
logH0 <- log(1e6); logF0=log(1)
hfrLL(lbeta, lnu, logH0, logF0)
# priors are all equally likely
# 
# minimize negative-log-likelihood
fit1 <- mle2(hfrLL, 
            start=list(lbeta=qlogis(1e-5), 
                       lnu=qlogis(.2), 
                       logH0=log(1e6), logF0=log(1) ),  
            method="Nelder-Mead",
            control=list(maxit=1E5,trace=0),
            trace=FALSE)
summary(fit1)
#
```

We can try to deviate from these finely tuned start values, but would do so very hesitantly. Our Vensim experience shows a very high sensitivity of ODE solutions to very small changes in $\beta$ and $\nu$. The most sensitive component of the model is the linkage between $H$ and $R$, namely, the non-linear integration of $F$.  

We calculate the ODE solution with estimated parameters as the $\mu$ of  this regression.

```{r}
par_mean_1 <- fit1@coef
par_varcov_1 <- fit1@vcov
#
samples_fit1 <- rmvnorm(1000, par_mean_1, par_varcov_1 ) 
summary(exp(samples_fit1))
#
hfrMU <- function(par_mean) {
  lbeta <- par_mean[[1]]; lnu <- par_mean[[2]] 
  logH0 <- par_mean[[3]]; logF0 <- par_mean[[2]]
  parms <- c(beta=plogis(lbeta), nu=plogis(lnu))
  y0 <- c(H=exp(logH0), F=exp(logF0), R=0)
  out <- ode(y=y0, weeks, hfr, parms)
  out[,4]
}
mu_1 <- hfrMU( par_mean_1 )
sigma_1 <- exp((par_varcov_1[3,3])^0.5)
#
samples_fit1 <- as_tibble(samples_fit1) |> 
  mutate(prediction = rnorm( 1000, mu_1, sigma_1) )
#
cumplant_pred <- samples_fit1$prediction
quantile(cumplant_pred, 0.045)
quantile(cumplant, 0.045)
#
quantile(cumplant_pred, 0.5)
quantile(cumplant, 0.5)
#
quantile(cumplant_pred, 0.945)
quantile(cumplant, 0.945)
```
The lower tail does not fitting as well as we might want, but the upper tail does. The median deviates by a wide margin. That divergence is a characteristic of the interactive relationship between hours of operation and forced outage, and non-commutative forced outage and restoration.

The plot might then possibly be misleading.

```{r}
#
theta <- as.numeric(c(plogis(coef(fit1)[1:2]),
                      exp(coef(fit1)[3:4])) )
parms <- c(beta=theta[1], nu = theta[2])
times <- seq(0,30,0.1)
y0 <- c(theta[3],theta[4],0)
out_HFR1 <- ode(y=y0, times, hfr, parms)
colnames(out_HFR1) <- c("time","H","F","R")
d_plot <- out_HFR1 |> as_tibble()

p <- ggplot( data=d_plot, aes( x=time, y=R )) +
  geom_line( color = "blue" ) +
  geom_point( data = d_plot_obs, aes( x=weeks, y=cumplant ), color = "red") +
  labs( title="Cumulative Restored Hours Forced Out", 
        x="Weeks", 
        y="Cumulative Restoration (hours)")
p
```


One more attempt at inference are the HPDIs (High Density Probability Intervals).

```{r}
library( tidybayes )
library( rethinking )
library( coda )
prob <- c(0.89)
#
# a helper function to make a table of lower and upper
# estimates of highest posterior density intervals
# for estimated parameters and predictions
#
hpdi_tbl <- function( samples, prob=0.89){
  column <- colnames(samples)
  x <- sapply(column, 
      function(c) HPDinterval(mcmc(eval(
        parse(text = concat("samples$", c)))), 
              prob = prob))
  rownames(x) <- c( concat("|lower ", prob), 
                    concat("upper ", prob, "|") )
  return( t(x) )
}
hpdi_samples <- hpdi_tbl( samples_fit1, prob=0.89)
hpdi_samples
# filter samples <= 5000
hpdi_samples <- samples_fit1 |> filter ( prediction > 5000 ) |>
  hpdi_tbl( prob=0.89)
```

That interval for _prediction_ is breath-takingly alarming! 

```{r}
set.seed(4242)
samples_fit1 |> filter( prediction > 5000 ) |>
  ggplot(aes(x = prediction)) +
  geom_histogram(color = "blue", alpha=0.7) +
  geom_vline( xintercept=hpdi_samples["prediction",1], 
              color="red") +
  geom_vline( xintercept=hpdi_samples["prediction",2], 
              color="red") +
  labs(title = "Posterior predictive distribution",
       subtitle = "89% High Posterior Density Interval",
       x = "predicted restored hours",
       y = "frequency of predicted outcomes"
       ) +
  theme(panel.grid = element_blank())
```

However we slice this data, we will see there are more than one set of behaviors embedded in this highly non-linear soup.

Ever being optimistic, we try to update our estimates using this routine.

```{r}
#
# Another fit piggy-backing on our first try `fit1`
# 
fit2 <- mle2(hfrLL, 
             start=as.list(coef(fit1)),
             fixed=list(logH0=coef(fit1)[3], 
                        logF0=coef(fit1)[4]), 
             method="Nelder-Mead",
             control=list(maxit=1E5,trace=2),
             trace=TRUE)
summary( fit2 )
par_mean_2 <- fit2@coef
par_varcov_2 <- fit2@vcov
library(mvtnorm)
samples_fit2 <- as.data.frame(rmvnorm(1000, par_mean_2, par_varcov_2 ) )
summary(exp(samples_fit2))

#
hfrMU_2 <- function(par_mean) {
  lbeta <- par_mean[[1]]; lnu <- par_mean[[2]] 
  logH0 <- coef(fit1)[3]; logF0 <- coef(fit1)[4]
  parms <- c(beta=plogis(lbeta), nu=plogis(lnu))
  y0 <- c(H=exp(logH0), F=exp(logF0), R=0)
  out <- ode(y=y0, weeks, hfr, parms)
  out[,4]
}
mu_2 <- hfrMU_2( par_mean_2 )
sigma_2 <- (par_varcov_1[3,3])^0.5
#
samples_fit2 <- as_tibble(samples_fit2) |> 
  mutate(prediction = rnorm( 1000, mu_1, sigma_1) )
#
# HPDI next
# 
hpdi_samples <- samples_fit2 |> filter ( prediction > 5000 ) |>
  hpdi_tbl( prob=0.89)
hpdi_samples
```

That interval for _prediction_ is still takes our statistical breath away! 

```{r}
set.seed(4242)
samples_fit1 |> filter( prediction > 5000 ) |>
  ggplot(aes(x = prediction)) +
  geom_histogram(color = "blue", alpha=0.7) +
  geom_vline( xintercept=hpdi_samples["prediction",1], 
              color="red") +
  geom_vline( xintercept=hpdi_samples["prediction",2], 
              color="red") +
  labs(title = "Posterior predictive distribution",
       subtitle = "beta, nu fit: 89% High Posterior Density Interval",
       x = "predicted restored hours",
       y = "frequency of predicted outcomes"
       ) +
  theme(panel.grid = element_blank())
```
#

```{r}
# and this initial plot
# 
theta2 <- as.numeric(c(plogis(coef(fit2)[1:2]),
                      exp(coef(fit1)[3:4])) ) # fit2 has only two parms 
parms <- c(beta=theta2[1], nu = theta2[2])
times <- seq(0,30,0.1)
y0 <- c(theta2[3],theta2[4],0)
out_HFR2 <- ode(y=y0, times, hfr, parms)
colnames(out_HFR2) <- c("time","H","F","R")
d_plot <- out_HFR2 |> as_tibble()
#
p <- ggplot( data = d_plot_obs, aes( x=weeks, y=cumplant)) +
  geom_quantile(data = d_plot_obs, method = "rqss", lambda = 0.6, color = "red", size=1.1) +
  geom_point( data = d_plot_obs, aes( x=weeks, y=cumplant ), color = "red") +
  geom_line( data = d_plot, aes( x=time, y=R), color = "gray30", size=1.5) +
  labs( title="Cumulative Restored Hours Forced Out", 
        subtitle="5%-95% quantile predictive intervals (red)",
        x="Weeks", 
        y="Cumulative Restoration (hours)")
p
#
```

That looks pretty tight!

We still miss the top, and skim the quantile interval at the high end for that matter. But when does the model predict that we are half-way through the restoration? It still seems that the answer lies somewhere between weeks 15 and 20 for both fits. So perhaps we have beaten our model with too large a hammer and have our answer already.

