dobson-cyclone-example.qmd

(Adapted from @dobson4e §1.6.5)

### Data {#sec-dobson-cyclones-data}

{{< include macros.qmd >}}

The `cyclones` dataset in the `dobson` package (@tbl-cyclones-data) records the number of tropical cyclones in Northeastern Australia during 13 November-to-April cyclone seasons (more details in @dobson4e §1.6.5 and `help(cyclones, package = "dobson")`). @fig-dobson-cyclone-time-series graphs the number of cyclones (y-axis) by season (x-axis). Let's use $Y_i$ to represent these counts, where $i$ is an indexing variable for the seasons and $Y_i$ is the number of cyclones in season $i$.

### Exploratory analysis

Suppose we want to learn about how many cyclones to expect per season.

```{r}
#| tbl-cap: "Number of tropical cyclones during a season from November to April in Northeastern Australia"
#| label: tbl-cyclones-data
library(dobson)
library(dplyr)
data(cyclones)
library(pander)
pander(cyclones |> relocate(season, .before = everything()))
```

```{r}
#| label: fig-dobson-cyclone-time-series
#| fig-cap: "Number of tropical cyclones per season in northeastern Australia, 1956-1969"
library(ggplot2)
library(dplyr)
cyclones |>
  mutate(years = years |> factor(levels = years)) |>
  ggplot(aes(x = years, y = number, group = 1)) +
  geom_point() +
  geom_line() +
  xlab("Season") +
  ylab("Number of cyclones") +
  expand_limits(y = 0) +
  theme(axis.text.x = element_text(vjust = .5, angle = 45))
```

There's no obvious correlation between adjacent seasons, so let's assume that each season is independent of the others.

Let's also assume that they are identically distributed; let's denote this distribution as $P(Y=y)$ (note that there's no index $i$ in this expression, since we are assuming the $Y_i$s are identically distributed). We can visualize the distribution using a bar plot (@fig-cyclones-bar-plot). @tbl-dobson-cyclones-sumstat provides summary statistics.

```{r}
#| label: fig-cyclones-bar-plot
#| fig-cap: "Bar plot of cyclones per season"
cyclones |>
  ggplot() +
  geom_histogram(aes(x = number)) +
  expand_limits(x = 0) +
  xlab("Number of cyclones") +
  ylab("Count (number of seasons)")
```


```{r}
#| label: tbl-dobson-cyclones-sumstat
#| tbl-cap: "Summary statistics for `cyclones` data"

n = nrow(cyclones)
sumx = cyclones |> pull(number) |> sum()
xbar =  cyclones |> pull(number) |> mean()

cyclones |> table1::table1(x = ~ number)

```

### Model

We want to estimate $P(Y=y)$; that is, $P(Y=y)$ is our [estimand](estimation.qmd#def-estimand).

We could estimate $P(Y=y)$ for each value of $y$ in $0:\infty$ separately ("nonparametrically") using the fraction of our data with $Y_i=y$, but then we would be estimating an infinitely large set of parameters, and we would have low precision. We will probably do better with a parametric model.

:::{#exr-cyclone-choose-dist}
What parametric probability distribution family might we use to model this empirical distribution?

:::

::::{.solution}
Let's use the Poisson. The Poisson distribution is appropriate for this data , because the data are counts that could theoretically take any integer value (discrete) in the range $0:\infty$. Visually, the plot of our data closely resembles a Poisson or binomial distribution. Since cyclones do not have an "upper limit" on the number of events we could potentially observe in one season, the Poisson distribution is more appropriate than the binomial.
::::

:::{#exr-def-poisson}

Write down the Poisson distribution's probability mass function.

:::

::::{.solution}

$$P(Y=y) = \frac{\lambda^{y} e^{-\lambda}}{y!}$${#eq-iid-model}
::::


### Estimating the model parameters using maximum likelihood

Now, we can estimate the parameter $\lambda$ for this distribution using maximum likelihood estimation.

What is the likelihood?

:::{#exr-poisson-likelihood}
Write down the likelihood (probability mass function or probability density function) of a single observation $x$, according to your model.
:::

::::{.solution}

$$
\ba
\mathcal{L}(\lambda; x) 
&= p(X=x|\Lambda = \lambda)\\
&= \frac{\lambda^x e^{-\lambda}}{x!}\\
\ea
$$

::::


::: {#exr-poisson-parameters}

Write down the vector of parameters in your model.
:::
::::{.solution}

There is only one parameter, $\lambda$:

$$\theta = (\lambda)$$
::::


::: {#exr-poisson-mean-and-variance}
Write down the population mean and variance of a single observation from your chosen probability model, as a function of the parameters (extra credit - derive them).

:::
::::{.solution}
* Population mean: $\text{E}[X]=\lambda$ 
* Population variance: $\text{Var}(X)=\lambda$

::::


::: {#exr-sample-likelihood}

Write down the likelihood of the full dataset.

:::

::::{.solution}

$$
\ba
\Lik(\lambda; \vx) 
&= \P(\vX = \vx) \\
&= \P(X_1 = x_1, X_2 = x_2, ..., X_{13} = x_{13}) \\
&= \prod_{i=1}^{13} \P(X_i = x_i) \\
&= \prod_{i=1}^{13} \frac{\lambda^{x_i} e^{-\lambda}}{x_i!}
\ea
$$

::::

::: {#exr-graph-likelihood}

Graph the likelihood as a function of $\lambda$.

:::: {.solution}

```{r}
#| label: fig-cyclone-lik
#| fig-cap: "Likelihood of Dobson cyclone data"

lik = function(lambda, y = cyclones$number, n = length(y)) 
{
lambda^sum(y) * exp(-n*lambda) / prod(factorial(y))
}

library(ggplot2)
lik_plot = 
ggplot() +
geom_function(fun = lik, n = 1001) +
xlim(min(cyclones$number), max(cyclones$number)) +
ylab("likelihood") +
xlab('lambda')

print(lik_plot)
```

::::

:::

::: {#exr-sample-llik}

Write down the log-likelihood of the full dataset.

::::{.solution}

$$
\begin{aligned}
\ell(\lambda; \vec x) &= \log{\mathcal{L}(\lambda;\vec{x})}\\
&= \log{\prod_{i = 1}^n\frac{\lambda^{x_i}\text{e}^{-\lambda}}{x_i!}}\\
&= \sum_{i = 1}^n\log{\frac{\lambda^{x_i}\text{e}^{-\lambda}}{x_i!}}\\
&= \sum_{i = 1}^n{\log{\lambda^{x_i}} +\log{\text{e}^{-\lambda}} - \log{x_i!}}\\
&= \sum_{i = 1}^n{x_i\log{\lambda} -\lambda - \log{x_i!}}\\
&= \sum_{i = 1}^nx_i\log{\lambda} - \sum_{i = 1}^n\lambda - \sum_{i = 1}^n\log{x_i!}\\
&= \sum_{i = 1}^nx_i\log{\lambda} - n\lambda - \sum_{i = 1}^n\log{x_i!}\\
\end{aligned}
$$

::::

:::

::: {#exr-graph-loglikelihood}

Graph the log-likelihood as a function of $\lambda$.

:::: {.solution}

```{r}
#| label: fig-cyclone-llik
#| fig-cap: "log-likelihood of Dobson cyclone data"


loglik = function(lambda, y = cyclones$number, n = length(y))
{
sum(y) * log(lambda) - n*lambda - sum(log(factorial(y)))
}

ll_plot = ggplot() +
geom_function(fun = loglik, n = 1001) +
xlim(min(cyclones$number), max(cyclones$number)) +
ylab("log-likelihood") +
xlab('lambda')
ll_plot

```

::::

:::

#### The score function

::: {#exr-cyclone-score-fn}

Derive the score function for the dataset.

::::{.solution}

The score function is the first derivative of the log-likelihood:

$$
\begin{aligned}
\ell'( \lambda; \vec x ) &= 
\deriv{\lambda}{\sum_{i = 1}^nx_i\log{\lambda} - n\lambda - \sum_{i = 1}^n\log{x_i!}}\\
&= \deriv{\lambda}\sum_{i = 1}^nx_i\log{\lambda} - \deriv{\lambda}n\lambda - \deriv{\lambda}\sum_{i = 1}^n\log{x_i!}\\
&= \sum_{i = 1}^nx_i\deriv{\lambda}\log{\lambda} - n\deriv{\lambda}\lambda - \sum_{i = 1}^n\deriv{\lambda}\log{x_i!}\\
&= \sum_{i = 1}^nx_i\frac{1}{\lambda} - n - 0\\
&= \frac{1}{\lambda} \sum_{i = 1}^nx_i- n
\\&= \paren{\frac{1}{\lambda} n \xbar} - n
\\&= \paren{\frac{1}{\lambda} `r sumx`} - `r n`
\end{aligned}
$$

::::

:::

:::{#exr-graph-score-function}
Graph the score function.


::::{.solution}

```{r}
#| label: fig-cyclone-score
#| fig-cap: "score function of Dobson cyclone data"


score = function(lambda, y = cyclones$number, n = length(y))
{
  (sum(y) / lambda) - n
}

ggplot() +
geom_function(fun = score, n = 1001) +
xlim(min(cyclones$number), max(cyclones$number)) +

ylab("l'(lambda)") +
xlab('lambda') +
geom_hline(yintercept = 0, col = 'red')

```

::::

:::

#### The Hessian matrix {#sec-hessian}

:::{#exr-hessian}

Derive the Hessian matrix.

::::{.solution}

The Hessian function for an iid sample is the 2nd derivative(s) of the log-likelihood:

$$
\begin{aligned}
\ell''( \lambda; \vec x ) &= \deriv{\lambda}\left(\frac{1}{\lambda} \sum_{i = 1}^nx_i- n\right)\\
&= \deriv{\lambda}\frac{1}{\lambda} \sum_{i = 1}^nx_i- \deriv{\lambda}n\\
&= -\frac{1}{\lambda^2} \sum_{i = 1}^nx_i\\
&= -\frac{1}{\lambda^2} n \bar x
\\&= -\frac{1}{\lambda^2} \cd `r sumx`
\end{aligned}
$$
::::

:::

:::{#exr-graph-hession}

Graph the Hessian.

::::{.solution}

```{r}
#| label: fig-cyclone-hessian
#| fig-cap: "Hessian function of Dobson cyclone data"


hessian = function(lambda, y = cyclones$number, n = length(y))
{
-sum(y)/(lambda^2)
}

ggplot() +
geom_function(fun = hessian, n = 1001) +
xlim(min(cyclones$number), max(cyclones$number)) +

ylab("l''(lambda)") +
xlab('lambda') +
geom_hline(yintercept = 0, col = 'red')


```

::::

:::

:::{#exr-score-equation}
Write the score equation (estimating equation).

::::{.solution}

$$\ell'( \lambda; \vec x ) = 0$$

::::

:::

:::{#exr-solve-score-equation}
Solve the estimating equation for $\lambda$:

::::{.solution}
$$
\begin{aligned}
0 &= \frac{1}{\lambda}\sum_{i = 1}^nx_i - n\\
n &= \frac{1}{\lambda}\sum_{i = 1}^nx_i\\
n\lambda &= \sum_{i = 1}^nx_i\\
\lambda &= 
\frac{1}{n}\sum_{i = 1}^nx_i\\
&=\bar x
\end{aligned}
$$
::::

:::

Let's call this solution of the estimating equation $\tilde \lambda$ for now:

$$\tilde \lambda \eqdef \bar x$$



:::{#exr-check-hessian}

Confirm that the Hessian $\ell''(\lambda; \vec x)$ is negative when evaluated at $\tilde \lambda$.

::::{.solution}

$$
\begin{aligned}
\ell''( \tilde\lambda; \vec x ) &= 
-\frac{1}{\tilde\lambda^2} n \bar x\\
&= -\frac{1}{\bar x^2} n\bar x\\
&= -\frac{n}{\bar x}\\
&<0\\
\end{aligned}
$$

::::

:::

::: {#exr-find-mle}

Find the MLE of $\lambda$.

:::: {.solution}

Since $\ell''(\tilde \lambda; \vec x)<0$, $\tilde \lambda$ is at least a local maximizer of the likelihood function $\mathcal L(\lambda)$. Since there is only one solution to the estimating equation and the Hessian is negative definite everywhere, $\tilde \lambda$ must also be the global maximizer of $\mathcal L(\lambda; \vec x)$:

```{r}
#| label: calc-mle
mle = mean(cyclones$number)
```

$$\hat{\lambda}_{\text{ML}} = \bar x = `r mle`$$

::::

:::

:::{#exr-graph-mle}

Graph the log-likelihood with the MLE superimposed.

::::{.solution}

```{r}
#| label: fig-cyclone-llik-mle
#| fig-cap: "log-likelihood of Dobson cyclone data with MLE"
library(dplyr)

mle_data = tibble(x = mle, y = loglik(mle))
ll_plot + geom_point(data = mle_data, aes(x = x, y = y), col = 'red')
```

::::
:::

#### Information matrices

```{r}
#| label: fig-obs-inf-matrix
#| fig-cap: "Observed information function of Dobson cyclone data"
obs_inf = function(...) -hessian(...)
ggplot() +
geom_function(fun = obs_inf, n = 1001) +
xlim(min(cyclones$number), max(cyclones$number)) +
ylab("I(lambda)") +
xlab('lambda') +
geom_hline(yintercept = 0, col = 'red') 
```

---

:::::::{#exm-dobson-cyclones-Newton}

##### Finding the MLE using the Newton-Raphson algorithm

::: notes

We found that the MLE was $\hat{\lambda} = \bar{x}$, 
by solving the score equation $\ell'(\lambda)=0$ for $\lambda$.

What if we hadn't been able to solve the score equation?

Then we could start with some initial guess $\esttmpl$, 
such as $\esttmpl = 3$, and use the [Newton-Raphson algorithm](#sec-newton-raphson).

:::

```{r}
#| label: initial guess
# specify initial guess:
cur_lambda_est = 3
```

:::::::

---

::: notes

In @exr-cyclone-score-fn, we found that the score function was:

:::

$$
\score( \lambda; \vec x ) = \paren{\frac{`r sumx`}{\lambda} } - n
$$

::: notes
In @exr-hessian, we found that the Hessian was:
:::

$$
\hessian( \lambda; \vec x ) = -\frac{`r sumx`}{\lambda^2} 
$$

---

::: notes
So we can approximate the the score function using the first-order [Taylor polynomial](https://en.wikipedia.org/wiki/Taylor%27s_theorem):
:::

$$
\ba
\score(\lambda) 
&\approx \score^*(\lambda)
\\ &\eqdef \score(\esttmpl) + \hessian(\esttmpl)(\lambda - \esttmpl)
\\ &= \paren{\frac{`r sumx`}{\esttmpl}  - n} + \paren{-\frac{`r sumx`}{\sqf{\esttmpl}}} (\lambda - \esttmpl)
\ea
$$

---

::: notes

@fig-cyclone-newton-step1 compares the true score function and the approximate score function at $\esttmpl = `r cur_lambda_est`$.

:::

```{r}
#| label: fig-cyclone-newton-step1
#| fig-cap: "score function of Dobson cyclone data and approximate score function"

approx_score = function(lambda, lhat, ...)
{
  score(lambda = lhat, ...) +
  hessian(lambda = lhat, ...) * (lambda - lhat)
}

point_size = 5

plot1 = ggplot() +
geom_function(
  fun = score, 
  aes(col = "true score function"), 
  n = 1001) +
geom_function(
  fun = approx_score, 
  aes(col = "approximate score function"),
  n = 1001, 
  args = list(lhat = cur_lambda_est)) +
geom_point(
  size = point_size,
  aes(x = cur_lambda_est, y = score(lambda = cur_lambda_est),
      col = "current estimate")
) +
geom_point(
size = point_size,
aes(
x = xbar,
y = 0,
col = "true MLE"
)
) +
xlim(min(cyclones$number), max(cyclones$number)) +
ylab("l'(lambda)") +
xlab('lambda') +
geom_hline(yintercept = 0)

print(plot1)
```

---

::: notes

This is equivalent to estimating the log-likelihood with a second-order Taylor polynomial:

:::

$$
\llik^*(\lambda) = 
\llik(\esttmpl) + 
(\lambda - \esttmpl) \score(\esttmpl) +
\frac{1}{2}\hessian(\esttmpl)(\lambda-\esttmpl)^2
$$

```{r}
#| label: fig-cyclone-newton-step1-loglik
#| fig-cap: "log-likelihood of Dobson cyclone data and approximate log-likelihood function"

approx_loglik = function(lambda, lhat, ...)
{
loglik(lambda = lhat, ...) +
score(lambda = lhat, ...) * (lambda - lhat) +
  1/2 * hessian(lambda = lhat, ...) * (lambda - lhat)^2
}

plot_loglik = ggplot() +
geom_function(
  fun = loglik, 
  aes(col = "true log-likelihood"), 
  n = 1001) +
geom_function(
  fun = approx_loglik, 
  aes(col = "approximate log-likelihood"),
  n = 1001, 
  args = list(lhat = cur_lambda_est)) +
geom_point(
  size = point_size,
  aes(x = cur_lambda_est, y = loglik(lambda = cur_lambda_est),
      col = "current estimate")
) +
geom_point(
size = point_size,
aes(
x = xbar,
y = loglik(xbar),
col = "true MLE"
)
) +
xlim(min(cyclones$number) - 1, max(cyclones$number)) +
ylab("l'(lambda)") +
xlab('lambda')

print(plot_loglik)
```

---

::: notes
The approximate score function, $\score^*(\lambda)$, is a linear function of $\lambda$, so it is easy to solve the corresponding approximate score equation, $\score^*(\lambda) = 0$, for $\lambda$:
:::

$$
\ba
\lambda 
&= \esttmpl - \score(\esttmpl) \cd \inv{\hessian(\esttmpl)}
\\ &= `r cur_lambda_est - score(cur_lambda_est) * hessian(cur_lambda_est)^-1`
\ea
$$

```{r}
#| label: NR-one-step
new_lambda_est <- 
   cur_lambda_est - 
   score(cur_lambda_est) * hessian(cur_lambda_est)^-1
```

---

```{r}
#| label: fig-cyclone-newton-step1-with-new-est
#| fig-cap: "score function of Dobson cyclone data and approximate score function"

plot2 = plot1 +
  geom_point(size = point_size,
             aes(x = new_lambda_est,
                 y = 0,
                 col = "new estimate")) +
  geom_segment(
    arrow = grid::arrow(),
    linewidth = 2,
    alpha = .7,
    aes(
      x = cur_lambda_est,
      y = approx_score(lhat = cur_lambda_est,
                       lambda = cur_lambda_est),
      xend = new_lambda_est,
      yend = 0,
      col = "update"
    )
  )
print(plot2)

```

---

::: notes

So we update $\esttmpl \leftarrow `r new_lambda_est`$ and repeat our estimation process:

:::

```{r}
#| label: fig-cyclone-newton-step2
#| fig-cap: "score function of Dobson cyclone data and approximate score function"

plot2 +
geom_function(
  fun = approx_score, 
  aes(col = "new approximate score function"),
  n = 1001, 
  args = list(lhat = new_lambda_est)) +
geom_point(
  size = point_size,
  aes(x = new_lambda_est, y = score(lambda = new_lambda_est),
      col = "new estimate")
)

```

---

::: notes
We repeat this process until the likelihood converges:
:::

```{r}
#| label: tbl-mle-converge
#| tbl-cap: "Convergence of Newton-Raphson Algorithm for finding MLE of `cyclone` data"


library(tibble)
cur_lambda_est = 3 # restarting
diff_loglik = Inf
tolerance = 10 ^ -4
max_iter = 100
NR_info = tibble(
  iteration = 0,
  lambda = cur_lambda_est |> num(digits = 4),
  likelihood = lik(cur_lambda_est),
  `log(likelihood)` = loglik(cur_lambda_est)  |> num(digits = 4),
  score = score(cur_lambda_est),
  hessian = hessian(cur_lambda_est)
)

for (cur_iter in 1:max_iter)
{
  new_lambda_est <-
    cur_lambda_est - score(cur_lambda_est) * hessian(cur_lambda_est) ^ -1
  
  diff_loglik = loglik(new_lambda_est) - loglik(cur_lambda_est)
  
  new_NR_info = tibble(
    iteration = cur_iter,
    lambda = new_lambda_est,
    likelihood = lik(new_lambda_est),
    `log(likelihood)` = loglik(new_lambda_est),
    score = score(new_lambda_est),
    hessian = hessian(new_lambda_est),
    `diff(loglik)` = diff_loglik
  )
  
  NR_info = NR_info |> bind_rows(new_NR_info)
  
  cur_lambda_est = new_lambda_est
  
  if (abs(diff_loglik) < tolerance)
    break
  
}

NR_info

```

::: notes
Compare with @exr-find-mle
:::

---

```{r}
#| label: fig-NR-converge
#| fig-cap: "Newton-Raphson algorithm for finding MLE of [model @eq-iid-model] for `cyclone` data"

ll_plot +
  geom_segment(
    data = NR_info,
    arrow = grid::arrow(),
    # linewidth = 2,
    alpha = .7,
    aes(
      x = lambda,
      xend = lead(lambda),
      y = `log(likelihood)`,
      yend = lead(`log(likelihood)`),
      col = factor(iteration)
    )
  )

```