R: Sample log posterior and gradient functions for select...

hmclearn-glm-posterior {hmclearn}

R Documentation

Sample log posterior and gradient functions for select generalized linear models and mixed effect models

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 `\sigma_\epsilon` in linear regression models
`b`	hyperparameter for the Inverse Gamma scale parameter for `\sigma_\epsilon` 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 `\nu` for the half-t prior of the log transformed error for linear mixed effects model `\gamma`
`nuxi`	hyperparameter `\nu` for the half-t prior of the random effects diagonal for all mixed effects models `\xi`
`Agamma`	hyperparameter `A` for the half-t prior of the log transformed error for linear mixed effects model `\gamma`
`Axi`	hyperparameter `A` for the half-t prior of the random effects diagonal for all mixed effects models`\xi`

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, \beta, \sigma) = \frac{1}{(2\pi\sigma^2)^{n/2}}\exp{\left(-\frac{1}{2\sigma^2} (y - X\beta)^T(y-X\beta) \right)}

with priors p(\sigma^2) \sim IG(a, b) and \beta \sim N(0, \sigma_\beta^2 I). The variance term is log transformed \gamma = \log\sigma The input parameter vector theta is of length k. The first k-1 parameters are for \beta, and the last parameter is \gamma 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, \beta, \sigma) = \frac{1}{(2\pi\sigma^2)^{n/2}}\exp{\left(-\frac{1}{2\sigma^2} (y - X\beta)^T(y-X\beta) \right)}

'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(\beta| X, y) = \prod_{i=1}^{n} \left(\frac{1}{1+e^{-X_i\beta}}\right)^{y_i} \left(\frac{e^{-X_i\beta}}{1+e^{-X_i\beta}}\right)^{1-y_i}

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

'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(\beta| X, y) = \prod_{i=1}^{n} \left(\frac{1}{1+e^{-X_i\beta}}\right)^{y_i} \left(\frac{e^{-X_i\beta}}{1+e^{-X_i\beta}}\right)^{1-y_i}

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

'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(\beta| y, X) = \prod_{i=1}^n \frac{e^{-e^{X_i\beta}}e^{y_iX_i\beta}}{y_i!}

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

'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(\beta| y, X) = \prod_{i=1}^n \frac{e^{-e^{X_i\beta}}e^{y_iX_i\beta}}{y_i!}

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

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 | \beta, u, \sigma_\epsilon) \propto (\sigma_\epsilon^2)^{-nd/2} e^{-\frac{1}{2\sigma_\epsilon^2}(y - X\beta - Zu)^T (y - X\beta - Zu)}

with priors \beta \sim N(0, \sigma_\beta^2 I), \sigma_\epsilon \sim half-t(A_\epsilon, nu_\epsilon), \lambda \sim half-t. The vector \xi is the diagonal of the covariance G log transformed hyperprior where u \sim N(0, G, \xi = \log\lambda and A_\xi, \nu_\xi are parameters for the transformed distribution The standard deviation of the error is log transformed, where \gamma = \log \sigma_\epsilon and \sigma_\epsilon \sim half-t. The parameters for \gamma are A_\gamma, \nu_\gamma The input parameter vector theta is of length k. The order of parameters for the vector is \beta, \tau, \gamma, \xi.

'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 | \beta, u, \sigma_\epsilon) \propto (\sigma_\epsilon^2)^{-n/2} e^{-\frac{1}{2\sigma_\epsilon^2}(y - X\beta - Zu)^T (y - X\beta - Zu)}

'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, \beta, u) = \prod_{i=1}^n\prod_{j=1}^d \left(\frac{1}{1 + e^{-X_{i}\beta - Z_{ij}u_i}}\right)^{y_{ij}} \left(\frac{e^{-X_i\beta - Z_{ij}u_i}}{1 + e^{-X_{i}\beta - Z_{ij}u_i}}\right)^{1-y_{ij}}

with priors \beta \sim N(0, \sigma_\beta^2 I), \sigma_\epsilon \sim half-t(A_\epsilon, nu_\epsilon), \lambda \sim half-t(A_\lambda, nu_\lambda ). The vector \lambda is the diagonal of the covariance G hyperprior where u \sim N(0, G, \xi = \log\lambda and A_\xi, \nu_\xi are parameters for the transformed distribution The input parameter vector theta is of length k. The order of parameters for the vector is \beta, \tau, \xi

'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, \beta, u) = \prod_{i=1}^n\prod_{j=1}^m \left(\frac{1}{1 + e^{-X_{i}\beta - Z_{ij}u_i}}\right)^{y_{ij}} \left(\frac{e^{-X_i\beta - Z_{ij}u_i}}{1 + e^{-X_{i}\beta - Z_{ij}u_i}}\right)^{1-y_{ij}}

'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, \beta, u) = \prod_{i=1}^n \prod_{j=1}^m \frac{e^{-e^{X_i\beta + Z_{ij}u_{ij}}}e^{y_i(X_i\beta + Z_{ij}u_{ij})}}{y_i!}

'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, \beta, u) = \prod_{i=1}^n \prod_{j=1}^m \frac{e^{-e^{X_i\beta + Z_{ij}u_{ij}}}e^{y_i(X_i\beta + Z_{ij}u_{ij})}}{y_i!}

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)

[Package hmclearn version 0.0.5 Index]