hmclearn-glm-posterior.RdThese 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.
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 )
| 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\) |
Numeric vector for the log posterior or gradient of the log posterior
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)
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)}$$ 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)
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\)
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\)
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\)
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\)
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\).
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)}$$ 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\)
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\)
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}} $$ 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\)
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!} $$ 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\)
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!} $$ 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\)
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.
# 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)#> Summary of MCMC simulation #>#> 2.5% 5% 25% 50% 75% 95% #> beta0 0.4778379 0.4847789 0.518149 0.5321988 0.5474873 0.5726653 #> beta1 -1.0435956 -1.0398451 -1.022090 -1.0099795 -0.9969146 -0.9695882 #> beta2 1.9749648 1.9834760 2.003620 2.0201366 2.0330953 2.0521270 #> beta3 -3.0212212 -3.0146309 -2.996771 -2.9813806 -2.9668917 -2.9442674 #> log_sigma_sq -3.2887695 -3.2522864 -3.150465 -3.0570739 -2.9441514 -2.8255475 #> 97.5% #> beta0 0.5781760 #> beta1 -0.9622133 #> beta2 2.0554530 #> beta3 -2.9308313 #> log_sigma_sq -2.7737642# 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)