---
title: Compartment Models
author: 36-467/667, Fall 2020
date: 12 November 2020 (Lecture 21)
output: slidy_presentation
bibliography: locusts.bib
---

\newcommand{\Expect}[1]{\mathbb{E}\left[ #1 \right]}
\newcommand{\Var}[1]{\mathrm{Var}\left[ #1 \right]}
\newcommand{\Binom}{\mathrm{Binom}}

```{r, include=FALSE}
# General set-up options
library(knitr)
opts_chunk$set(size="small", background="white",
               cache=TRUE, autodep=TRUE,
               tidy=FALSE, warning=FALSE, message=FALSE,
               echo=FALSE)
```

## Previously

- [Markov models](http://www.stat.cmu.edu/~cshalizi/dst/20/lectures/18/lecture-18.html): The future of the process is independent of its past, _given_ its present state
    + Everything boils down to the probabilities or rates of transitions between states
- [Inference for Markov models](http://www.stat.cmu.edu/~cshalizi/dst/20/lectures/19/lecture-19.html): Based on the likelihood, _or_ on other ways of assessing fit to the data
- [Epidemic models](http://www.stat.cmu.edu/~cshalizi/dst/20/lectures/20/lecture-20.html): Every individual in the population is in one of a limited number of states of health; transitions between states happen at rates which depend on how many individuals are in each state

## Compartment Models

- Many sciences deal with populations, where each individual is in one of a limited number of states
    + Epidemiology (as we saw)
	+ Ecology
	+ Sociology (think: occupations, classes, education, racial/ethnic groups, religion, ...)
	+ Demography (age, sex, ethnicity, religion, ...)
	+ Chemical kinetics (molecular "species")
- All these use **compartment models**
    + Compartments = status for individuals
	+ Each individual is in 1 and only 1 compartment at any time
        * Compartments are often intersections, e.g., "Black males aged 61--65"
	+ Rates of transitions from one compartment to another, often depending on occupancy of compartments

## Basic Definitions and Assumptions {.smaller}

- Start with a "closed" population, constant size $n$
- Population is divided into $r$ **compartments**
- State of population at time $t$,  $X(t)$, is a vector of counts,
\begin{equation}
  X(t) \equiv [X_1(t), X_2(t), \ldots X_r(t)]
\end{equation}
    + $X_i(t)=$ how many individuals are in compartment $i$ at time $t$
- $X(t)$ changes because individuals move from compartment to compartment, randomly
- $p_{ij}=$ probability that individual in compartment $i$ will move to compartment $j$ at the next time step, the $i\rightarrow j$ **transition rate**
	+ If $p_{ij} =0$ the $i\rightarrow j$ transition is **forbidden** (or **suppressed**)
    + $p_{ii}=$ probability of staying in compartment $i$
- _Assumption_: $p_{ij}$ is a function of the _current_ state $X(t)$ and not the past or future states
    + $p_{ij}$ might be constant for some $i, j$
- _Assumption_: All individuals in the population make _independent_ transitions, given $X(t)$
    + Transitions from compartment $i$ are independent and identically distributed, given $X(t)$
- These assumptions could be wrong!

## Consequences of the Assumptions {.smaller}

- $Y_{ij}(t+1) =$ number of individuals transitioning from compartment $i$ to compartment $j$ at time $t+1$, the **flow** or **flux** from $i$ to $j$
- By the assumptions,
\begin{eqnarray}
  Y_{ij}(t+1) & \sim & \Binom(X_i(t), p_{ij}(X(t)))\\
  \Expect{Y_{ij}(t+1)|X(t)} & = & X_i(t) p_{ij}(X(t))\\
  \Var{Y_{ij}(t+1)|X(t)} & = & X_i(t) p_{ij}(X(t))(1-p_{ij}(X(t)))
\end{eqnarray}
- **Stock-flow consistency**:
\begin{eqnarray}
  \label{eqn:stock-flow} X_i(t+1) & = & X_i(t) + \sum_{j\neq i}{Y_{ji}(t+1) - Y_{ij}(t+1)}\\
  \text{(new level)} & = & \text{(old level)} + \sum{\text{(in flow)} - \text{(out flow)}}
\end{eqnarray}
- Implication: $X(t)$ is a Markov process
    + Because $Y_{ij}$s only depend on current $X(t)$, not past
    + The state space is finite, but very large


## Large-Population ($n$) Asymptotics

- We want to consider very large populations ($n\rightarrow\infty$) making lots of small steps over short time intervals $h$ ($h\rightarrow 0$) (so we go from time $t$ to $t+h$, not $t$ to $t+1$)
- Assume $p_{ij}(X) = h \rho_{ij}(X/n)$
    + Meaning that the _distribution_ over compartments matters, not the raw _numbers_ in each compartment
- Define $x(t) = X(t)/n$
- As $n\rightarrow\infty$ and $h\rightarrow 0$, we go to deterministic differential equations,
\[
\frac{dx_i}{dt} = \sum_{j \neq i}{\rho_{ji}(x(t)) - \rho_{ij}(x(t))}
\]
- There are some caveats --- see the handout

## Epidemic Models

- SIR model from [last time](http://www.stat.cmu.edu/~cshalizi/dst/20/lectures/20/lecture-20.html): 3 compartments, transition rates depend on counts in each compartment
    + We'd now write $X_S(t), X_I(t), X_R(t)$ instead of $S(t), I(t), R(t)$
	+ $Y_{SI}(t+1) = C(t+1)$ from last time, $Y_{IR}(t+1) = D(t+1)$
	+ $p_{SI} \propto X_I(t)$, $p_{IR} =$ constant
	+ All other transitions are forbidden
- SIS model: 2 compartments, all transitions both ways
- SIRS model: 3 compartments, allow an $R\rightarrow S$ transition but no $R\rightarrow I$ or $I \rightarrow S$
    + Often used for influenza

## Sociological Application: Modeling the Spread of Novelties

- Including: ideas, practices, technologies, use of products, beliefs, ...
- Some practices spread by personal contact
    + E.g., to learn how to make pots in a certain way, you (probably) need to learn from someone who already knows how to make that kind of pot
    + Others can be spread by artifacts (you see a new kind of pot and can guess how to make it), or get re-invented
- _Assume_ that everyone either has adopted the new practice (and can/will spread it), _or_ that they haven't
    + $\Rightarrow$ 2 compartments, call them $A$ (adopters) and $N$ (non-adopters)
- _Assume_ that the new practice spreads by personal contact
    + So $p_{NA}$ depends on $X_A$,  say $\propto X_A$
- _Assume_ that no one gives up the practice, so $p_{AN} = 0$


## Comments on the Assumptions

- Assuming that the practice has to spread by personal contact:
    + Some practices can be spread by _artifacts_
	    * Maybe you see a new-to-you type of pot and you can figure out how to make it
		* Ideas can be spread by books, videos, etc., as well as by people
	+ Some practices get re-invented independently
- Assuming that every adopter is a potential spreader:
    + Not necessarily true by any means
	    * Maybe you just want to make pots and aren't interested in teaching anyone
- Assuming the rate of adoption $\propto X_A$:
    + Presumes some sort of random mixing within the population
	+ Presumes no independent adoptions (could add a constant term to the rate)
	+ Presumes no "diminishing returns" to encountering lots of adopters
	+ Presumes presumes there aren't other influences on the decision to adopt or not
	    * Maybe un-common practices are _always_ less attractive, because it helps to coordinate with others
- Assuming no one gives up the new practice, $p_{AN} = 0$:
    + Plenty of people try new technologies, products, beliefs, etc., and then give them up!

$\Rightarrow$ Whether the assumptions are good matches to reality is something you need to investigate case by case, not presume!

## Comment: Analogy to Epidemic Models; What to Call This

- This is like an SIR model:
    + Adopters $A$ $\Leftrightarrow$ Infectious $I$
	+ Non-adopters $N$ $\Leftrightarrow$ Susceptibles $S$
	+ But set $p_{IR} = 0$ (i.e., in terms of last time, $\gamma=0$)
	+ $\therefore$ The infectious stay infectious forever
- _Some_ kind of analogy between the transmission of ideas/behavior and the transmission of disease is very old and very common
    + E.g., @Siegfried-germs-and-ideas is a careful working-out of the analogy by a French geographer in 1960
- The _mathematical_ models, and the realization that it matches epidemiological _models_, comes from the 1950s
- Talk about "diffusion of innovations" [@Rogers-diffusion-of-innovations] and it sounds good
- Talk about "viruses of the mind" [@Dawkins-viruses-of-the-mind] and it sounds bad
- It's hard come up with a really neutral name
    + The anthropologist Dan Sperber suggested "epidemiology of representations" [@Sperber-explaining], which I like but is a bit of a mouthful
	

## Consequences of the Model Assumptions

\begin{eqnarray}
  Y_{NA}(t+1) & \sim & \Binom(X_N(t), \alpha X_A(t))\\
  X_A(t+1) & = & X_A(t) + Y_{NA}(t+1)\\
  X_N(t+1) & = & X_N(t) - Y_{NA}(t+1)\\
  \Expect{X_A(t+1)|X_A(t)}  - X_A(t) & = & X_N(t)\alpha X_A(t)\\
  & = & \alpha X_A(t)(n-X_A(t))
\end{eqnarray}

- Notice: expected _increment_ to $X_A$ is quadratic in $X_A$, small if $X_A$ is small, small if $X_A$ is big, maximized when $X_A = n/2$
- Notice: if $X_A \ll n$, then the increment is roughly $\propto X_A$, so initially there's exponential growth
- "Slowly, and then all at once" but also "Nothing grows forever"

## Consequences of the assumptions

```{r, echo=TRUE}
alpha <- 1e-6
n <- 1e5
xa <- vector(length=200)
xa[1] <- 1
for (t in 2:length(xa)) {
    xa[t] <- xa[t-1] + alpha*xa[t-1]*(n-xa[t-1])
}
plot(1:length(xa), xa, xlab="t", ylab=expression(X[A](t)), type="l")
```

- This looks like the logistic function we use in logistic regression, $\frac{e^t}{1+e^t}$
- This is, in fact, called a "logistic growth curve"
    + Actually, the name "logistic regression" comes from this sort of logistic growth curve, not the other way around
	+ "Logistic" here from the sense of logistics, of supplying something [@Bacaer-on-population-dynamics]

## Consequences of the assumptions

```{r, echo=TRUE}
sim.diffusion <- function(n, alpha, xa.1=1, T) {
    xa <- vector(length=T)
    xa[1] <- xa.1
    for (t in 2:length(xa)) {
        xa[t] <- xa[t-1] + rbinom(n=1, size=n-xa[t-1], prob=alpha*xa[t-1])
    }
    return(xa)
}
```

```{r}
SI.sims <- replicate(100, sim.diffusion(1e5, 1e-6, 1, 200))
plot(1:length(xa), xa, xlab="t", ylab=expression(X[I](t)), type="l")
for (sim in 1:ncol(SI.sims)) {
    lines(1:nrow(SI.sims), SI.sims[,sim], col="grey")
}
lines(1:length(xa), xa)
legend("topleft", legend=c("simulations", "deterministic approximation"),
       col=c("grey", "black"), lty="solid")
```
  
## "Conversion to Islam in the Medieval Period" {.smaller}

- An interesting, important, and analytically-clean example from
  @Bulliet-conversion
- Islam is a "missionary" religion; its adherents (often) go out to seek converts
    + Some other religions aren't missionary ones and conversion is rare or even impossible (Judaism, Shinto, etc.)
- Conversion requires personal interaction; you are converted _by_ someone
    + In some religions you might get converted just by reading (and believing) a book, but in all religions that's rare
- During the initial expansion, say +630 to 900, _lots_ of people converted to Islam
    + Previously they were polytheists, Christians, Jews, Zoroastrians, Buddhists, Hindus, Manicheans, etc.
	+ There were many reasons for conversion --- there were certainly practical advantages if living under Islamic rule --- but part of it was definitely sincere religious conversions
	     * @Hodgson-venture, @Bulliet-view, @Ansary-destiny all try to present _why_ people might've wanted to convert to Islam
- @Bulliet-conversion used medieval biographical dictionaries to determine _when_ many different conversions happened, in different countries
    + The biographies give the subject's dates, and their formal name, in the style "A, son of B, son of C, son of D, ... son of Z"
	+ Converts typically assumed a Muslim name
	+ Tracking back to the latest non-Muslim name in a genealogy gives the number of generations to conversion
	    * Especially clear for Persian-speaking countries (modern Iran, Afghanistan, Uzbekistan, etc.) where the pre-Islamic names were linguistically very different from Arabic names, and also there weren't many Jews and Christians with Biblical names (which are also common in Islam)
	+ Assuming so many years per generation says (roughly) when the conversion happened
	+ Bulliet grouped by 25-year intervals
	


## "Conversion to Islam in the Medieval Period"

```{r}
bulliet <- read.csv("bulliet.csv")
```

```{r}
plot(Conversions ~ Period.End, data=bulliet,
     xlab="Years since introduction of Islam",
     ylab="Number of conversions per 25 year period",
     main="Persia/Iran", type="l")
```

## "Conversion to Islam in the Medieval Period"

```{r}
plot(cumsum(Conversions) ~ Period.End, data=bulliet,
     xlab="Years since introduction of Islam",
     ylab="Cumulative number of conversions",
     main="Persia/Iran", type="l")
```

- How could we fit the model to this data?

1. Regress the increments $X_A(t+1) - X_A(t)$ on $X_A(t)(n-X_A(t))$
2. Fit the curve $X_A(t)$ to the deterministic curve implied by the model
3. Use maximum likelihood

## "Conversion to Islam in the Medieval Period"

![](bulliet-fig-4.png)

## Summary

- Divide the population into mutually exclusive and jointly exhaustive compartments
- Assume both that transition rates between compartments only depend on the present state of the population, _and_ that individuals transition (conditionally) independently 
- Result: Markov model for the state of the population, with stable large-$n$ behavior
- Applications: epidemiology, spread of new ideas/practices, lots of other social processes, demography, ...


## References
