1########################################################################## 2## MCMChpoisson.R 3## 4## MCMChpoisson() samples from the posterior distribution of a 5## Poisson hierarchical linear regression model in R using linked C++ 6## code in Scythe 7## 8## The code uses Algorithm 2 of Chib & Carlin (1999) for efficient 9## inference of (\beta | Y, sigma^2, Vb). 10## 11## Chib, S. & Carlin, B. P. (1999) On MCMC sampling in hierarchical 12## longitudinal models, Statistics and Computing, 9, 17-26 13## 14#################################################################### 15## 16## Original code by Ghislain Vieilledent, may 2011 17## CIRAD UR B&SEF 18## ghislain.vieilledent@cirad.fr / ghislainv@gmail.com 19## 20#################################################################### 21## 22## The initial version of this file was generated by the 23## auto.Scythe.call() function in the MCMCpack R package 24## written by: 25## 26## Andrew D. Martin 27## Dept. of Political Science 28## Washington University in St. Louis 29## admartin@wustl.edu 30## 31## Kevin M. Quinn 32## Dept. of Government 33## Harvard University 34## kevin_quinn@harvard.edu 35## 36## This software is distributed under the terms of the GNU GENERAL 37## PUBLIC LICENSE Version 2, June 1991. See the package LICENSE 38## file for more information. 39## 40## Copyright (C) 2011 Andrew D. Martin and Kevin M. Quinn 41## 42#################################################################### 43## 44## Revisions: 45## - G. Vieilledent, May 9 2011 [initial file] 46## 47#################################################################### 48 49 50#' Markov Chain Monte Carlo for the Hierarchical Poisson Linear Regression 51#' Model using the log link function 52#' 53#' MCMChpoisson generates a sample from the posterior distribution of a 54#' Hierarchical Poisson Linear Regression Model using the log link function and 55#' Algorithm 2 of Chib and Carlin (1999). This model uses a multivariate Normal 56#' prior for the fixed effects parameters, an Inverse-Wishart prior on the 57#' random effects variance matrix, and an Inverse-Gamma prior on the variance 58#' modelling over-dispersion. The user supplies data and priors, and a sample 59#' from the posterior distribution is returned as an mcmc object, which can be 60#' subsequently analyzed with functions provided in the coda package. 61#' 62#' 63#' \code{MCMChpoisson} simulates from the posterior distribution sample using 64#' the blocked Gibbs sampler of Chib and Carlin (1999), Algorithm 2. The 65#' simulation is done in compiled C++ code to maximize efficiency. Please 66#' consult the coda documentation for a comprehensive list of functions that 67#' can be used to analyze the posterior sample. 68#' 69#' The model takes the following form: 70#' 71#' \deqn{y_i \sim \mathcal{P}oisson(\lambda_i)} 72#' 73#' With latent variables \eqn{\phi(\lambda_i)}, \eqn{\phi} being the 74#' log link function: 75#' 76#' \deqn{\phi(\lambda_i) = X_i \beta + W_i b_i + \varepsilon_i} 77#' 78#' Where each group \eqn{i} have \eqn{k_i} observations. 79#' 80#' Where the random effects: 81#' 82#' \deqn{b_i \sim \mathcal{N}_q(0,V_b)} 83#' 84#' And the over-dispersion terms: 85#' 86#' \deqn{\varepsilon_i \sim \mathcal{N}(0, \sigma^2 I_{k_i})} 87#' 88#' We assume standard, conjugate priors: 89#' 90#' \deqn{\beta \sim \mathcal{N}_p(\mu_{\beta},V_{\beta})} 91#' 92#' And: 93#' 94#' \deqn{\sigma^{2} \sim \mathcal{IG}amma(\nu, 1/\delta)} 95#' 96#' And: 97#' 98#' \deqn{V_b \sim \mathcal{IW}ishart(r, rR)} See Chib and Carlin (1999) for more 99#' details. 100#' 101#' \emph{NOTE:} We do not provide default parameters for the priors on the 102#' precision matrix for the random effects. When fitting one of these models, 103#' it is of utmost importance to choose a prior that reflects your prior 104#' beliefs about the random effects. Using the \code{dwish} and \code{rwish} 105#' functions might be useful in choosing these values. 106#' 107#' @param fixed A two-sided linear formula of the form 'y~x1+...+xp' describing 108#' the fixed-effects part of the model, with the response on the left of a '~' 109#' operator and the p fixed terms, separated by '+' operators, on the right. 110#' Response variable y must be 0 or 1 (Binomial process). 111#' 112#' @param random A one-sided formula of the form '~x1+...+xq' specifying the 113#' model for the random effects part of the model, with the q random terms, 114#' separated by '+' operators. 115#' 116#' @param group String indicating the name of the grouping variable in 117#' \code{data}, defining the hierarchical structure of the model. 118#' 119#' @param data A data frame containing the variables in the model. 120#' 121#' @param burnin The number of burnin iterations for the sampler. 122#' 123#' @param mcmc The number of Gibbs iterations for the sampler. Total number of 124#' Gibbs iterations is equal to \code{burnin+mcmc}. \code{burnin+mcmc} must be 125#' divisible by 10 and superior or equal to 100 so that the progress bar can be 126#' displayed. 127#' 128#' @param thin The thinning interval used in the simulation. The number of mcmc 129#' iterations must be divisible by this value. 130#' 131#' @param seed The seed for the random number generator. If NA, the Mersenne 132#' Twister generator is used with default seed 12345; if an integer is passed 133#' it is used to seed the Mersenne twister. 134#' 135#' @param verbose A switch (0,1) which determines whether or not the progress 136#' of the sampler is printed to the screen. Default is 1: a progress bar is 137#' printed, indicating the step (in \%) reached by the Gibbs sampler. 138#' 139#' @param beta.start The starting values for the \eqn{\beta} vector. This 140#' can either be a scalar or a p-length vector. The default value of NA will 141#' use the OLS \eqn{\beta} estimate of the corresponding Gaussian Linear 142#' Regression without random effects. If this is a scalar, that value will 143#' serve as the starting value mean for all of the betas. 144#' 145#' @param sigma2.start Scalar for the starting value of the residual error 146#' variance. The default value of NA will use the OLS estimates of the 147#' corresponding Gaussian Linear Regression without random effects. 148#' 149#' @param Vb.start The starting value for variance matrix of the random 150#' effects. This must be a square q-dimension matrix. Default value of NA uses 151#' an identity matrix. 152#' 153#' @param mubeta The prior mean of \eqn{\beta}. This can either be a 154#' scalar or a p-length vector. If this takes a scalar value, then that value 155#' will serve as the prior mean for all of the betas. The default value of 0 156#' will use a vector of zeros for an uninformative prior. 157#' 158#' @param Vbeta The prior variance of \eqn{\beta}. This can either be a 159#' scalar or a square p-dimension matrix. If this takes a scalar value, then 160#' that value times an identity matrix serves as the prior variance of beta. 161#' Default value of 1.0E6 will use a diagonal matrix with very large variance 162#' for an uninformative flat prior. 163#' 164#' @param r The shape parameter for the Inverse-Wishart prior on variance 165#' matrix for the random effects. r must be superior or equal to q. Set r=q for 166#' an uninformative prior. See the NOTE for more details. 167#' 168#' @param R The scale matrix for the Inverse-Wishart prior on variance matrix 169#' for the random effects. This must be a square q-dimension matrix. Use 170#' plausible variance regarding random effects for the diagonal of R. See the 171#' NOTE for more details. 172#' 173#' @param nu The shape parameter for the Inverse-Gamma prior on the residual 174#' error variance. Default value is \code{nu=delta=0.001} for uninformative 175#' prior. 176#' 177#' @param delta The rate (1/scale) parameter for the Inverse-Gamma prior on the 178#' residual error variance. Default value is \code{nu=delta=0.001} for 179#' uninformative prior. 180#' 181#' @param FixOD A switch (0,1) which determines whether or not the variance for 182#' over-dispersion (sigma2) should be fixed (1) or not (0). Default is 0, 183#' parameter sigma2 is estimated. If FixOD=1, sigma2 is fixed to the value 184#' provided for \code{sigma2.start}. 185#' 186#' @param ... further arguments to be passed 187#' 188#' @return 189#' 190#' \item{mcmc}{An mcmc object that contains the posterior sample. This 191#' object can be summarized by functions provided by the coda 192#' package. The posterior sample of the deviance \eqn{D}, with 193#' \eqn{D=-2\log(\prod_i P(y_i|\lambda_i))}, is also provided.} 194#' 195#' \item{lambda.pred}{Predictive posterior mean for the exponential of 196#' the latent variables. The approximation of Diggle et al. (2004) is 197#' used to marginalized with respect to over-dispersion terms: 198#' 199#' \deqn{E[\lambda_i|\beta,b_i,\sigma^2]=\phi^{-1}((X_i \beta+W_i b_i)+0.5 \sigma^2)}} 200#' 201#' @export 202#' 203#' @author Ghislain Vieilledent <ghislain.vieilledent@@cirad.fr> 204#' 205#' @seealso \code{\link[coda]{plot.mcmc}}, \code{\link[coda]{summary.mcmc}} 206#' 207#' @references Siddhartha Chib and Bradley P. Carlin. 1999. ``On MCMC Sampling 208#' in Hierarchical Longitudinal Models.'' \emph{Statistics and Computing.} 9: 209#' 17-26. 210#' 211#' Daniel Pemstein, Kevin M. Quinn, and Andrew D. Martin. 2007. \emph{Scythe 212#' Statistical Library 1.0.} \url{http://scythe.lsa.umich.edu}. 213#' 214#' Andrew D. Martin and Kyle L. Saunders. 2002. ``Bayesian Inference for 215#' Political Science Panel Data.'' Paper presented at the 2002 Annual Meeting 216#' of the American Political Science Association. 217#' 218#' Martyn Plummer, Nicky Best, Kate Cowles, and Karen Vines. 2006. ``Output 219#' Analysis and Diagnostics for MCMC (CODA)'', \emph{R News}. 6(1): 7-11. 220#' \url{https://CRAN.R-project.org/doc/Rnews/Rnews_2006-1.pdf}. 221#' 222#' @keywords models hierarchical models mixed models glmm Poisson MCMC bayesian 223#' 224#' @examples 225#' 226#' \dontrun{ 227#' #======================================== 228#' # Hierarchical Poisson Linear Regression 229#' #======================================== 230#' 231#' #== Generating data 232#' 233#' # Constants 234#' nobs <- 1000 235#' nspecies <- 20 236#' species <- c(1:nspecies,sample(c(1:nspecies),(nobs-nspecies),replace=TRUE)) 237#' 238#' # Covariates 239#' X1 <- runif(n=nobs,min=-1,max=1) 240#' X2 <- runif(n=nobs,min=-1,max=1) 241#' X <- cbind(rep(1,nobs),X1,X2) 242#' W <- X 243#' 244#' # Target parameters 245#' # beta 246#' beta.target <- matrix(c(0.1,0.1,0.1),ncol=1) 247#' # Vb 248#' Vb.target <- c(0.05,0.05,0.05) 249#' # b 250#' b.target <- cbind(rnorm(nspecies,mean=0,sd=sqrt(Vb.target[1])), 251#' rnorm(nspecies,mean=0,sd=sqrt(Vb.target[2])), 252#' rnorm(nspecies,mean=0,sd=sqrt(Vb.target[3]))) 253#' 254#' # Response 255#' lambda <- vector() 256#' Y <- vector() 257#' for (n in 1:nobs) { 258#' lambda[n] <- exp(X[n,]%*%beta.target+W[n,]%*%b.target[species[n],]) 259#' Y[n] <- rpois(1,lambda[n]) 260#' } 261#' 262#' # Data-set 263#' Data <- as.data.frame(cbind(Y,lambda,X1,X2,species)) 264#' plot(Data$X1,Data$lambda) 265#' 266#' #== Call to MCMChpoisson 267#' model <- MCMChpoisson(fixed=Y~X1+X2, random=~X1+X2, group="species", 268#' data=Data, burnin=5000, mcmc=1000, thin=1,verbose=1, 269#' seed=NA, beta.start=0, sigma2.start=1, 270#' Vb.start=1, mubeta=0, Vbeta=1.0E6, 271#' r=3, R=diag(c(0.1,0.1,0.1)), nu=0.001, delta=0.001, FixOD=1) 272#' 273#' #== MCMC analysis 274#' 275#' # Graphics 276#' pdf("Posteriors-MCMChpoisson.pdf") 277#' plot(model$mcmc) 278#' dev.off() 279#' 280#' # Summary 281#' summary(model$mcmc) 282#' 283#' # Predictive posterior mean for each observation 284#' model$lambda.pred 285#' 286#' # Predicted-Observed 287#' plot(Data$lambda,model$lambda.pred) 288#' abline(a=0,b=1) 289#' 290#' ## #Not run 291#' ## #You can also compare with lme4 results 292#' ## #== lme4 resolution 293#' ## library(lme4) 294#' ## model.lme4 <- lmer(Y~X1+X2+(1+X1+X2|species),data=Data,family="poisson") 295#' ## summary(model.lme4) 296#' ## plot(fitted(model.lme4),model$lambda.pred,main="MCMChpoisson/lme4") 297#' ## abline(a=0,b=1) 298#' } 299#' 300MCMChpoisson <- function (fixed, random, group, data, burnin=5000, 301 mcmc=10000, thin=10, 302 verbose=1, seed=NA, beta.start=NA, sigma2.start=NA, 303 Vb.start=NA, mubeta=0, Vbeta=1.0E6, r, R, 304 nu=0.001, delta=0.001, FixOD=0, ...) 305{ 306 307 #======== 308 # Basic checks 309 #======== 310 check.group.hmodels(group, data) 311 check.mcmc.parameters.hmodels(burnin, mcmc, thin) 312 check.verbose.hmodels(verbose) 313 check.FixOD.hmodels(FixOD) 314 check.offset(list(...)) 315 316 #======== 317 # Seed 318 #======== 319 seed <- form.seeds.hmodels(seed) 320 321 #======== 322 # Form response and model matrices 323 #======== 324 mf.fixed <- model.frame(formula=fixed,data=data) 325 X <- model.matrix(attr(mf.fixed,"terms"),data=mf.fixed) 326 Y <- model.response(mf.fixed) 327 check.Y.Poisson.hmodels(Y) 328 mf.random <- model.frame(formula=random,data=data) 329 W <- model.matrix(attr(mf.random,"terms"),data=mf.random) 330 331 #======== 332 # Model parameters 333 #======== 334 nobs <- nrow(X) 335 IdentGroup <- as.numeric(as.factor(as.character(data[,names(data)==as.character(group)])))-1 336 LevelsGroup <- sort(unique(IdentGroup+1)) 337 LevelsGroup.Name <- sort(unique(as.character(data[,names(data)==as.character(group)]))) 338 ngroup <- length(LevelsGroup) 339 np <- ncol(X) 340 nq <- ncol(W) 341 ngibbs <- mcmc+burnin 342 nthin <- thin 343 nburn <- burnin 344 nsamp <- mcmc/thin 345 346 #======== 347 # Form and check starting parameters 348 #======== 349 beta.start <- form.beta.start.hmodels(fixed,data,beta.start,np,family="poisson",defaults=NA) 350 sigma2.start <- form.sigma2.start.hmodels(fixed,data,sigma2.start,family="poisson") 351 Vb.start <- form.Vb.start.hmodels(Vb.start,nq) 352 353 #======== 354 # Form priors 355 #======== 356 mvn.prior <- form.mvn.prior.hmodels(mubeta,Vbeta,np) 357 mubeta <- mvn.prior[[1]] 358 Vbeta <- mvn.prior[[2]] 359 wishart.prior <- form.wishart.prior.hmodels(r,R,nq) 360 r <- wishart.prior[[1]] 361 R <- wishart.prior[[2]] 362 check.ig.prior.hmodels(nu,delta) 363 s1 <- nu 364 s2 <- delta 365 366 #======== 367 # Parameters to save 368 #======== 369 beta_vect <- rep(c(beta.start),each=nsamp) 370 Vb_vect <- rep(c(Vb.start),each=nsamp) 371 b_vect <- rep(0,nq*ngroup*nsamp) 372 V <- rep(sigma2.start,nsamp) 373 lambda_pred <- rep(0.5,nobs) 374 Deviance <- rep(0,nsamp) 375 376 #======== 377 # call C++ code to draw sample 378 #======== 379 Sample <- .C("cMCMChpoisson", 380 #= Constants and data 381 ngibbs=as.integer(ngibbs), nthin=as.integer(nthin), nburn=as.integer(nburn),## Number of iterations, burning and samples 382 nobs=as.integer(nobs), ngroup=as.integer(ngroup), ## Constants 383 np=as.integer(np), nq=as.integer(nq), ## Constants 384 IdentGroup=as.integer(IdentGroup), 385 Y_vect=as.double(c(Y)), ## Response variable 386 X_vect=as.double(c(X)), ## Covariates 387 W_vect=as.double(c(W)), ## Covariates 388 #= Parameters to save 389 beta_vect.nonconst=as.double(beta_vect), ## Fixed parameters of the regression 390 b_vect.nonconst=as.double(b_vect), ## Random effects on intercept and slope 391 Vb_vect.nonconst=as.double(Vb_vect), ## Variance-covariance of random effects 392 V.nonconst=as.double(V), ## Variance of residuals 393 #= Defining priors 394 mubeta_vect=as.double(c(mubeta)), Vbeta_vect=as.double(c(Vbeta)), 395 r=as.double(r), R_vect=as.double(c(R)), 396 s1_V=as.double(s1), s2_V=as.double(s2), 397 #= Diagnostic 398 Deviance.nonconst=as.double(Deviance), 399 lambda_pred.nonconst=as.double(lambda_pred), ## Predictive posterior mean 400 #= Seeds 401 seed=as.integer(seed), 402 #= Verbose 403 verbose=as.integer(verbose), 404 #= Overdispersion 405 FixOD=as.integer(FixOD), 406 PACKAGE="MCMCpack") 407 408 #= Matrix of MCMC samples 409 Matrix <- matrix(NA,nrow=nsamp,ncol=np+nq*ngroup+nq*nq+2) 410 names.fixed <- paste("beta.",colnames(X),sep="") 411 names.random <- paste("b.",rep(colnames(W),each=ngroup),".",rep(LevelsGroup.Name,nq),sep="") 412 names.variances <- c(paste("VCV.",colnames(W),".",rep(colnames(W),each=nq),sep=""),"sigma2") 413 colnames(Matrix) <- c(names.fixed,names.random,names.variances,"Deviance") 414 415 #= Filling-in the matrix 416 Matrix[,c(1:np)] <- matrix(Sample[[12]],ncol=np) 417 Matrix[,c((np+1):(np+nq*ngroup))] <- matrix(Sample[[13]],ncol=nq*ngroup) 418 Matrix[,c((np+nq*ngroup+1):(np+nq*ngroup+nq*nq))] <- matrix(Sample[[14]],ncol=nq*nq) 419 Matrix[,ncol(Matrix)-1] <- Sample[[15]] 420 Matrix[,ncol(Matrix)] <- Sample[[22]] 421 422 #= Transform Sample list in an MCMC object 423 MCMC <- mcmc(Matrix,start=nburn+1,end=ngibbs,thin=nthin) 424 425 #= Output 426 return (list(mcmc=MCMC,lambda.pred=Sample[[23]])) 427} 428 429#################################################################### 430## END 431#################################################################### 432