hmclearn-glm-posterior: Sample log posterior and gradient functions for select...

In hmclearn: Fit Statistical Models Using Hamiltonian Monte Carlo

Description

These functions can be used to fit common generalized linear models and mixed effect models. See the accompanying vignettes for details on the derivations of the log posterior and gradient. In addition, these functions can be used as templates to build custom models to fit using HMC.

Usage

linear_posterior(theta, y, X, a = 1e-04, b = 1e-04, sig2beta = 1000)

g_linear_posterior(theta, y, X, a = 1e-04, b = 1e-04, sig2beta = 1000)

logistic_posterior(theta, y, X, sig2beta = 1000)

g_logistic_posterior(theta, y, X, sig2beta = 1000)

poisson_posterior(theta, y, X, sig2beta = 1000)

g_poisson_posterior(theta, y, X, sig2beta = 1000)

lmm_posterior(
  theta,
  y,
  X,
  Z,
  n,
  d,
  nrandom = 1,
  nugamma = 1,
  nuxi = 1,
  Agamma = 25,
  Axi = 25,
  sig2beta = 1000
)

g_lmm_posterior(
  theta,
  y,
  X,
  Z,
  n,
  d,
  nrandom = 1,
  nugamma = 1,
  nuxi = 1,
  Agamma = 25,
  Axi = 25,
  sig2beta = 1000
)

glmm_bin_posterior(
  theta,
  y,
  X,
  Z,
  n,
  nrandom = 1,
  nuxi = 1,
  Axi = 25,
  sig2beta = 1000
)

g_glmm_bin_posterior(
  theta,
  y,
  X,
  Z,
  n,
  nrandom = 1,
  nuxi = 1,
  Axi = 25,
  sig2beta = 1000
)

glmm_poisson_posterior(
  theta,
  y,
  X,
  Z,
  n,
  nrandom = 1,
  nuxi = 1,
  Axi = 25,
  sig2beta = 1000
)

g_glmm_poisson_posterior(
  theta,
  y,
  X,
  Z,
  n,
  nrandom = 1,
  nuxi = 1,
  Axi = 25,
  sig2beta = 1000
)

Arguments

theta	vector of parameters. See details below for the order of parameters for each model
y	numeric vector for the dependent variable for all models
X	numeric design matrix of fixed effect parameters for all models
a	hyperparameter for the Inverse Gamma shape parameter for σ_ε in linear regression models
b	hyperparameter for the Inverse Gamma scale parameter for σ_ε in linear regression models
sig2beta	diagonal covariance of prior for linear predictors is multivariate normal with mean 0 for linear regression and linear mixed effect models.
Z	numeric design matrix of random effect parameters for all mixed effects models
n	number of observations for standard glm models, or number of subjects for all mixed effect models
d	number of observations per subject for mixed effects models, but an input for linear mixed effect models only.
nrandom	number of random effects covariance parameters for all mixed effects models
nugamma	hyperparameter ν for the half-t prior of the log transformed error for linear mixed effects model γ
nuxi	hyperparameter ν for the half-t prior of the random effects diagonal for all mixed effects models ξ
Agamma	hyperparameter A for the half-t prior of the log transformed error for linear mixed effects model γ
Axi	hyperparameter A for the half-t prior of the random effects diagonal for all mixed effects modelsξ

Value

Numeric vector for the log posterior or gradient of the log posterior

Generalized Linear Models with available posterior and gradient functions

'linear_posterior(theta, y, X, a=1e-4, b=1e-4, sig2beta = 1e3)'

The log posterior function for linear regression

f(y | X, β, σ) = \frac{1}{(2πσ^2)^{n/2}}\exp{≤ft(-\frac{1}{2σ^2} (y - Xβ)^T(y-Xβ) \right)}

with priors p(σ^2) \sim IG(a, b) and β \sim N(0, σ_β^2 I). The variance term is log transformed γ = \logσ The input parameter vector theta is of length k. The first k-1 parameters are for β, and the last parameter is γ Note that the Inverse Gamma prior can be problematic for certain applications with low variance, such as hierarchical models. See Gelman (2006)

'g_linear_posterior(theta, y, X, a = 1e-04, b = 1e-04, sig2beta=1e3)'

Gradient of the log posterior for a linear regression model with Normal prior for the linear parameters and Inverse Gamma for the error term.

f(y | X, β, σ) = \frac{1}{(2πσ^2)^{n/2}}\exp{≤ft(-\frac{1}{2σ^2} (y - Xβ)^T(y-Xβ) \right)}

with priors p(σ^2) \sim IG(a, b) and β \sim N(0, σ_β^2 I). The variance term is log transformed γ = \logσ The input parameter vector theta is of length k. The first k-1 parameters are for β, and the last parameter is γ Note that the Inverse Gamma prior can be problematic for certain applications with low variance, such as hierarchical models. See Gelman (2006)

'logistic_posterior(theta, y, X, sig2beta=1e3) '

Log posterior for a logistic regression model with Normal prior for the linear parameters. The likelihood function for logistic regression

f(β| X, y) = ∏_{i=1}^{n} ≤ft(\frac{1}{1+e^{-X_iβ}}\right)^{y_i} ≤ft(\frac{e^{-X_iβ}}{1+e^{-X_iβ}}\right)^{1-y_i}

with priors β \sim N(0, σ_β^2 I). The input parameter vector theta is of length k, containing parameter values for β

'g_logistic_posterior(theta, y, X, sig2beta=1e3) '

Gradient of the log posterior for a logistic regression model with Normal prior for the linear parameters. The likelihood function for logistic regression

f(β| X, y) = ∏_{i=1}^{n} ≤ft(\frac{1}{1+e^{-X_iβ}}\right)^{y_i} ≤ft(\frac{e^{-X_iβ}}{1+e^{-X_iβ}}\right)^{1-y_i}

with priors β \sim N(0, σ_β^2 I). The input parameter vector theta is of length k, containing parameter values for β

'poisson_posterior(theta, y, X, sig2beta=1e3) '

Log posterior for a Poisson regression model with Normal prior for the linear parameters. The likelihood function for poisson regression

f(β| y, X) = ∏_{i=1}^n \frac{e^{-e^{X_iβ}}e^{y_iX_iβ}}{y_i!}

with priors β \sim N(0, σ_β^2 I). The input parameter vector theta is of length k, containing parameter values for β

'g_poisson_posterior(theta, y, X, sig2beta=1e3) '

Gradient of the log posterior for a Poisson regression model with Normal prior for the linear parameters. The likelihood function for poisson regression

f(β| y, X) = ∏_{i=1}^n \frac{e^{-e^{X_iβ}}e^{y_iX_iβ}}{y_i!}

with priors β \sim N(0, σ_β^2 I). The input parameter vector theta is of length k, containing parameter values for β

Generalized Linear Mixed Effect with available posterior and gradient functions

'lmm_posterior(theta, y, X, Z, n, d, nrandom = 1, nueps = 1, nuxi = 1, Aeps = 25, Axi = 25, sig2beta = 1e3) '

The log posterior function for linear mixed effects regression

f(y | β, u, σ_ε) \propto (σ_ε^2)^{-nd/2} e^{-\frac{1}{2σ_ε^2}(y - Xβ - Zu)^T (y - Xβ - Zu)}

with priors β \sim N(0, σ_β^2 I), σ_ε \sim half-t(A_ε, nu_ε), λ \sim half-t. The vector ξ is the diagonal of the covariance G log transformed hyperprior where u \sim N(0, G, ξ = \logλ and A_ξ, ν_ξ are parameters for the transformed distribution The standard deviation of the error is log transformed, where γ = \log σ_ε and σ_ε \sim half-t. The parameters for γ are A_γ, ν_γ The input parameter vector theta is of length k. The order of parameters for the vector is β, τ, γ, ξ.

'g_lmm_posterior(theta, y, X, Z, n, d, nrandom = 1, nueps = 1, nuxi = 1, Aeps = 25, Axi = 25, sig2beta = 1e3)'

Gradient of the log posterior for a linear mixed effects regression model

f(y | β, u, σ_ε) \propto (σ_ε^2)^{-n/2} e^{-\frac{1}{2σ_ε^2}(y - Xβ - Zu)^T (y - Xβ - Zu)}

with priors β \sim N(0, σ_β^2 I), σ_ε \sim half-t(A_ε, nu_ε), λ \sim half-t. The vector ξ is the diagonal of the covariance G log transformed hyperprior where u \sim N(0, G, ξ = \logλ and A_ξ, ν_ξ are parameters for the transformed distribution The standard deviation of the error is log transformed, where γ = \log σ_ε and σ_ε \sim half-t. The parameters for γ are A_γ, ν_γ The input parameter vector theta is of length k. The order of parameters for the vector is β, τ, γ, ξ

'glmm_bin_posterior(theta, y, X, Z, n, nrandom = 1, nuxi = 1, Axi = 25, sig2beta=1e3)'

The log posterior function for logistic mixed effects regression

f(y | X, Z, β, u) = ∏_{i=1}^n∏_{j=1}^d ≤ft(\frac{1}{1 + e^{-X_{i}β - Z_{ij}u_i}}\right)^{y_{ij}} ≤ft(\frac{e^{-X_iβ - Z_{ij}u_i}}{1 + e^{-X_{i}β - Z_{ij}u_i}}\right)^{1-y_{ij}}

with priors β \sim N(0, σ_β^2 I), σ_ε \sim half-t(A_ε, nu_ε), λ \sim half-t(A_λ, nu_λ ). The vector λ is the diagonal of the covariance G hyperprior where u \sim N(0, G, ξ = \logλ and A_ξ, ν_ξ are parameters for the transformed distribution The input parameter vector theta is of length k. The order of parameters for the vector is β, τ, ξ

'g_glmm_bin_posterior(theta, y, X, Z, n, nrandom = 1, nuxi = 1, Axi = 25, sig2beta = 1e3) '

Gradient of the log posterior function for logistic mixed effects regression

f(y | X, Z, β, u) = ∏_{i=1}^n∏_{j=1}^m ≤ft(\frac{1}{1 + e^{-X_{i}β - Z_{ij}u_i}}\right)^{y_{ij}} ≤ft(\frac{e^{-X_iβ - Z_{ij}u_i}}{1 + e^{-X_{i}β - Z_{ij}u_i}}\right)^{1-y_{ij}}

with priors β \sim N(0, σ_β^2 I), σ_ε \sim half-t(A_ε, nu_ε), λ \sim half-t(A_λ, nu_λ ). The vector λ is the diagonal of the covariance G hyperprior where u \sim N(0, G, ξ = \logλ and A_ξ, ν_ξ are parameters for the transformed distribution The input parameter vector theta is of length k. The order of parameters for the vector is β, τ, ξ

'glmm_poisson_posterior(theta, y, X, Z, n, nrandom = 1, nuxi = 1, Axi = 25, sig2beta = 1e3) '

Log posterior for a Poisson mixed effect regression

f(y | X, Z, β, u) = ∏_{i=1}^n ∏_{j=1}^m \frac{e^{-e^{X_iβ + Z_{ij}u_{ij}}}e^{y_i(X_iβ + Z_{ij}u_{ij})}}{y_i!}

with priors β \sim N(0, σ_β^2 I), σ_ε \sim half-t(A_ε, nu_ε), λ \sim half-t(A_λ, nu_λ ). The vector λ is the diagonal of the covariance G hyperprior where u \sim N(0, G, ξ = \logλ and A_ξ, ν_ξ are parameters for the transformed distribution The input parameter vector theta is of length k. The order of parameters for the vector is β, τ, ξ

'g_glmm_poisson_posterior(theta, y, X, Z, n, nrandom = 1, nuxi = 1, Axi = 25, sig2beta = 1e3) '

Gradient of the log posterior for a Poisson mixed effect regression

f(y | X, Z, β, u) = ∏_{i=1}^n ∏_{j=1}^m \frac{e^{-e^{X_iβ + Z_{ij}u_{ij}}}e^{y_i(X_iβ + Z_{ij}u_{ij})}}{y_i!}

with priors β \sim N(0, σ_β^2 I), σ_ε \sim half-t(A_ε, nu_ε), λ \sim half-t(A_λ, nu_λ ). The vector λ is the diagonal of the covariance G hyperprior where u \sim N(0, G, ξ = \logλ and A_ξ, ν_ξ are parameters for the transformed distribution The input parameter vector theta is of length k. The order of parameters for the vector is β, τ, ξ

References

Gelman, A. (2006). Prior distributions for variance parameters in hierarchical models (comment on article by Browne and Draper). Bayesian analysis, 1(3), 515-534.

Chan, J. C. C., & Jeliazkov, I. (2009). MCMC estimation of restricted covariance matrices. Journal of Computational and Graphical Statistics, 18(2), 457-480.

Betancourt, M., & Girolami, M. (2015). Hamiltonian Monte Carlo for hierarchical models. Current trends in Bayesian methodology with applications, 79, 30.

Examples

# Linear regression example
set.seed(521)
X <- cbind(1, matrix(rnorm(300), ncol=3))
betavals <- c(0.5, -1, 2, -3)
y <- X%*%betavals + rnorm(100, sd=.2)

f1_hmc <- hmc(N = 500,
          theta.init = c(rep(0, 4), 1),
          epsilon = 0.01,
          L = 10,
          logPOSTERIOR = linear_posterior,
          glogPOSTERIOR = g_linear_posterior,
          varnames = c(paste0("beta", 0:3), "log_sigma_sq"),
          param=list(y=y, X=X), parallel=FALSE, chains=1)

summary(f1_hmc, burnin=100)


# poisson regression example
set.seed(7363)
X <- cbind(1, matrix(rnorm(40), ncol=2))
betavals <- c(0.8, -0.5, 1.1)
lmu <- X %*% betavals
y <- sapply(exp(lmu), FUN = rpois, n=1)

f2_hmc <- hmc(N = 500,
          theta.init = rep(0, 3),
          epsilon = 0.01,
          L = 10,
          logPOSTERIOR = poisson_posterior,
          glogPOSTERIOR = g_poisson_posterior,
          varnames = paste0("beta", 0:2),
          param = list(y=y, X=X),
          parallel=FALSE, chains=1)

hmclearn documentation built on Oct. 23, 2020, 8:04 p.m.