[r-cran-mcmcpack] 58/90: Imported Upstream version 1.1-1
Andreas Tille
tille at debian.org
Fri Dec 16 09:07:47 UTC 2016
This is an automated email from the git hooks/post-receive script.
tille pushed a commit to branch master
in repository r-cran-mcmcpack.
commit 1c54db78ef10b96a25db54890844220fdc032aa8
Author: Andreas Tille <tille at debian.org>
Date: Fri Dec 16 08:07:20 2016 +0100
Imported Upstream version 1.1-1
---
DESCRIPTION | 8 +-
NAMESPACE | 4 +
R/MCMChierBetaBinom.R | 252 +++++++++++++++++++++++++
R/MCMChlogit.R | 181 ++++++++++++++++++
R/MCMChpoisson.R | 181 ++++++++++++++++++
R/MCMChregress.R | 177 ++++++++++++++++++
R/hidden-hmodels.R | 348 ++++++++++++++++++++++++++++++++++
data/PErisk.rda | Bin 1966 -> 1968 bytes
data/SupremeCourt.rda | Bin 471 -> 471 bytes
man/MCMCbinaryChange.Rd | 0
man/MCMChlogit.Rd | 292 +++++++++++++++++++++++++++++
man/MCMChpoisson.Rd | 279 ++++++++++++++++++++++++++++
man/MCMChregress.Rd | 254 +++++++++++++++++++++++++
src/MCMChierBetaBinom.cc | 458 +++++++++++++++++++++++++++++++++++++++++++++
src/MCMChlogit.cc | 475 +++++++++++++++++++++++++++++++++++++++++++++++
src/MCMChpoisson.cc | 463 +++++++++++++++++++++++++++++++++++++++++++++
src/MCMChregress.cc | 381 +++++++++++++++++++++++++++++++++++++
17 files changed, 3749 insertions(+), 4 deletions(-)
diff --git a/DESCRIPTION b/DESCRIPTION
index 8599f39..87286f6 100644
--- a/DESCRIPTION
+++ b/DESCRIPTION
@@ -1,6 +1,6 @@
Package: MCMCpack
-Version: 1.0-11
-Date: 2011-4-26
+Version: 1.1-1
+Date: 2011-7-8
Title: Markov chain Monte Carlo (MCMC) Package
Author: Andrew D. Martin, Kevin M. Quinn, and Jong Hee Park
Maintainer: Andrew D. Martin <admartin at wustl.edu>
@@ -18,6 +18,6 @@ Description: This package contains functions to perform Bayesian
License: GPL-2
SystemRequirements: gcc (>= 4.0)
URL: http://mcmcpack.wustl.edu
-Packaged: 2011-04-26 14:48:33 UTC; adm
+Packaged: 2011-07-08 17:10:47 UTC; adm
Repository: CRAN
-Date/Publication: 2011-04-26 15:12:35
+Date/Publication: 2011-07-23 13:02:20
diff --git a/NAMESPACE b/NAMESPACE
index 811315c..6ce78ee 100644
--- a/NAMESPACE
+++ b/NAMESPACE
@@ -20,7 +20,11 @@ export(
MCMCdynamicEI,
MCMCdynamicIRT1d,
MCMCfactanal,
+## MCMChierBetaBinom,
MCMChierEI,
+ MCMChlogit,
+ MCMChpoisson,
+ MCMChregress,
MCMCirt1d,
MCMCirtHier1d,
MCMCirtKd,
diff --git a/R/MCMChierBetaBinom.R b/R/MCMChierBetaBinom.R
new file mode 100644
index 0000000..d474586
--- /dev/null
+++ b/R/MCMChierBetaBinom.R
@@ -0,0 +1,252 @@
+###########################################################################
+## sample from the posterior distribution of a hierarchical beta binomial
+## model in R using linked C++ code in Scythe
+##
+## y_{ij} ~ Binomial(s_{ij}, theta_{ij})
+## theta_{ij} ~ Beta(alpha_j, beta_j)
+## alpha_j ~ Pareto(1, a)
+## beta_j ~ Pareto(1, b)
+##
+## KQ 5/24/2011 - 6/25/2011
+##
+## This software is distributed under the terms of the GNU GENERAL
+## PUBLIC LICENSE Version 2, June 1991. See the package LICENSE
+## file for more information.
+##
+##
+## Copyright (C) 2003-2007 Andrew D. Martin and Kevin M. Quinn
+## Copyright (C) 2007-present Andrew D. Martin, Kevin M. Quinn,
+## and Jong Hee Park
+##########################################################################
+
+
+"MCMChierBetaBinom" <- function(y, s, i.labels, j.labels, burnin=1000,
+ mcmc=10000, thin=1,
+ verbose=0, seed=NA, theta.start=NA,
+ alpha.start=NA, beta.start=NA,
+ a=0, b=0){
+
+ ## checks
+ check.mcmc.parameters(burnin, mcmc, thin)
+ max.length <- max(length(y), length(s), length(i.labels), length(j.labels))
+ min.length <- min(length(y), length(s), length(i.labels), length(j.labels))
+ if (max.length != min.length){
+ stop("y, s, i.labels, j.labels not all of same length\n")
+ }
+
+ na.indic <- is.na(y)
+ na.indic <- na.indic + is.na(s)
+ na.indic <- na.indic + is.na(i.labels)
+ na.indic <- na.indic + is.na(j.labels)
+ na.indic <- na.indic != 0
+ y <- y[!na.indic]
+ s <- s[!na.indic]
+ i.labels <- i.labels[!na.indic]
+ j.labels <- j.labels[!na.indic]
+
+ if (min(y) < 0){
+ stop("y cannot have an element less than 0\n")
+ }
+ if (max(y > s) > 0){
+ stop("y[i] cannot be greater than s[i] for any i\n")
+ }
+
+ ## more checking needed
+
+
+ ## seeds
+ seeds <- form.seeds(seed)
+ lecuyer <- seeds[[1]]
+ seed.array <- seeds[[2]]
+ lecuyer.stream <- seeds[[3]]
+
+
+
+ i.labels <- as.character(i.labels)
+ j.labels <- as.character(j.labels)
+ i.labels.unique <- unique(i.labels)
+ j.labels.unique <- unique(j.labels)
+ i.labels.num <- rep(NA, length(i.labels))
+ j.labels.num <- rep(NA, length(j.labels))
+ for (i in 1:length(i.labels)){
+ i.labels.num[i] <- which(i.labels.unique == i.labels[i])
+ }
+ for (j in 1:length(j.labels)){
+ j.labels.num[j] <- which(j.labels.unique == j.labels[j])
+ }
+ i.labels.num.unique <- unique(i.labels.num)
+ j.labels.num.unique <- unique(j.labels.num)
+
+
+
+ ## starting values
+ if (length(theta.start) == 1){
+ if (is.numeric(theta.start)){
+ theta.start <- rep(theta.start, length(y))
+ }
+ else if (is.na(theta.start)){
+ theta.start <- (y + 0.01) / (s + 0.02)
+ }
+ else{
+ theta.start <- rep(0.5, length(y))
+ }
+ }
+ if (length(theta.start) != length(y)){
+ stop("length(theta.start) != length(y)\n")
+ }
+ if (max(theta.start) >= 1.0){
+ stop("elements of theta.start must be less than 1.0\n")
+ }
+ if (min(theta.start) <= 0.0){
+ stop("elements of theta.start must be greater than 0.0\n")
+ }
+
+ if (length(alpha.start) == 1){
+ if (is.na(alpha.start)){
+ alpha.start <- rep(1.001, length(j.labels.unique))
+ }
+ }
+ if (length(alpha.start) != length(j.labels.unique)){
+ stop("length(alpha.start) != length(unique(j.labels))\n")
+ }
+ if (length(beta.start) == 1){
+ if (is.na(beta.start)){
+ beta.start <- rep(1.001, length(j.labels.unique))
+ }
+ }
+ if (length(beta.start) != length(j.labels.unique)){
+ stop("length(beta.start) != length(unique(j.labels))\n")
+ }
+
+
+ accepts <- rep(0, length(j.labels.unique))
+
+
+
+
+
+ ## get reasonable values for base.sigma
+ base.sigma <- rep(1, length(j.labels.unique))
+
+
+
+ ## call C++ code to draw the sample
+ sample <- matrix(data=0, mcmc/thin, (length(y) +
+ 2*length(j.labels.unique)) )
+
+
+ posterior <- .C("hierBetaBinom",
+ samdata = as.double(sample),
+ samrow = as.integer(nrow(sample)),
+ samcol = as.integer(ncol(sample)),
+ y = as.integer(y),
+ s = as.integer(s),
+ theta = as.double(theta.start),
+ alpha = as.double(alpha.start),
+ beta = as.double(beta.start),
+ a = as.double(a),
+ b = as.double(b),
+ ilabels = as.integer(i.labels.num),
+ jlabels = as.integer(j.labels.num),
+ ilabelsunique = as.integer(i.labels.num.unique),
+ jlabelsunique = as.integer(j.labels.num.unique),
+ n = as.integer(length(y)),
+ ni = as.integer(length(i.labels.unique)),
+ nj = as.integer(length(j.labels.unique)),
+ burnin = as.integer(burnin),
+ mcmc = as.integer(mcmc),
+ thin = as.integer(thin),
+ lecuyer = as.integer(lecuyer),
+ seedarray = as.integer(seed.array),
+ lecuyerstream = as.integer(lecuyer.stream),
+ verbose = as.integer(verbose),
+ accepts = as.integer(accepts),
+ basesigma = as.double(base.sigma),
+ PACKAGE="MCMCpack")
+
+
+ sample <- matrix(posterior$samdata,
+ posterior$samrow,
+ posterior$samcol,
+ byrow=FALSE)
+
+ output <- mcmc(data=sample, start=burnin+1, end=burnin+mcmc, thin=thin)
+
+ theta.names <- paste("theta.", i.labels, ".", j.labels, sep="")
+ alpha.names <- paste("alpha.", j.labels.unique, sep="")
+ beta.names <- paste("beta.", j.labels.unique, sep="")
+
+ varnames(output) <- c(theta.names, alpha.names, beta.names)
+
+ attr(output, "title") <- "MCMChierBetaBinom Posterior Sample"
+ attr(output, "acceptance.rates") <- posterior$accepts / (posterior$mcmc + posterior$burnin)
+ return(output)
+
+
+} ## end MCMChierBetaBinom
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
diff --git a/R/MCMChlogit.R b/R/MCMChlogit.R
new file mode 100644
index 0000000..5dda99f
--- /dev/null
+++ b/R/MCMChlogit.R
@@ -0,0 +1,181 @@
+##########################################################################
+## MCMChlogit.R
+##
+## MCMChlogit() samples from the posterior distribution of a
+## Binomial hierarchical linear regression model in R using linked C++
+## code in Scythe
+##
+## The code uses Algorithm 2 of Chib & Carlin (1999) for efficient
+## inference of (\beta | Y, sigma^2, Vb).
+##
+## Chib, S. & Carlin, B. P. (1999) On MCMC sampling in hierarchical
+## longitudinal models, Statistics and Computing, 9, 17-26
+##
+####################################################################
+##
+## Original code by Ghislain Vieilledent, may 2011
+## CIRAD UR B&SEF
+## ghislain.vieilledent at cirad.fr / ghislainv at gmail.com
+##
+####################################################################
+##
+## The initial version of this file was generated by the
+## auto.Scythe.call() function in the MCMCpack R package
+## written by:
+##
+## Andrew D. Martin
+## Dept. of Political Science
+## Washington University in St. Louis
+## admartin at wustl.edu
+##
+## Kevin M. Quinn
+## Dept. of Government
+## Harvard University
+## kevin_quinn at harvard.edu
+##
+## This software is distributed under the terms of the GNU GENERAL
+## PUBLIC LICENSE Version 2, June 1991. See the package LICENSE
+## file for more information.
+##
+## Copyright (C) 2011 Andrew D. Martin and Kevin M. Quinn
+##
+####################################################################
+##
+## Revisions:
+## - G. Vieilledent, May 9 2011 [initial file]
+##
+####################################################################
+
+
+MCMChlogit <- function (fixed, random, group, data, burnin=5000,
+ mcmc=10000, thin=10,
+ verbose=1, seed=NA, beta.start=NA, sigma2.start=NA,
+ Vb.start=NA, mubeta=0, Vbeta=1.0E6, r, R,
+ nu=0.001, delta=0.001, FixOD=0, ...)
+{
+
+ #========
+ # Basic checks
+ #========
+ check.group.hmodels(group, data)
+ check.mcmc.parameters.hmodels(burnin, mcmc, thin)
+ check.verbose.hmodels(verbose)
+ check.FixOD.hmodels(FixOD)
+ check.offset(list(...))
+
+ #========
+ # Seed
+ #========
+ seed <- form.seeds.hmodels(seed)
+
+ #========
+ # Form response and model matrices
+ #========
+ mf.fixed <- model.frame(formula=fixed,data=data)
+ X <- model.matrix(attr(mf.fixed,"terms"),data=mf.fixed)
+ Y <- model.response(mf.fixed)
+ check.Y.Binomial.hmodels(Y)
+ mf.random <- model.frame(formula=random,data=data)
+ W <- model.matrix(attr(mf.random,"terms"),data=mf.random)
+
+ #========
+ # Model parameters
+ #========
+ nobs <- nrow(X)
+ IdentGroup <- as.numeric(as.factor(as.character(data[,names(data)==as.character(group)])))-1
+ LevelsGroup <- sort(unique(IdentGroup+1))
+ LevelsGroup.Name <- sort(unique(as.character(data[,names(data)==as.character(group)])))
+ ngroup <- length(LevelsGroup)
+ np <- ncol(X)
+ nq <- ncol(W)
+ ngibbs <- mcmc+burnin
+ nthin <- thin
+ nburn <- burnin
+ nsamp <- mcmc/thin
+
+ #========
+ # Form and check starting parameters
+ #========
+ beta.start <- form.beta.start.hmodels(fixed,data,beta.start,np,family="binomial",defaults=NA)
+ sigma2.start <- form.sigma2.start.hmodels(fixed,data,sigma2.start,family="binomial")
+ Vb.start <- form.Vb.start.hmodels(Vb.start,nq)
+
+ #========
+ # Form priors
+ #========
+ mvn.prior <- form.mvn.prior.hmodels(mubeta,Vbeta,np)
+ mubeta <- mvn.prior[[1]]
+ Vbeta <- mvn.prior[[2]]
+ wishart.prior <- form.wishart.prior.hmodels(r,R,nq)
+ r <- wishart.prior[[1]]
+ R <- wishart.prior[[2]]
+ check.ig.prior.hmodels(nu,delta)
+ s1 <- nu
+ s2 <- delta
+
+ #========
+ # Parameters to save
+ #========
+ beta_vect <- rep(c(beta.start),each=nsamp)
+ Vb_vect <- rep(c(Vb.start),each=nsamp)
+ b_vect <- rep(0,nq*ngroup*nsamp)
+ V <- rep(sigma2.start,nsamp)
+ theta_pred <- rep(0.5,nobs)
+ Deviance <- rep(0,nsamp)
+
+ #========
+ # call C++ code to draw sample
+ #========
+ Sample <- .C("MCMChlogit",
+ #= Constants and data
+ ngibbs=as.integer(ngibbs), nthin=as.integer(nthin), nburn=as.integer(nburn),## Number of iterations, burning and samples
+ nobs=as.integer(nobs), ngroup=as.integer(ngroup), ## Constants
+ np=as.integer(np), nq=as.integer(nq), ## Constants
+ IdentGroup=as.integer(IdentGroup),
+ Y_vect=as.double(c(Y)), ## Response variable
+ X_vect=as.double(c(X)), ## Covariates
+ W_vect=as.double(c(W)), ## Covariates
+ #= Parameters to save
+ beta_vect.nonconst=as.double(beta_vect), ## Fixed parameters of the regression
+ b_vect.nonconst=as.double(b_vect), ## Random effects on intercept and slope
+ Vb_vect.nonconst=as.double(Vb_vect), ## Variance-covariance of random effects
+ V.nonconst=as.double(V), ## Variance of residuals
+ #= Defining priors
+ mubeta_vect=as.double(c(mubeta)), Vbeta_vect=as.double(c(Vbeta)),
+ r=as.double(r), R_vect=as.double(c(R)),
+ s1_V=as.double(s1), s2_V=as.double(s2),
+ #= Diagnostic
+ Deviance.nonconst=as.double(Deviance),
+ theta_pred.nonconst=as.double(theta_pred), ## Predictive posterior mean
+ #= Seeds
+ seed=as.integer(seed),
+ #= Verbose
+ verbose=as.integer(verbose),
+ #= Overdispersion
+ FixOD=as.integer(FixOD),
+ PACKAGE="MCMCpack")
+
+ #= Matrix of MCMC samples
+ Matrix <- matrix(NA,nrow=nsamp,ncol=np+nq*ngroup+nq*nq+2)
+ names.fixed <- paste("beta.",colnames(X),sep="")
+ names.random <- paste("b.",rep(colnames(W),each=ngroup),".",rep(LevelsGroup.Name,nq),sep="")
+ names.variances <- c(paste("VCV.",colnames(W),".",rep(colnames(W),each=nq),sep=""),"sigma2")
+ colnames(Matrix) <- c(names.fixed,names.random,names.variances,"Deviance")
+
+ #= Filling-in the matrix
+ Matrix[,c(1:np)] <- matrix(Sample[[12]],ncol=np)
+ Matrix[,c((np+1):(np+nq*ngroup))] <- matrix(Sample[[13]],ncol=nq*ngroup)
+ Matrix[,c((np+nq*ngroup+1):(np+nq*ngroup+nq*nq))] <- matrix(Sample[[14]],ncol=nq*nq)
+ Matrix[,ncol(Matrix)-1] <- Sample[[15]]
+ Matrix[,ncol(Matrix)] <- Sample[[22]]
+
+ #= Transform Sample list in an MCMC object
+ MCMC <- mcmc(Matrix,start=nburn+1,end=ngibbs,thin=nthin)
+
+ #= Output
+ return (list(mcmc=MCMC,theta.pred=Sample[[23]]))
+}
+
+####################################################################
+## END
+####################################################################
diff --git a/R/MCMChpoisson.R b/R/MCMChpoisson.R
new file mode 100644
index 0000000..796be7e
--- /dev/null
+++ b/R/MCMChpoisson.R
@@ -0,0 +1,181 @@
+##########################################################################
+## MCMChpoisson.R
+##
+## MCMChpoisson() samples from the posterior distribution of a
+## Poisson hierarchical linear regression model in R using linked C++
+## code in Scythe
+##
+## The code uses Algorithm 2 of Chib & Carlin (1999) for efficient
+## inference of (\beta | Y, sigma^2, Vb).
+##
+## Chib, S. & Carlin, B. P. (1999) On MCMC sampling in hierarchical
+## longitudinal models, Statistics and Computing, 9, 17-26
+##
+####################################################################
+##
+## Original code by Ghislain Vieilledent, may 2011
+## CIRAD UR B&SEF
+## ghislain.vieilledent at cirad.fr / ghislainv at gmail.com
+##
+####################################################################
+##
+## The initial version of this file was generated by the
+## auto.Scythe.call() function in the MCMCpack R package
+## written by:
+##
+## Andrew D. Martin
+## Dept. of Political Science
+## Washington University in St. Louis
+## admartin at wustl.edu
+##
+## Kevin M. Quinn
+## Dept. of Government
+## Harvard University
+## kevin_quinn at harvard.edu
+##
+## This software is distributed under the terms of the GNU GENERAL
+## PUBLIC LICENSE Version 2, June 1991. See the package LICENSE
+## file for more information.
+##
+## Copyright (C) 2011 Andrew D. Martin and Kevin M. Quinn
+##
+####################################################################
+##
+## Revisions:
+## - G. Vieilledent, May 9 2011 [initial file]
+##
+####################################################################
+
+
+MCMChpoisson <- function (fixed, random, group, data, burnin=5000,
+ mcmc=10000, thin=10,
+ verbose=1, seed=NA, beta.start=NA, sigma2.start=NA,
+ Vb.start=NA, mubeta=0, Vbeta=1.0E6, r, R,
+ nu=0.001, delta=0.001, FixOD=0, ...)
+{
+
+ #========
+ # Basic checks
+ #========
+ check.group.hmodels(group, data)
+ check.mcmc.parameters.hmodels(burnin, mcmc, thin)
+ check.verbose.hmodels(verbose)
+ check.FixOD.hmodels(FixOD)
+ check.offset(list(...))
+
+ #========
+ # Seed
+ #========
+ seed <- form.seeds.hmodels(seed)
+
+ #========
+ # Form response and model matrices
+ #========
+ mf.fixed <- model.frame(formula=fixed,data=data)
+ X <- model.matrix(attr(mf.fixed,"terms"),data=mf.fixed)
+ Y <- model.response(mf.fixed)
+ check.Y.Poisson.hmodels(Y)
+ mf.random <- model.frame(formula=random,data=data)
+ W <- model.matrix(attr(mf.random,"terms"),data=mf.random)
+
+ #========
+ # Model parameters
+ #========
+ nobs <- nrow(X)
+ IdentGroup <- as.numeric(as.factor(as.character(data[,names(data)==as.character(group)])))-1
+ LevelsGroup <- sort(unique(IdentGroup+1))
+ LevelsGroup.Name <- sort(unique(as.character(data[,names(data)==as.character(group)])))
+ ngroup <- length(LevelsGroup)
+ np <- ncol(X)
+ nq <- ncol(W)
+ ngibbs <- mcmc+burnin
+ nthin <- thin
+ nburn <- burnin
+ nsamp <- mcmc/thin
+
+ #========
+ # Form and check starting parameters
+ #========
+ beta.start <- form.beta.start.hmodels(fixed,data,beta.start,np,family="poisson",defaults=NA)
+ sigma2.start <- form.sigma2.start.hmodels(fixed,data,sigma2.start,family="poisson")
+ Vb.start <- form.Vb.start.hmodels(Vb.start,nq)
+
+ #========
+ # Form priors
+ #========
+ mvn.prior <- form.mvn.prior.hmodels(mubeta,Vbeta,np)
+ mubeta <- mvn.prior[[1]]
+ Vbeta <- mvn.prior[[2]]
+ wishart.prior <- form.wishart.prior.hmodels(r,R,nq)
+ r <- wishart.prior[[1]]
+ R <- wishart.prior[[2]]
+ check.ig.prior.hmodels(nu,delta)
+ s1 <- nu
+ s2 <- delta
+
+ #========
+ # Parameters to save
+ #========
+ beta_vect <- rep(c(beta.start),each=nsamp)
+ Vb_vect <- rep(c(Vb.start),each=nsamp)
+ b_vect <- rep(0,nq*ngroup*nsamp)
+ V <- rep(sigma2.start,nsamp)
+ lambda_pred <- rep(0.5,nobs)
+ Deviance <- rep(0,nsamp)
+
+ #========
+ # call C++ code to draw sample
+ #========
+ Sample <- .C("MCMChpoisson",
+ #= Constants and data
+ ngibbs=as.integer(ngibbs), nthin=as.integer(nthin), nburn=as.integer(nburn),## Number of iterations, burning and samples
+ nobs=as.integer(nobs), ngroup=as.integer(ngroup), ## Constants
+ np=as.integer(np), nq=as.integer(nq), ## Constants
+ IdentGroup=as.integer(IdentGroup),
+ Y_vect=as.double(c(Y)), ## Response variable
+ X_vect=as.double(c(X)), ## Covariates
+ W_vect=as.double(c(W)), ## Covariates
+ #= Parameters to save
+ beta_vect.nonconst=as.double(beta_vect), ## Fixed parameters of the regression
+ b_vect.nonconst=as.double(b_vect), ## Random effects on intercept and slope
+ Vb_vect.nonconst=as.double(Vb_vect), ## Variance-covariance of random effects
+ V.nonconst=as.double(V), ## Variance of residuals
+ #= Defining priors
+ mubeta_vect=as.double(c(mubeta)), Vbeta_vect=as.double(c(Vbeta)),
+ r=as.double(r), R_vect=as.double(c(R)),
+ s1_V=as.double(s1), s2_V=as.double(s2),
+ #= Diagnostic
+ Deviance.nonconst=as.double(Deviance),
+ lambda_pred.nonconst=as.double(lambda_pred), ## Predictive posterior mean
+ #= Seeds
+ seed=as.integer(seed),
+ #= Verbose
+ verbose=as.integer(verbose),
+ #= Overdispersion
+ FixOD=as.integer(FixOD),
+ PACKAGE="MCMCpack")
+
+ #= Matrix of MCMC samples
+ Matrix <- matrix(NA,nrow=nsamp,ncol=np+nq*ngroup+nq*nq+2)
+ names.fixed <- paste("beta.",colnames(X),sep="")
+ names.random <- paste("b.",rep(colnames(W),each=ngroup),".",rep(LevelsGroup.Name,nq),sep="")
+ names.variances <- c(paste("VCV.",colnames(W),".",rep(colnames(W),each=nq),sep=""),"sigma2")
+ colnames(Matrix) <- c(names.fixed,names.random,names.variances,"Deviance")
+
+ #= Filling-in the matrix
+ Matrix[,c(1:np)] <- matrix(Sample[[12]],ncol=np)
+ Matrix[,c((np+1):(np+nq*ngroup))] <- matrix(Sample[[13]],ncol=nq*ngroup)
+ Matrix[,c((np+nq*ngroup+1):(np+nq*ngroup+nq*nq))] <- matrix(Sample[[14]],ncol=nq*nq)
+ Matrix[,ncol(Matrix)-1] <- Sample[[15]]
+ Matrix[,ncol(Matrix)] <- Sample[[22]]
+
+ #= Transform Sample list in an MCMC object
+ MCMC <- mcmc(Matrix,start=nburn+1,end=ngibbs,thin=nthin)
+
+ #= Output
+ return (list(mcmc=MCMC,lambda.pred=Sample[[23]]))
+}
+
+####################################################################
+## END
+####################################################################
diff --git a/R/MCMChregress.R b/R/MCMChregress.R
new file mode 100644
index 0000000..ed7bd96
--- /dev/null
+++ b/R/MCMChregress.R
@@ -0,0 +1,177 @@
+##########################################################################
+## MCMChregress.R
+##
+## MCMChregress() samples from the posterior distribution of a
+## Gaussian hierarchical linear regression model in R using linked C++
+## code in Scythe
+##
+## The code uses Algorithm 2 of Chib & Carlin (1999) for efficient
+## inference of (\beta | Y, sigma^2, Vb).
+##
+## Chib, S. & Carlin, B. P. (1999) On MCMC sampling in hierarchical
+## longitudinal models, Statistics and Computing, 9, 17-26
+##
+####################################################################
+##
+## Original code by Ghislain Vieilledent, may 2011
+## CIRAD UR B&SEF
+## ghislain.vieilledent at cirad.fr / ghislainv at gmail.com
+##
+####################################################################
+##
+## The initial version of this file was generated by the
+## auto.Scythe.call() function in the MCMCpack R package
+## written by:
+##
+## Andrew D. Martin
+## Dept. of Political Science
+## Washington University in St. Louis
+## admartin at wustl.edu
+##
+## Kevin M. Quinn
+## Dept. of Government
+## Harvard University
+## kevin_quinn at harvard.edu
+##
+## This software is distributed under the terms of the GNU GENERAL
+## PUBLIC LICENSE Version 2, June 1991. See the package LICENSE
+## file for more information.
+##
+## Copyright (C) 2011 Andrew D. Martin and Kevin M. Quinn
+##
+####################################################################
+##
+## Revisions:
+## - G. Vieilledent, May 4 2011 [initial file]
+##
+####################################################################
+
+
+MCMChregress <- function (fixed, random, group, data, burnin=1000,
+ mcmc=10000, thin=10,
+ verbose=1, seed=NA, beta.start=NA, sigma2.start=NA,
+ Vb.start=NA, mubeta=0, Vbeta=1.0E6, r, R,
+ nu=0.001, delta=0.001, ...)
+{
+
+ #========
+ # Basic checks
+ #========
+ check.group.hmodels(group, data)
+ check.mcmc.parameters.hmodels(burnin, mcmc, thin)
+ check.verbose.hmodels(verbose)
+ check.offset(list(...))
+
+ #========
+ # Seed
+ #========
+ seed <- form.seeds.hmodels(seed)
+
+ #========
+ # Form response and model matrices
+ #========
+ mf.fixed <- model.frame(formula=fixed,data=data)
+ X <- model.matrix(attr(mf.fixed,"terms"),data=mf.fixed)
+ Y <- model.response(mf.fixed)
+ mf.random <- model.frame(formula=random,data=data)
+ W <- model.matrix(attr(mf.random,"terms"),data=mf.random)
+
+ #========
+ # Model parameters
+ #========
+ nobs <- nrow(X)
+ IdentGroup <- as.numeric(as.factor(as.character(data[,names(data)==as.character(group)])))-1
+ LevelsGroup <- sort(unique(IdentGroup+1))
+ LevelsGroup.Name <- sort(unique(as.character(data[,names(data)==as.character(group)])))
+ ngroup <- length(LevelsGroup)
+ np <- ncol(X)
+ nq <- ncol(W)
+ ngibbs <- mcmc+burnin
+ nthin <- thin
+ nburn <- burnin
+ nsamp <- mcmc/thin
+
+ #========
+ # Form and check starting parameters
+ #========
+ beta.start <- form.beta.start.hmodels(fixed,data,beta.start,np,family="gaussian",defaults=NA)
+ sigma2.start <- form.sigma2.start.hmodels(fixed,data,sigma2.start,family="gaussian")
+ Vb.start <- form.Vb.start.hmodels(Vb.start,nq)
+
+ #========
+ # Form priors
+ #========
+ mvn.prior <- form.mvn.prior.hmodels(mubeta,Vbeta,np)
+ mubeta <- mvn.prior[[1]]
+ Vbeta <- mvn.prior[[2]]
+ wishart.prior <- form.wishart.prior.hmodels(r,R,nq)
+ r <- wishart.prior[[1]]
+ R <- wishart.prior[[2]]
+ check.ig.prior.hmodels(nu,delta)
+ s1 <- nu
+ s2 <- delta
+
+ #========
+ # Parameters to save
+ #========
+ beta_vect <- rep(c(beta.start),each=nsamp)
+ Vb_vect <- rep(c(Vb.start),each=nsamp)
+ b_vect <- rep(0,nq*ngroup*nsamp)
+ V <- rep(sigma2.start,nsamp)
+ Y_pred <- rep(0,nobs)
+ Deviance <- rep(0,nsamp)
+
+ #========
+ # call C++ code to draw sample
+ #========
+ Sample <- .C("MCMChregress",
+ #= Constants and data
+ ngibbs=as.integer(ngibbs), nthin=as.integer(nthin), nburn=as.integer(nburn),## Number of iterations, burning and samples
+ nobs=as.integer(nobs), ngroup=as.integer(ngroup), ## Constants
+ np=as.integer(np), nq=as.integer(nq), ## Constants
+ IdentGroup=as.integer(IdentGroup),
+ Y_vect=as.double(c(Y)), ## Response variable
+ X_vect=as.double(c(X)), ## Covariates
+ W_vect=as.double(c(W)), ## Covariates
+ #= Parameters to save
+ beta_vect.nonconst=as.double(beta_vect), ## Fixed parameters of the regression
+ b_vect.nonconst=as.double(b_vect), ## Random effects on intercept and slope
+ Vb_vect.nonconst=as.double(Vb_vect), ## Variance-covariance of random effects
+ V.nonconst=as.double(V), ## Variance of residuals
+ #= Defining priors
+ mubeta_vect=as.double(c(mubeta)), Vbeta_vect=as.double(c(Vbeta)),
+ r=as.double(r), R_vect=as.double(c(R)),
+ s1_V=as.double(s1), s2_V=as.double(s2),
+ #= Diagnostic
+ Deviance.nonconst=as.double(Deviance),
+ Y_pred.nonconst=as.double(Y_pred), ## Predictive posterior mean
+ #= Seeds
+ seed=as.integer(seed),
+ #= Verbose
+ verbose=as.integer(verbose),
+ PACKAGE="MCMCpack")
+
+ #= Matrix of MCMC samples
+ Matrix <- matrix(NA,nrow=nsamp,ncol=np+nq*ngroup+nq*nq+2)
+ names.fixed <- paste("beta.",colnames(X),sep="")
+ names.random <- paste("b.",rep(colnames(W),each=ngroup),".",rep(LevelsGroup.Name,nq),sep="")
+ names.variances <- c(paste("VCV.",colnames(W),".",rep(colnames(W),each=nq),sep=""),"sigma2")
+ colnames(Matrix) <- c(names.fixed,names.random,names.variances,"Deviance")
+
+ #= Filling-in the matrix
+ Matrix[,c(1:np)] <- matrix(Sample[[12]],ncol=np)
+ Matrix[,c((np+1):(np+nq*ngroup))] <- matrix(Sample[[13]],ncol=nq*ngroup)
+ Matrix[,c((np+nq*ngroup+1):(np+nq*ngroup+nq*nq))] <- matrix(Sample[[14]],ncol=nq*nq)
+ Matrix[,ncol(Matrix)-1] <- Sample[[15]]
+ Matrix[,ncol(Matrix)] <- Sample[[22]]
+
+ #= Transform Sample list in an MCMC object
+ MCMC <- mcmc(Matrix,start=nburn+1,end=ngibbs,thin=nthin)
+
+ #= Output
+ return (list(mcmc=MCMC,Y.pred=Sample[[23]]))
+}
+
+####################################################################
+## END
+####################################################################
diff --git a/R/hidden-hmodels.R b/R/hidden-hmodels.R
new file mode 100644
index 0000000..e1aa929
--- /dev/null
+++ b/R/hidden-hmodels.R
@@ -0,0 +1,348 @@
+#########################################################################
+##
+## Following functions perform checks for MCMCpack hierarchical models
+## (denoted MCMCh...).
+## ghislain.vieilledent at cirad.fr, May 5 2011
+##
+#########################################################################
+
+##=======================================================================
+##
+## Check group
+##
+##=======================================================================
+
+check.group.hmodels <- function(group,data) {
+
+ if (!(group %in% colnames(data))) {
+ cat("Error: group must be a string which gives the name of the grouping variable in data \n")
+ stop("Please respecify and call ", calling.function(), " again.",
+ call.=FALSE)
+ }
+ return(0)
+}
+
+##=======================================================================
+##
+## Check mcmc parameters
+##
+##=======================================================================
+
+check.mcmc.parameters.hmodels <- function(burnin, mcmc, thin) {
+
+ if(mcmc %% thin != 0) {
+ cat("Error: MCMC iterations not evenly divisible by thinning interval.\n")
+ stop("Please respecify and call ", calling.function(), " again.",
+ call.=FALSE)
+ }
+ if(mcmc <= 0) {
+ cat("Error: MCMC iterations must be strictly positive.\n")
+ stop("Please respecify and call ", calling.function(), " again.",
+ call.=FALSE)
+ }
+ if(burnin < 0) {
+ cat("Error: Burnin iterations negative.\n")
+ stop("Please respecify and call ", calling.function(), " again.",
+ call.=FALSE)
+ }
+ if(((burnin+mcmc) %% 10 != 0) || (burnin+mcmc)<100) {
+ cat("Error: Value 'burnin+mcmc' should be divisible by 10 and >= 100.\n")
+ stop("Please respecify and call ", calling.function(), " again.",
+ call.=FALSE)
+ }
+ if(thin < 1) {
+ cat("Error: Thinning interval must be superior or equal to 1.\n")
+ stop("Please respecify and call ", calling.function(), " again.",
+ call.=FALSE)
+ }
+ return(0)
+}
+
+##=======================================================================
+##
+## Check verbose
+##
+##=======================================================================
+
+check.verbose.hmodels <- function (verbose) {
+ if (!(verbose%in%c(0,1))) {
+ cat("Error: verbose must take value 0 or 1.\n")
+ stop("Please respecify and call ", calling.function(), " again.",
+ call.=FALSE)
+ }
+ return(0)
+}
+
+##=======================================================================
+##
+## Check FixOD
+##
+##=======================================================================
+
+check.FixOD.hmodels <- function (FixOD) {
+ if (!(FixOD%in%c(0,1))) {
+ cat("Error: FixOD must take value 0 or 1.\n")
+ stop("Please respecify and call ", calling.function(), " again.",
+ call.=FALSE)
+ }
+ return(0)
+}
+
+##=======================================================================
+##
+## Check Y in c(0,1) for Binomial process
+##
+##=======================================================================
+
+check.Y.Binomial.hmodels <- function (Y) {
+ if (sum(!(c(Y)%in%c(0,1)))>0) {
+ cat("Error: Response variable must take value 0 or 1.\n")
+ stop("Please respecify and call ", calling.function(), " again.",
+ call.=FALSE)
+ }
+ return(0)
+}
+
+##=======================================================================
+##
+## Check Y is a positive integer for Poisson process
+##
+##=======================================================================
+
+check.Y.Poisson.hmodels <- function (Y) {
+ if (sum(!(c(Y)>=0 && c(Y)%%1==0)>0)) {
+ cat("Error: Response variable must be a positive integer.\n")
+ stop("Please respecify and call ", calling.function(), " again.",
+ call.=FALSE)
+ }
+ return(0)
+}
+
+##=======================================================================
+##
+## Check and form seed
+##
+##=======================================================================
+
+form.seeds.hmodels <- function(seed) {
+ if(length(seed)!=1) {
+ cat("Error: Mersenne seed should be of length 1.\n")
+ stop("Please respecify and call ", calling.function(), " again.", call.=FALSE)
+ }
+ if(length(seed)==1) {
+ if(is.na(seed)) {
+ seed <- 12345
+ }
+ else {
+ seed <- as.integer(seed)
+ }
+ if(seed < 0) {
+ cat("Error: Mersenne seed negative.\n")
+ stop("Please respecify and call ", calling.function(), " again.", call.=FALSE)
+ }
+ }
+ return(seed)
+}
+
+##=======================================================================
+##
+## Check and form starting values
+##
+##=======================================================================
+
+#===============================
+# form beta.start
+#===============================
+
+form.beta.start.hmodels <- function (fixed,data,beta.start,np,family,defaults=NA) {
+
+ if (is.na(beta.start)[1] & is.na(defaults)[1]){ # use GLM estimates
+ beta.start <- matrix(coef(glm(fixed, family=family, data=data)), np, 1)
+ }
+ else if(is.na(beta.start)[1] & !is.na(defaults)[1]){ # use passed values
+ beta.start <- matrix(defaults,np,1)
+ }
+ else if(is.null(dim(beta.start))) {
+ beta.start <- beta.start * matrix(1,np,1)
+ }
+ else if(!all(dim(beta.start) == c(np,1))) {
+ cat("Error: beta.start not conformable.\n")
+ stop("Please respecify and call ", calling.function(), " again.\n",
+ call.=FALSE)
+ }
+ return(beta.start)
+}
+
+#===============================
+# form sigma2.start
+#===============================
+
+form.sigma2.start.hmodels <- function (fixed,data,sigma2.start,family) {
+
+ if (is.na(sigma2.start)[1]){ # use GLM estimates
+ sigma2.start <- var(residuals(glm(fixed, family=family, data=data)))
+ }
+ else {
+ sigma2.start <- as.integer(sigma2.start[1])
+ }
+ if (sigma2.start<=0) {
+ cat("Error: Starting value for sigma2 negative.\n")
+ stop("Please respecify and call ", calling.function(), " again.",
+ call.=FALSE)
+ }
+ return(sigma2.start)
+}
+
+#===============================
+# form Vb.start
+#===============================
+
+form.Vb.start.hmodels <- function (Vb.start,nq) {
+
+ if (is.na(Vb.start)[1]) {
+ Vb.start <- diag(1,nq)
+ }
+ else if(is.null(dim(Vb.start))) {
+ Vb.start <- Vb.start * diag(nq)
+ }
+ if ((dim(Vb.start)[1] != nq) || (dim(Vb.start)[2] != nq)) {
+ cat("Error: Vb.start not conformable.\n")
+ stop("Please respecify and call ", calling.function(), " again.",
+ call.=FALSE)
+ }
+ if (sum(diag(Vb.start)>0)!=nq) {
+ cat("Error: Vb.start should have positive values on the diagonal.\n")
+ stop("Please respecify and call ", calling.function(), " again.",
+ call.=FALSE)
+ }
+ return(Vb.start)
+}
+
+##=======================================================================
+##
+## Check and form priors
+##
+##=======================================================================
+
+#===============================
+# form multivariate Normal prior
+#===============================
+
+form.mvn.prior.hmodels <- function(mubeta, Vbeta, np) {
+
+ # prior mean
+ if(is.null(dim(mubeta))) {
+ mubeta <- mubeta * matrix(1,nrow=np,ncol=1)
+ }
+ if((dim(mubeta)[1] != np) || (dim(mubeta)[2] != 1)) {
+ cat("Error: in N(mubeta,Vbeta) prior, mubeta not conformable.\n")
+ stop("Please respecify and call ", calling.function(), " again.",
+ call.=FALSE)
+ }
+
+ # prior variance
+ if(is.null(dim(Vbeta))) {
+ if (length(Vbeta) > np){
+ cat("Vbeta was passed as a vector longer than np.
+ Vbeta must be either a scalar or a matrix.")
+ stop("Please respecify and call ", calling.function(), " again.\n", call.=FALSE)
+ }
+ Vbeta <- Vbeta * diag(np)
+ }
+ if((dim(Vbeta)[1] != np) || (dim(Vbeta)[2] != np)) {
+ cat("Error: in N(mubeta,Vbeta), prior Vbeta not conformable [p times p].\n")
+ stop("Please respecify and call ", calling.function(), " again.\n",
+ call.=FALSE)
+ }
+
+ # check Vbeta for symmetry
+ symproblem <- FALSE
+ for (i in 1:np){
+ for (j in i:np){
+ if (Vbeta[i,j] != Vbeta[j,i]){
+ symproblem <- TRUE
+ }
+ }
+ }
+ if (symproblem){
+ cat("Vbeta is not symmetric.\n")
+ stop("Please respecify and call ", calling.function(), " again.\n",
+ call.=FALSE)
+ }
+
+ # check Vbeta for positive number on the diagonal
+ if (sum(diag(Vbeta)>0)!=np) {
+ cat("Error: Vbeta should have positive values on the diagonal.\n")
+ stop("Please respecify and call ", calling.function(), " again.\n",
+ call.=FALSE)
+ }
+ return(list(mubeta,Vbeta))
+}
+
+#===================
+# form Wishart prior
+#===================
+
+form.wishart.prior.hmodels <- function(r, R, nq) {
+
+ # check to see if degrees of freedom produces proper prior
+ if(r < nq) {
+ cat("Error: in Wishart(r,rR) prior, r less than q.\n")
+ stop("Please respecify and call ", calling.function(), " again.\n")
+ }
+
+ # form the prior scale matrix
+ if(is.null(dim(R))) {
+ R <- R * diag(nq)
+ }
+ if((dim(R)[1] != nq) | (dim(R)[2] != nq)) {
+ cat("Error: in Wishart(r,rR) prior, R not comformable [q times q].\n")
+ stop("Please respecify and call ", calling.function(), " again.\n")
+ }
+
+ # check R for symmetry
+ symproblem <- FALSE
+ for (i in 1:nq){
+ for (j in i:nq){
+ if (R[i,j] != R[j,i]){
+ symproblem <- TRUE
+ }
+ }
+ }
+ if (symproblem){
+ cat("R is not symmetric.\n")
+ stop("Please respecify and call ", calling.function(), " again.\n",
+ call.=FALSE)
+ }
+
+ # check R for positive number on the diagonal
+ if (sum(diag(R)>0)!=nq) {
+ cat("Error: R should have positive values on the diagonal.\n")
+ stop("Please respecify and call ", calling.function(), " again.\n",
+ call.=FALSE)
+ }
+ return(list(r,R))
+}
+
+#==========================
+# check inverse Gamma prior
+#==========================
+
+check.ig.prior.hmodels <- function(nu, delta) {
+
+ if(nu <= 0) {
+ cat("Error: in IG(nu,delta) prior, nu less than or equal to zero.\n")
+ stop("Please respecify and call ", calling.function(), " again.\n",
+ call.=FALSE)
+ }
+ if(delta <= 0) {
+ cat("Error: in IG(nu,delta) prior, delta less than or equal to zero.\n")
+ stop("Please respecify and call ", calling.function(), " again.\n",
+ call.=FALSE)
+ }
+ return(0)
+}
+
+#==========================
+# END
+#==========================
diff --git a/data/PErisk.rda b/data/PErisk.rda
index f272314..d60ae32 100644
Binary files a/data/PErisk.rda and b/data/PErisk.rda differ
diff --git a/data/SupremeCourt.rda b/data/SupremeCourt.rda
index 7866968..3b07c6c 100644
Binary files a/data/SupremeCourt.rda and b/data/SupremeCourt.rda differ
diff --git a/man/MCMCbinaryChange.Rd b/man/MCMCbinaryChange.Rd
old mode 100644
new mode 100755
diff --git a/man/MCMChlogit.Rd b/man/MCMChlogit.Rd
new file mode 100644
index 0000000..5265153
--- /dev/null
+++ b/man/MCMChlogit.Rd
@@ -0,0 +1,292 @@
+\name{MCMChlogit}
+
+\alias{MCMChlogit}
+
+\title{Markov Chain Monte Carlo for the Hierarchical Binomial Linear
+ Regression Model using the logit link function}
+
+\description{MCMChlogit generates a sample from the posterior
+ distribution of a Hierarchical Binomial Linear Regression Model using
+ the logit link function and Algorithm 2 of Chib and Carlin
+ (1999). This model uses a multivariate Normal prior for the fixed
+ effects parameters, an Inverse-Wishart prior on the random effects
+ variance matrix, and an Inverse-Gamma prior on the variance modelling
+ over-dispersion. The user supplies data and priors, and a sample from
+ the posterior distribution is returned as an mcmc object, which can be
+ subsequently analyzed with functions provided in the coda package.}
+
+\usage{MCMChlogit(fixed, random, group, data, burnin=5000,
+mcmc=10000, thin=10, verbose=1, seed=NA, beta.start=NA,
+sigma2.start=NA, Vb.start=NA, mubeta=0, Vbeta=1.0E6, r, R,
+nu=0.001, delta=0.001, FixOD=0, ...)}
+
+\arguments{
+
+ \item{fixed}{A two-sided linear formula of the form 'y~x1+...+xp'
+ describing the fixed-effects part of the model, with the response on
+ the left of a '~' operator and the p fixed terms, separated by '+'
+ operators, on the right. Response variable y must be 0 or 1
+ (Binomial process).}
+
+ \item{random}{A one-sided formula of the form '~x1+...+xq' specifying
+ the model for the random effects part of the model, with the q
+ random terms, separated by '+' operators.}
+
+ \item{group}{String indicating the name of the grouping variable
+ in \code{data}, defining the hierarchical structure of the model.}
+
+ \item{data}{A data frame containing the variables in the model.}
+
+ \item{burnin}{The number of burnin iterations for the sampler.}
+
+ \item{mcmc}{The number of Gibbs iterations for the sampler. Total
+ number of Gibbs iterations is equal to
+ \code{burnin+mcmc}. \code{burnin+mcmc} must be divisible by 10 and
+ superior or equal to 100 so that the progress bar can be displayed.}
+
+ \item{thin}{The thinning interval used in the simulation. The number of
+ mcmc iterations must be divisible by this value.}
+
+ \item{seed}{The seed for the random number generator. If NA, the
+ Mersenne Twister generator is used with default seed 12345; if an
+ integer is passed it is used to seed the Mersenne twister.}
+
+ \item{verbose}{A switch (0,1) which determines whether or not the
+ progress of the sampler is printed to the screen. Default is 1: a
+ progress bar is printed, indicating the step (in \%) reached by the
+ Gibbs sampler.}
+
+ \item{beta.start}{The starting values for the \eqn{\beta}{beta}
+ vector. This can either be a scalar or a p-length vector. The
+ default value of NA will use the OLS \eqn{\beta}{beta} estimate of
+ the corresponding Gaussian Linear Regression without random
+ effects. If this is a scalar, that value will serve as the starting
+ value mean for all of the betas.}
+
+ \item{sigma2.start}{Scalar for the starting value of the residual
+ error variance. The default value of NA will use the OLS estimates
+ of the corresponding Gaussian Linear Regression without random
+ effects.}
+
+ \item{Vb.start}{The starting value for variance matrix of the random
+ effects. This must be a square q-dimension matrix. Default value of NA
+ uses an identity matrix.}
+
+ \item{mubeta}{The prior mean of \eqn{\beta}{beta}. This can either be
+ a scalar or a p-length vector. If this takes a scalar value, then
+ that value will serve as the prior mean for all of the betas. The
+ default value of 0 will use a vector of zeros for an uninformative
+ prior.}
+
+ \item{Vbeta}{The prior variance of \eqn{\beta}{beta}. This can either
+ be a scalar or a square p-dimension matrix. If this takes a scalar value,
+ then that value times an identity matrix serves as the prior
+ variance of beta. Default value of 1.0E6 will use a diagonal matrix
+ with very large variance for an uninformative flat prior.}
+
+ \item{r}{The shape parameter for the Inverse-Wishart prior on variance
+ matrix for the random effects. r must be superior or equal to
+ q. Set r=q for an uninformative prior. See the NOTE for more details}
+
+ \item{R}{The scale matrix for the Inverse-Wishart prior on variance matrix
+ for the random effects. This must be a square q-dimension
+ matrix. Use plausible variance regarding random effects for the
+ diagonal of R. See the NOTE for more details}
+
+ \item{nu}{The shape parameter for the Inverse-Gamma prior on the
+ residual error variance. Default value is \code{nu=delta=0.001} for
+ uninformative prior.}
+
+ \item{delta}{The rate (1/scale) parameter for the Inverse-Gamma prior on the
+ residual error variance. Default value is \code{nu=delta=0.001} for
+ uninformative prior.}
+
+ \item{FixOD}{A switch (0,1) which determines whether or not the
+ variance for over-dispersion (sigma2) should be fixed (1) or not
+ (0). Default is 0, parameter sigma2 is estimated. If FixOD=1, sigma2
+ is fixed to the value provided for \code{sigma2.start}.}
+
+ \item{...}{further arguments to be passed}
+
+}
+
+\details{
+
+ \code{MCMChlogit} simulates from the posterior distribution sample
+ using the blocked Gibbs sampler of Chib and Carlin (1999), Algorithm
+ 2. The simulation is done in compiled C++ code to maximize
+ efficiency. Please consult the coda documentation for a comprehensive
+ list of functions that can be used to analyze the posterior sample.
+
+ The model takes the following form:
+
+ \deqn{y_i \sim \mathcal{B}ernoulli(\theta_i)}{y_i ~ Bernoulli(theta_i)}
+
+ With latent variables \eqn{\phi(\theta_i)}{phi(theta)},
+ \eqn{\phi}{phi} being the logit link function:
+ \deqn{\phi(\theta_i) = X_i \beta + W_i b_i + \varepsilon_i}{phi(theta_i) = X_i * beta + W_i *
+ b_i + epsilon_i}
+
+ Where each group \eqn{i} have \eqn{k_i} observations.
+
+ Where the random effects:
+ \deqn{b_i \sim \mathcal{N}_q(0,V_b)}{b_i ~ N_q(0,Vb)}
+ And the over-dispersion terms:
+ \deqn{\varepsilon_i \sim \mathcal{N}(0, \sigma^2 I_{k_i})}{epsilon_i ~ N(0,
+ sigma^2 I_{k_i})}
+ We assume standard, conjugate priors:
+ \deqn{\beta \sim \mathcal{N}_p(\mu_{\beta},V_{\beta})}{beta ~ N_p(mubeta,Vbeta)}
+ And:
+ \deqn{\sigma^{2} \sim \mathcal{IG}amma(\nu, 1/\delta)}{sigma^2 ~
+ IGamma(nu, 1/delta)}
+ And:
+ \deqn{V_b \sim \mathcal{IW}ishart(r, rR)}{Vb ~ IWishart(r,
+ rR)}
+ See Chib and Carlin (1999) for more details.
+
+ \emph{NOTE:} We do not provide default parameters for the priors on
+ the precision matrix for the random effects. When fitting one of
+ these models, it is of utmost importance to choose a prior that
+ reflects your prior beliefs about the random effects. Using the
+ \code{dwish} and \code{rwish} functions might be useful in choosing
+ these values.
+}
+
+\value{
+
+ \item{mcmc}{An mcmc object that contains the posterior sample. This
+ object can be summarized by functions provided by the coda
+ package. The posterior sample of the deviance \eqn{D}{D}, with
+ \eqn{D=-2\log(\prod_i P(y_i|\theta_i))}{%
+ D=-2log(prod_i P(y_i|theta_i))}, is also provided.}
+
+ \item{theta.pred}{Predictive posterior mean for the inverse-logit of
+ the latent variables. The approximation of Diggle et
+ al. (2004) is used to marginalized with respect to over-dispersion terms:
+ \deqn{E[\theta_i|\beta,b_i,\sigma^2]=\phi^{-1}((X_i
+ \beta+W_i b_i)/\sqrt{(16\sqrt{3}/15\pi)^2\sigma^2+1})}{%
+ E[theta_i|beta,b_i,sigma^2]=phi^(-1)((X_i*beta+W_i*b_i) /%
+ sqrt((16*sqrt(3)/15*pi)^2*sigma^2+1))}
+ }
+}
+
+\references{
+ Siddhartha Chib and Bradley P. Carlin. 1999. ``On MCMC Sampling in
+ Hierarchical Longitudinal Models.'' \emph{Statistics and Computing.} 9:
+ 17-26.
+
+ Daniel Pemstein, Kevin M. Quinn, and Andrew D. Martin. 2007.
+ \emph{Scythe Statistical Library 1.0.} \url{http://scythe.wustl.edu}.
+
+ Andrew D. Martin and Kyle L. Saunders. 2002. ``Bayesian Inference for
+ Political Science Panel Data.'' Paper presented at the 2002 Annual Meeting
+ of the American Political Science Association.
+
+ Martyn Plummer, Nicky Best, Kate Cowles, and Karen Vines. 2002.
+ \emph{Output Analysis and Diagnostics for MCMC (CODA)}.
+ \url{http://www-fis.iarc.fr/coda/}.
+
+ Diggle P., Heagerty P., Liang K., and Zeger S. 2004. ``Analysis of Longitudinal Data.''
+ \emph{Oxford University Press}, 2sd Edition.
+
+}
+
+\author{
+ Ghislain Vieilledent <ghislain.vieilledent at cirad.fr>
+}
+
+\seealso{
+ \code{\link[coda]{plot.mcmc}}, \code{\link[coda]{summary.mcmc}}
+}
+
+\examples{
+\dontrun{
+#========================================
+# Hierarchical Binomial Linear Regression
+#========================================
+
+#== inv.logit function
+inv.logit <- function(x, min=0, max=1) {
+ p <- exp(x)/(1+exp(x))
+ p <- ifelse( is.na(p) & !is.na(x), 1, p ) # fix problems with +Inf
+ return(p*(max-min)+min)
+}
+
+#== Generating data
+
+# Constants
+nobs <- 1000
+nspecies <- 20
+species <- c(1:nspecies,sample(c(1:nspecies),(nobs-nspecies),replace=TRUE))
+
+# Covariates
+X1 <- runif(n=nobs,min=-10,max=10)
+X2 <- runif(n=nobs,min=-10,max=10)
+X <- cbind(rep(1,nobs),X1,X2)
+W <- X
+
+# Target parameters
+# beta
+beta.target <- matrix(c(0.3,0.2,0.1),ncol=1)
+# Vb
+Vb.target <- c(0.5,0.05,0.05)
+# b
+b.target <- cbind(rnorm(nspecies,mean=0,sd=sqrt(Vb.target[1])),
+ rnorm(nspecies,mean=0,sd=sqrt(Vb.target[2])),
+ rnorm(nspecies,mean=0,sd=sqrt(Vb.target[3])))
+
+# Response
+theta <- vector()
+Y <- vector()
+for (n in 1:nobs) {
+ theta[n] <- inv.logit(X[n,]\%*\%beta.target+W[n,]\%*\%b.target[species[n],])
+ Y[n] <- rbinom(n=1,size=1,prob=theta[n])
+}
+
+# Data-set
+Data <- as.data.frame(cbind(Y,theta,X1,X2,species))
+plot(Data$X1,Data$theta)
+
+#== Call to MCMChlogit
+model <- MCMChlogit(fixed=Y~X1+X2, random=~X1+X2, group="species",
+ data=Data, burnin=5000, mcmc=1000, thin=1,verbose=1,
+ seed=NA, beta.start=0, sigma2.start=1,
+ Vb.start=1, mubeta=0, Vbeta=1.0E6,
+ r=3, R=diag(c(1,0.1,0.1)), nu=0.001, delta=0.001, FixOD=1)
+
+#== MCMC analysis
+
+# Graphics
+pdf("Posteriors-MCMChlogit.pdf")
+plot(model$mcmc)
+dev.off()
+
+# Summary
+summary(model$mcmc)
+
+# Predictive posterior mean for each observation
+model$theta.pred
+
+# Predicted-Observed
+plot(Data$theta,model$theta.pred)
+abline(a=0,b=1)
+
+## #Not run
+## #You can also compare with lme4 results
+## #== lme4 resolution
+## library(lme4)
+## model.lme4 <- lmer(Y~X1+X2+(1+X1+X2|species),data=Data,family="binomial")
+## summary(model.lme4)
+## plot(fitted(model.lme4),model$theta.pred,main="MCMChlogit/lme4")
+## abline(a=0,b=1)
+}
+}
+
+\keyword{models}
+\keyword{hierarchical models}
+\keyword{mixed models}
+\keyword{glmm}
+\keyword{logit}
+\keyword{Binomial}
+\keyword{MCMC}
+\keyword{bayesian}
diff --git a/man/MCMChpoisson.Rd b/man/MCMChpoisson.Rd
new file mode 100644
index 0000000..cc63275
--- /dev/null
+++ b/man/MCMChpoisson.Rd
@@ -0,0 +1,279 @@
+\name{MCMChpoisson}
+
+\alias{MCMChpoisson}
+
+\title{Markov Chain Monte Carlo for the Hierarchical Poisson Linear
+ Regression Model using the log link function}
+
+\description{MCMChpoisson generates a sample from the posterior
+ distribution of a Hierarchical Poisson Linear Regression Model using
+ the log link function and Algorithm 2 of Chib and Carlin
+ (1999). This model uses a multivariate Normal prior for the fixed
+ effects parameters, an Inverse-Wishart prior on the random effects
+ variance matrix, and an Inverse-Gamma prior on the variance modelling
+ over-dispersion. The user supplies data and priors, and a sample from
+ the posterior distribution is returned as an mcmc object, which can be
+ subsequently analyzed with functions provided in the coda package.}
+
+\usage{MCMChpoisson(fixed, random, group, data, burnin=5000,
+mcmc=10000, thin=10, verbose=1, seed=NA, beta.start=NA,
+sigma2.start=NA, Vb.start=NA, mubeta=0, Vbeta=1.0E6, r, R,
+nu=0.001, delta=0.001, FixOD=0, ...)}
+
+\arguments{
+
+ \item{fixed}{A two-sided linear formula of the form 'y~x1+...+xp'
+ describing the fixed-effects part of the model, with the response on
+ the left of a '~' operator and the p fixed terms, separated by '+'
+ operators, on the right. Response variable y must be 0 or 1
+ (Binomial process).}
+
+ \item{random}{A one-sided formula of the form '~x1+...+xq' specifying
+ the model for the random effects part of the model, with the q
+ random terms, separated by '+' operators.}
+
+ \item{group}{String indicating the name of the grouping variable
+ in \code{data}, defining the hierarchical structure of the model.}
+
+ \item{data}{A data frame containing the variables in the model.}
+
+ \item{burnin}{The number of burnin iterations for the sampler.}
+
+ \item{mcmc}{The number of Gibbs iterations for the sampler. Total
+ number of Gibbs iterations is equal to
+ \code{burnin+mcmc}. \code{burnin+mcmc} must be divisible by 10 and
+ superior or equal to 100 so that the progress bar can be displayed.}
+
+ \item{thin}{The thinning interval used in the simulation. The number of
+ mcmc iterations must be divisible by this value.}
+
+ \item{seed}{The seed for the random number generator. If NA, the
+ Mersenne Twister generator is used with default seed 12345; if an
+ integer is passed it is used to seed the Mersenne twister.}
+
+ \item{verbose}{A switch (0,1) which determines whether or not the
+ progress of the sampler is printed to the screen. Default is 1: a
+ progress bar is printed, indicating the step (in \%) reached by the
+ Gibbs sampler.}
+
+ \item{beta.start}{The starting values for the \eqn{\beta}{beta}
+ vector. This can either be a scalar or a p-length vector. The
+ default value of NA will use the OLS \eqn{\beta}{beta} estimate of
+ the corresponding Gaussian Linear Regression without random
+ effects. If this is a scalar, that value will serve as the starting
+ value mean for all of the betas.}
+
+ \item{sigma2.start}{Scalar for the starting value of the residual
+ error variance. The default value of NA will use the OLS estimates
+ of the corresponding Gaussian Linear Regression without random
+ effects.}
+
+ \item{Vb.start}{The starting value for variance matrix of the random
+ effects. This must be a square q-dimension matrix. Default value of NA
+ uses an identity matrix.}
+
+ \item{mubeta}{The prior mean of \eqn{\beta}{beta}. This can either be
+ a scalar or a p-length vector. If this takes a scalar value, then
+ that value will serve as the prior mean for all of the betas. The
+ default value of 0 will use a vector of zeros for an uninformative
+ prior.}
+
+ \item{Vbeta}{The prior variance of \eqn{\beta}{beta}. This can either
+ be a scalar or a square p-dimension matrix. If this takes a scalar value,
+ then that value times an identity matrix serves as the prior
+ variance of beta. Default value of 1.0E6 will use a diagonal matrix
+ with very large variance for an uninformative flat prior.}
+
+ \item{r}{The shape parameter for the Inverse-Wishart prior on variance
+ matrix for the random effects. r must be superior or equal to
+ q. Set r=q for an uninformative prior. See the NOTE for more details}
+
+ \item{R}{The scale matrix for the Inverse-Wishart prior on variance matrix
+ for the random effects. This must be a square q-dimension
+ matrix. Use plausible variance regarding random effects for the
+ diagonal of R. See the NOTE for more details}
+
+ \item{nu}{The shape parameter for the Inverse-Gamma prior on the
+ residual error variance. Default value is \code{nu=delta=0.001} for
+ uninformative prior.}
+
+ \item{delta}{The rate (1/scale) parameter for the Inverse-Gamma prior on the
+ residual error variance. Default value is \code{nu=delta=0.001} for
+ uninformative prior.}
+
+ \item{FixOD}{A switch (0,1) which determines whether or not the
+ variance for over-dispersion (sigma2) should be fixed (1) or not
+ (0). Default is 0, parameter sigma2 is estimated. If FixOD=1, sigma2
+ is fixed to the value provided for \code{sigma2.start}.}
+
+ \item{...}{further arguments to be passed}
+
+}
+
+\details{
+
+ \code{MCMChpoisson} simulates from the posterior distribution sample
+ using the blocked Gibbs sampler of Chib and Carlin (1999), Algorithm
+ 2. The simulation is done in compiled C++ code to maximize
+ efficiency. Please consult the coda documentation for a comprehensive
+ list of functions that can be used to analyze the posterior sample.
+
+ The model takes the following form:
+
+ \deqn{y_i \sim \mathcal{P}oisson(\lambda_i)}{y_i ~ Poisson(lambda_i)}
+
+ With latent variables \eqn{\phi(\lambda_i)}{phi(lambda)},
+ \eqn{\phi}{phi} being the log link function:
+ \deqn{\phi(\lambda_i) = X_i \beta + W_i b_i + \varepsilon_i}{phi(lambda_i) = X_i * beta + W_i *
+ b_i + epsilon_i}
+
+ Where each group \eqn{i} have \eqn{k_i} observations.
+
+ Where the random effects:
+ \deqn{b_i \sim \mathcal{N}_q(0,V_b)}{b_i ~ N_q(0,Vb)}
+ And the over-dispersion terms:
+ \deqn{\varepsilon_i \sim \mathcal{N}(0, \sigma^2 I_{k_i})}{epsilon_i ~ N(0,
+ sigma^2 I_{k_i})}
+ We assume standard, conjugate priors:
+ \deqn{\beta \sim \mathcal{N}_p(\mu_{\beta},V_{\beta})}{beta ~ N_p(mubeta,Vbeta)}
+ And:
+ \deqn{\sigma^{2} \sim \mathcal{IG}amma(\nu, 1/\delta)}{sigma^2 ~
+ IGamma(nu, 1/delta)}
+ And:
+ \deqn{V_b \sim \mathcal{IW}ishart(r, rR)}{Vb ~ IWishart(r,
+ rR)}
+ See Chib and Carlin (1999) for more details.
+
+ \emph{NOTE:} We do not provide default parameters for the priors on
+ the precision matrix for the random effects. When fitting one of
+ these models, it is of utmost importance to choose a prior that
+ reflects your prior beliefs about the random effects. Using the
+ \code{dwish} and \code{rwish} functions might be useful in choosing
+ these values.
+}
+
+\value{
+
+ \item{mcmc}{An mcmc object that contains the posterior sample. This
+ object can be summarized by functions provided by the coda
+ package. The posterior sample of the deviance \eqn{D}{D}, with
+ \eqn{D=-2\log(\prod_i P(y_i|\lambda_i))}{%
+ D=-2log(prod_i P(y_i|lambda_i))}, is also provided.}
+
+ \item{lambda.pred}{Predictive posterior mean for the exponential of
+ the latent variables. The approximation of Diggle et
+ al. (2004) is used to marginalized with respect to over-dispersion terms:
+ \deqn{E[\lambda_i|\beta,b_i,\sigma^2]=\phi^{-1}((X_i
+ \beta+W_i b_i)+0.5 \sigma^2)}{%
+ E[lambda_i|beta,b_i,sigma^2]=phi^(-1)((X_i*beta+W_i*b_i)+0.5*sigma^2)}
+ }
+}
+
+\references{
+ Siddhartha Chib and Bradley P. Carlin. 1999. ``On MCMC Sampling in
+ Hierarchical Longitudinal Models.'' \emph{Statistics and Computing.} 9:
+ 17-26.
+
+ Daniel Pemstein, Kevin M. Quinn, and Andrew D. Martin. 2007.
+ \emph{Scythe Statistical Library 1.0.} \url{http://scythe.wustl.edu}.
+
+ Andrew D. Martin and Kyle L. Saunders. 2002. ``Bayesian Inference for
+ Political Science Panel Data.'' Paper presented at the 2002 Annual Meeting
+ of the American Political Science Association.
+
+ Martyn Plummer, Nicky Best, Kate Cowles, and Karen Vines. 2002.
+ \emph{Output Analysis and Diagnostics for MCMC (CODA)}.
+ \url{http://www-fis.iarc.fr/coda/}.
+}
+
+\author{
+ Ghislain Vieilledent <ghislain.vieilledent at cirad.fr>
+}
+
+\seealso{
+ \code{\link[coda]{plot.mcmc}}, \code{\link[coda]{summary.mcmc}}
+}
+
+\examples{
+\dontrun{
+#========================================
+# Hierarchical Poisson Linear Regression
+#========================================
+
+#== Generating data
+
+# Constants
+nobs <- 1000
+nspecies <- 20
+species <- c(1:nspecies,sample(c(1:nspecies),(nobs-nspecies),replace=TRUE))
+
+# Covariates
+X1 <- runif(n=nobs,min=-1,max=1)
+X2 <- runif(n=nobs,min=-1,max=1)
+X <- cbind(rep(1,nobs),X1,X2)
+W <- X
+
+# Target parameters
+# beta
+beta.target <- matrix(c(0.1,0.1,0.1),ncol=1)
+# Vb
+Vb.target <- c(0.05,0.05,0.05)
+# b
+b.target <- cbind(rnorm(nspecies,mean=0,sd=sqrt(Vb.target[1])),
+ rnorm(nspecies,mean=0,sd=sqrt(Vb.target[2])),
+ rnorm(nspecies,mean=0,sd=sqrt(Vb.target[3])))
+
+# Response
+lambda <- vector()
+Y <- vector()
+for (n in 1:nobs) {
+ lambda[n] <- exp(X[n,]\%*\%beta.target+W[n,]\%*\%b.target[species[n],])
+ Y[n] <- rpois(1,lambda[n])
+}
+
+# Data-set
+Data <- as.data.frame(cbind(Y,lambda,X1,X2,species))
+plot(Data$X1,Data$lambda)
+
+#== Call to MCMChpoisson
+model <- MCMChpoisson(fixed=Y~X1+X2, random=~X1+X2, group="species",
+ data=Data, burnin=5000, mcmc=1000, thin=1,verbose=1,
+ seed=NA, beta.start=0, sigma2.start=1,
+ Vb.start=1, mubeta=0, Vbeta=1.0E6,
+ r=3, R=diag(c(0.1,0.1,0.1)), nu=0.001, delta=0.001, FixOD=1)
+
+#== MCMC analysis
+
+# Graphics
+pdf("Posteriors-MCMChpoisson.pdf")
+plot(model$mcmc)
+dev.off()
+
+# Summary
+summary(model$mcmc)
+
+# Predictive posterior mean for each observation
+model$lambda.pred
+
+# Predicted-Observed
+plot(Data$lambda,model$lambda.pred)
+abline(a=0,b=1)
+
+## #Not run
+## #You can also compare with lme4 results
+## #== lme4 resolution
+## library(lme4)
+## model.lme4 <- lmer(Y~X1+X2+(1+X1+X2|species),data=Data,family="poisson")
+## summary(model.lme4)
+## plot(fitted(model.lme4),model$lambda.pred,main="MCMChpoisson/lme4")
+## abline(a=0,b=1)
+}
+}
+
+\keyword{models}
+\keyword{hierarchical models}
+\keyword{mixed models}
+\keyword{glmm}
+\keyword{Poisson}
+\keyword{MCMC}
+\keyword{bayesian}
diff --git a/man/MCMChregress.Rd b/man/MCMChregress.Rd
new file mode 100644
index 0000000..84fe37e
--- /dev/null
+++ b/man/MCMChregress.Rd
@@ -0,0 +1,254 @@
+\name{MCMChregress}
+
+\alias{MCMChregress}
+
+\title{Markov Chain Monte Carlo for the Hierarchical Gaussian Linear
+ Regression Model}
+
+\description{MCMChregress generates a sample from the posterior
+ distribution of a Hierarchical Gaussian Linear Regression Model using
+ Algorithm 2 of Chib and Carlin (1999). This model uses a multivariate
+ Normal prior for the fixed effects parameters, an Inverse-Wishart
+ prior on the random effects variance matrix, and an Inverse-Gamma
+ prior on the residual error variance. The user supplies data and
+ priors, and a sample from the posterior distribution is returned as an
+ mcmc object, which can be subsequently analyzed with functions
+ provided in the coda package.}
+
+\usage{MCMChregress(fixed, random, group, data, burnin=1000,
+mcmc=10000, thin=10, verbose=1, seed=NA, beta.start=NA,
+sigma2.start=NA, Vb.start=NA, mubeta=0, Vbeta=1.0E6, r, R,
+nu=0.001, delta=0.001, ...)}
+
+\arguments{
+
+ \item{fixed}{A two-sided linear formula of the form 'y~x1+...+xp'
+ describing the fixed-effects part of the model, with the response on
+ the left of a '~' operator and the p fixed terms, separated by '+'
+ operators, on the right.}
+
+ \item{random}{A one-sided formula of the form '~x1+...+xq' specifying
+ the model for the random effects part of the model, with the q
+ random terms, separated by '+' operators.}
+
+ \item{group}{String indicating the name of the grouping variable
+ in \code{data}, defining the hierarchical structure of the model.}
+
+ \item{data}{A data frame containing the variables in the model.}
+
+ \item{burnin}{The number of burnin iterations for the sampler.}
+
+ \item{mcmc}{The number of Gibbs iterations for the sampler. Total
+ number of Gibbs iterations is equal to
+ \code{burnin+mcmc}. \code{burnin+mcmc} must be divisible by 10 and
+ superior or equal to 100 so that the progress bar can be displayed.}
+
+ \item{thin}{The thinning interval used in the simulation. The number of
+ mcmc iterations must be divisible by this value.}
+
+ \item{seed}{The seed for the random number generator. If NA, the
+ Mersenne Twister generator is used with default seed 12345; if an
+ integer is passed it is used to seed the Mersenne twister.}
+
+ \item{verbose}{A switch (0,1) which determines whether or not the
+ progress of the sampler is printed to the screen. Default is 1: a
+ progress bar is printed, indicating the step (in \%) reached by the
+ Gibbs sampler.}
+
+ \item{beta.start}{The starting values for the \eqn{\beta}{beta}
+ vector. This can either be a scalar or a p-length vector. The
+ default value of NA will use the OLS \eqn{\beta}{beta} estimate of
+ the corresponding Gaussian Linear Regression without random
+ effects. If this is a scalar, that value will serve as the starting
+ value mean for all of the betas.}
+
+ \item{sigma2.start}{Scalar for the starting value of the residual
+ error variance. The default value of NA will use the OLS estimates
+ of the corresponding Gaussian Linear Regression without random
+ effects.}
+
+ \item{Vb.start}{The starting value for variance matrix of the random
+ effects. This must be a square q-dimension matrix. Default value of NA
+ uses an identity matrix.}
+
+ \item{mubeta}{The prior mean of \eqn{\beta}{beta}. This can either be
+ a scalar or a p-length vector. If this takes a scalar value, then
+ that value will serve as the prior mean for all of the betas. The
+ default value of 0 will use a vector of zeros for an uninformative
+ prior.}
+
+ \item{Vbeta}{The prior variance of \eqn{\beta}{beta}. This can either
+ be a scalar or a square p-dimension matrix. If this takes a scalar value,
+ then that value times an identity matrix serves as the prior
+ variance of beta. Default value of 1.0E6 will use a diagonal matrix
+ with very large variance for an uninformative flat prior.}
+
+ \item{r}{The shape parameter for the Inverse-Wishart prior on variance
+ matrix for the random effects. r must be superior or equal to
+ q. Set r=q for an uninformative prior. See the NOTE for more details}
+
+ \item{R}{The scale matrix for the Inverse-Wishart prior on variance matrix
+ for the random effects. This must be a square q-dimension
+ matrix. Use plausible variance regarding random effects for the
+ diagonal of R. See the NOTE for more details}
+
+ \item{nu}{The shape parameter for the Inverse-Gamma prior on the
+ residual error variance. Default value is \code{nu=delta=0.001} for
+ uninformative prior.}
+
+ \item{delta}{The rate (1/scale) parameter for the Inverse-Gamma prior on the
+ residual error variance. Default value is \code{nu=delta=0.001} for
+ uninformative prior.}
+
+ \item{...}{further arguments to be passed}
+
+}
+
+\details{
+
+ \code{MCMChregress} simulates from the posterior distribution sample
+ using the blocked Gibbs sampler of Chib and Carlin (1999), Algorithm
+ 2. The simulation is done in compiled C++ code to maximize
+ efficiency. Please consult the coda documentation for a comprehensive
+ list of functions that can be used to analyze the posterior sample.
+
+ The model takes the following form:
+ \deqn{y_i = X_i \beta + W_i b_i + \varepsilon_i}{y_i = X_i * beta + W_i *
+ b_i + epsilon_i}
+
+ Where each group \eqn{i} have \eqn{k_i} observations.
+
+ Where the random effects:
+ \deqn{b_i \sim \mathcal{N}_q(0,V_b)}{b_i ~ N_q(0,Vb)}
+ And the errors:
+ \deqn{\varepsilon_i \sim \mathcal{N}(0, \sigma^2 I_{k_i})}{epsilon_i ~ N(0,
+ sigma^2 I_{k_i})}
+ We assume standard, conjugate priors:
+ \deqn{\beta \sim \mathcal{N}_p(\mu_{\beta},V_{\beta})}{beta ~ N_p(mubeta,Vbeta)}
+ And:
+ \deqn{\sigma^{2} \sim \mathcal{IG}amma(\nu, 1/\delta)}{sigma^2 ~
+ IGamma(nu, 1/delta)}
+ And:
+ \deqn{V_b \sim \mathcal{IW}ishart(r, rR)}{Vb ~ IWishart(r,
+ rR)}
+ See Chib and Carlin (1999) for more details.
+
+ \emph{NOTE:} We do not provide default parameters for the priors on
+ the precision matrix for the random effects. When fitting one of
+ these models, it is of utmost importance to choose a prior that
+ reflects your prior beliefs about the random effects. Using the
+ \code{dwish} and \code{rwish} functions might be useful in choosing
+ these values.
+}
+
+\value{
+
+ \item{mcmc}{An mcmc object that contains the posterior sample. This
+ object can be summarized by functions provided by the coda
+ package. The posterior sample of the deviance \eqn{D}{D}, with
+ \eqn{D=-2\log(\prod_i P(y_i|\beta,b_i,\sigma^2))}{%
+ D=-2log(prod_i P(y_i|beta,b_i,sigma^2))}, is also provided.}
+
+ \item{Y.pred}{Predictive posterior mean for each observation.}
+}
+
+\references{
+ Siddhartha Chib and Bradley P. Carlin. 1999. ``On MCMC Sampling in
+ Hierarchical Longitudinal Models.'' \emph{Statistics and Computing.} 9:
+ 17-26.
+
+ Daniel Pemstein, Kevin M. Quinn, and Andrew D. Martin. 2007.
+ \emph{Scythe Statistical Library 1.0.} \url{http://scythe.wustl.edu}.
+
+ Andrew D. Martin and Kyle L. Saunders. 2002. ``Bayesian Inference for
+ Political Science Panel Data.'' Paper presented at the 2002 Annual Meeting
+ of the American Political Science Association.
+
+ Martyn Plummer, Nicky Best, Kate Cowles, and Karen Vines. 2002.
+ \emph{Output Analysis and Diagnostics for MCMC (CODA)}.
+ \url{http://www-fis.iarc.fr/coda/}.
+}
+
+\author{
+ Ghislain Vieilledent <ghislain.vieilledent at cirad.fr>
+}
+
+\seealso{
+ \code{\link[coda]{plot.mcmc}}, \code{\link[coda]{summary.mcmc}}
+}
+
+\examples{
+\dontrun{
+#========================================
+# Hierarchical Gaussian Linear Regression
+#========================================
+
+#== Generating data
+
+# Constants
+nobs <- 1000
+nspecies <- 20
+species <- c(1:nspecies,sample(c(1:nspecies),(nobs-nspecies),replace=TRUE))
+
+# Covariates
+X1 <- runif(n=nobs,min=0,max=10)
+X2 <- runif(n=nobs,min=0,max=10)
+X <- cbind(rep(1,nobs),X1,X2)
+W <- X
+
+# Target parameters
+# beta
+beta.target <- matrix(c(0.1,0.3,0.2),ncol=1)
+# Vb
+Vb.target <- c(0.5,0.2,0.1)
+# b
+b.target <- cbind(rnorm(nspecies,mean=0,sd=sqrt(Vb.target[1])),
+ rnorm(nspecies,mean=0,sd=sqrt(Vb.target[2])),
+ rnorm(nspecies,mean=0,sd=sqrt(Vb.target[3])))
+# sigma2
+sigma2.target <- 0.02
+
+# Response
+Y <- vector()
+for (n in 1:nobs) {
+ Y[n] <- rnorm(n=1,
+ mean=X[n,]\%*\%beta.target+W[n,]\%*\%b.target[species[n],],
+ sd=sqrt(sigma2.target))
+}
+
+# Data-set
+Data <- as.data.frame(cbind(Y,X1,X2,species))
+plot(Data$X1,Data$Y)
+
+#== Call to MCMChregress
+model <- MCMChregress(fixed=Y~X1+X2, random=~X1+X2, group="species",
+ data=Data, burnin=1000, mcmc=1000, thin=1,verbose=1,
+ seed=NA, beta.start=0, sigma2.start=1,
+ Vb.start=1, mubeta=0, Vbeta=1.0E6,
+ r=3, R=diag(c(1,0.1,0.1)), nu=0.001, delta=0.001)
+
+#== MCMC analysis
+
+# Graphics
+pdf("Posteriors-MCMChregress.pdf")
+plot(model$mcmc)
+dev.off()
+
+# Summary
+summary(model$mcmc)
+
+# Predictive posterior mean for each observation
+model$Y.pred
+
+# Predicted-Observed
+plot(Data$Y,model$Y.pred)
+abline(a=0,b=1)
+}
+}
+
+\keyword{models}
+\keyword{hierarchical models}
+\keyword{mixed models}
+\keyword{Gaussian}
+\keyword{MCMC}
+\keyword{bayesian}
diff --git a/src/MCMChierBetaBinom.cc b/src/MCMChierBetaBinom.cc
new file mode 100644
index 0000000..6a16f25
--- /dev/null
+++ b/src/MCMChierBetaBinom.cc
@@ -0,0 +1,458 @@
+//////////////////////////////////////////////////////////////////////////
+// MCMChierBetaBinom.cc is C++ code to fit a hierarchical beta binomial
+// model
+//
+// y_{ij} ~ Binomial(s_{ij}, theta_{ij})
+// theta_{ij} ~ Beta(alpha_j, beta_j)
+// alpha_j ~ Pareto(1, a)
+// beta_j ~ Pareto(1, b)
+//
+//
+// uses adaptive Metropolis scheme similar to that of Haario et al. (2001)
+// as implemented by Roberts and Rosenthal (2006/2008)
+// "Examples of Adaptive MCMC"
+//
+// This software is distributed under the terms of the GNU GENERAL
+// PUBLIC LICENSE Version 2, June 1991. See the package LICENSE
+// file for more information.
+//
+// 5/28/2011 KQ
+//
+// Copyright (C) 2003-2007 Andrew D. Martin and Kevin M. Quinn
+// Copyright (C) 2007-present Andrew D. Martin, Kevin M. Quinn,
+// and Jong Hee Park
+//////////////////////////////////////////////////////////////////////////
+
+
+
+
+
+#ifndef MCMCHIERBETABINOM_CC
+#define MCMCHIERBETABINOM_CC
+
+#include <iostream>
+
+#include "matrix.h"
+#include "algorithm.h"
+#include "distributions.h"
+#include "stat.h"
+#include "la.h"
+#include "ide.h"
+#include "smath.h"
+#include "MCMCrng.h"
+#include "MCMCfcds.h"
+
+#include <R.h> // needed to use Rprintf()
+#include <R_ext/Utils.h> // needed to allow user interrupts
+
+typedef Matrix<double,Row,View> rmview;
+
+using namespace std;
+using namespace scythe;
+
+
+// used to access Ydata like a 2d array
+#define M(ROW,COL,NROWS) (COL*NROWS+ROW)
+
+
+
+
+// log of the pareto density
+double logdpareto(const double& x, const double& xm, const double& a){
+ double logfunval;
+ if (x > xm && a > 0){
+ logfunval = log(a) + a*log(xm) - (a+1)*log(x);
+ }
+ else{
+ logfunval = -numeric_limits<double>::infinity();
+ }
+ return logfunval;
+}
+
+
+
+
+// log of the full conditional density for alpha_j, beta_j
+double logABfcd(const double& alpha, const double& beta,
+ const vector<const double*>& theta,
+ const double& a, const double& b){
+
+ double term1 = 0.0;
+ double term2 = 0.0;
+ const int len_theta = theta.size();
+ if (alpha > 1.0 && beta > 1.0){
+ for (int i=0; i<len_theta; ++i){
+ term1 += lndbeta1(*theta[i], alpha, beta);
+ }
+ }
+ else{
+ term1 + -numeric_limits<double>::infinity();
+ }
+
+ // a and/or b <= 0 is treated as improper uniform prior
+ if (a > 0){
+ term2 += logdpareto(alpha, 1.0, a);
+ }
+ if (b > 0){
+ term2 += logdpareto(beta, 1.0, b);
+ }
+ double logfcd = term1 + term2;
+
+ return logfcd;
+
+}
+
+
+// candidate generating function
+template <typename RNGTYPE>
+Matrix<double> mixcangen(const double& alpha, const double& beta,
+ const double& mixfrac, const double& base_sigma,
+ const double& alpha_var_n, const double& beta_var_n,
+ const double& ab_cov_n,
+ rng<RNGTYPE>& stream){
+
+ Matrix<double> ab_V(2, 2, false);
+ ab_V = alpha_var_n, ab_cov_n, ab_cov_n, beta_var_n;
+ ab_V = (5.6644 * ab_V) / 2.0;
+
+ Matrix<double> ab_mean(2, 1, false);
+ ab_mean = alpha, beta;
+
+ Matrix<double> ab_can(2, 1, false);
+ double u = stream.runif();
+ if (u < mixfrac){
+ double alpha_can = stream.rnorm(alpha, base_sigma);
+ double beta_can = stream.rnorm(beta, base_sigma);
+ ab_can = alpha_can, beta_can;
+ }
+ else{
+ ab_can = stream.rmvnorm(ab_mean, ab_V);
+ }
+
+ return ab_can;
+
+}
+
+
+
+
+template <typename RNGTYPE>
+void hierBetaBinom_impl(rng<RNGTYPE>& stream,
+ double* sampledata, const int samplerow,
+ const int samplecol,
+ const int* y,
+ const int* s,
+ const double* theta_start,
+ const double* alpha_start,
+ const double* beta_start,
+ const double a,
+ const double b,
+ const int* ilabels,
+ const int* jlabels,
+ const int* ilabelsunique,
+ const int* jlabelsunique,
+ const int n,
+ const int ni,
+ const int nj,
+ const int burnin,
+ const int mcmc,
+ const int thin,
+ const int verbose,
+ int * accepts,
+ const double* base_sigma){
+
+ const int tot_iter = burnin + mcmc;
+ const int nstore = mcmc/thin;
+
+ // these probably should not be hard coded
+ const double mixfrac = 0.05;
+
+
+
+ double* theta;
+ theta = new double[n];
+ for (int i=0; i<n; ++i){
+ theta[i] = theta_start[i];
+ }
+
+ double* alpha;
+ alpha = new double[nj];
+ for (int i=0; i<nj; ++i){
+ alpha[i] = alpha_start[i];
+ }
+
+ double* beta;
+ beta = new double[nj];
+ for (int i=0; i<nj; ++i){
+ beta[i] = beta_start[i];
+ }
+
+ double* alpha_storage;
+ alpha_storage = new double[nj * tot_iter];
+ for (int i=0; i<(nj*tot_iter); ++i){
+ alpha_storage[i] = -999.0;
+ }
+
+ double* beta_storage;
+ beta_storage = new double[nj * tot_iter];
+ for (int i=0; i<(nj*tot_iter); ++i){
+ beta_storage[i] = -999.0;
+ }
+
+ double* sums_alpha_n;
+ sums_alpha_n = new double[nj];
+ for (int i=0; i<nj; ++i){
+ sums_alpha_n[i] = 0.0;
+ }
+
+ double* sums_beta_n;
+ sums_beta_n = new double[nj];
+ for (int i=0; i<nj; ++i){
+ sums_beta_n[i] = 0.0;
+ }
+
+
+
+
+
+
+ vector< vector<const double*> > js_thetas_ptr;
+ js_thetas_ptr.reserve(nj);
+ for (int j=0; j<nj; ++j){
+ vector<const double*> holder;
+ holder.reserve(n);
+ for (int i=0; i<n; ++i){
+ if (jlabels[i] == (j+1)){
+ holder.push_back( (const double*)&theta[i] );
+ }
+ }
+ js_thetas_ptr.push_back(holder);
+ }
+
+
+
+ int count = 0; // counter for output storage
+ // MCMC iterations
+ for (int iter=0; iter<tot_iter; ++iter){
+
+ // sample [theta_{ij} | alpha_j, beta_j, y, s]
+ for (int i=0; i<n; ++i){
+ int cluster_id = jlabels[i];
+ theta[i] = stream.rbeta(static_cast<double>(y[i]) + alpha[cluster_id-1],
+ static_cast<double>(s[i] - y[i]) +
+ beta[cluster_id - 1]);
+ }
+
+
+
+
+
+
+
+ // sample [alpha_j, beta_j | theta, y, s, a, b]
+ for (int j=0; j<nj; ++j){
+ const int len_thetaj = js_thetas_ptr[j].size();
+ double logfcd_cur = logABfcd(alpha[j], beta[j],
+ js_thetas_ptr[j], a, b);
+ if (iter < 10){
+ double alpha_can = stream.rnorm(alpha[j], base_sigma[j]);
+ double beta_can = stream.rnorm(beta[j], base_sigma[j]);
+ double accept_numer;
+ double accept_denom;
+ accept_numer = logABfcd(alpha_can, beta_can,
+ js_thetas_ptr[j], a, b);
+ accept_denom = logABfcd(alpha[j], beta[j],
+ js_thetas_ptr[j], a, b);
+ const double ratio = exp(accept_numer - accept_denom);
+
+ if (stream.runif() < ratio) {
+ alpha[j] = alpha_can;
+ beta[j] = beta_can;
+ ++accepts[j];
+ }
+ } // end iter < 10
+ else{ // adaptive MCMC after the first 10 iterations
+ // use iter and not iter -1 in 2 lines below b/c of start w/ 0
+ double alpha_mean_n = sums_alpha_n[j] / static_cast<double>(iter);
+ double beta_mean_n = sums_beta_n[j] / static_cast<double>(iter);
+ double alpha_var_n = 0.0;
+ double beta_var_n = 0.0;
+ double ab_cov_n = 0.0;
+ for (int i=0; i<iter; ++i){
+ const double alpha_dev = (alpha_mean_n - alpha_storage[M(i, j, tot_iter)]);
+ const double beta_dev = (beta_mean_n - beta_storage[M(i, j, tot_iter)]);
+ alpha_var_n += std::pow(alpha_dev, 2);
+ beta_var_n += std::pow(beta_dev, 2);
+ ab_cov_n += alpha_dev * beta_dev;
+ }
+ alpha_var_n = alpha_var_n / static_cast<double>(iter);
+ beta_var_n = beta_var_n / static_cast<double>(iter);
+ ab_cov_n = ab_cov_n / static_cast<double>(iter);
+ if (alpha_var_n <= 0.0){
+ alpha_var_n = std::pow(base_sigma[j], 2);
+ ab_cov_n = 0.0;
+ }
+ if (beta_var_n <= 0.0){
+ beta_var_n = std::pow(base_sigma[j], 2);
+ ab_cov_n = 0.0;
+ }
+
+ Matrix<double> ab_can = mixcangen(alpha[j], beta[j], mixfrac,
+ base_sigma[j], alpha_var_n,
+ beta_var_n, ab_cov_n, stream);
+ double alpha_can = ab_can(0);
+ double beta_can = ab_can(1);
+ double accept_numer;
+ double accept_denom;
+ accept_numer = logABfcd(alpha_can, beta_can,
+ js_thetas_ptr[j], a, b);
+ accept_denom = logABfcd(alpha[j], beta[j],
+ js_thetas_ptr[j], a, b);
+ const double ratio = exp(accept_numer - accept_denom);
+
+ if (stream.runif() < ratio) {
+ alpha[j] = alpha_can;
+ beta[j] = beta_can;
+ ++accepts[j];
+ }
+
+ } // end else iter not less than 10
+ sums_alpha_n[j] += alpha[j];
+ sums_beta_n[j] += beta[j];
+ alpha_storage[M(iter, j, tot_iter)] = alpha[j];
+ beta_storage[M(iter, j, tot_iter)] = beta[j];
+ } // end j loop
+
+
+
+
+
+
+
+
+
+ // store values
+ if (iter >= burnin && (iter % thin == 0)) {
+ for (int i=0; i<n; ++i){
+ sampledata[M(count, i, samplerow)] = theta[i];
+ }
+ for (int i=0; i<nj; ++i){
+ sampledata[M(count, (i+n), samplerow)] = alpha[i];
+ }
+ for (int i=0; i<nj; ++i){
+ sampledata[M(count, (i+n+nj), samplerow)] = beta[i];
+ }
+ ++count;
+ }
+
+
+
+
+
+
+
+ // print output to screen
+ if(verbose > 0 && iter % verbose == 0){
+ Rprintf("\n\niteration %i of %i \n", (iter+1), tot_iter);
+ }
+
+
+
+
+
+
+
+
+ // allow user interrupts
+ R_CheckUserInterrupt();
+
+ } // end MCMC iterations
+
+
+
+
+
+ // clear memory
+ delete [] theta;
+ delete [] alpha;
+ delete [] beta;
+ delete [] alpha_storage;
+ delete [] beta_storage;
+ delete [] sums_alpha_n;
+ delete [] sums_beta_n;
+
+
+} // end hierBetaBinom_impl
+
+
+
+
+
+
+
+
+
+
+
+
+
+extern "C"{
+
+ // function called by R to fit model
+ void
+ hierBetaBinom(double* sampledata,
+ const int* samplerow,
+ const int* samplecol,
+ const int* y,
+ const int* s,
+ const double* theta_start,
+ const double* alpha_start,
+ const double* beta_start,
+ const double* a,
+ const double* b,
+ const int* ilabels,
+ const int* jlabels,
+ const int* ilabelsunique,
+ const int* jlabelsunique,
+ const int* n,
+ const int* ni,
+ const int* nj,
+ const int* burnin, const int* mcmc, const int* thin,
+ const int *uselecuyer,
+ const int *seedarray,
+ const int *lecuyerstream,
+ const int* verbose,
+ int *accepts,
+ const double* base_sigma){
+
+
+ MCMCPACK_PASSRNG2MODEL(hierBetaBinom_impl,
+ sampledata,
+ *samplerow,
+ *samplecol,
+ y,
+ s,
+ theta_start,
+ alpha_start,
+ beta_start,
+ *a,
+ *b,
+ ilabels,
+ jlabels,
+ ilabelsunique,
+ jlabelsunique,
+ *n,
+ *ni,
+ *nj,
+ *burnin, *mcmc, *thin,
+ *verbose,
+ accepts,
+ base_sigma);
+
+ } // end hierBetaBinom
+
+} // end extern "C"
+
+
+
+#endif
diff --git a/src/MCMChlogit.cc b/src/MCMChlogit.cc
new file mode 100644
index 0000000..49f2391
--- /dev/null
+++ b/src/MCMChlogit.cc
@@ -0,0 +1,475 @@
+////////////////////////////////////////////////////////////////////
+// MCMChlogit.cc
+//
+// MCMChlogit samples from the posterior distribution of a
+// Binomial hierarchical linear regression model with a logit link function
+//
+// The code uses Algorithm 2 of Chib & Carlin (1999) for efficient
+// inference of (\beta | Y, sigma^2, Vb).
+//
+// Chib, S. & Carlin, B. P. (1999) On MCMC sampling in hierarchical
+// longitudinal models, Statistics and Computing, 9, 17-26
+//
+////////////////////////////////////////////////////////////////////
+//
+// Original code by Ghislain Vieilledent, may 2011
+// CIRAD UR B&SEF
+// ghislain.vieilledent at cirad.fr / ghislainv at gmail.com
+//
+////////////////////////////////////////////////////////////////////
+//
+// The initial version of this file was generated by the
+// auto.Scythe.call() function in the MCMCpack R package
+// written by:
+//
+// Andrew D. Martin
+// Dept. of Political Science
+// Washington University in St. Louis
+// admartin at wustl.edu
+//
+// Kevin M. Quinn
+// Dept. of Government
+// Harvard University
+// kevin_quinn at harvard.edu
+//
+// This software is distributed under the terms of the GNU GENERAL
+// PUBLIC LICENSE Version 2, June 1991. See the package LICENSE
+// file for more information.
+//
+// Copyright (C) 2011 Andrew D. Martin and Kevin M. Quinn
+//
+////////////////////////////////////////////////////////////////////
+//
+// Revisions:
+// - This file was initially generated on Wed May 4 10:42:50 2011
+// - G. Vieilledent, on May 9 2011
+//
+////////////////////////////////////////////////////////////////////
+
+#include "matrix.h"
+#include "distributions.h"
+#include "stat.h"
+#include "la.h"
+#include "ide.h"
+#include "smath.h"
+#include "MCMCrng.h"
+#include "MCMCfcds.h"
+
+#include <R.h> // needed to use Rprintf()
+#include <R_ext/Utils.h> // needed to allow user interrupts
+
+using namespace scythe;
+using namespace std;
+
+/* Fonction logit */
+double logit (double x) {
+ double Result=log(x)-log(1-x);
+ return Result;
+}
+/* Fonction invlogit */
+double invlogit (double x) {
+ double Result=1/(1+exp(-x));
+ return Result;
+}
+
+extern "C"{
+
+ /* Gibbs sampler function */
+ void MCMChlogit (
+
+ // Constants and data
+ const int *ngibbs, const int *nthin, const int *nburn, // Number of iterations, burning and samples
+ const int *nobs, const int *ngroup, // Constants
+ const int *np, const int *nq, // Number of fixed and random covariates
+ const int *IdentGroup, // Vector of group
+ const double *Y_vect, // Observed response variable
+ const double *X_vect, // Covariate for fixed effects
+ const double *W_vect, // Covariate for random effects
+ // Parameters to save
+ double *beta_vect, // Fixed effects
+ double *b_vect, // Random effects
+ double *Vb_vect, // Variance of random effects
+ double *V, // Variance of residuals
+ // Defining priors
+ const double *mubeta_vect, const double *Vbeta_vect,
+ const double *r, const double *R_vect,
+ const double *s1_V, const double *s2_V,
+ // Diagnostic
+ double *Deviance,
+ double *theta_pred, // Annual mortality rate
+ // Seeds
+ const int *seed,
+ // Verbose
+ const int *verbose,
+ // Overdispersion
+ const int *FixOD
+
+ ) {
+
+ ////////////////////////////////////////////////////////////////////////////////
+ //%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+ // Defining and initializing objects
+
+ ///////////////////////////
+ // Redefining constants //
+ const int NGIBBS=ngibbs[0];
+ const int NTHIN=nthin[0];
+ const int NBURN=nburn[0];
+ const int NSAMP=(NGIBBS-NBURN)/NTHIN;
+ const int NOBS=nobs[0];
+ const int NGROUP=ngroup[0];
+ const int NP=np[0];
+ const int NQ=nq[0];
+
+ ///////////////
+ // Constants //
+ // Number of observations by group k
+ int *nobsk = new int[NGROUP];
+ for (int k=0; k<NGROUP; k++) {
+ nobsk[k]=0;
+ for (int n=0; n<NOBS; n++) {
+ if (IdentGroup[n]==k) {
+ nobsk[k]+=1;
+ }
+ }
+ }
+ // Position of each group in the data-set
+ int **posk_arr = new int*[NGROUP];
+ for (int k=0; k<NGROUP; k++) {
+ posk_arr[k] = new int[nobsk[k]];
+ int repk=0;
+ for (int n=0; n<NOBS; n++) {
+ if (IdentGroup[n]==k) {
+ posk_arr[k][repk]=n;
+ repk++;
+ }
+ }
+ }
+ // Small fixed matrices indexed on k for data access
+ Matrix<double> *Xk_arr = new Matrix<double>[NGROUP];
+ Matrix<double> *Wk_arr = new Matrix<double>[NGROUP];
+ for(int k=0; k<NGROUP; k++) {
+ Xk_arr[k] = Matrix<double>(nobsk[k],NP);
+ Wk_arr[k] = Matrix<double>(nobsk[k],NQ);
+ for (int m=0; m<nobsk[k]; m++) {
+ for (int p=0; p<NP; p++) {
+ Xk_arr[k](m,p)=X_vect[p*NOBS+posk_arr[k][m]];
+ }
+ for (int q=0; q<NQ; q++) {
+ Wk_arr[k](m,q)=W_vect[q*NOBS+posk_arr[k][m]];
+ }
+ }
+ }
+ Matrix<double> *tXk_arr = new Matrix<double>[NGROUP];
+ Matrix<double> *tWk_arr = new Matrix<double>[NGROUP];
+ Matrix<double> *cpXk_arr = new Matrix<double>[NGROUP];
+ Matrix<double> *tXWk_arr = new Matrix<double>[NGROUP];
+ Matrix<double> *tWXk_arr = new Matrix<double>[NGROUP];
+ for(int k=0; k<NGROUP; k++) {
+ tXk_arr[k] = t(Xk_arr[k]);
+ tWk_arr[k] = t(Wk_arr[k]);
+ cpXk_arr[k] = crossprod(Xk_arr[k]);
+ tXWk_arr[k] = t(Xk_arr[k])*Wk_arr[k];
+ tWXk_arr[k] = t(Wk_arr[k])*Xk_arr[k];
+ }
+
+ ////////////////////
+ // Priors objects //
+ Matrix<double> mubeta(NP,1,mubeta_vect);
+ Matrix<double> Vbeta(NP,NP,Vbeta_vect);
+ Matrix<double> R(NQ,NQ,R_vect);
+
+ /////////////////////////////////////
+ // Initializing running parameters //
+ Matrix<double> *bk_run = new Matrix<double>[NGROUP]; // Random effects
+ for (int k=0;k<NGROUP;k++) {
+ bk_run[k] = Matrix<double>(NQ,1);
+ for (int q=0; q<NQ; q++) {
+ bk_run[k](q)=b_vect[q*NGROUP*NSAMP+k*NSAMP];
+ }
+ }
+ Matrix<double> beta_run(NP,1,false); // Unicolumn matrix of fixed effects
+ for (int p=0; p<NP; p++) {
+ beta_run(p)=beta_vect[p*NSAMP];
+ }
+ Matrix<double> Vb_run(NQ,NQ,true,0.0);
+ for (int q=0; q<NQ; q++) {
+ for (int qprim=0; qprim<NQ; qprim++) {
+ Vb_run(q,qprim)=Vb_vect[qprim*NQ*NSAMP+q*NSAMP];
+ }
+ }
+ double V_run=V[0];
+ Matrix<double> *logit_thetak_run = new Matrix<double>[NGROUP];
+ for (int k=0;k<NGROUP;k++) {
+ logit_thetak_run[k] = Matrix<double>(nobsk[k],1);
+ for (int m=0; m<nobsk[k]; m++) {
+ logit_thetak_run[k](m)=logit(theta_pred[posk_arr[k][m]]);
+ theta_pred[posk_arr[k][m]]=0; // We reinitialize theta_pred to zero to compute the posterior mean
+ }
+ }
+ double Deviance_run=Deviance[0];
+
+ ////////////////////////////////////////////////////////////
+ // Proposal variance and acceptance for adaptive sampling //
+ double *sigmap = new double[NOBS];
+ double *nA = new double[NOBS];
+ for (int n=0; n<NOBS; n++) {
+ nA[n]=0;
+ sigmap[n]=1;
+ }
+ double *Ar = new double[NOBS]; // Acceptance rate
+ for (int n=0; n<NOBS; n++) {
+ Ar[n]=0.0;
+ }
+
+ ////////////
+ // Message//
+ Rprintf("\nRunning the Gibbs sampler. It may be long, keep cool :)\n\n");
+ R_FlushConsole();
+ //R_ProcessEvents(); for windows
+
+ /////////////////////
+ // Set random seed //
+ mersenne myrng;
+ myrng.initialize(*seed);
+
+ ///////////////////////////////////////////////////////////////////////////////////////
+ //%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+ // Gibbs sampler (see Chib et Carlin, 1999 p.4 "Blocking for Gaussian mixed models")
+
+ for (int g=0;g<NGIBBS;g++) {
+
+ ////////////////////////////////////////////////
+ // vector logit_theta : Metropolis algorithm //
+
+ for (int k=0; k<NGROUP; k++) {
+
+ for (int m=0; m<nobsk[k]; m++) {
+
+ // Which one ?
+ int w=posk_arr[k][m];
+
+ // Mean of the prior
+ Matrix<double> logit_theta_hat=Xk_arr[k](m,_)*beta_run+Wk_arr[k](m,_)*bk_run[k];
+
+ // Proposal
+ double logit_theta_prop=myrng.rnorm(logit_thetak_run[k](m),sigmap[w]);
+
+ // theta
+ double theta_prop=invlogit(logit_theta_prop);
+ double theta_run=invlogit(logit_thetak_run[k](m));
+
+ // Ratio of probabilities
+ double p_prop=log(dbinom(Y_vect[w],1,theta_prop))+
+ log(dnorm(logit_theta_prop,logit_theta_hat(0),sqrt(V_run)));
+
+ double p_now=log(dbinom(Y_vect[w],1,theta_run))+
+ log(dnorm(logit_thetak_run[k](m),logit_theta_hat(0),sqrt(V_run)));
+
+ double r=exp(p_prop-p_now); // ratio
+ double z=myrng.runif();
+
+ // Actualization
+ if (z < r) {
+ logit_thetak_run[k](m)=logit_theta_prop;
+ nA[w]++;
+ }
+ }
+ }
+
+
+ //////////////////////////////////
+ // vector beta: Gibbs algorithm //
+
+ // invVi, sum_V and sum_v
+ Matrix<double> sum_V(NP,NP);
+ Matrix<double> sum_v(NP,1);
+ for (int k=0; k<NGROUP; k++) {
+ sum_V += (1/V_run)*cpXk_arr[k]-pow(1/V_run,2)*tXWk_arr[k]*invpd(invpd(Vb_run)+tWk_arr[k]*(1/V_run)*Wk_arr[k])*tWXk_arr[k];
+ sum_v += (1/V_run)*tXk_arr[k]*logit_thetak_run[k]-pow(1/V_run,2)*tXWk_arr[k]*invpd(invpd(Vb_run)+tWk_arr[k]*(1/V_run)*Wk_arr[k])*tWk_arr[k]*logit_thetak_run[k];
+ }
+
+ // big_V
+ Matrix<double> big_V=invpd(invpd(Vbeta)+sum_V/V_run);
+
+ // small_v
+ Matrix<double> small_v=invpd(Vbeta)*mubeta+sum_v/V_run;
+
+ // Draw in the posterior distribution
+ beta_run=myrng.rmvnorm(big_V*small_v,big_V);
+
+
+ ///////////////////////////////
+ // vector b: Gibbs algorithm //
+
+ // Loop on group
+ for (int k=0; k<NGROUP; k++) {
+
+ // big_Vk
+ Matrix<double> big_Vk=invpd(invpd(Vb_run)+crossprod(Wk_arr[k])/V_run);
+
+ // small_vk
+ Matrix<double> small_vk=(t(Wk_arr[k])*(logit_thetak_run[k]-Xk_arr[k]*beta_run))/V_run;
+
+ // Draw in the posterior distribution
+ bk_run[k]=myrng.rmvnorm(big_Vk*small_vk,big_Vk);
+ }
+
+
+ ////////////////////////////////////////////
+ // vector of variance Vb: Gibbs algorithm //
+ Matrix<double> SSBb(NQ,NQ,true,0.0);
+ for(int k=0; k<NGROUP; k++) {
+ SSBb+=bk_run[k]*t(bk_run[k]);
+ }
+ int Vb_dof=(*r)+(NGROUP);
+ Matrix<double> Vb_scale = invpd(SSBb+(*r)*R);
+ Vb_run=invpd(myrng.rwish(Vb_dof,Vb_scale));
+
+
+ ////////////////
+ // variance V //
+
+ // e
+ Matrix<double> e(1,1,true,0.0);
+ for (int k=0; k<NGROUP; k++) {
+ e+=crossprod(logit_thetak_run[k]-Xk_arr[k]*beta_run-Wk_arr[k]*bk_run[k]);
+ }
+
+ // Parameters
+ double S1=*s1_V+(NOBS/2); //shape
+ double S2=*s2_V+0.5*e(0); //rate
+
+ // Draw in the posterior distribution
+ if (FixOD[0]==1) {
+ V_run=V[0];
+ }
+ else {
+ V_run=1/myrng.rgamma(S1,S2);
+ }
+
+
+ //////////////////////////////////////////////////
+ //// Deviance
+
+ // logLikelihood
+ double logLk=0;
+ for (int k=0; k<NGROUP; k++) {
+ for (int m=0; m<nobsk[k]; m++) {
+ // Which one ?
+ int w=posk_arr[k][m];
+ // theta_run
+ double theta_run=invlogit(logit_thetak_run[k](m));
+ // L
+ logLk+=log(dbinom(Y_vect[w],1,theta_run));
+ }
+ }
+
+ // Deviance
+ Deviance_run=-2*logLk;
+
+
+ //////////////////////////////////////////////////
+ // Output
+ if(((g+1)>NBURN) && (((g+1)%(NTHIN))==0)){
+ int isamp=((g+1)-NBURN)/(NTHIN);
+ for (int p=0; p<NP; p++) {
+ beta_vect[p*NSAMP+(isamp-1)]=beta_run(p);
+ }
+ for (int k=0; k<NGROUP; k++) {
+ for (int q=0; q<NQ; q++) {
+ b_vect[q*NGROUP*NSAMP+k*NSAMP+(isamp-1)]=bk_run[k](q);
+ }
+ }
+ for (int q=0; q<NQ; q++) {
+ for (int qprim=0; qprim<NQ; qprim++) {
+ Vb_vect[qprim*NQ*NSAMP+q*NSAMP+(isamp-1)]=Vb_run(q,qprim);
+ }
+ }
+ V[isamp-1]=V_run;
+ Deviance[isamp-1]=Deviance_run;
+ for (int k=0;k<NGROUP;k++){
+ for (int m=0; m<nobsk[k]; m++) {
+ int w=posk_arr[k][m];
+ const double C=pow(16*sqrt(3)/(15*M_PI),2);
+ Matrix<double> logit_theta_hat=Xk_arr[k](m,_)*beta_run+Wk_arr[k](m,_)*bk_run[k];
+ theta_pred[w]+=invlogit(logit_theta_hat(0)/sqrt(1+C*V_run))/(NSAMP); // We compute the mean of NSAMP values
+ }
+ }
+ }
+
+
+ ///////////////////////////////////////////////////////
+ // Adaptive sampling (on the burnin period)
+ int DIV=0;
+ if (NGIBBS >=1000) DIV=100;
+ else DIV=NGIBBS/10;
+ if((g+1)%DIV==0 && (g+1)<=NBURN){
+ for (int n=0; n<NOBS; n++) {
+ const double ropt=0.44;
+ Ar[n]=nA[n]/(DIV*1.0);
+ if(Ar[n]>=ropt) sigmap[n]=sigmap[n]*(2-(1-Ar[n])/(1-ropt));
+ else sigmap[n]=sigmap[n]/(2-Ar[n]/ropt);
+ nA[n]=0.0; // We reinitialize the number of acceptance to zero
+ }
+ }
+ if((g+1)%DIV==0 && (g+1)>NBURN){
+ for (int n=0; n<NOBS; n++) {
+ Ar[n]=nA[n]/(DIV*1.0);
+ nA[n]=0.0; // We reinitialize the number of acceptance to zero
+ }
+ }
+
+
+ //////////////////////////////////////////////////
+ // Progress bar
+ double Perc=100*(g+1)/(NGIBBS);
+ if(((g+1)%(NGIBBS/100))==0 && (*verbose==1)){
+ Rprintf("*");
+ R_FlushConsole();
+ //R_ProcessEvents(); for windows
+ if(((g+1)%(NGIBBS/10))==0){
+ double mAr=0; // Mean acceptance rate
+ for (int n=0; n<NOBS; n++) {
+ mAr+=Ar[n]/(NOBS);
+ }
+ Rprintf(":%.1f%%, mean accept. rate=%.3f\n",Perc,mAr);
+ R_FlushConsole();
+ //R_ProcessEvents(); for windows
+ }
+ }
+
+
+ //////////////////////////////////////////////////
+ // User interrupt
+ void R_CheckUserInterrupt(void); // allow user interrupts
+
+ } // end MCMC loop
+
+ ///////////////
+ // Delete memory allocation
+ delete[] nobsk;
+ for(int k=0; k<NGROUP; k++) {
+ delete[] posk_arr[k];
+ }
+ delete[] sigmap;
+ delete[] nA;
+ delete[] posk_arr;
+ delete[] Xk_arr;
+ delete[] Wk_arr;
+ delete[] tXk_arr;
+ delete[] tWk_arr;
+ delete[] cpXk_arr;
+ delete[] tXWk_arr;
+ delete[] tWXk_arr;
+ delete[] bk_run;
+ delete[] logit_thetak_run;
+ delete[] Ar;
+
+ } // end MCMChlogit
+
+} // end extern "C"
+
+////////////////////////////////////////////////////////////////////
+// END
+////////////////////////////////////////////////////////////////////
diff --git a/src/MCMChpoisson.cc b/src/MCMChpoisson.cc
new file mode 100644
index 0000000..fb4ad62
--- /dev/null
+++ b/src/MCMChpoisson.cc
@@ -0,0 +1,463 @@
+////////////////////////////////////////////////////////////////////
+// MCMChpoisson.cc
+//
+// MCMChpoisson samples from the posterior distribution of a
+// Poisson hierarchical linear regression model with a log link function
+//
+// The code uses Algorithm 2 of Chib & Carlin (1999) for efficient
+// inference of (\beta | Y, sigma^2, Vb).
+//
+// Chib, S. & Carlin, B. P. (1999) On MCMC sampling in hierarchical
+// longitudinal models, Statistics and Computing, 9, 17-26
+//
+////////////////////////////////////////////////////////////////////
+//
+// Original code by Ghislain Vieilledent, may 2011
+// CIRAD UR B&SEF
+// ghislain.vieilledent at cirad.fr / ghislainv at gmail.com
+//
+////////////////////////////////////////////////////////////////////
+//
+// The initial version of this file was generated by the
+// auto.Scythe.call() function in the MCMCpack R package
+// written by:
+//
+// Andrew D. Martin
+// Dept. of Political Science
+// Washington University in St. Louis
+// admartin at wustl.edu
+//
+// Kevin M. Quinn
+// Dept. of Government
+// Harvard University
+// kevin_quinn at harvard.edu
+//
+// This software is distributed under the terms of the GNU GENERAL
+// PUBLIC LICENSE Version 2, June 1991. See the package LICENSE
+// file for more information.
+//
+// Copyright (C) 2011 Andrew D. Martin and Kevin M. Quinn
+//
+////////////////////////////////////////////////////////////////////
+//
+// Revisions:
+// - This file was initially generated on Wed May 4 10:42:50 2011
+// - G. Vieilledent, on May 9 2011
+//
+////////////////////////////////////////////////////////////////////
+
+#include "matrix.h"
+#include "distributions.h"
+#include "stat.h"
+#include "la.h"
+#include "ide.h"
+#include "smath.h"
+#include "MCMCrng.h"
+#include "MCMCfcds.h"
+
+#include <R.h> // needed to use Rprintf()
+#include <R_ext/Utils.h> // needed to allow user interrupts
+
+using namespace scythe;
+using namespace std;
+
+extern "C"{
+
+ /* Gibbs sampler function */
+ void MCMChpoisson (
+
+ // Constants and data
+ const int *ngibbs, const int *nthin, const int *nburn, // Number of iterations, burning and samples
+ const int *nobs, const int *ngroup, // Constants
+ const int *np, const int *nq, // Number of fixed and random covariates
+ const int *IdentGroup, // Vector of group
+ const double *Y_vect, // Observed response variable
+ const double *X_vect, // Covariate for fixed effects
+ const double *W_vect, // Covariate for random effects
+ // Parameters to save
+ double *beta_vect, // Fixed effects
+ double *b_vect, // Random effects
+ double *Vb_vect, // Variance of random effects
+ double *V, // Variance of residuals
+ // Defining priors
+ const double *mubeta_vect, const double *Vbeta_vect,
+ const double *r, const double *R_vect,
+ const double *s1_V, const double *s2_V,
+ // Diagnostic
+ double *Deviance,
+ double *lambda_pred, // Annual mortality rate
+ // Seeds
+ const int *seed,
+ // Verbose
+ const int *verbose,
+ // Overdispersion
+ const int *FixOD
+
+ ) {
+
+ ////////////////////////////////////////////////////////////////////////////////
+ //%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+ // Defining and initializing objects
+
+ ///////////////////////////
+ // Redefining constants //
+ const int NGIBBS=ngibbs[0];
+ const int NTHIN=nthin[0];
+ const int NBURN=nburn[0];
+ const int NSAMP=(NGIBBS-NBURN)/NTHIN;
+ const int NOBS=nobs[0];
+ const int NGROUP=ngroup[0];
+ const int NP=np[0];
+ const int NQ=nq[0];
+
+ ///////////////
+ // Constants //
+ // Number of observations by group k
+ int *nobsk = new int[NGROUP];
+ for (int k=0; k<NGROUP; k++) {
+ nobsk[k]=0;
+ for (int n=0; n<NOBS; n++) {
+ if (IdentGroup[n]==k) {
+ nobsk[k]+=1;
+ }
+ }
+ }
+ // Position of each group in the data-set
+ int **posk_arr = new int*[NGROUP];
+ for (int k=0; k<NGROUP; k++) {
+ posk_arr[k] = new int[nobsk[k]];
+ int repk=0;
+ for (int n=0; n<NOBS; n++) {
+ if (IdentGroup[n]==k) {
+ posk_arr[k][repk]=n;
+ repk++;
+ }
+ }
+ }
+ // Small fixed matrices indexed on k for data access
+ Matrix<double> *Xk_arr = new Matrix<double>[NGROUP];
+ Matrix<double> *Wk_arr = new Matrix<double>[NGROUP];
+ for(int k=0; k<NGROUP; k++) {
+ Xk_arr[k] = Matrix<double>(nobsk[k],NP);
+ Wk_arr[k] = Matrix<double>(nobsk[k],NQ);
+ for (int m=0; m<nobsk[k]; m++) {
+ for (int p=0; p<NP; p++) {
+ Xk_arr[k](m,p)=X_vect[p*NOBS+posk_arr[k][m]];
+ }
+ for (int q=0; q<NQ; q++) {
+ Wk_arr[k](m,q)=W_vect[q*NOBS+posk_arr[k][m]];
+ }
+ }
+ }
+ Matrix<double> *tXk_arr = new Matrix<double>[NGROUP];
+ Matrix<double> *tWk_arr = new Matrix<double>[NGROUP];
+ Matrix<double> *cpXk_arr = new Matrix<double>[NGROUP];
+ Matrix<double> *tXWk_arr = new Matrix<double>[NGROUP];
+ Matrix<double> *tWXk_arr = new Matrix<double>[NGROUP];
+ for(int k=0; k<NGROUP; k++) {
+ tXk_arr[k] = t(Xk_arr[k]);
+ tWk_arr[k] = t(Wk_arr[k]);
+ cpXk_arr[k] = crossprod(Xk_arr[k]);
+ tXWk_arr[k] = t(Xk_arr[k])*Wk_arr[k];
+ tWXk_arr[k] = t(Wk_arr[k])*Xk_arr[k];
+ }
+
+ ////////////////////
+ // Priors objects //
+ Matrix<double> mubeta(NP,1,mubeta_vect);
+ Matrix<double> Vbeta(NP,NP,Vbeta_vect);
+ Matrix<double> R(NQ,NQ,R_vect);
+
+ /////////////////////////////////////
+ // Initializing running parameters //
+ Matrix<double> *bk_run = new Matrix<double>[NGROUP]; // Random effects
+ for (int k=0;k<NGROUP;k++) {
+ bk_run[k] = Matrix<double>(NQ,1);
+ for (int q=0; q<NQ; q++) {
+ bk_run[k](q)=b_vect[q*NGROUP*NSAMP+k*NSAMP];
+ }
+ }
+ Matrix<double> beta_run(NP,1,false); // Unicolumn matrix of fixed effects
+ for (int p=0; p<NP; p++) {
+ beta_run(p)=beta_vect[p*NSAMP];
+ }
+ Matrix<double> Vb_run(NQ,NQ,true,0.0);
+ for (int q=0; q<NQ; q++) {
+ for (int qprim=0; qprim<NQ; qprim++) {
+ Vb_run(q,qprim)=Vb_vect[qprim*NQ*NSAMP+q*NSAMP];
+ }
+ }
+ double V_run=V[0];
+ Matrix<double> *log_lambdak_run = new Matrix<double>[NGROUP];
+ for (int k=0;k<NGROUP;k++) {
+ log_lambdak_run[k] = Matrix<double>(nobsk[k],1);
+ for (int m=0; m<nobsk[k]; m++) {
+ log_lambdak_run[k](m)=log(lambda_pred[posk_arr[k][m]]);
+ lambda_pred[posk_arr[k][m]]=0; // We reinitialize lambda_pred to zero to compute the posterior mean
+ }
+ }
+ double Deviance_run=Deviance[0];
+
+ ////////////////////////////////////////////////////////////
+ // Proposal variance and acceptance for adaptive sampling //
+ double *sigmap = new double[NOBS];
+ double *nA = new double[NOBS];
+ for (int n=0; n<NOBS; n++) {
+ nA[n]=0;
+ sigmap[n]=1;
+ }
+ double *Ar = new double[NOBS]; // Acceptance rate
+ for (int n=0; n<NOBS; n++) {
+ Ar[n]=0.0;
+ }
+
+ ////////////
+ // Message//
+ Rprintf("\nRunning the Gibbs sampler. It may be long, keep cool :)\n\n");
+ R_FlushConsole();
+ //R_ProcessEvents(); for windows
+
+ /////////////////////
+ // Set random seed //
+ mersenne myrng;
+ myrng.initialize(*seed);
+
+ ///////////////////////////////////////////////////////////////////////////////////////
+ //%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+ // Gibbs sampler (see Chib et Carlin, 1999 p.4 "Blocking for Gaussian mixed models")
+
+ for (int g=0;g<NGIBBS;g++) {
+
+ ////////////////////////////////////////////////
+ // vector log_lambda : Metropolis algorithm //
+
+ for (int k=0; k<NGROUP; k++) {
+
+ for (int m=0; m<nobsk[k]; m++) {
+
+ // Which one ?
+ int w=posk_arr[k][m];
+
+ // Mean of the prior
+ Matrix<double> log_lambda_hat=Xk_arr[k](m,_)*beta_run+Wk_arr[k](m,_)*bk_run[k];
+
+ // Proposal
+ double log_lambda_prop=myrng.rnorm(log_lambdak_run[k](m),sigmap[w]);
+
+ // lambda
+ double lambda_prop=exp(log_lambda_prop);
+ double lambda_run=exp(log_lambdak_run[k](m));
+
+ // Ratio of probabilities
+ double p_prop=log(dpois(Y_vect[w],lambda_prop))+
+ log(dnorm(log_lambda_prop,log_lambda_hat(0),sqrt(V_run)));
+
+ double p_now=log(dpois(Y_vect[w],lambda_run))+
+ log(dnorm(log_lambdak_run[k](m),log_lambda_hat(0),sqrt(V_run)));
+
+ double r=exp(p_prop-p_now); // ratio
+ double z=myrng.runif();
+
+ // Actualization
+ if (z < r) {
+ log_lambdak_run[k](m)=log_lambda_prop;
+ nA[w]++;
+ }
+ }
+ }
+
+
+ //////////////////////////////////
+ // vector beta: Gibbs algorithm //
+
+ // invVi, sum_V and sum_v
+ Matrix<double> sum_V(NP,NP);
+ Matrix<double> sum_v(NP,1);
+ for (int k=0; k<NGROUP; k++) {
+ sum_V += (1/V_run)*cpXk_arr[k]-pow(1/V_run,2)*tXWk_arr[k]*invpd(invpd(Vb_run)+tWk_arr[k]*(1/V_run)*Wk_arr[k])*tWXk_arr[k];
+ sum_v += (1/V_run)*tXk_arr[k]*log_lambdak_run[k]-pow(1/V_run,2)*tXWk_arr[k]*invpd(invpd(Vb_run)+tWk_arr[k]*(1/V_run)*Wk_arr[k])*tWk_arr[k]*log_lambdak_run[k];
+ }
+
+ // big_V
+ Matrix<double> big_V=invpd(invpd(Vbeta)+sum_V/V_run);
+
+ // small_v
+ Matrix<double> small_v=invpd(Vbeta)*mubeta+sum_v/V_run;
+
+ // Draw in the posterior distribution
+ beta_run=myrng.rmvnorm(big_V*small_v,big_V);
+
+
+ ///////////////////////////////
+ // vector b: Gibbs algorithm //
+
+ // Loop on group
+ for (int k=0; k<NGROUP; k++) {
+
+ // big_Vk
+ Matrix<double> big_Vk=invpd(invpd(Vb_run)+crossprod(Wk_arr[k])/V_run);
+
+ // small_vk
+ Matrix<double> small_vk=(t(Wk_arr[k])*(log_lambdak_run[k]-Xk_arr[k]*beta_run))/V_run;
+
+ // Draw in the posterior distribution
+ bk_run[k]=myrng.rmvnorm(big_Vk*small_vk,big_Vk);
+ }
+
+
+ ////////////////////////////////////////////
+ // vector of variance Vb: Gibbs algorithm //
+ Matrix<double> SSBb(NQ,NQ,true,0.0);
+ for(int k=0; k<NGROUP; k++) {
+ SSBb+=bk_run[k]*t(bk_run[k]);
+ }
+ int Vb_dof=(*r)+(NGROUP);
+ Matrix<double> Vb_scale = invpd(SSBb+(*r)*R);
+ Vb_run=invpd(myrng.rwish(Vb_dof,Vb_scale));
+
+
+ ////////////////
+ // variance V //
+
+ // e
+ Matrix<double> e(1,1,true,0.0);
+ for (int k=0; k<NGROUP; k++) {
+ e+=crossprod(log_lambdak_run[k]-Xk_arr[k]*beta_run-Wk_arr[k]*bk_run[k]);
+ }
+
+ // Parameters
+ double S1=*s1_V+(NOBS/2); //shape
+ double S2=*s2_V+0.5*e(0); //rate
+
+ // Draw in the posterior distribution
+ if (FixOD[0]==1) {
+ V_run=V[0];
+ }
+ else {
+ V_run=1/myrng.rgamma(S1,S2);
+ }
+
+
+ //////////////////////////////////////////////////
+ //// Deviance
+
+ // logLikelihood
+ double logLk=0;
+ for (int k=0; k<NGROUP; k++) {
+ for (int m=0; m<nobsk[k]; m++) {
+ // Which one ?
+ int w=posk_arr[k][m];
+ // lambda_run
+ double lambda_run=exp(log_lambdak_run[k](m));
+ // L
+ logLk+=log(dpois(Y_vect[w],lambda_run));
+ }
+ }
+
+ // Deviance
+ Deviance_run=-2*logLk;
+
+
+ //////////////////////////////////////////////////
+ // Output
+ if(((g+1)>NBURN) && (((g+1)%(NTHIN))==0)){
+ int isamp=((g+1)-NBURN)/(NTHIN);
+ for (int p=0; p<NP; p++) {
+ beta_vect[p*NSAMP+(isamp-1)]=beta_run(p);
+ }
+ for (int k=0; k<NGROUP; k++) {
+ for (int q=0; q<NQ; q++) {
+ b_vect[q*NGROUP*NSAMP+k*NSAMP+(isamp-1)]=bk_run[k](q);
+ }
+ }
+ for (int q=0; q<NQ; q++) {
+ for (int qprim=0; qprim<NQ; qprim++) {
+ Vb_vect[qprim*NQ*NSAMP+q*NSAMP+(isamp-1)]=Vb_run(q,qprim);
+ }
+ }
+ V[isamp-1]=V_run;
+ Deviance[isamp-1]=Deviance_run;
+ for (int k=0;k<NGROUP;k++){
+ for (int m=0; m<nobsk[k]; m++) {
+ int w=posk_arr[k][m];
+ Matrix<double> log_lambda_hat=Xk_arr[k](m,_)*beta_run+Wk_arr[k](m,_)*bk_run[k];
+ lambda_pred[w]+=exp(log_lambda_hat(0)+0.5*V_run)/(NSAMP); // We compute the mean of NSAMP values
+ }
+ }
+ }
+
+
+ ///////////////////////////////////////////////////////
+ // Adaptive sampling (on the burnin period)
+ int DIV=0;
+ if (NGIBBS >=1000) DIV=100;
+ else DIV=NGIBBS/10;
+ if((g+1)%DIV==0 && (g+1)<=NBURN){
+ for (int n=0; n<NOBS; n++) {
+ const double ropt=0.44;
+ Ar[n]=nA[n]/(DIV*1.0);
+ if(Ar[n]>=ropt) sigmap[n]=sigmap[n]*(2-(1-Ar[n])/(1-ropt));
+ else sigmap[n]=sigmap[n]/(2-Ar[n]/ropt);
+ nA[n]=0.0; // We reinitialize the number of acceptance to zero
+ }
+ }
+ if((g+1)%DIV==0 && (g+1)>NBURN){
+ for (int n=0; n<NOBS; n++) {
+ Ar[n]=nA[n]/(DIV*1.0);
+ nA[n]=0.0; // We reinitialize the number of acceptance to zero
+ }
+ }
+
+
+ //////////////////////////////////////////////////
+ // Progress bar
+ double Perc=100*(g+1)/(NGIBBS);
+ if(((g+1)%(NGIBBS/100))==0 && (*verbose==1)){
+ Rprintf("*");
+ R_FlushConsole();
+ //R_ProcessEvents(); for windows
+ if(((g+1)%(NGIBBS/10))==0){
+ double mAr=0; // Mean acceptance rate
+ for (int n=0; n<NOBS; n++) {
+ mAr+=Ar[n]/(NOBS);
+ }
+ Rprintf(":%.1f%%, mean accept. rate=%.3f\n",Perc,mAr);
+ R_FlushConsole();
+ //R_ProcessEvents(); for windows
+ }
+ }
+
+
+ //////////////////////////////////////////////////
+ // User interrupt
+ void R_CheckUserInterrupt(void); // allow user interrupts
+
+ } // end MCMC loop
+
+ ///////////////
+ // Delete memory allocation
+ delete[] nobsk;
+ for(int k=0; k<NGROUP; k++) {
+ delete[] posk_arr[k];
+ }
+ delete[] sigmap;
+ delete[] nA;
+ delete[] posk_arr;
+ delete[] Xk_arr;
+ delete[] Wk_arr;
+ delete[] tXk_arr;
+ delete[] tWk_arr;
+ delete[] cpXk_arr;
+ delete[] tXWk_arr;
+ delete[] tWXk_arr;
+ delete[] bk_run;
+ delete[] log_lambdak_run;
+ delete[] Ar;
+
+ } // end MCMChpoisson
+
+} // end extern "C"
+
+////////////////////////////////////////////////////////////////////
+// END
+////////////////////////////////////////////////////////////////////
diff --git a/src/MCMChregress.cc b/src/MCMChregress.cc
new file mode 100644
index 0000000..dfe42f9
--- /dev/null
+++ b/src/MCMChregress.cc
@@ -0,0 +1,381 @@
+////////////////////////////////////////////////////////////////////
+// MCMChregress.cc
+//
+// MCMChregress samples from the posterior distribution of a
+// Gaussian hierarchical linear regression model
+//
+// The code uses Algorithm 2 of Chib & Carlin (1999) for efficient
+// inference of (\beta | Y, sigma^2, Vb).
+//
+// Chib, S. & Carlin, B. P. (1999) On MCMC sampling in hierarchical
+// longitudinal models, Statistics and Computing, 9, 17-26
+//
+////////////////////////////////////////////////////////////////////
+//
+// Original code by Ghislain Vieilledent, may 2011
+// CIRAD UR B&SEF
+// ghislain.vieilledent at cirad.fr / ghislainv at gmail.com
+//
+////////////////////////////////////////////////////////////////////
+//
+// The initial version of this file was generated by the
+// auto.Scythe.call() function in the MCMCpack R package
+// written by:
+//
+// Andrew D. Martin
+// Dept. of Political Science
+// Washington University in St. Louis
+// admartin at wustl.edu
+//
+// Kevin M. Quinn
+// Dept. of Government
+// Harvard University
+// kevin_quinn at harvard.edu
+//
+// This software is distributed under the terms of the GNU GENERAL
+// PUBLIC LICENSE Version 2, June 1991. See the package LICENSE
+// file for more information.
+//
+// Copyright (C) 2011 Andrew D. Martin and Kevin M. Quinn
+//
+////////////////////////////////////////////////////////////////////
+//
+// Revisions:
+// - This file was initially generated on Wed May 4 10:42:50 2011
+// - G. Vieilledent, on May 4 2011
+//
+////////////////////////////////////////////////////////////////////
+
+#include "matrix.h"
+#include "distributions.h"
+#include "stat.h"
+#include "la.h"
+#include "ide.h"
+#include "smath.h"
+#include "MCMCrng.h"
+#include "MCMCfcds.h"
+
+#include <R.h> // needed to use Rprintf()
+#include <R_ext/Utils.h> // needed to allow user interrupts
+
+using namespace scythe;
+using namespace std;
+
+extern "C"{
+
+ /* Gibbs sampler function */
+ void MCMChregress (
+
+ // Constants and data
+ const int *ngibbs, const int *nthin, const int *nburn, // Number of iterations, burning and samples
+ const int *nobs, const int *ngroup, // Constants
+ const int *np, const int *nq, // Number of fixed and random covariates
+ const int *IdentGroup, // Vector of group
+ const double *Y_vect, // Observed response variable
+ const double *X_vect, // Covariate for fixed effects
+ const double *W_vect, // Covariate for random effects
+ // Parameters to save
+ double *beta_vect, // Fixed effects
+ double *b_vect, // Random effects
+ double *Vb_vect, // Variance of random effects
+ double *V, // Variance of residuals
+ // Defining priors
+ const double *mubeta_vect, const double *Vbeta_vect,
+ const double *r, const double *R_vect,
+ const double *s1_V, const double *s2_V,
+ // Diagnostic
+ double *Deviance,
+ double *Y_pred, // Fitted values (predictive posterior mean)
+ // Seeds
+ const int *seed,
+ // Verbose
+ const int *verbose
+
+ ) {
+
+ ////////////////////////////////////////////////////////////////////////////////
+ //%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+ // Defining and initializing objects
+
+ ///////////////////////////
+ // Redefining constants //
+ const int NGIBBS=ngibbs[0];
+ const int NTHIN=nthin[0];
+ const int NBURN=nburn[0];
+ const int NSAMP=(NGIBBS-NBURN)/NTHIN;
+ const int NOBS=nobs[0];
+ const int NGROUP=ngroup[0];
+ const int NP=np[0];
+ const int NQ=nq[0];
+
+ ///////////////
+ // Constants //
+ // Number of observations by group k
+ int *nobsk = new int[NGROUP];
+ for (int k=0; k<NGROUP; k++) {
+ nobsk[k]=0;
+ for (int n=0; n<NOBS; n++) {
+ if (IdentGroup[n]==k) {
+ nobsk[k]+=1;
+ }
+ }
+ }
+ // Position of each group in the data-set
+ int **posk_arr = new int*[NGROUP];
+ for (int k=0; k<NGROUP; k++) {
+ posk_arr[k] = new int[nobsk[k]];
+ int repk=0;
+ for (int n=0; n<NOBS; n++) {
+ if (IdentGroup[n]==k) {
+ posk_arr[k][repk]=n;
+ repk++;
+ }
+ }
+ }
+ // Small fixed matrices indexed on k for data access
+ Matrix<double> *Yk_arr = new Matrix<double>[NGROUP];
+ Matrix<double> *Xk_arr = new Matrix<double>[NGROUP];
+ Matrix<double> *Wk_arr = new Matrix<double>[NGROUP];
+ for(int k=0; k<NGROUP; k++) {
+ Xk_arr[k] = Matrix<double>(nobsk[k],NP);
+ Wk_arr[k] = Matrix<double>(nobsk[k],NQ);
+ Yk_arr[k] = Matrix<double>(nobsk[k],1);
+ for (int m=0; m<nobsk[k]; m++) {
+ for (int p=0; p<NP; p++) {
+ Xk_arr[k](m,p)=X_vect[p*NOBS+posk_arr[k][m]];
+ }
+ for (int q=0; q<NQ; q++) {
+ Wk_arr[k](m,q)=W_vect[q*NOBS+posk_arr[k][m]];
+ }
+ Yk_arr[k](m,0)=Y_vect[posk_arr[k][m]];
+ Y_pred[posk_arr[k][m]]=0; // We initialize Y_pred to zero to compute the predictive posterior mean
+ }
+ }
+ Matrix<double> *tXk_arr = new Matrix<double>[NGROUP];
+ Matrix<double> *tWk_arr = new Matrix<double>[NGROUP];
+ Matrix<double> *cpXk_arr = new Matrix<double>[NGROUP];
+ Matrix<double> *tXWk_arr = new Matrix<double>[NGROUP];
+ Matrix<double> *tWXk_arr = new Matrix<double>[NGROUP];
+ Matrix<double> *tXYk_arr = new Matrix<double>[NGROUP];
+ Matrix<double> *tWYk_arr = new Matrix<double>[NGROUP];
+ for(int k=0; k<NGROUP; k++) {
+ tXk_arr[k] = t(Xk_arr[k]);
+ tWk_arr[k] = t(Wk_arr[k]);
+ cpXk_arr[k] = crossprod(Xk_arr[k]);
+ tXWk_arr[k] = t(Xk_arr[k])*Wk_arr[k];
+ tWXk_arr[k] = t(Wk_arr[k])*Xk_arr[k];
+ tXYk_arr[k] = t(Xk_arr[k])*Yk_arr[k];
+ tWYk_arr[k] = t(Wk_arr[k])*Yk_arr[k];
+ }
+
+ ////////////////////
+ // Priors objects //
+ Matrix<double> mubeta(NP,1,mubeta_vect);
+ Matrix<double> Vbeta(NP,NP,Vbeta_vect);
+ Matrix<double> R(NQ,NQ,R_vect);
+
+ /////////////////////////////////////
+ // Initializing running parameters //
+ Matrix<double> *bk_run = new Matrix<double>[NGROUP]; // Random effects
+ for (int k=0;k<NGROUP;k++) {
+ bk_run[k] = Matrix<double>(NQ,1);
+ for (int q=0; q<NQ; q++) {
+ bk_run[k](q)=b_vect[q*NGROUP*NSAMP+k*NSAMP];
+ }
+ }
+ Matrix<double> beta_run(NP,1,false); // Unicolumn matrix of fixed effects
+ for (int p=0; p<NP; p++) {
+ beta_run(p)=beta_vect[p*NSAMP];
+ }
+ Matrix<double> Vb_run(NQ,NQ,true,0.0);
+ for (int q=0; q<NQ; q++) {
+ for (int qprim=0; qprim<NQ; qprim++) {
+ Vb_run(q,qprim)=Vb_vect[qprim*NQ*NSAMP+q*NSAMP];
+ }
+ }
+ double V_run=V[0];
+ double Deviance_run=Deviance[0];
+
+ ////////////
+ // Message//
+ Rprintf("\nRunning the Gibbs sampler. It may be long, keep cool :)\n\n");
+ R_FlushConsole();
+ //R_ProcessEvents(); for windows
+
+ /////////////////////
+ // Set random seed //
+ mersenne myrng;
+ myrng.initialize(*seed);
+
+
+ ///////////////////////////////////////////////////////////////////////////////////////
+ //%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+ // Gibbs sampler (see Chib et Carlin, 1999 p.4 "Blocking for Gaussian mixed models")
+
+ for (int g=0;g<NGIBBS;g++) {
+
+ //////////////////////////////////
+ // vector beta: Gibbs algorithm //
+
+ // invVi, sum_V and sum_v
+ // see http://en.wikipedia.org/wiki/Woodbury_matrix_identity with A=(1/V_run)*Id
+ Matrix<double> sum_V(NQ,NQ);
+ Matrix<double> sum_v(NQ,1);
+ for (int k=0; k<NGROUP; k++) {
+ sum_V += (1/V_run)*cpXk_arr[k]-pow(1/V_run,2)*tXWk_arr[k]*invpd(invpd(Vb_run)+tWk_arr[k]*(1/V_run)*Wk_arr[k])*tWXk_arr[k];
+ sum_v += (1/V_run)*tXYk_arr[k]-pow(1/V_run,2)*tXWk_arr[k]*invpd(invpd(Vb_run)+tWk_arr[k]*(1/V_run)*Wk_arr[k])*tWYk_arr[k];
+ }
+
+ // big_V
+ Matrix<double> big_V=invpd(invpd(Vbeta)+sum_V/V_run);
+
+ // small_v
+ Matrix<double> small_v=invpd(Vbeta)*mubeta+sum_v/V_run;
+
+ // Draw in the posterior distribution
+ beta_run=myrng.rmvnorm(big_V*small_v,big_V);
+
+
+ ///////////////////////////////
+ // vector b: Gibbs algorithm //
+
+ // Loop on group
+ for (int k=0; k<NGROUP; k++) {
+
+ // big_Vk
+ Matrix<double> big_Vk=invpd(invpd(Vb_run)+crossprod(Wk_arr[k])/V_run);
+
+ // small_vk
+ Matrix<double> small_vk=(t(Wk_arr[k])*(Yk_arr[k]-Xk_arr[k]*beta_run))/V_run;
+
+ // Draw in the posterior distribution
+ bk_run[k]=myrng.rmvnorm(big_Vk*small_vk,big_Vk);
+ }
+
+
+ ////////////////////////////////////////////
+ // vector of variance Vb: Gibbs algorithm //
+ Matrix<double> SSBb(NQ,NQ,true,0.0);
+ for(int k=0; k<NGROUP; k++) {
+ SSBb+=bk_run[k]*t(bk_run[k]);
+ }
+ int Vb_dof=(*r)+(NGROUP);
+ Matrix<double> Vb_scale = invpd(SSBb+(*r)*R);
+ Vb_run=invpd(myrng.rwish(Vb_dof,Vb_scale));
+
+
+ ////////////////
+ // variance V //
+
+ // e
+ Matrix<double> e(1,1,true,0.0);
+ for (int k=0; k<NGROUP; k++) {
+ e+=crossprod(Yk_arr[k]-Xk_arr[k]*beta_run-Wk_arr[k]*bk_run[k]);
+ }
+
+ // Parameters
+ double S1=*s1_V+(NOBS/2); //shape
+ double S2=*s2_V+0.5*e(0); //rate
+
+ // Draw in the posterior distribution
+ V_run=1/myrng.rgamma(S1,S2);
+
+ //////////////////////////////////////////////////
+ //// Deviance
+
+ // logLikelihood
+ double logLk=0;
+ for (int k=0; k<NGROUP; k++) {
+ for (int m=0; m<nobsk[k]; m++) {
+ // Which one ?
+ int w=posk_arr[k][m];
+ // Y_hat
+ Matrix<double> Y_hat=Xk_arr[k](m,_)*beta_run+Wk_arr[k](m,_)*bk_run[k];
+ // L
+ logLk+=log(dnorm(Y_vect[w],Y_hat(0),sqrt(V_run)));
+ }
+ }
+
+ // Deviance
+ Deviance_run=-2*logLk;
+
+
+ //////////////////////////////////////////////////
+ // Output
+ if(((g+1)>NBURN) && (((g+1)%(NTHIN))==0)) {
+ int isamp=((g+1)-NBURN)/(NTHIN);
+ for (int p=0; p<NP; p++) {
+ beta_vect[p*NSAMP+(isamp-1)]=beta_run(p);
+ }
+ for (int k=0; k<NGROUP; k++) {
+ for (int q=0; q<NQ; q++) {
+ b_vect[q*NGROUP*NSAMP+k*NSAMP+(isamp-1)]=bk_run[k](q);
+ }
+ }
+ for (int q=0; q<NQ; q++) {
+ for (int qprim=0; qprim<NQ; qprim++) {
+ Vb_vect[qprim*NQ*NSAMP+q*NSAMP+(isamp-1)]=Vb_run(q,qprim);
+ }
+ }
+ V[isamp-1]=V_run;
+ Deviance[isamp-1]=Deviance_run;
+ for (int k=0; k<NGROUP; k++) {
+ for (int m=0; m<nobsk[k]; m++) {
+ // Which one ?
+ int w=posk_arr[k][m];
+ // Y_hat
+ Matrix<double> Y_hat=Xk_arr[k](m,_)*beta_run+Wk_arr[k](m,_)*bk_run[k];
+ // Y_pred
+ Y_pred[w]+=Y_hat(0)/NSAMP;
+ }
+ }
+ }
+
+
+ //////////////////////////////////////////////////
+ // Progress bar
+ double Perc=100*(g+1)/(NGIBBS);
+ if(((g+1)%(NGIBBS/100))==0 && (*verbose==1)){
+ Rprintf("*");
+ R_FlushConsole();
+ //R_ProcessEvents(); for windows
+ if(((g+1)%(NGIBBS/10))==0){
+ Rprintf(":%.1f%%\n",Perc);
+ R_FlushConsole();
+ //R_ProcessEvents(); for windows
+ }
+ }
+
+
+ //////////////////////////////////////////////////
+ // User interrupt
+ void R_CheckUserInterrupt(void); // allow user interrupts
+
+ } // end MCMC loop
+
+ ///////////////
+ // Delete memory allocation
+ delete[] nobsk;
+ for(int k=0; k<NGROUP; k++) {
+ delete[] posk_arr[k];
+ }
+ delete[] posk_arr;
+ delete[] Yk_arr;
+ delete[] Xk_arr;
+ delete[] Wk_arr;
+ delete[] tXk_arr;
+ delete[] tWk_arr;
+ delete[] cpXk_arr;
+ delete[] tXWk_arr;
+ delete[] tWXk_arr;
+ delete[] tXYk_arr;
+ delete[] tWYk_arr;
+ delete[] bk_run;
+
+ } // end MCMChregress
+
+} // end extern "C"
+
+////////////////////////////////////////////////////////////////////
+// END
+////////////////////////////////////////////////////////////////////
--
Alioth's /usr/local/bin/git-commit-notice on /srv/git.debian.org/git/debian-science/packages/r-cran-mcmcpack.git
More information about the debian-science-commits
mailing list