glmbayes: Bayesian Generalized Linear Models (iid Samples)
In knygren/glmbayes: Bayesian Generalized Linear Models (iid Samples)

Description Details Author(s) References Examples

Generates iid samples for Bayesian Generalized Linear Models.

The glmbayes package produces iid samples for Bayesian Genereralized Linear Models and is intended as a Bayesian version of the glm function for classical models. Estimation can be performed using three main functions. For models with fixed dispersion parameters, the rglmb function is the workhorse function and comes with a minimialistic interface for the input and output. It is also suitable for use as part of block Gibbs sampling procedures. The glmb function is essentially a wrapper function for the rglmb function that provides an interface closer to that of the glm function. The rGamma_reg function can be leveraged in order to produce samples for the dispersion parameters associated with the gaussian and Gamma link functions. Most methods defined for the output of the glm function are also defined for the glmb, rglmb, and rGamma_reg functions (see their respective documentation for details).

For the regression parameters, multivariate normal priors are assumed. Simulation for the gaussian family with the identify link function is performed using standard procedures for multivariate normal densities. For all other families and link functions, simulation is performed using the likelihood subgradient approach of Nygren and Nygren (2006). This approach involves the construction of an enveloping function for the full posterior density followed by accept-reject sampling. For models that are approximately multivariate normal, the expected number of draws required per acceptance are bounded from above as noted in Nygren and Nygren (2006).

Currently implemented models include the gaussian (identity link), poisson/quasipoisson (log link), binomial/quasibinomial (logit, probit, and cloglog links), and Gamma (log link) families. These models all have log-concave likelihood functions that allow us to leverage the likelihood-subgradient approach for the iid sampling. Models that fail to have log-concave likelihood functions are not implemented. Our demos (viewable by entering the demo() command) provides examples of each of these families and links.

The current implementation requires separate use of the rGamma_reg function in order to generate samples for dispersion parameters (gaussian, Gammma, quasipoisson, quasi-binomial families). Our demos include examples of the joint use of the rglmb and rGamma_reg to produce samples for both regression and dispersion parameters using two-block Gibbs samplers. As these two-block Gibbs samplers likely are geometrically ergodic, future implementations may incorporate these two-block Gibbs samplers into the rglmb and glmb functions by leveraging theoretical bounds om convergence rates derived using Rosenthal (1996) type drift and minorization conditions.

The rglmb function can also be used in Block-Gibbs sampling implementations for Hierarchical Bayesian models. The demos associated with this package contains examples of such models.

Kjell Nygren

Dobson, A. J. (1990) An Introduction to Generalized Linear Models. London: Chapman and Hall.

Hastie, T. J. and Pregibon, D. (1992) Generalized linear models. Chapter 6 of Statistical Models in S eds J. M. Chambers and T. J. Hastie, Wadsworth & Brooks/Cole. McCullagh P. and Nelder, J. A. (1989) Generalized Linear Models. London: Chapman and Hall.

Nygren, K.N. and Nygren, L.M (2006) Likelihood Subgradient Densities. Journal of the American Statistical Association. vol.101, no.475, pp 1144-1156.

Raiffa, Howard and Schlaifer, R (1961) Applied Statistical Decision Theory. Boston: Clinton Press, Inc.

Venables, W. N. and Ripley, B. D. (2002) Modern Applied Statistics with S. New York: Springer.

## Dobson (1990) Page 93: Randomized Controlled Trial :
counts <- c(18,17,15,20,10,20,25,13,12)
outcome <- gl(3,1,9)
treatment <- gl(3,3)
print(d.AD <- data.frame(treatment, outcome, counts))

## Call to glm
glm.D93 <- glm(counts ~ outcome + treatment, 
               family = poisson())

## Using glmb
## Step 1: Set up Prior
ps=Prior_Setup(counts ~ outcome + treatment)
mu=ps$mu
V=ps$Sigma
# Step2A: Check the Prior
Prior_Check(counts ~ outcome + treatment,family = poisson(),
            pfamily=dNormal(mu=mu,Sigma=V))
# Step2B: Update and Re-Check the Prior
mu[1,1]=log(mean(counts))
Prior_Check(counts ~ outcome + treatment,family = poisson(),
            pfamily=dNormal(mu=mu,Sigma=V))
# Step 3: Call the glmb function
glmb.D93<-glmb(counts ~ outcome + treatment, family=poisson(), 
               pfamily=dNormal(mu=mu,Sigma=V))

## ----Printed_Views------------------------------------------------------------
## Printed view of the output from the glm function 
print(glm.D93)
## Printed view of the output from the glmb function 
print(glmb.D93)

## ----Methods---------------------------------------------------------------
## Methods for class "lm"
methods(class="lm")

## Methods for class "glm"
methods(class="glm")

## Methods for class "glmb"
methods(class="glmb")

## ----summary--------------------------------------------------------------
## summary output for the "glm" class
summary(glm.D93)

## summary output for the "glm" class
summary(glmb.D93)

## ----fitted outputs-------------------------------------------------------
## fitted outputs for the glm function
fitted(glm.D93)

## ----glmb fitted outputs------------------------------------------------------
## mean of fitted outputs for the glm function
colMeans(fitted(glmb.D93))

## ----predictions----------------------------------------------------------
## predictions for the glm function
predict(glm.D93)

## predictions for the glmb function
colMeans(predict(glmb.D93)) 

## ----residuals------------------------------------------------------------
## residuals for the glm function
residuals(glm.D93)

## residuals for the glmb function
colMeans(residuals(glmb.D93))

## ----vcov-----------------------------------------------------------------
## vcov for the glm function
vcov(glm.D93)

## vcov for the glm function
vcov(glmb.D93)

## ----confint--------------------------------------------------------------
## confint for the glm function
confint(glm.D93)

## confint for the glm function
confint(glmb.D93)

## ----AIC/DIC------------------------------------------------------------------
## AIC for the glm function (equivalent degrees of freedom and the AIC)
extractAIC(glm.D93)

## DIC for the glmb function (estimated effective number of parameters and the DIC)
extractAIC(glmb.D93)

## ----Deviance-------------------------------------------------------------
## Deviance for the glm function
deviance(glm.D93)

## Deviance for the glmb function
mean(deviance(glmb.D93))

## ----logLik---------------------------------------------------------------
## Deviance for the glm function
logLik(glm.D93)

## Deviance for the glmb function
mean(logLik(glmb.D93))

## ----Model Frame----------------------------------------------------------
## Model Frame for the glm function
model.frame(glm.D93)

## Model Frame for the glmb function
model.frame(glmb.D93$glm)

## ----formula--------------------------------------------------------------
## formula for the glm function
formula(glm.D93)

## ----formula-------------------------------------------------------------
## formula for the glmb function
formula(glmb.D93)

## ----family--------------------------------------------------------------
## family for the glm function
family(glm.D93)

## family for the glmb function
family(glmb.D93$glm)

## ----nobs-----------------------------------------------------------------
## nobs for the glm function
nobs(glm.D93)

## nobs for the glmb function
nobs(glmb.D93)

## ----show-----------------------------------------------------------------
## show for the glm function
show(glm.D93)

## show for the glmb function
show(glmb.D93)