Skip to contents

Response-family constructors used by gllvmTMB

Source: R/families.R
families.Rd

Build family descriptors for likelihoods that are not supplied directly by stats::family(). The descriptor records the family name, link functions, and any fixed shape control consumed by gllvmTMB().

Usage

Beta(link = "logit")

lognormal(link = "log")

gengamma(link = "log")

gamma_mix(link = "log", p_extreme = NULL)

lognormal_mix(link = "log", p_extreme = NULL)

nbinom2_mix(link = "log", p_extreme = NULL)

nbinom2(link = "log")

nbinom1(link = "log")

truncated_poisson(link = "log")

truncated_nbinom2(link = "log")

truncated_nbinom1(link = "log")

student(link = "identity", df = NULL)

tweedie(link = "log", p = NULL)

censored_poisson(link = "log")

zi_poisson(link = "log")

zi_nbinom2(link = "log")

zi_binomial(link = "logit")

delta_gamma(link1, link2 = "log", type = c("standard", "poisson-link"))

delta_gamma_mix(link1 = "logit", link2 = "log", p_extreme = NULL)

delta_gengamma(link1, link2 = "log", type = c("standard", "poisson-link"))

delta_lognormal(link1, link2 = "log", type = c("standard", "poisson-link"))

delta_lognormal_mix(
  link1,
  link2 = "log",
  type = c("standard", "poisson-link"),
  p_extreme = NULL
)

delta_truncated_nbinom2(link1 = "logit", link2 = "log")

delta_truncated_nbinom1(link1 = "logit", link2 = "log")

delta_poisson_link_gamma(link1 = "log", link2 = "log")

delta_poisson_link_lognormal(link1 = "log", link2 = "log")

betabinomial(link = "logit")

delta_beta(link1 = "logit", link2 = "logit")

Arguments

A single link name. The fitting engine applies the narrower family-specific link restrictions documented in gllvmTMB().

p_extreme

NULL to estimate the mixture probability, or one finite number strictly between zero and one to fix it.

df

NULL to estimate Student-t degrees of freedom, or one finite number greater than one to fix it.

p

NULL to estimate the Tweedie power, or one finite number strictly between one and two to fix it.

link1

Link for the occurrence component of a delta model.

link2

Link for the positive component of a delta model.

type

Delta construction: "standard" uses a binary occurrence component, while "poisson-link" requests the alternative Poisson-link parameterisation.

Value

A family object. One-part descriptors contain family, link, linkfun, and linkinv, plus constructor-specific metadata where applicable. Delta descriptors contain their two component families and delta = TRUE.

Details

These functions construct descriptors; successful construction does not by itself mean that a route is in the dependable release core. Ordinary Gaussian, Poisson, NB2, and binomial Laplace routes have the strongest first-release evidence. Other families and delta, mixture, truncated, ordinal, or mixed-family routes have narrower evidence. Read vignette("current-limits", package = "gllvmTMB") before relying on one of those routes.

Mixed-family input

Supply a list of descriptors when response traits use different families. The data column family selects an element for each row. Set attr(families, "family_var") to select a differently named column.

Ordering contract. Name the list to match the selector column's values, e.g. list(gaussian = gaussian(), student = student()) – matching by name is unambiguous and always accepted. An unnamed list is matched positionally to the selector column's levels (levels() for a factor, alphabetical sort(unique(...)) otherwise); this is validated, not merely assumed, by checking that reading against each family's own name (e.g. a "student" level naturally means student()) agrees with where the list put it. When the two readings agree – which is true whenever the list is already written in level order – the fit proceeds silently. When they disagree, that disagreement is precisely the silent-swap failure mode: list(student(), gaussian()) against family = rep(c("student","gaussian"), each = n) previously paired the student rows with gaussian() and vice versa, with no warning and a converged fit; the list is refused with an error showing both readings. A selector column whose values carry no name evidence at all (arbitrary labels such as "count"/"binary") still uses list order, but the resolved pairing is reported once so it stays auditable.

Ordinarily the family may vary between traits but not within one: all rows of a given trait must share one family and link, and a fit that mixes them stops with an error. Poisson-log with binomial-logit, for instance, is not a coherent common scale merely because it converges.

Integrated multi-source models (experimental)

To integrate more than two sources – a portal stream, digitised literature records, checklists, a structured survey – declare every source and its observation law with isdm_sources(). The two-source contract below is the two-source case of the same rule and keeps working unchanged.

Integrated two-source models (experimental)

There is one admitted exception, for combining opportunistic presence-only records with a structured detection/non-detection survey of the same species. Here the family genuinely does vary within a trait: the portal rows are Poisson-log counts and the survey rows are Bernoulli with a complementary log-log link, which is what makes both arms consistent with one shared underlying intensity. Poisson and Bernoulli carry no dispersion parameter, so nothing per-trait becomes ambiguous.

To reach it, the fit must match this contract exactly:

  • family = list(gbif = poisson(), survey_pa = binomial("cloglog")), with attr(family, "family_var") <- "isdm_family";

  • a source column holding only "gbif" and "survey", with both present;

  • an isdm_family column equal to "gbif" on portal rows and "survey_pa" on survey rows.

Within that contract an offset() is admitted on the cloglog arm, where it is a known change-of-support (sampled area or visit effort) rather than an ordinary log-rate exposure. Anything short of the full contract keeps the ordinary one-family-per-trait refusal.

Everything such a fit reports is relative intensity. Presence-only data cannot identify absolute abundance, occupancy, or detectability, so those quantities are not estimated and must not be read off the output. This interface is experimental and may change.

Delta defaults and compatibility wrappers

When link1 is omitted, standard delta constructors use "logit" and Poisson-link constructors use "log". The older delta_poisson_link_gamma() and delta_poisson_link_lognormal() names remain as deprecated wrappers.

gengamma() uses the Prentice parameterisation. Its fitted gengamma_Q value controls shape; the limiting case at zero is lognormal.

Several constructors remain exported for compatibility although the current multivariate fitter rejects them. In particular, this applies to generalized-gamma, mixture, censored Poisson, truncated NB1, and most delta variants. A rejected constructor is not silently mapped to another family.

Names ending in _mix() describe two mean components that share a scale parameter. p_extreme is the probability assigned to the component with the larger mean; leave it NULL to estimate that probability.

NB2 has variance \(\mu + \mu^2 / \phi\); its extra variance therefore grows quadratically with the mean.

NB1 has variance \(\mu(1 + \phi)\), so its extra variance grows linearly with the mean.

student() estimates degrees of freedom unless df is supplied.

censored_poisson() is a right-censored Poisson family (log link only): each row is either an ordinary observed count, or is only known to satisfy Y >= C for some censoring limit C (e.g. a detection ceiling). Supply the response as cbind(y, censored) ~ ..., where column 1 is the observed count (uncensored rows) or the censoring limit C (censored rows), and column 2 is a strict {0, 1} right-censoring indicator per row. A plain y ~ ... response is also accepted and treats every row as uncensored (identical to poisson()). Left-censoring and interval-censoring are not supported – right-censoring only. fitted()/predict(type = "response") report the mean of the underlying uncensored Poisson process (exp(eta)), not a censoring-adjusted expectation, matching how truncated_poisson() reports its untruncated mean.

zi_poisson(), zi_nbinom2(), and zi_binomial() are TRUE zero-inflation mixtures: the count process (log-link mean, or logit-link probability for zi_binomial()) is active at every observation, including y = 0, on top of a per-trait, intercept-only structural-zero probability zi (logit link, no covariates and no random effects on the zero part). This is different from delta_lognormal()/delta_gamma(), which are hurdle models – a strictly positive component conditional on presence, with no second zero-generating process (see vignette("current-limits", package = "gllvmTMB") for guidance on choosing between these models). Boundary: the count process (its fixed effects, latent() structure, and every correlation gllvmTMB reports) is conditional on the non-structural component; zi itself carries no covariates, no random effects, and no reported interval. zi_nbinom2() reuses the ordinary nbinom2() per-trait dispersion convention (not a shared scalar across traits). zi_binomial() requires multi-trial data (cbind(successes, failures) or a trials column) with at least one row per trait carrying n_trials >= 2; single-trial (0/1) responses do not identify the mixture and are refused with binomial() named as the working alternative.

References

Prentice, R.L. 1974. A log gamma model and its maximum likelihood estimation. Biometrika 61(3): 539–544. doi:10.1093/biomet/61.3.539

Thorson, J.T., Stewart, I.J., and Punt, A.E. 2011. Accounting for fish shoals in single- and multi-species survey data using mixture distribution models. Can. J. Fish. Aquat. Sci. 68(9): 1681–1693. doi:10.1139/f2011-086 .

Hilbe, J. M. 2011. Negative binomial regression. Cambridge University Press.

Thorson, J.T. 2018. Three problems with the conventional delta-model for biomass sampling data, and a computationally efficient alternative. Canadian Journal of Fisheries and Aquatic Sciences, 75(9), 1369-1382. doi:10.1139/cjfas-2017-0266

Examples

Beta(link = "logit")
#> 
#> Family: Beta 
#> Link function: logit 
#> 
## Mixed-family example: per-row family chosen by a `family` column. Trait 1
## stays Gaussian; trait 2 is recorded as presence/absence.
sim <- simulate_site_trait(n_sites = 40, n_species = 10, n_traits = 2,
                           mean_species_per_site = 5,
                           Lambda_B = matrix(c(0.8, 0.6), 2, 1), seed = 1)
dat <- sim$data
is_bin <- dat$trait == levels(dat$trait)[2]
dat$value[is_bin] <- as.integer(dat$value[is_bin] > 0)
dat$family <- factor(ifelse(is_bin, "binomial", "gaussian"),
                     levels = c("gaussian", "binomial"))
fit <- gllvmTMB(
  value ~ 0 + trait + latent(0 + trait | site, d = 1, unique = TRUE),
  data   = dat,
  family = list(gaussian = gaussian(), binomial = binomial()),
  trait  = "trait",
  unit   = "site"
)
#> ℹ Skipping the default between-unit Psi for 1 binary / categorical-contrast /
#>   single-observation ordinal trait under the family-specific identifiability
#>   gate.
#> ℹ Affected trait: trait_2.
#> ℹ Single-trial binomial traits trait_2 have one 0/1 trial per (trait, unit)
#>   cell. Multi-trial data (`cbind(successes, failures)` or `weights = n_trials`)
#>   can identify the diagonal; an explicit `indep()` term is a separate
#>   deliberate model choice.
#> • Mapped `theta_diag_B[t]` and the corresponding `s_B` row off.
#> This message is displayed once per session.
fit
#> Stacked-trait gllvmTMB fit
#>   Traits = 2, site = 40 
#>   Covstructs: latent_unit, indep_unit 
#>   Fixed effects (b_fix): 2
#>   Per-trait link (mixed-family fit):
#>    trait     link
#>  trait_1 identity
#>  trait_2    logit
#>   ML log L = -361.465   convergence = 0   engine = Laplace
#>   Run gllvmTMB_diagnose(fit) for a full health check, or summary(fit) for parameter estimates.

## Integrated two-source model: portal counts + survey detections, one
## shared ecological predictor. `dat` is long, one row per
## (cell, species, source); `log_support` is log sampled area or effort,
## `bias` a portal-only recording-effort covariate.
set.seed(1)
dat <- expand.grid(
  cell_id = factor(paste0("c", 1:60)),
  trait   = factor(c("sp1", "sp2")),
  source  = c("gbif", "survey")
)
env <- rnorm(60)  # one value per cell
u <- rnorm(60)    # unmeasured cell effect shared by both species
dat$env <- env[as.integer(dat$cell_id)]
dat$isdm_gbif <- as.integer(dat$source == "gbif")
dat$bias <- ifelse(dat$source == "gbif", rnorm(nrow(dat)), 0)
dat$log_support <- ifelse(dat$source == "gbif", log(2), log(0.5))
eta <- c(0.2, -0.1)[dat$trait] + c(0.5, -0.4)[dat$trait] * dat$env +
  c(0.6, 0.4)[dat$trait] * u[as.integer(dat$cell_id)] +
  dat$log_support + dat$isdm_gbif * (-0.3 + 0.4 * dat$bias)
dat$value <- ifelse(dat$source == "gbif",
                    rpois(nrow(dat), exp(eta)),
                    rbinom(nrow(dat), 1, 1 - exp(-exp(eta))))
dat$isdm_family <- factor(
  ifelse(dat$source == "gbif", "gbif", "survey_pa"),
  levels = c("gbif", "survey_pa")
)
isdm_fam <- list(gbif = poisson(), survey_pa = binomial(link = "cloglog"))
attr(isdm_fam, "family_var") <- "isdm_family"
fit_isdm <- gllvmTMB(
  value ~ 0 + trait + trait:env + trait:isdm_gbif + trait:bias +
    offset(log_support) + latent(0 + trait | cell_id, d = 1, unique = FALSE),
  data   = dat,
  family = isdm_fam,
  trait  = "trait",
  unit   = "cell_id"
)
#> The integrated multi-source route is experimental.
#> ℹ It combines presence-only count streams with structured
#>   detection/non-detection data under one shared ecological linear predictor.
#> ℹ Everything it reports is relative intensity: presence-only data cannot
#>   identify absolute abundance, occupancy, or detectability, and this fit does
#>   not estimate them.
#> ℹ Source-specific spatial structure is only weakly identified on small designs;
#>   treat a source-only field as a nuisance adjustment unless your design is
#>   large enough to support it.
#> ! Give every presence-only arm its own reporting-rate term (an interaction with
#>   a source indicator). Without one, the arms share an absolute intercept and
#>   the fit implicitly claims the absolute intensity that presence-only data
#>   cannot identify.
#> → Check convergence and positive-definiteness on every fit, and expect this
#>   interface to change.
#> This message is displayed once per session.
lognormal(link = "log")
#> 
#> Family: lognormal 
#> Link function: log 
#> 
gamma_mix(link = "log")
#> 
#> Family: gamma_mix 
#> Link function: log 
#> 
lognormal_mix(link = "log")
#> 
#> Family: lognormal_mix 
#> Link function: log 
#> 
nbinom2_mix(link = "log")
#> 
#> Family: nbinom2_mix 
#> Link function: log 
#> 
nbinom2(link = "log")
#> 
#> Family: nbinom2 
#> Link function: log 
#> 
nbinom1(link = "log")
#> 
#> Family: nbinom1 
#> Link function: log 
#> 
truncated_poisson(link = "log")
#> 
#> Family: truncated_poisson 
#> Link function: log 
#> 
truncated_nbinom2(link = "log")
#> 
#> Family: truncated_nbinom2 
#> Link function: log 
#> 
truncated_nbinom1(link = "log")
#> 
#> Family: truncated_nbinom1 
#> Link function: log 
#> 
student(link = "identity") # estimate df
#> Student-t degrees of freedom parameter will be estimated. This used to be fixed
#> at 3 by default. To fix it, supply a value to `df` (e.g., `df = 3`).
#> 
#> Family: student 
#> Link function: identity 
#> 
student(link = "identity", df = 3) # fix df at 3
#> Student-t degrees of freedom parameter fixed at 3. To estimate it, set `df =
#> NULL`.
#> 
#> Family: student 
#> Link function: identity 
#> 
tweedie(link = "log")
#> 
#> Family: tweedie 
#> Link function: log 
#> 
censored_poisson(link = "log")
#> 
#> Family: censored_poisson 
#> Link function: log 
#> 
## A small simulated rank-1 GLLVM fit with one right-censored trait
## (censoring limit C = 6).
set.seed(1)
n_site <- 100L
n_trait <- 3L
u <- rnorm(n_site)
lambda <- c(0.6, -0.5, 0.4)
beta <- c(1.2, 0.9, 1.4)
C <- 6
eta <- outer(u, lambda) + matrix(beta, n_site, n_trait, byrow = TRUE)
mu <- exp(eta)
y_true <- matrix(rpois(n_site * n_trait, mu), n_site, n_trait)
censored <- y_true >= C
y_obs <- ifelse(censored, C, y_true)
dat <- data.frame(
  site     = factor(rep(seq_len(n_site), n_trait)),
  trait    = factor(rep(seq_len(n_trait), each = n_site)),
  y        = as.vector(y_obs),
  censored = as.integer(as.vector(censored))
)
fit <- gllvmTMB(
  cbind(y, censored) ~ 0 + trait + latent(0 + trait | site, d = 1, unique = FALSE),
  data = dat, family = censored_poisson(), trait = "trait", unit = "site"
)
fit
#> Stacked-trait gllvmTMB fit
#>   Traits = 3, site = 100 
#>   Covstructs: latent_unit 
#>   Fixed effects (b_fix): 3
#>   ML log L = -526.908   convergence = 0   engine = Laplace
#>   Run gllvmTMB_diagnose(fit) for a full health check, or summary(fit) for parameter estimates.
## A small simulated rank-1 GLLVM fit in which a quarter of the counts are
## structural zeros.
set.seed(1)
n_site <- 100L
n_trait <- 3L
u <- rnorm(n_site)
lambda <- c(0.6, -0.5, 0.4)
beta <- c(1.2, 0.9, 1.4)
pi_true <- 0.25
eta <- outer(u, lambda) + matrix(beta, n_site, n_trait, byrow = TRUE)
mu <- exp(eta)
z <- matrix(rbinom(n_site * n_trait, 1L, 1 - pi_true), n_site, n_trait)
y <- matrix(rpois(n_site * n_trait, mu), n_site, n_trait) * z
dat <- data.frame(
  site  = factor(rep(seq_len(n_site), n_trait)),
  trait = factor(rep(seq_len(n_trait), each = n_site)),
  y     = as.vector(y)
)
fit <- gllvmTMB(
  y ~ 0 + trait + latent(0 + trait | site, d = 1, unique = FALSE),
  data = dat, family = zi_poisson(), trait = "trait", unit = "site"
)
fit$report$zi  # per-trait fitted structural-zero probability
#> [1] 0.2215543 0.2794186 0.2616444
delta_gamma()
#> 
#> Family: binomial Gamma 
#> Link function: logit log 
#> 
delta_gamma_mix()
#> 
#> Family: binomial gamma_mix 
#> Link function: logit log 
#> 
delta_gengamma()
#> 
#> Family: binomial gengamma 
#> Link function: logit log 
#> 
delta_lognormal()
#> 
#> Family: binomial lognormal 
#> Link function: logit log 
#> 
delta_lognormal_mix()
#> 
#> Family: binomial lognormal_mix 
#> Link function: logit log 
#> 
delta_truncated_nbinom2()
#> 
#> Family: binomial truncated_nbinom2 
#> Link function: logit log 
#> 
delta_truncated_nbinom1()
#> 
#> Family: binomial truncated_nbinom1 
#> Link function: logit log 
#> 
betabinomial(link = "logit")
#> 
#> Family: betabinomial 
#> Link function: logit 
#> 
delta_beta()
#> 
#> Family: binomial Beta 
#> Link function: logit logit 
#>