Introduction to likelihood.contr

The Problem

Suppose you observe failure times from a Weibull distribution, but some units are right-censored: you only know they survived past a cutoff. Exact observations contribute log f(t) to the log-likelihood while censored ones contribute log S(t). Both share the same shape and scale parameters, but the math is different for each row.

likelihood.contr handles this by letting you declare each observation type’s contribution separately, then composing them into a single model.

Declaring Contributions

contr_name() generates a contribution from any R distribution with d<name> and p<name> functions. Four censoring types are supported:

Type	Log-likelihood contribution
`"exact"`	`log f(x; theta)` via `d<name>`
`"right"`	`log S(x; theta)` via `p<name>(lower.tail=FALSE)`
`"left"`	`log F(x; theta)` via `p<name>`
`"interval"`	`log [F(b; theta) - F(a; theta)]` via `p<name>`

library(likelihood.contr)
library(likelihood.model)

exact <- contr_name("weibull", "exact", ob_col = "t")
right <- contr_name("weibull", "right", ob_col = "t")

Composing a Model

likelihood_contr() combines contributions and specifies how to dispatch rows to the right contribution. Here obs_type = "status" means the status column determines which contribution handles each row:

model <- likelihood_contr(
  obs_type = "status",
  exact = exact,
  right = right
)
model
#> Likelihood Contribution Model
#> -----------------------------
#> Observation types: exact, right 
#> Dispatch method: column 'status'
#> Assumptions:
#>  - iid

The contribution names (exact, right) must match the values in the dispatch column.

Simulating Censored Data

set.seed(42)
true_shape <- 2
true_scale <- 5
censor_time <- 4

raw_times <- rweibull(300, shape = true_shape, scale = true_scale)
df <- data.frame(
  t      = pmin(raw_times, censor_time),
  status = ifelse(raw_times <= censor_time, "exact", "right")
)
table(df$status)
#> 
#> exact right 
#>   142   158

Fitting the Model

The fit() generic returns a solver function. Call it with data and starting values:

result <- suppressWarnings(
  fit(model)(df, par = c(shape = 1.5, scale = 4))
)
summary(result)
#> Maximum Likelihood Estimate (Fisherian)
#> ----------------------------------------
#> 
#> Coefficients:
#>       Estimate Std. Error   2.5% 97.5%
#> shape   2.0672     0.1611 1.7514 2.383
#> scale   4.9571     0.2309 4.5046 5.410
#> 
#> Log-likelihood: -383.9 
#> AIC: 771.8 
#> Number of observations: 300

The result is a fisher_mle object with the usual methods:

coef(result)
#>    shape    scale 
#> 2.067153 4.957124
confint(result)
#>           2.5%    97.5%
#> shape 1.751370 2.382936
#> scale 4.504556 5.409692

Evaluating the Log-Likelihood Directly

Every likelihood.model generic returns a closure. This is useful for profiling or plotting:

ll_fn <- loglik(model)

# Evaluate at two different parameter vectors
ll_fn(df, par = c(shape = 2, scale = 5))
#> [1] -383.973
ll_fn(df, par = c(shape = 1, scale = 3))
#> [1] -485.7462

Function-Based Dispatch

When the observation type depends on multiple columns or a computed condition, pass a function instead of a column name:

model_fn <- likelihood_contr(
  obs_type = function(df) ifelse(df$delta == 1, "exact", "right"),
  exact = contr_name("exp", "exact", ob_col = "t"),
  right = contr_name("exp", "right", ob_col = "t")
)

df_fn <- data.frame(t = c(0.5, 1.0, 2.0), delta = c(1, 0, 1))
loglik(model_fn)(df_fn, par = c(rate = 1.5))
#> [1] -4.43907