Likelihood-Based Statistical Inference in the Fisherian Tradition

Facilitates building likelihood models in the Fisherian tradition following Richard Royall (1997, ISBN:978-0412044113) "Statistical Evidence: A Likelihood Paradigm". Defines generic methods for working with likelihoods (loglik(), score(), hess_loglik(), fim()) and provides functions for pure likelihood-based inference (support(), relative_likelihood(), likelihood_interval(), profile_loglik()).

likelihood.model

Likelihood-based statistical inference in the Fisherian tradition for R.

Why likelihood?

Classical statistics often jumps from data to p-values or posterior distributions. The likelihood approach takes a different path: the likelihood function is the complete summary of what the data say about the parameters. No priors, no long-run frequency guarantees -- just the evidence in the data.

This package makes it easy to:

• Specify likelihood models for standard distributions or custom observation types
• Estimate parameters via maximum likelihood, with analytical or numerical derivatives
• Quantify evidence using the likelihood function directly -- support functions, relative likelihood, and likelihood intervals
• Handle censored data (left, right, interval) as a natural part of the likelihood framework

Installation

# From CRAN
install.packages("likelihood.model")

# Development version from GitHub
devtools::install_github("queelius/likelihood.model")

Quick start

Fit a distribution

library(likelihood.model)

# Any R distribution works: "norm", "weibull", "exp", "gamma", ...
model <- likelihood_name("weibull", ob_col = "time", censor_col = "status")
mle <- fit(model)(df, par = c(shape = 1, scale = 1))

coef(mle)      # parameter estimates
se(mle)        # standard errors
confint(mle)   # confidence intervals

Handle censored data properly

Ignoring censoring biases your estimates. likelihood.model handles it correctly:

df <- data.frame(
  time   = c(1.2, 3.4, 5.0, 2.1, 5.0),
  status = c("exact", "exact", "right", "exact", "right")
)

model <- likelihood_name("weibull", ob_col = "time", censor_col = "status")
mle <- fit(model)(df, par = c(shape = 1, scale = 1))

Right-censored observations contribute through the survival function rather than the density -- the likelihood framework handles this seamlessly.

Fisherian inference

Instead of confidence intervals that rely on repeated-sampling arguments, likelihood intervals identify parameter values well-supported by the data:

# Likelihood interval: {theta : L(theta)/L(theta_hat) >= 1/8}
li <- likelihood_interval(mle, df, model, k = 8)

# How well does the data support theta1 over theta2?
evidence(model, df, theta1 = c(2, 3), theta2 = c(1, 5))

Three ways to build a model

likelihood_name() -- the quick path

Wraps any R distribution that has d<name> and p<name> functions. Supports exact, left-censored, right-censored, and interval-censored observations automatically.

model <- likelihood_name("norm", ob_col = "x", censor_col = "censor")

Specialized models -- the fast path

When you have analytical derivatives, estimation is 10-100x faster than numerical differentiation:

# Weibull with hand-derived score and Hessian
model <- weibull_uncensored("x")

# Exponential with closed-form MLE (bypasses optim entirely)
model <- exponential_lifetime("t", censor_col = "status")
mle <- fit(model)(df)  # no initial guess needed

likelihood_contr_model -- the flexible path

For complex models with different observation types, define type-specific log-likelihood contributions:

model <- likelihood_contr_model$new(
  obs_type = function(df) df$type,
  logliks = list(exact = loglik_exact, censored = loglik_censored)
)

The likelihood calculus

Every likelihood model provides a consistent set of operations:

loglik(model)       # log-likelihood function
score(model)        # score function (gradient)
hess_loglik(model)  # Hessian matrix
fim(model)          # Fisher information matrix

If you provide only loglik(), the package computes score() and hess_loglik() via numerical differentiation (Richardson extrapolation). Provide analytical versions for better speed and accuracy.

Inference toolkit

Maximum likelihood estimation

mle <- fit(model)(df, par = c(1, 1))
coef(mle)           # point estimates
vcov(mle)           # variance-covariance matrix
confint(mle)        # Wald confidence intervals
summary(mle)        # full summary

Bootstrap

boot <- sampler(model, df = df, par = c(1, 1))(n = 1000)
se(boot)            # bootstrap standard errors
confint(boot)       # bootstrap confidence intervals (BCa)
bias(boot)          # bootstrap bias estimate

Model comparison

lrt(null_model, alt_model, df,            # likelihood ratio test
    null_par = p0, alt_par = p1, dof = 1)
aic(mle)                                   # Akaike information criterion
bic(mle)                                   # Bayesian information criterion

Fisherian likelihood inference

support(mle, theta, df, model)             # S(theta) = logL(theta) - logL(theta_hat)
relative_likelihood(mle, theta, df, model) # R(theta) = L(theta) / L(theta_hat)
likelihood_interval(mle, df, model, k = 8) # {theta : R(theta) >= 1/k}
profile_loglik(mle, df, model, param = 1)  # profile out nuisance parameters
evidence(model, df, theta1, theta2)        # log-likelihood ratio

Ecosystem

likelihood.model interoperates with two companion packages:

• algebraic.mle -- MLE result objects with a unified interface (params, rmap, marginal)
• algebraic.dist -- probability distribution objects for simulation and analysis

Learn more

• Getting started -- worked examples from basic MLE to Fisherian inference
• Named distributions -- censoring, interval data, model comparison
• Likelihood contributions -- custom models for series systems and mixed observation types
• Exponential lifetime -- closed-form MLE, analytical derivatives, Monte Carlo validation
• algebraic.mle integration -- the three-package ecosystem
• Reference -- full API documentation