Constructing Likelihood Functions with Existing R Tools

R has no shortage of machinery for building likelihoods: optim() and optimize() in base R, maxLik, bbmle::mle2(), ProfileLikelihood, and cond for exact conditional GLM inference.¹ They look different, but for concurve they all reduce to the same thing: a grid of parameter values and a profile log-likelihood evaluated on it. This vignette shows several of those tools feeding one adapter.

Note: maxLik, bbmle, and cond are not concurve dependencies, so the chunks that use them run only if you have them installed. The optim() route is base R and always runs.

One adapter for all of them

Whatever produces the log-likelihood, this helper normalizes it and packs it into the five columns concurve uses (values, likelihood, loglikelihood, support, deviancestat) — the object curve_lik() returns and ggcurve() draws.

as_concurve_lik <- function(values, loglik) {
  loglik_rel <- loglik - max(loglik)
  support    <- exp(loglik_rel)
  df <- data.frame(
    values        = values,
    likelihood    = support,
    loglikelihood = loglik_rel,
    support       = support,
    deviancestat  = -loglik_rel
  )
  class(df) <- c("data.frame", "concurve")
  df
}

The only rule that matters for validity: when the model has more than one parameter, the log-likelihood you evaluate at each grid point must be profiled over the others (maximized over the nuisance parameters), not merely sliced at their estimates.² We return to that at the end, because it is the single most common way a “likelihood” plot goes wrong.

Some data to reuse:

set.seed(123)
x <- rnorm(50, mean = 9.3, sd = 3.2)
c(mean = mean(x), sd = sd(x))
#>     mean       sd 
#> 9.410091 2.962784

Route 1: base R `optim()` (+ observed-information SEs)

Minimize the negative log-likelihood to get the MLE and, from the Hessian, the observed-information standard errors — exactly the pattern in the pasted script.

nll_normal <- function(theta, x) {
  mu <- theta[1]; sigma <- theta[2]
  -sum(dnorm(x, mean = mu, sd = sigma, log = TRUE))
}

fit <- optim(c(mu = 8, sigma = 3), nll_normal, x = x,
             method = "BFGS", hessian = TRUE)
fit$par
#>       mu    sigma 
#> 9.410091 2.933006
sqrt(diag(solve(fit$hessian)))   # SE(mu), SE(sigma)
#>        mu     sigma 
#> 0.4147897 0.2933005

To turn this into a likelihood function for the mean, profile out sigma at each candidate mu — and here profiling is free, since the conditional MLE of sigma given mu is closed-form.

n <- length(x)
mu_grid <- seq(fit$par["mu"] - 3, fit$par["mu"] + 3, length.out = 4000)
loglik_mu <- vapply(mu_grid, function(m) {
  sigma2_hat <- mean((x - m)^2)          # profile: MLE of sigma^2 given mu
  -(n / 2) * (log(2 * pi) + log(sigma2_hat) + 1)
}, numeric(1))

lik_optim <- as_concurve_lik(mu_grid, loglik_mu)

ggcurve(lik_optim, type = "l1", nullvalue = round(fit$par["mu"], 2),
        xaxis = expression(mu),
        title = "Likelihood for the Mean (via optim profiling)")
#> Warning: Using `size` aesthetic for lines was deprecated in ggplot2 3.4.0.
#> ℹ Please use `linewidth` instead.
#> ℹ The deprecated feature was likely used in the concurve package.
#>   Please report the issue at <https://github.com/zadrafi/concurve/issues>.
#> This warning is displayed once per session.
#> Call `lifecycle::last_lifecycle_warnings()` to see where this warning was
#> generated.

Route 2: `maxLik`

maxLik maximizes a log-likelihood you supply (returning the value, not its negative) and gives a tidy summary() with SEs. To get the function rather than just the point, evaluate the profiled log-likelihood on a grid and hand it to the adapter.

library(maxLik)

loglik_normal <- function(theta, x) {
  sum(dnorm(x, mean = theta[1], sd = sqrt(theta[2]), log = TRUE))
}
start <- c(mu = 8, sig2 = 9)
theta_mle <- maxLik(loglik_normal, start = start, x = x)
summary(theta_mle)

# profile likelihood for mu, reusing the closed-form sigma^2 profile
loglik_mu2 <- vapply(mu_grid, function(m) {
  loglik_normal(c(m, mean((x - m)^2)), x = x)
}, numeric(1))

ggcurve(as_concurve_lik(mu_grid, loglik_mu2), type = "l1",
        nullvalue = round(coef(theta_mle)[1], 2),
        xaxis = expression(mu),
        title = "Likelihood for the Mean (via maxLik)")

maxLik is not installed, so this route is not evaluated. Install it with install.packages("maxLik") to run it.

Route 3: `bbmle::mle2()` with built-in profiling

bbmle is purpose-built for this: mle2() fits by maximum likelihood and profile() returns the profile likelihood for each parameter, which is exactly what we want. bbmle reports the profile on the signed square-root deviance ( $z$ ) scale; squaring recovers the deviance, and $\exp(-\tfrac12 z^2)$ recovers the relative likelihood.

library(bbmle)

nll <- function(mu, sigma) -sum(dnorm(x, mean = mu, sd = sigma, log = TRUE))
m <- mle2(nll, start = list(mu = 8, sigma = 3), data = list(x = x))

pr <- profile(m)
pr_mu <- as.data.frame(pr)                      # columns: param, z, par.vals...
pr_mu <- pr_mu[pr_mu$param == "mu", ]

mu_vals <- pr_mu[["par.vals.mu"]]
if (is.null(mu_vals)) mu_vals <- pr_mu[["mu"]]  # column name varies by version

lik_bbmle <- as_concurve_lik(mu_vals, -0.5 * pr_mu$z^2)  # loglik_rel = -z^2/2

ggcurve(lik_bbmle, type = "l1", nullvalue = round(coef(m)["mu"], 2),
        xaxis = expression(mu),
        title = "Profile Likelihood for the Mean (via bbmle)")

bbmle is not installed, so this route is not evaluated. Install it with install.packages("bbmle") to run it.

Route 4: exact conditional GLM inference with `cond`

For a logistic model, cond computes higher-order exact conditional likelihood inference for a coefficient of interest, eliminating the nuisance parameters by conditioning rather than by asymptotic approximation.² This is the same spirit as the exact conditional odds-ratio likelihood in the Constructing Valid Likelihood Functions vignette, but for a full regression.

library(cond)
data(babies)

mod <- glm(cbind(r1, r2) ~ day + lull - 1, family = binomial, data = babies)
cond_fit <- cond(mod, offset = lullyes)   # exact conditional inference for 'lull'
summary(cond_fit)
plot(cond_fit, which = 2)                 # conditional (modified profile) likelihood

cond is not installed, so this route is not evaluated. Install it with install.packages("cond") to run it. Its plot() draws the conditional and modified-profile likelihoods directly; to bring them into concurve, pass the object’s grid of parameter values and conditional log-likelihood to as_concurve_lik().

cond plots the conditional and modified-profile likelihoods itself; to render them in concurve, extract its grid of values and conditional log-likelihood and pass them to as_concurve_lik() as above.

The pitfall: slicing is not profiling

The pasted script contains a likelihood.normal.mu(mu, sig2 = 1, x) that fixes $\sigma^2 = 1$ . That is a slice through the likelihood at an assumed variance, not a profile. It peaks in the right place but carries the wrong curvature — so every interval read off it is wrong unless $\sigma^2$ really is 1.

loglik_slice <- vapply(mu_grid, function(m) {
  -0.5 * sum((x - m)^2)                  # fixes sigma^2 = 1
}, numeric(1))

lik_slice <- as_concurve_lik(mu_grid, loglik_slice)

ggcurve(lik_slice, type = "l1", nullvalue = round(fit$par["mu"], 2),
        xaxis = expression(mu),
        title = "WRONG: mean likelihood with sigma^2 fixed at 1",
        subtitle = "Too narrow -- the real sigma is about 3.2, not 1")

Compare its width to the profiled version from Route 1: the slice is far too tight because it pretends the variance is known and small. When you build a likelihood with any of the tools above, profile the nuisance parameters — optim() over them, bbmle::profile(), maxLik on a profiled objective, or a package like cond that conditions them away exactly.^2,3

Likelihood meets the P-value function

The payoff of building likelihoods properly is that they line up with the consonance (P-value) function for the same parameter — the two are coordinate views of one inferential object. Here we make that visible with plot_compare() using the native constructors curve_lik_glm() and curve_analytic() (this section evaluates when those functions are available in your installed version of concurve).

mod <- glm(am ~ mpg, family = binomial, data = mtcars)
est <- coef(mod)["mpg"]
se  <- sqrt(diag(vcov(mod)))["mpg"]

# profile likelihood: nuisance intercept maximized out at every grid point
lik_prof <- curve_lik_glm(mod, "mpg", steps = 200)

# Wald likelihood: the normal approximation implied by (est, se) --
# exactly the likelihood the consonance function from curve_analytic encodes
grid <- lik_prof[[1]]$values
lik_wald <- as_curve_lik(grid, dnorm(grid, mean = est, sd = se, log = TRUE))

plot_compare(lik_prof[[1]], lik_wald[[1]], type = "l1",
             title = "Profile vs Wald Likelihood",
             subtitle = "Same coefficient; the gap between them is the skew the Wald approximation ignores",
             xaxis = expression(beta[mpg]))

The two curves agree near the maximum and part company in the shoulders — which is precisely where interval limits live. The numeric version of the same statement: the 1/6.83 support interval from the profile likelihood reproduces the 95% consonance interval from the profile-based function, while the Wald consonance interval matches the Wald likelihood.

# 95% consonance interval, Wald construction
cons <- curve_analytic(estimate = est, se = se, dist = "z")
ci95 <- cons[[1]][which.min(abs(cons[[1]]$intrvl.level - 0.95)),
                  c("lower.limit", "upper.limit")]

rbind(
  `profile 1/6.83 support` = unlist(curve_support(lik_prof[[1]], 6.83)[c("lower.limit","upper.limit")]),
  `Wald    1/6.83 support` = unlist(curve_support(lik_wald[[1]], 6.83)[c("lower.limit","upper.limit")]),
  `Wald    95% consonance` = unlist(ci95)
)
#>                        lower.limit upper.limit
#> profile 1/6.83 support  0.12184097   0.5875499
#> Wald    1/6.83 support  0.08183696   0.5322194
#> Wald    95% consonance  0.08194283   0.5321136

The Wald support interval and the Wald consonance interval are identical by construction ( $e^{-z^2/2}$ at $z = 1.96$ is $1/6.83$ ); the profile versions differ from both exactly as much as the likelihood is skewed. Report the profile-based ones.

Adding these tools to your package

If you want any of the guarded routes to run when the vignettes build on CRAN, add the packages to Suggests in DESCRIPTION (maxLik, bbmle, cond) and keep the requireNamespace() guards — that is the pattern CRAN expects for optional, illustrative dependencies.

Cite R Packages

Please remember to cite the packages that you use.

citation("concurve")
#> To cite package 'concurve' in publications use:
#> 
#>   Rafi Z, Vigotsky A (2026). _concurve: Computes and Plots
#>   Compatibility (Confidence) Intervals, P-Values, S-Values, &
#>   Likelihood Intervals to Form Consonance, Surprisal, & Likelihood
#>   Functions_. R package version 3.0.0,
#>   <https://CRAN.R-project.org/package=concurve>.
#> 
#>   Rafi Z, Greenland S (2020). "Semantic and Cognitive Tools to Aid
#>   Statistical Science: Replace Confidence and Significance by
#>   Compatibility and Surprise." _BMC Medical Research Methodology_,
#>   *20*, 244. ISSN 1471-2288. doi:10.1186/s12874-020-01105-9
#>   <https://doi.org/10.1186/s12874-020-01105-9>.
#>   <https://doi.org/10.1186/s12874-020-01105-9>.
#> 
#> To see these entries in BibTeX format, use 'print(<citation>,
#> bibtex=TRUE)', 'toBibtex(.)', or set
#> 'options(citation.bibtex.max=999)'.

References

Choi L. ProfileLikelihood: Profile likelihood for a parameter in commonly used statistical models. 2011. https://CRAN.R-project.org/package=ProfileLikelihood.

Cox DR, Hinkley DV. Theoretical Statistics. Chapman and Hall/CRC; 1974. doi:10.1201/b14832

Royall R. Statistical Evidence: A Likelihood Paradigm. CRC Press; 1997.