Extract ordinal cutpoints from a fitted gllvmTMB model

For traits fitted with ordinal_probit() or ordinal_logit(), returns a tidy data frame with the K - 2 estimated cutpoints \(\tau_2, \ldots, \tau_{K-1}\) per trait, with optional standard errors from the joint sdreport.

Usage

extract_cutpoints(fit, quiet = FALSE)

Arguments

fit: A fit returned by gllvmTMB() with at least one ordinal_probit() or ordinal_logit() trait.
quiet: Suppress the note emitted when tau_se cannot be filled in because the fit has no sd_report. Internal callers that display only the cutpoint estimates (print(), tidy()) pass TRUE: a note explaining an absent standard-error column is noise when no such column is shown.

Value

A data frame with columns

trait: Trait label (factor level from data[[trait]]).
cutpoint_index: Integer index \(k \in \{2, \ldots, K-1\}\).
cutpoint_label: Character label "cutpoint_<k>".
tau_estimate: Estimated \(\tau_k\) on the trait's own latent scale (probit for an ordinal_probit() trait, logit for an ordinal_logit() trait).
tau_se: Standard error from the joint sdreport, or NA if the report is unavailable.

If the fit contains no ordinal_probit() or ordinal_logit() traits, returns a zero-row data frame with the same columns.

Details

Convention: gllvmTMB follows Hadfield (2015) — \(\tau_1 = 0\) is fixed for identifiability and the K - 2 free cutpoints are reported as cutpoint_2, cutpoint_3, etc. This differs from brms, which reports K - 1 cutpoints as Intercept[1..K-1].

This extractor includes Wald standard errors, so it refuses fits made with non-unit likelihood weights. The fitted cutpoint parameters remain point estimates, but weighted-objective uncertainty is not calibrated.

References

Hadfield, J. D. (2015). Increasing the efficiency of MCMC for hierarchical phylogenetic models of categorical traits using reduced mixed models. Methods Ecol. Evol. 6:706-714. doi:10.1111/2041-210X.12354

Examples

## Three ordinal traits scored on 5 plants at each of 30 sites: vigour
## (3 categories), leaf damage (4) and flowering (3). Each score comes
## from a hidden continuous "liability" -- the trait mean, plus a site
## effect, plus noise with standard deviation 1 -- cut at thresholds
## (stored in cuts, starting at tau_1 = 0): score 1 if the liability is
## below tau_1, score 2 if it lies between tau_1 and tau_2, and so on.
set.seed(1)
cuts <- list(vigour = c(0, 0.8), damage = c(0, 0.6, 1.4),
             flower = c(0, 1))
mu <- c(vigour = 0.3, damage = 0.2, flower = -0.2)
n_site <- 30
dat <- expand.grid(plant = 1:5, site = seq_len(n_site),
                   trait = names(cuts))
site_eff <- matrix(rnorm(n_site * 3, sd = 0.8), n_site, 3)
liability <- mu[as.character(dat$trait)] +
  site_eff[cbind(dat$site, as.integer(dat$trait))] + rnorm(nrow(dat))
dat$value <- NA_integer_
for (tr in names(cuts)) {
  rows <- dat$trait == tr
  dat$value[rows] <- findInterval(liability[rows], cuts[[tr]]) + 1L
}
dat$site <- factor(dat$site)

## Ordinal responses are coded 1, 2, ..., K. indep() gives each trait
## its own between-site variance; the 5 plants per site make it
## estimable.
fit <- gllvmTMB(
  value ~ 0 + trait + indep(0 + trait | site),
  data   = dat,
  trait  = "trait",
  unit   = "site",
  family = ordinal_probit()
)

## The cutpoints are the thresholds on the liability (probit) scale.
## The first threshold is fixed at tau_1 = 0 to identify the model, so
## a trait with K categories has K - 2 estimated cutpoints: one each
## for vigour and flowering, two for damage. tau_se is a large-sample
## (Wald) standard error.
cp <- extract_cutpoints(fit)
## Add the thresholds used to simulate the data, for comparison
cp$tau_true <- c(cuts$vigour[2], cuts$damage[2:3], cuts$flower[2])
cp
#>    trait cutpoint_index cutpoint_label tau_estimate     tau_se tau_true
#> 1 vigour              2     cutpoint_2    0.7857366 0.12111462      0.8
#> 2 damage              2     cutpoint_2    0.5256471 0.09667978      0.6
#> 3 damage              3     cutpoint_3    1.3259814 0.15039542      1.4
#> 4 flower              2     cutpoint_2    0.8211681 0.12714513      1.0

## Standard errors are optional. A fit made with se = FALSE skips them
## and gives the same cutpoints, with tau_se = NA. extract_cutpoints()
## then prints a note saying so; quiet = TRUE turns the note off.
fit_no_se <- gllvmTMB(
  value ~ 0 + trait + indep(0 + trait | site),
  data    = dat,
  trait   = "trait",
  unit    = "site",
  family  = ordinal_probit(),
  control = gllvmTMBcontrol(se = FALSE)
)
extract_cutpoints(fit_no_se, quiet = TRUE)
#>    trait cutpoint_index cutpoint_label tau_estimate tau_se
#> 1 vigour              2     cutpoint_2    0.7857366     NA
#> 2 damage              2     cutpoint_2    0.5256471     NA
#> 3 damage              3     cutpoint_3    1.3259814     NA
#> 4 flower              2     cutpoint_2    0.8211681     NA