---
title: "The MLE Ecosystem"
author: "Alexander Towell"
date: "`r Sys.Date()`"
output:
  rmarkdown::html_vignette:
    toc: true
    toc_depth: 3
vignette: >
  %\VignetteIndexEntry{The MLE Ecosystem}
  %\VignetteEngine{knitr::rmarkdown}
  %\VignetteEncoding{UTF-8}
---

```{r setup, include=FALSE}
knitr::opts_chunk$set(
  collapse = TRUE,
  comment = "#>",
  fig.width = 6,
  fig.height = 4,
  eval = FALSE
)
```

## Introduction

`compositional.mle` is part of an ecosystem of R packages that share a
single design principle from *Structure and Interpretation of Computer
Programs* (SICP): **closure**---the result of combining things is itself
the same kind of thing, ready to be combined again.

This principle appears at three levels:

| Level | Package | Closure property |
|-------|---------|------------------|
| **Solvers** | `compositional.mle` | Composing solvers yields a solver |
| **Estimates** | `algebraic.mle` | Combining MLEs yields an MLE |
| **Tests** | `hypothesize` | Combining hypothesis tests yields a test |

This vignette walks through that story. All packages referenced below
are available on CRAN.

## Composing Solvers

The core idea of `compositional.mle` is that solvers are first-class
functions that compose via operators:

- `%>>%` (sequential chaining) --- pipe the result of one solver into
  the next
- `%|%` (parallel racing) --- run solvers simultaneously, keep the best
- `with_restarts()` --- retry a solver from multiple starting points

Because each operator returns a solver, you can compose freely:

```{r}
library(compositional.mle)

# Generate data
set.seed(42)
y <- rnorm(200, mean = 5, sd = 2)

# Define the log-likelihood (numDeriv handles derivatives by default)
ll <- function(theta) {
  mu <- theta[1]; sigma <- theta[2]
  -length(y) * log(sigma) - 0.5 * sum((y - mu)^2) / sigma^2
}

problem <- mle_problem(loglike = ll)

# Coarse-to-fine: grid search for a starting region, then L-BFGS-B
strategy <- grid_search(lower = c(0, 0.5), upper = c(10, 5), n = 5) %>>%
  lbfgsb(lower = c(-Inf, 1e-4), upper = c(Inf, Inf))

result <- strategy(problem, theta0 = c(0, 1))
result
```

```
#> MLE Solver Result
#> -----------------
#> Method:     grid -> lbfgsb
#> Converged:  TRUE
#> Iterations: 10
#> Log-likelihood: -424.84
#>
#> Estimate:
#> [1] 5.0609 1.9498
```

You can race different approaches and let the best log-likelihood win:

```{r}
race <- bfgs() %|% nelder_mead() %|% gradient_ascent()
result_race <- race(problem, theta0 = c(0, 1))
```

Or combine restarts with chaining:

```{r}
robust <- with_restarts(
  gradient_ascent() %>>% bfgs(),
  n = 10,
  sampler = uniform_sampler(lower = c(-10, 0.1), upper = c(10, 10))
)
result_robust <- robust(problem, theta0 = c(0, 1))
```

See `vignette("strategy-design")` for a thorough treatment of solver
composition.

## The `mle` Result Object

Every solver returns an `mle` object (defined by `algebraic.mle`). This
is the shared currency of the ecosystem---any package that produces or
consumes `mle` objects works with every other package.

The `mle` class provides a uniform interface:

```{r}
library(algebraic.mle)

# Extract estimates and uncertainty
params(result)       # parameter estimates (theta-hat)
se(result)           # standard errors
vcov(result)         # variance-covariance matrix
as.numeric(logLik(result))   # log-likelihood at the MLE
confint(result)      # 95% confidence intervals
```

```
#> params(result)
#> [1] 5.0609 1.9498
#>
#> se(result)
#> [1] 0.1379 0.0975
#>
#> vcov(result)
#>              [,1]         [,2]
#> [1,]  0.019009 -0.000007
#> [2,] -0.000007  0.009506
#>
#> as.numeric(logLik(result))
#> [1] -424.84
#>
#> confint(result)
#>          2.5 %   97.5 %
#> [1,] 4.7906  5.3312
#> [2,] 1.7587  2.1410
```

The key point: **results from any solver have the same interface**.
Whether you used Newton-Raphson, a grid search, or a composed pipeline,
`params()`, `se()`, and `confint()` work the same way.

## Combining Independent MLEs

Suppose three independent labs each estimate the mean of a
Normal$(\mu, \sigma)$ population from their own data. Each produces an
`mle` object. `algebraic.mle` provides `mle_weighted()` to combine them
via **inverse-variance weighting** using Fisher information:

```{r}
# Three independent datasets from the same population
set.seed(123)
y1 <- rnorm(50,  mean = 10, sd = 3)
y2 <- rnorm(100, mean = 10, sd = 3)
y3 <- rnorm(75,  mean = 10, sd = 3)

# Fit each independently
make_problem <- function(data) {
  mle_problem(loglike = function(theta) {
    mu <- theta[1]; sigma <- theta[2]
    -length(data) * log(sigma) - 0.5 * sum((data - mu)^2) / sigma^2
  })
}

solver <- lbfgsb(lower = c(-Inf, 1e-4), upper = c(Inf, Inf))

fit1 <- solver(make_problem(y1), theta0 = c(0, 1))
fit2 <- solver(make_problem(y2), theta0 = c(0, 1))
fit3 <- solver(make_problem(y3), theta0 = c(0, 1))

# Combine: inverse-variance weighting via Fisher information
combined <- mle_weighted(list(fit1, fit2, fit3))
```

```{r}
# Individual SEs vs combined
data.frame(
  Source   = c("Lab 1 (n=50)", "Lab 2 (n=100)", "Lab 3 (n=75)", "Combined"),
  Estimate = c(params(fit1)[1], params(fit2)[1],
               params(fit3)[1], params(combined)[1]),
  SE       = c(se(fit1)[1], se(fit2)[1], se(fit3)[1], se(combined)[1])
)
```

```
#>            Source Estimate     SE
#> 1  Lab 1 (n=50)   10.207 0.4243
#> 2 Lab 2 (n=100)    9.869 0.2998
#> 3  Lab 3 (n=75)    9.753 0.3464
#> 4      Combined    9.934 0.1995
```

The combined estimate has a smaller standard error than any individual
estimate---this is the power of information pooling. And because
`mle_weighted()` returns an `mle` object, you can pass it to any
function that accepts MLEs: `confint()`, `se()`, further weighting, or
hypothesis testing. **Combining MLEs yields an MLE** (closure).

## Hypothesis Testing with `hypothesize`

The [`hypothesize`](https://CRAN.R-project.org/package=hypothesize)
package provides composable hypothesis tests. Like solvers and
estimates, **combining tests yields a test**.

### Creating tests from MLE results

`wald_test()` tests whether a single parameter equals a hypothesised
value. `lrt()` compares nested models via the likelihood ratio:

```{r}
library(hypothesize)

# Wald test: is mu = 10?
w <- wald_test(estimate = params(combined)[1],
               se = se(combined)[1],
               null_value = 10)
w
```

```
#> Hypothesis test ( wald_test )
#> -----------------------------
#> Test statistic: 0.1094
#> P-value:       0.7409
#> Degrees of freedom: 1
#> Significant at 5% level: FALSE
```

```{r}
# Likelihood ratio test: does sigma differ from 3?
# Compare full model (mu, sigma free) vs restricted (mu free, sigma = 3)
ll_full <- as.numeric(logLik(fit2))  # full model log-likelihood
ll_null <- sum(dnorm(y2, mean = params(fit2)[1], sd = 3, log = TRUE))

lr <- lrt(null_loglik = ll_null, alt_loglik = ll_full, dof = 1)
lr
```

```
#> Hypothesis test ( likelihood_ratio_test )
#> ------------------------------------------
#> Test statistic: 0.2851
#> P-value:       0.5934
#> Degrees of freedom: 1
#> Significant at 5% level: FALSE
```

### Combining independent tests

When you have independent tests (e.g., each lab tests $H_0: \mu = 10$),
`fisher_combine()` merges them using Fisher's method:

```{r}
# Each lab tests H0: mu = 10
w1 <- wald_test(params(fit1)[1], se(fit1)[1], null_value = 10)
w2 <- wald_test(params(fit2)[1], se(fit2)[1], null_value = 10)
w3 <- wald_test(params(fit3)[1], se(fit3)[1], null_value = 10)

# Combine independent tests
combined_test <- fisher_combine(w1, w2, w3)
combined_test
```

```
#> Hypothesis test ( fisher_combined_test )
#> -----------------------------------------
#> Test statistic: 3.4287
#> P-value:       0.7534
#> Degrees of freedom: 6
#> Significant at 5% level: FALSE
```

```{r}
# Multiple testing correction
adjusted <- adjust_pval(list(w1, w2, w3), method = "BH")
```

The combined test is itself a `hypothesis_test` object---you can
extract `pval()`, check `is_significant_at()`, or combine it further.
This is the closure property at work.

### The testing interface

All `hypothesis_test` objects share a uniform interface, mirroring how
all `mle` objects share `params()` / `se()`:

```{r}
pval(w)                    # extract p-value
test_stat(w)               # extract test statistic
dof(w)                     # degrees of freedom
is_significant_at(w, 0.05) # check significance at alpha = 0.05
confint(w)                 # confidence interval (for Wald tests)
```

## Automatic Differentiation

Derivatives are pluggable in `mle_problem()`. The solver never knows
(or cares) where they come from:

```{r}
# Strategy 1: numDeriv (default) --- zero effort
p1 <- mle_problem(loglike = ll)

# Strategy 2: hand-coded analytic derivatives
p2 <- mle_problem(
  loglike = ll,
  score = function(theta) {
    mu <- theta[1]; sigma <- theta[2]; n <- length(y)
    c(sum(y - mu) / sigma^2,
      -n / sigma + sum((y - mu)^2) / sigma^3)
  },
  fisher = function(theta) {
    n <- length(y); sigma <- theta[2]
    matrix(c(n / sigma^2, 0, 0, 2 * n / sigma^2), 2, 2)
  }
)

# Strategy 3: automatic differentiation
# Any AD package that provides score(f, theta) and hessian(f, theta)
# can be plugged in here, e.g.:
# p3 <- mle_problem(
#   loglike = ll,
#   score   = function(theta) ad_pkg::score(ll, theta),
#   fisher  = function(theta) ad_pkg::hessian(ll, theta)
# )
```

Each strategy has trade-offs:

| Approach | Accuracy | Effort | Speed |
|----------|----------|--------|-------|
| Numerical (`numDeriv`) | $O(h^2)$ truncation error | None (default) | Slow (extra evaluations) |
| Hand-coded analytic | Exact | High (error-prone) | Fast |
| Automatic differentiation (AD package) | Exact (machine precision) | Low | Moderate |

The general recommendation: **start with the default** (no explicit
derivatives) during model development, then **switch to an AD package**
or hand-coded derivatives when you need accuracy or performance. For
derivative-free problems, use `nelder_mead()`.

| Situation | Recommended | Reason |
|-----------|-------------|--------|
| Quick prototyping | `numDeriv` (default) | Zero effort, just write the log-likelihood |
| Production / accuracy-critical | AD package or hand-coded | Exact derivatives |
| Complex models (mixtures, hierarchical) | AD package | Hand-coding is error-prone and tedious |
| Non-differentiable objective | `nelder_mead()` | Derivative-free simplex method |

## The `likelihood.model` Bridge

The [`likelihood.model`](https://github.com/queelius/likelihood.model) package
provides a Fisherian likelihood framework with generics (`loglik`, `score`,
`hess_loglik`, `fim`) and contribution-based models for heterogeneous observation
types.

`compositional.mle` provides a bridge: pass a `likelihood_model` object and
your data directly to `mle_problem()`, and it will auto-extract the loglik,
score, and Fisher information functions:

```{r likelihood-model-bridge, eval=FALSE}
library(likelihood.model)

# model is a likelihood_model object from likelihood.model
problem <- mle_problem(model, data = my_data)

# Now solve with any composed strategy
solver <- lbfgsb(lower = c(-Inf, 1e-4), upper = c(Inf, Inf))
result <- solver(problem, theta0 = c(0, 1))
```

This means any model in the `likelihood.model` ecosystem can immediately be
optimized with the full suite of `compositional.mle` strategies---no glue
code required.

## The Ecosystem at a Glance

The packages form a pipeline: **define** a model, **solve** for the MLE,
**combine** estimates, and **test** hypotheses---with closure at every
level.

| Package | CRAN | Purpose |
|---------|------|---------|
| [`algebraic.mle`](https://CRAN.R-project.org/package=algebraic.mle) | Yes | MLE algebra: `params()`, `se()`, `confint()`, `mle_weighted()` |
| [`algebraic.dist`](https://CRAN.R-project.org/package=algebraic.dist) | Yes | Distribution algebra: convolution, mixture, truncation |
| [`compositional.mle`](https://github.com/queelius/compositional.mle) | --- | Solver composition: `%>>%`, `%|%`, `with_restarts()` |
| [`hypothesize`](https://CRAN.R-project.org/package=hypothesize) | Yes | Test composition: `wald_test()`, `lrt()`, `fisher_combine()` |
| [`flexhaz`](https://github.com/queelius/flexhaz) | --- | Distributions via arbitrary hazard functions |
| [`likelihood.model`](https://github.com/queelius/likelihood.model) | --- | Fisherian likelihood framework; bridge via `mle_problem(model, data)` |

The typical workflow:

1. **Define** your model as a log-likelihood and pass it to
   `mle_problem()`
2. **Solve** with a composed solver strategy from `compositional.mle`
   --- get an `mle` object
3. **Combine** independent estimates with `mle_weighted()` from
   `algebraic.mle` --- get a (better) `mle` object
4. **Test** hypotheses with `hypothesize` --- get a `hypothesis_test`
   object
5. **Combine** tests with `fisher_combine()` --- get a (stronger)
   `hypothesis_test` object

At every step, the result is the same type as the input. That is the
SICP closure property, and it is what makes the ecosystem composable.
