04_maximum_likelihood.qmd

# Maximum likelihood estimation

------------------------------------------------------------------------

## Maximum likelihood: the probabilistic approach

[**Core Idea:**]{style="color: #FFFF00;"} Find parameter values that make the observed data most probable

------------------------------------------------------------------------

**Mathematical Formulation:** $$\hat{\theta} = \arg\max_{\theta} L(\theta) = \arg\max_{\theta} \prod_{i=1}^{n} f(y_i | \theta)$$

Where:

- $L(\theta)$ = likelihood function
- $f(y_i | \theta)$ = probability density of observation $i$

------------------------------------------------------------------------

**In Practice:** Maximize log-likelihood $$\hat{\theta} = \arg\max_{\theta} \ell(\theta) = \arg\max_{\theta} \sum_{i=1}^{n} \log f(y_i | \theta)$$

------------------------------------------------------------------------

## Why maximum likelihood?

### Theoretical advantages:

-   Principled [statistical framework]{style="color: #FFFF00;"} 
-   Provides [uncertainty quantification]{style="color: #FFFF00;"}
-   Enables [model comparison (AIC, BIC)]{style="color: #FFFF00;"}
-   [Asymptotically]{style="color: #FFFF00;"} optimal properties

------------------------------------------------------------------------

### Practical benefits

-   Confidence intervals
-   Hypothesis testing
-   Model selection
-   Incorporates different error structures

------------------------------------------------------------------------

## Choosing a probability distribution

**For [Count Data]{style="color: #FFFF00;"} (e.g., Cases):**

-   **Poisson**: $Y_i \sim \text{Poisson}(\lambda_i)$
-   **Negative Binomial**: $Y_i \sim \text{NB}(\mu_i, \phi)$

------------------------------------------------------------------------

## Choosing a probability distribution

**For [Continuous Data:]{style="color: #FFFF00;"}**

-   **Normal**: $Y_i \sim N(\mu_i, \sigma^2)$
-   **Log-normal**: $\log Y_i \sim N(\log \mu_i, \sigma^2)$

::: {.callout-note}
[**For Our SIR Example:**]{style="color:yellow"} We'll use [Poisson]{style="yellow"} since we're modeling case counts
:::

------------------------------------------------------------------------

## MLE implementation: Poisson likelihood

```{r mle-implementation, echo=TRUE}
#| code-line-numbers: "1-2|3-6|7-8|10-13|16-17|18|18-24|25-26|27" 
# Define negative log-likelihood function
nll_function <- function(beta, gamma, data) {
  # Simulate model
  out <- ode(y = init_conds, times = data$time,
             func = sir_model, parms = c(beta = beta, gamma = gamma))

  # Model predictions (scaled to cases)
  predicted <- out[,"I"] * 1000

  # Poisson negative log-likelihood
  nll <- -sum(dpois(data$observed, lambda = predicted, log = TRUE))

  return(nll)
}
```

------------------------------------------------------------------------

```{r mle-fitting, echo=TRUE}
#| code-line-numbers: "1-2|3|4|6|7-8|1-8|10-11|12" 
# Use optimization to find MLE
library(bbmle)
fit_mle <- mle2(nll_function,
                start = list(beta = 0.2, gamma = 0.1),
                data = list(data = data),
                method = "L-BFGS-B",
                lower = c(0.01, 0.01),
                upper = c(1.0, 0.5))

# Extract results
mle_params <- coef(fit_mle)
mle_se <- sqrt(diag(vcov(fit_mle)))
```

------------------------------------------------------------------------

### MLE results

```{r mle-results, echo=FALSE}
cat("Beta:", round(mle_params[1], 3), "±", round(mle_se[1], 3), "\n")
cat("Gamma:", round(mle_params[2], 3), "±", round(mle_se[2], 3), "\n")
cat("R0:", round(mle_params[1]/mle_params[2], 2), "\n")
```

------------------------------------------------------------------------

## Uncertainty quantification with MLE

```{r mle-uncertainty, echo=TRUE, fig.width=10, fig.height=6}
# Profile likelihood for uncertainty
prof <- profile(fit_mle)
plot(prof, absVal = TRUE, main = "Profile Likelihood")
```

------------------------------------------------------------------------

```{r mle-ci, echo=TRUE}
# Confidence intervals
confint(fit_mle, level = 0.95)
```

------------------------------------------------------------------------

### MLE estimates vs true values
```{r mle-comparison, echo=FALSE}
# Compare with true values
cat("True Beta:", true_params[1], "\n")
cat("MLE Beta:", round(mle_params[1], 3), "\n")
cat("True Gamma:", true_params[2], "\n")
cat("MLE Gamma:", round(mle_params[2], 3), "\n")
```

------------------------------------------------------------------------

## Model comparison with MLE

```{r model-comparison, echo=TRUE}
#| code-line-numbers: "1|2|3|5-6|7-8|10-11|13-16" 
# Fit different models and compare
# Model 1: SIR with Poisson
# Model 2: SIR with negative binomial

# Negative binomial likelihood
nll_nb <- function(beta, gamma, phi, data) {
  out <- ode(y = init_conds, times = data$time,
             func = sir_model, parms = c(beta = beta, gamma = gamma))

  # Extract and scale predictions to cases
  predicted <- out[,"I"] * 1000

  # Negative Binomial negative log-likelihood
  nll <- -sum(dnbinom(data$observed, mu = predicted, size = phi, log = TRUE))

  return(nll)
}
```

------------------------------------------------------------------------

```{r model-fit, echo=TRUE}
#| code-line-numbers: "1-2|3|5|6-7|1-7" 
# Fit NB model
fit_nb <- mle2(nll_nb,
               start = list(beta = 0.2, gamma = 0.1, phi = 10),
               data = list(data = data),
               method = "L-BFGS-B",
               lower = c(0.01, 0.01, 0.1),
               upper = c(1.0, 0.5, 100))
```

------------------------------------------------------------------------

### Model comparison
```{r model-comparison-results, echo=FALSE}
# Compare models
cat("Poisson AIC:", AIC(fit_mle), "\n")
cat("Negative Binomial AIC:", AIC(fit_nb), "\n")
cat("Delta AIC:", AIC(fit_nb) - AIC(fit_mle), "\n")
```

------------------------------------------------------------------------

## Strengths of maximum likelihood

**Statistical Rigor:**

-   Principled probabilistic framework
-   Asymptotic optimality properties
-   Natural uncertainty quantification
-   Enables formal hypothesis testing

------------------------------------------------------------------------

## Practical benefits

-   Confidence intervals and standard errors
-   Model comparison via AIC/BIC
-   Handles different error structures
-   Extensible to complex models

------------------------------------------------------------------------

## Limitations of maximum likelihood

**Computational Challenges:**

-   More complex than least squares
-   Requires optimization algorithms
-   Can get stuck in local minima
-   Sensitive to starting values

------------------------------------------------------------------------

## Statistical assumptions

-   Requires specification of error distribution
-   Assumes model structure is correct
-   Asymptotic properties may not hold
-   Can be sensitive to outliers

------------------------------------------------------------------------

## Identifiability issues

-   Still suffers from parameter identifiability
-   Profile likelihood can be computationally expensive
-   May not converge for complex models

------------------------------------------------------------------------

## R implementation practicals {background-color="#447099" transition="fade-in"}

-   Let's turn to the [tutorials](https://github.com/jamesmbaazam/mppr_intro_to_fitting_practicals)