---
title: Model-Based Nonlinear Prediction
author: 36-467/667, Fall 2020
date: 3 December 2020 (Lecture 25)
output: slidy_presentation
bibliography: locusts.bib
---

```{r, include=FALSE}
library(knitr)
opts_chunk$set(size="small",background="white", highlight=FALSE,
               cache=TRUE, autodep=TRUE,
               tidy=TRUE, tidy.opts=list(comment=FALSE), warning=FALSE, message=FALSE,
               echo=FALSE)
set.seed(2020-12-03)
```


\newcommand{\Expect}[1]{\mathbb{E}\left[ #1 \right]}
\newcommand{\Prob}[1]{\mathbb{P}\left( #1 \right)}
\newcommand{\Var}[1]{\mathrm{Var}\left[ #1 \right]}


## Housekeeping

> - HW12 is now due _tomorrow_, Friday 4 December, at 6 pm
> - HW13 will be due Friday 11 December at 6 pm
> - I promise I will in fact grade all the after-class exercises

## Model-Based Forecasting

- What?
    + Build a parametric, or at least highly structured, model of the dynamics and the measurement process
	    * There's an internal state $S(t)$ which evolves according to some model (usually a Markov model)
		* The model for the dynamics of $S(t)$ is where we try to build in the "physics" (or biology or economics or ecology...)
		* There're observables $X(t)$ which are functions of $S(t)$ plus noise
		* The observation model is where we try to build in our knowledge of the measurement process
		* Both $S(t)$ and $X(t)$ might be big vectors, perhaps with spatial structure, etc.
	+ Estimate the parameters which show up in the process and observation models
	+ Estimate the state based on the history of observations
	+ Run the model it forward
	+ Run the model forward many times to get at uncertainty
- Why?
    + Non-stationarity and trends
	+ Efficient use of information in the data
	+ Checking whether a model is right

## Fixed Parameters, Fixed State

- Say we know the current state $S(t)$ and the parameters $\theta$
- We can then run the model forward (i.e., simulate it)
    + Get a sequence of future states $S(t+1), S(t+2), \ldots S(t+T)$
	+ With the states get a sequence of observations $X(t+1), X(t+2), \ldots X(t+T)$
- But even with fixed state and parameters, the model will usually
produce multiple possible futures
- $\Rightarrow$ Run the model many times to see the distribution of possible
trajectories
    + Just plotting all the $X$ trajectories can be very informative, but it may also be just a "spaghetti plot"


## Fixed Parameters, Unknown State

- The state is often not known exactly
    + What does the monkey have in mind for its hand?
	+ How many lynxes and hares are there?
	+ How many people are really in compartments $S$, $I$, $R$?
- We can estimate the state
    + E.g., $\Expect{S(t)|X(1), \ldots X(t)}$
- But there's some uncertainty around this estimate
- Ideally we'd use $\Prob{S(t)|X(1), \ldots X(t)}$
- **Ensemble forecasting**: Draw _samples_ from the state distribution and run each one forward
    + This is what weather forecasting does
	+ Today's predicted probability of snow in Pittsburgh on 21 December 2020 = what fraction of runs in the ensemble, started today, have evolved to snow on 21 December 2020?
    + This should remind you of the particle filter...
	
## Learning About the State as We Go

- Start with an ensemble forecast:
    + Many guesses initially about $S(t)$
	+ Extrapolate forward to guesses about $S(t+1), S(t+2), \ldots$
	+ Leading to guesses about $X(t+1), X(t+2), \ldots$, which are the predictions
- Now we actually observe $X(t+1)$
- **Data assimilation**: Use $X(t+1)$ to update our guesses for $S(t+1)$ and re-do the forecasting going forward
    + This is _exactly_ what the particle filter does
	+ Lots of numerical schemes for data assimilation to try to run it faster at some approximation cost [@Big-Laplace-Gauss-filter-paper]


## Unknown Parameters, Unknown State

- The parameters are usually not known exactly
- Different parameters $\Rightarrow$ different predictions
- How to deal with this?
	1. "Propagation of error"
    2. Bootstrap-type approaches


## "Propagation of error"

- If we have a parameter estimate $\hat{\theta}$ and standard error $\hat{se}$,
then estimate $f(\theta)$ as $f(\hat{\theta})$
- An approximate standard error for $f(\theta)$ would be $\left|\frac{d f}{d\theta}\right| \hat{se}$
    + Approximation becomes exact as $\hat{se} \rightarrow 0$
- Extends to multiple arguments: we estimate $f(\theta_1, \theta_2, \ldots \theta_p)$ by $f(\hat{\theta}_1, \hat{\theta}_2, \ldots \hat{\theta}_p)$, with
standard error
\[
\approx \sqrt{\sum_{i=1}^{p}{\left(\frac{\partial f}{\partial \theta_i}\right)^2 \hat{se}_i^2} + 2\sum_{i=1}^{p}{\sum_{j=i+1}^{p-1}{\frac{\partial f}{\partial \theta_i}\frac{\partial f}{\partial \theta_j} \mathrm{Cov}[\hat{\theta}_i, \hat{\theta}_j]}}}
\]
    + Just a Taylor expansion; see [handout on website](http://www.stat.cmu.edu/~cshalizi/dst/20/lectures/25/propagation-of-error.pdf)
- Lets you propagate uncertainty from parameters to functions of parameters (like predictions) just by taking derivatives
    + Taking derivatives in a big, complicated model can be hard!
- If you're willing to believe $\hat{\theta}$ has a Gaussian distribution, then $f(\hat{\theta})$ does too (approximately)
    + This is called the **delta method**


## Bootstrap approaches to unknown parameters

- Bootstrap the training data
    + Model-based bootstraps recommended when dealing with non-stationary data
- Estimate parameters on each bootstrap run
- For each parameter estimate, run your forecast as if those parameter values were fixed
    + May want to run an ensemble forecast for each parameter value
- Aggregate over all simulation runs




## Model Checking and Diagnostics

- Before relying on the model's predictions, check the model
    + Usually a good idea before wasting time on things like confidence sets of parameters
- You remember a lot of "diagnostics" for linear regression models from your regression class
- For dynamic models, the basic diagnostic is that
\[
\frac{X(t) - \Expect{X(t)|X(1), \ldots X(t-1)}}{\sqrt{\Var{X(t)|X(1), \ldots X(t-1)}}}
\]
should be mean 0 with variance 1
- If $X(t)$ always has a continuous distribution given the past, say with
CDF $F_t(r) \equiv \Prob{X(t) \leq r| X(1), \ldots X(t-1)}$, then the sequence $Q(t) = F_t(X(t))$ should be uniformly distributed on $[0,1]$ and IID
    + Similar but more complicated rules for discrete or discontinuous distributions
- The basic idea: once you've removed or adjusted for everything the model predicts, what's left should be pure noise

## Model Checking by Comparisons

- There may be multiple models (perhaps variations on a theme)
- One way to chose would be a pure prediction contest, i.e., try to get the best forecast of a common set of variables
- May be more illuminating to look for _qualitative_ distinctions which you can compare to the data (are there cycles? a hard lower or upper bound on some process?)
- May be even more illuminating to look for new phenomena/variables
- This gets us in to philosophy of science [@Glymour-theory-and-evidence;@Mayo-error] but sometimes that can't be helped
- Alternatively: present results from all the models which you can't rule out



## Summing Up

> 1. With fixed state and fixed parameters, run the model forward multiple times to get multiple simulated trajectories
> 2. With unknown state but fixed parameters, get an ensemble of states and go to (1)
>     + Updating the state as we see more = data assimilation, often through the particle filter
> 3. With unknown state and unknown parameters, get an ensemble of parameters and go to (2)
> 4. Try to do model checking _before_ investing lots of time into model-based forecasting


## Backup: Visualizing Many Simulation Runs

- Every run of the simulation gives us a trajectory, which is a function of time;
that's why I write it $X(t)$
- The easiest way to display all these at once is just to plot them one atop the other
- We've been doing this without remarking on it, as with the SIR model
we looked at in [Lecture 20](http://www.stat.cmu.edu/~cshalizi/dst/20/lectures/20/lecture-20.html)

```{r}
# Taken from Lecture 20 code [http://www.stat.cmu.edu/~cshalizi/dst/20/lectures/20/lecture-20.Rmd]

sim.SIR <- function(n, beta, gamma, s.initial=n-1, i.initial=1, r.initial=0, T) {
    stopifnot(s.initial+i.initial+r.initial == n)
    states <- matrix(NA, nrow=3, ncol=T)
    rownames(states) <- c("S", "I", "R")
    states[,1] <- c(s.initial, i.initial, r.initial)
    for (t in 2:T) {
        contagions <- rbinom(n=1, size=states["S",t-1],
                             prob=beta*states["I",t-1]/n)
        removals <- rbinom(n=1, size=states["I",t-1],
                           prob=gamma)
        states["S",t] <- states["S",t-1] - contagions
        states["I",t] <- states["I",t-1] + contagions - removals
        states["R",t] <- states["R",t-1] + removals
    }
    return(states)
}

# Run the simulation b times, with the given settings, and plot them all
# together
# Inputs: number of replicates (b)
  # population size (n)
  # infection parameter (beta)
  # removal parameter (gamma)
  # duration of simulation (T)
  # optional further arguments to sim.SIR (...)
# Output: the 3D array of simulation runs, invisibly
# Side-effect: the plot is made
# Calls: sim.SIR()
plot.many.SIR.runs <- function(b, n, beta, gamma, T, ...) {
    # Run the simulation b times, with the same parameters
    many.runs <- replicate(b, sim.SIR(n=n, beta=beta, gamma=gamma, T=T))
    # Create an empty plot with the right limits, labels, etc.
    plot(0, type="n",
         xlab="t", ylab="People (count)",
         main=paste("n=", n, ", beta=", beta, ", gamma=", gamma, sep=""),
         ylim=c(0,n),
         xlim=c(1,T))
    # Now plot all 3 state variables for each simulation run, distinguishing
      # variables by line-type and the runs by color
    for (run in 1:b) {
        lines(many.runs["S",,run], lty="solid", col=run)
        lines(many.runs["I",,run], lty="dashed", col=run)
        lines(many.runs["R",,run], lty="dotted", col=run)
    }
    legend("right", legend=c(expression(S(t)),
                        expression(I(t)),
                        expression(R(t))),
       lty=c("solid", "dashed", "dotted"),
       cex=0.5)
    invisible(many.runs)
}
```

```{r}
plot.many.SIR.runs(b=10, n=1000, beta=1e-1, gamma=1e-2, T=800)
```


## Backup: Visualizing Many Simulation Runs

- With many runs, this can get kind of messy ("spaghetti plots"):
```{r}
sir.1000.runs <- plot.many.SIR.runs(b=1000, n=1000, beta=1e-1, gamma=1e-2, T=800)
```

- There are many options for summarizing things like this in more comprehensible and visually-pleasing ways.  Common ones:
    1. Mean $\pm 2$ standard deviations, computed at each value of $t$
	2. Median plus high and low quantiles, computed at each value of $t$
	3. "Functional boxplots"
	

## Backup: Visualizing Many Simulation Runs

#### Mean $\pm 2$ standard deviations, at each value of $t$

```{r}
mean.I <- apply(sir.1000.runs["I",,],1,mean)
sd.I <- apply(sir.1000.runs["I",,],1,sd)
plot(1:length(mean.I), mean.I,
     xlim=c(1,length(mean.I)),
     ylim=c(0,1000),
     ylab="I(t)",
     xlab="t",
     main="Summarizing I(t) across simulations",
     sub="Mean +- 2 SD",
     type="l")
lines(1:length(mean.I), mean.I-2*sd.I, lty="dashed")
lines(1:length(mean.I), mean.I+2*sd.I, lty="dashed")
legend("topright", legend=c("Mean", "+-2 SD"), lty=c("solid", "dashed"))
```

Comments:

- Fine idea if everything's Gaussian
- Still informative in many other circumstances but may also give un-"physical" (or un-"biological") limits, like the negative values here
    + OTOH this does tell us the variance at those times is really big


## Backup: Visualizing Many Simulation Runs

#### Median and high and low quantiles, at each value of $t$

```{r}
median.I <- apply(sir.1000.runs["I",,],1,median)
quantiles.I <-  apply(sir.1000.runs["I",,],1,quantile,probs=c(0.05,0.25,0.75,0.95))
plot(1:length(median.I), median.I,
     xlim=c(1,length(median.I)),
     ylim=c(0,1000),
     ylab="I(t)",
     xlab="t",
     main="Summarizing I(t) across simulations",
     sub="Median and high/low quantiles",
     type="l")
lines(1:length(median.I), quantiles.I["5%",], lty="dotted")
lines(1:length(median.I), quantiles.I["25%",], lty="dashed")
lines(1:length(median.I), quantiles.I["75%",], lty="dashed")
lines(1:length(median.I), quantiles.I["95%",], lty="dotted")
legend("topright", legend=c("Median", "50% interval", "90% interval"),
       lty=c("solid", "dashed", "dotted"))
```

- Comments:
    + Common to pick say 25% and 75% quantiles, getting a 50% probability interval
    + Or 5% and 95% to get a 90% probability interval
    + Or, as in the plot, both (using shading or color to reinforce the sense of moving away from the median)
- Drawbacks of both the time-by-time quantile and time-by-time mean approach:
    + The typical trajectory isn't, necessarily, near the typical value at all times!
	+ These "fixed-time" summaries can give a misleading impression of the size and timing of extremes, if the timing of extremes is itself variable
	+ E.g., in an epidemic model like this, where the time of peak infection varies, the fixed-time quantiles can miss the peak [@Juul-et-al-fixed-time-descriptive-statistics]
	


## Backup: Visualizing Many Simulation Runs

#### "Functional boxplots" [@Sun-Genton-functional-boxplots] or "curve boxplots" [@curve-boxplots]

- Look at the distance between entire curves
- Try to find the "deepest" curve, the one most nearly at the center of the distribution
- Order curves outward by depth from that deepest curve
- Stop when you have, say, 90% of the curves
- Draw the resulting tube around the deepest curve, and highlight the outliers falling beyond it (as in a regular boxplot)


```{r}
library(fda) # Functional data analysis package, for functional boxplots command
invisible(fbplot(sir.1000.runs["I",,], xlab="t", ylab="I(t)", ylim=c(0,1000)))
```





## References
