Profile likelihood: routes and fallbacks • gllvmTMB

This reference answers one practical question: when you request a profile interval, does gllvmTMB actually compute one for your target?

Direct parameters and simple linear contrasts can use likelihood profiles. Nonlinear derived quantities do not currently have profile intervals because the constrained fits cannot yet be checked reliably in every case. Empirical coverage has also not been established across all supported models, so this page does not call any interval method universally best or publication-ready.

The key rule is simple: choose a route that is actually available for the target, then inspect the returned bounds and any method or status columns. A documented numerical-Wald fallback must retain its label, while a withdrawn nonlinear route must stop rather than being presented as a profile result.

library(gllvmTMB)
set.seed(42)

Current route map

The interval API is target-specific. The table below shows what a reader can request now; it is not a claim about empirical coverage.

Target	Start here	Current behavior
Direct fixed effect, dispersion, or variance component	`profile_targets(fit)` followed by `confint(fit, parm = ..., method = "profile")`	Direct likelihood profile when the target is marked ready.
Repeatability	`extract_repeatability()`	Wald is the default; requesting profile stops with a typed withdrawal because no full-covariance profile route is available.
Trait correlation	`extract_correlations()`	Point-only is the default. Fisher-z and bootstrap bounds are opt-in; a profile interval is not currently available.
Communality	`extract_communality()`	Point estimates are the default; Wald and bootstrap bounds are opt-in; a profile interval is not currently available.
Phylogenetic signal	`extract_phylo_signal()`	A simple two-component decomposition can use a profile contrast. Richer decompositions return labelled numerical-Wald bounds; a general profile constraint is not implemented.
Predictor-informed latent effect $B_{lv}$	`extract_lv_effects()`	Point estimates and experimental delta-method Wald summaries are available; profile and bootstrap intervals are not currently available.
Full covariance matrix	`confint(fit, parm = "Sigma_unit", ...)`	Route depends on the fitted structure. Read the returned method and any fallback message rather than assuming every matrix element was directly profiled.

Bootstrap is not an automatic remedy for an unusable profile. It relies on simulation and successful refits, and every active random-effect tier must be redrawn rather than silently conditioned on its fitted mode. Use enough replicates to estimate the requested tail probabilities, report failed refits, and treat small runs only as software smoke tests.

Fit one model from long or wide data

The same Gaussian model can be fitted through the canonical long interface or the compact wide data-frame interface. The rest of the article uses the long fit. Here unit = "site" names the between-site tier and unit_obs = "site_species" names the within-site sampling cell used by the second covariance term. No third grouping role is present, so this model does not pass cluster or cluster2.

s <- simulate_site_trait(
  n_sites               = 50,
  n_species             = 6,
  n_traits              = 3,
  mean_species_per_site = 4,
  Lambda_B              = matrix(c(0.9, 0.4, -0.3), 3, 1),
  psi_B                 = c(0.40, 0.30, 0.50),
  psi_W                 = c(0.30, 0.40, 0.30),
  beta                  = matrix(0, 3, 2),
  seed                  = 42
)
df_long <- s$data

library(tidyr)
df_wide <- tidyr::pivot_wider(
  df_long,
  id_cols = c(site, site_species),
  names_from = trait,
  values_from = value
)

# Long format (canonical)
fit_long <- gllvmTMB(
  value ~ 0 + trait +
    latent(0 + trait | site, d = 1, unique = TRUE) +
    indep(0 + trait | site_species),
  data = df_long,
  trait = "trait",
  unit = "site",
  unit_obs = "site_species"
)

# Wide data-frame format
fit_wide <- gllvmTMB(
  traits(trait_1, trait_2, trait_3) ~ 1 +
    latent(1 | site, d = 1, unique = TRUE) +
    indep(1 | site_species),
  data = df_wide,
  unit = "site",
  unit_obs = "site_species"
)

data.frame(
  fit = c("long", "wide"),
  logLik = c(as.numeric(logLik(fit_long)), as.numeric(logLik(fit_wide))),
  converged = c(fit_long$fit_health$converged,
                fit_wide$fit_health$converged)
)
#>    fit    logLik converged
#> 1 long -754.2025      TRUE
#> 2 wide -754.2025      TRUE

fit <- fit_long

Do not interpret an interval until the point fit is stable. If fit$fit_health$converged is false, repair the point surface first. A failed Hessian blocks or weakens Wald inference, but it does not turn an unstable fit into a valid profile problem. See Fit diagnostics and Convergence and start values.

Direct targets: inventory first, then profile

profile_targets() lists the fitted quantities that can be sent directly to confint(..., method = "profile"). This public inventory avoids making readers guess internal parameter names.

targets <- profile_targets(fit, ready_only = TRUE)
direct_summary <- subset(
  targets,
  target_class %in% c("fixed_effect", "variance")
)
head(direct_summary[c("parm", "target_class", "estimate", "profile_note")], 12)
#>        parm target_class     estimate profile_note
#> 1  b_fix[1] fixed_effect  1.136393713        ready
#> 2  b_fix[2] fixed_effect -0.679912771        ready
#> 3  b_fix[3] fixed_effect  0.341051849        ready
#> 7   sd_B[1]     variance  0.003004326        ready
#> 8   sd_B[2]     variance  0.604893187        ready
#> 9   sd_B[3]     variance  0.676798959        ready
#> 10  sd_W[1]     variance  0.793567166        ready
#> 11  sd_W[2]     variance  0.909648633        ready
#> 12  sd_W[3]     variance  0.878504383        ready

The first between-unit diagonal standard deviation is near its natural boundary in this fixture. Request its profile by the public parm label:

ci_sd <- confint(
  fit,
  parm = "sd_B[1]",
  method = "profile",
  level = 0.95
)
ci_sd
#>         2.5 %    97.5 %
#> sd_B[1]     0 0.1120251

The lower endpoint is 0, not NA: the profile reached the natural variance boundary. Current profile helpers use these conventions:

a finite number is an estimated threshold crossing;
a natural boundary such as 0, 1, or -1 is returned when that is the corresponding endpoint;
an infinite endpoint means the threshold was not crossed before the parameter limit;
NA means the profile failed or did not provide enough usable evaluations.

The direct wrapper uses the standard one-degree-of-freedom likelihood-ratio cutoff. Variance-component boundaries are non-regular, and this implementation has not been broadly calibrated for every boundary regime.

Why this page does not offer a t-based profile cutoff

The released direct-profile routes use the standard one-degree-of-freedom likelihood-ratio reference. A t-based cutoff would require defensible, target-specific degrees of freedom and repeated-sampling evidence. The number of units minus the latent rank is not a generally justified substitute for a Satterthwaite or Kenward–Roger derivation. gllvmTMB therefore does not teach or expose an automatic t-based profile correction in this release.

This is an evidence boundary, not a claim that t-based inference is unimportant. A future route needs to define the target-specific degrees of freedom, show how they are estimated, and demonstrate coverage across the intended sample-size regimes before it becomes reader-facing.

Derived targets: read the method that came back

Repeatability uses Wald; profile is withdrawn

For the full covariance definition of repeatability, the default is a delta-method Wald interval. method = "profile" does not perform a constrained likelihood profile and deliberately stops with gllvmTMB_repeatability_profile_withdrawn; the former profile token targeted only a diagonal-companion ratio and omitted shared latent variance. Use the default Wald route or the parametric-bootstrap route, and report its returned method and limitations.

rep_wald <- extract_repeatability(fit, level = 0.95)
rep_wald
#>     trait         R     lower     upper method
#> 1 trait_1 0.5906338 0.4488507 0.7187941   wald
#> 2 trait_2 0.3515303 0.2160166 0.5160929   wald
#> 3 trait_3 0.3734982 0.2347709 0.5367058   wald

The bounds are formed on a transformed scale and back-transformed to $[0,1]$ ; they are not symmetric on the natural repeatability scale. They are Wald bounds, not profile intervals, and this article does not establish their repeated-sampling calibration.

Correlation intervals are opt-in

cors <- extract_correlations(
  fit,
  tier = "unit",
  level = 0.95,
  method = "fisher-z"
)
cors
#>   tier trait_i trait_j correlation      lower     upper   method
#> 1    B trait_1 trait_2  0.42928528  0.1714207 0.6321027 fisher-z
#> 2    B trait_1 trait_3 -0.06666904 -0.3387309 0.2156811 fisher-z
#> 3    B trait_2 trait_3 -0.02862032 -0.3045419 0.2517328 fisher-z
#>         interval_status
#> 1 heuristic_unvalidated
#> 2 heuristic_unvalidated
#> 3 heuristic_unvalidated

Fisher-z keeps correlation bounds inside $[-1,1]$ , but its effective-sample size approximation can still be optimistic for complicated latent or non-Gaussian models. A profile interval is not currently available because the required constrained refits cannot yet be checked reliably. Here interval_status = "heuristic_unvalidated" makes the Fisher-z limitation explicit: the bounds use the requested nominal level, but repeated-sampling coverage has not been established for this model class.

Phylogenetic signal is decomposition-specific

For a simple two-component decomposition, the package can profile a contrast of the two log-standard deviations. For richer phylogenetic plus non-phylogenetic decompositions, extract_phylo_signal() currently uses labelled numerical-Wald bounds. A general profile constraint and a working general bootstrap interval are not available for those richer decompositions. Do not relabel the fallback.

What profile likelihood means

For a scalar target $\theta$ and nuisance parameters $\eta$ , the profile log-likelihood is

$\ell_p(\theta_0) = \max_{\eta}\ell(\theta_0,\eta).$

The implemented direct route finds the values where the likelihood-ratio deviance crosses the standard $\chi^2_1$ threshold. Derived quantities such as correlations or variance proportions require a different constrained-refit problem because they are nonlinear functions of several fitted parameters. The package withholds those nonlinear penalty-profile routes until it can enforce the constraint tightly, verify every constrained optimizer result, and expose failure diagnostics. That distinction explains why one target may profile directly while another currently offers only point, Wald, or bootstrap routes.

Troubleshoot the interval that came back

Before diagnosing a profile, confirm that a profile actually ran. Read the returned method and status first. For withdrawn nonlinear targets such as repeatability and cross-family correlations, a profile request stops with a typed error instead of silently substituting a Wald result.

Observed result	First check	Interpretation and safe next action
Endpoint equals a natural limit such as `0`, `1`, or `-1`	Confirm the target transformation and returned method.	This can be a valid one-sided result. Report the boundary rather than replacing it with `NA`.
`NA` endpoint	Check fit health, route messages, and whether finite profile evaluations exist on that side.	The endpoint is unavailable or the numerical profile failed. Do not reinterpret it as a natural boundary.
`Inf` or `-Inf`	Inspect the parameter scale, search range, and profile shape.	The threshold may not have been crossed, the transformed natural limit may be infinite, or the search may have failed. Rule out numerical truncation before claiming the parameter is substantively unbounded.
Equal lower and upper bounds	Inspect the status for a pinned value, structural zero, or unavailable route.	This need not be an interval estimate. Report why the value is fixed.
Profile much wider than Wald	Inspect curve shape and fit health.	Disagreement is a diagnostic signal, not proof that either method is correct. Flatness can reflect weak information; jaggedness can reflect failed constrained refits.
Finite bootstrap percentile after an unbounded or failed profile	Check unconditional simulation support, failed-refit count, successful-draw distribution, and Monte Carlo resolution.	A finite percentile from finite successful draws does not prove identification or repair a broken generative simulation.

For a direct target, confint_inspect() returns the profile curve, endpoint comparison, and diagnostic flags without requiring an internal engine parameter name:

inspection <- confint_inspect(fit, parm = "sd_B[1]")
inspection$bounds
#>                 parm estimate_natural lower_natural upper_natural
#> theta_diag_B sd_B[1]      0.003004326             0     0.1120251
#>              wald_lower_natural wald_upper_natural wald_profile_disagree_lower
#> theta_diag_B                  0                Inf                       FALSE
#>              wald_profile_disagree_upper
#> theta_diag_B                          NA
inspection$diagnostics
#> [1] "hits_lower_bound"

confint_inspect() is for direct targets listed by profile_targets(). Derived correlations, communalities, proportions, and latent effects do not currently have public nonlinear profile curves; a raw direct-parameter profile does not diagnose a nonlinear derived-target constraint.

When covariance components remain confounded, simplify the model or collect more informative data. Loading constraints can choose an orientation or encode a prespecified scientific structure, but they do not create information that separates redundant variance components. Bootstrap likewise depends on a credible fitted generative model and successful refits; it is not a universal rescue for multimodality, weak identification, or an unsupported simulation route.

References

Pawitan, Y. (2001). In All Likelihood: Statistical Modelling and Inference Using Likelihood. Oxford University Press, chapter 9.
Venzon, D. J. & Moolgavkar, S. H. (1988). A method for computing profile-likelihood-based confidence intervals. Applied Statistics 37, 87–94. https://doi.org/10.2307/2347496
Self, S. G. & Liang, K.-Y. (1987). Asymptotic properties of maximum likelihood estimators and likelihood-ratio tests under nonstandard conditions. Journal of the American Statistical Association 82, 605–610. https://doi.org/10.2307/2289471
Kristensen, K., Nielsen, A., Berg, C. W., Skaug, H., & Bell, B. M. (2016). TMB: Automatic Differentiation and Laplace Approximation. Journal of Statistical Software 70(5). https://doi.org/10.18637/jss.v070.i05