Penalty / prior specification for a phylogenetic location-scale fit
Source:R/penalty.R
drm_phylo_penalty.RdBuilds an optional penalty (a weakly-informative prior) for the phylogenetic
standard deviations and, optionally, the phylogenetic cross-parameter
correlation of a drmTMB() fit. Passing the result to the penalty argument
of drmTMB() switches the estimator from plain maximum likelihood to a
penalized / maximum-a-posteriori (MAP) estimator.
Arguments
- sd_u, sd_alpha
Penalised-complexity prior scale and tail probability for each phylogenetic SD: a priori
P(sd > sd_u) = sd_alpha.sd_umust be positive andsd_alphamust lie in(0, 1).- cor_sd
Optional standard deviation of a mean-zero normal penalty on the phylogenetic cross-parameter correlation parameter.
NULL(the default) applies no correlation penalty. A non-NULLcor_sdrequires a coupled phylogenetic model with at least two phylogenetic SDs (for example a coupled location-scale or bivariate fit); a location-only phylogenetic model has a single phylogenetic SD and no correlation parameter to penalize, sodrmTMB()errors rather than silently ignoringcor_sd.
Details
The standard-deviation penalty is a penalised-complexity (PC) prior (Simpson
et al. 2017): an exponential prior on the SD scale with mass at zero, which
regularises a weakly-identified phylogenetic SD (for example a scale-side
phylogenetic field at about one observation per tip) toward the simpler
"no phylogenetic variance" model in the sense that the prior places its
mass there — the reported point estimate itself never reaches it; see
below. The rate is lambda = -log(sd_alpha) / sd_u
so that, a priori, P(sd > sd_u) = sd_alpha holds exactly for this prior.
The optional correlation penalty is a mean-zero normal on the unconstrained
phylogenetic correlation parameter.
The optimizer works on log(sd), not sd, so the penalty actually added to
the objective is the negative log-density of the exponential prior expressed
in log(sd), including the |d sd / d log(sd)| = sd change-of-variables
Jacobian: rate * sd - log(sd) - log(rate). Because a maximum-a-posteriori
(MAP) estimate is not invariant to reparameterisation, the fitted penalised
sd_phylo is the mode of this log-scale expression, not of the SD-scale
prior, and that mode is never zero: the -log(sd) Jacobian term diverges as
sd -> 0, so under a flat likelihood the penalty alone pulls the estimate to
sd = 1/rate (0.334 at the defaults sd_u = 1, sd_alpha = 0.05, where
rate = -log(0.05) = 2.996). The SD-scale prior still has mass at zero as
documented; only the reported point estimate does not sit there. A
penalised sd_phylo must therefore not be used to test a null hypothesis of
no phylogenetic signal: a likelihood-ratio or Wald test against zero on a
MAP fit is not valid, because the estimator cannot report zero regardless of
the data.
A penalized fit is a MAP point estimate, not a maximum-likelihood fit: its
standard errors are credible-interval-shaped, and likelihood-ratio tests or
AIC across penalized fits are not standard. logLik() returns the
unpenalized data log-likelihood; the penalty contribution is stored
separately on the fit as fit$phylo_penalty.
References
Simpson, D., Rue, H., Riebler, A., Martins, T. G., & Sorbye, S. H. (2017). Penalising model component complexity: a principled, practical approach to constructing priors. Statistical Science, 32(1), 1-28.
Chung, Y., Rabe-Hesketh, S., Dorie, V., Gelman, A., & Liu, J. (2013). A nondegenerate penalized likelihood estimator for variance parameters in multilevel models. Psychometrika, 78(4), 685-709.
Examples
# Penalised-complexity prior: a priori P(phylogenetic SD > 1) = 0.05.
pen <- drm_phylo_penalty(sd_u = 1, sd_alpha = 0.05)
pen$rate
#> [1] 2.995732
# Also penalize the phylogenetic correlation in a coupled location-scale or
# bivariate phylogenetic model.
pen_cor <- drm_phylo_penalty(sd_u = 1, sd_alpha = 0.05, cor_sd = 0.5)
pen_cor$cor_sd
#> [1] 0.5