rlmb: The Bayesian Linear Model Distribution
In knygren/glmbayes: Bayesian Generalized Linear Models (iid Samples)

Description Usage Arguments Details Value References See Also Examples

rlmb is used to generate iid samples from Bayesian Linear Models with multivariate normal priors. The model is specified by providing a data vector, a design matrix, and a pfamily (determining the prior distribution).

rlmb(n = 1, y, x, pfamily, offset = rep(0, nobs), weights = NULL)

## S3 method for class 'rlmb'
print(x, digits = max(3, getOption("digits") - 3), ...)

`n`	number of draws to generate. If `length(n) > 1`, the length is taken to be the number required.
`y`	a vector of observations of length `m`.
`x`	for `rlmb` a design matrix of dimension `m * p` and for `print.rlmb` the object to be printed.
`pfamily`	a description of the prior distribution and associated constants to be used in the model. This should be a pfamily function (see `pfamily` for details of pfamily functions.)
`offset`	this can be used to specify an a priori known component to be included in the linear predictor during fitting. This should be `NULL` or a numeric vector or matrix of extents matching those of the response. One or more `offset` terms can be included in the formula instead or as well, and if more than one are specified their sum is used. See `model.offset`.
`weights`	an optional vector of weights to be used in the fitting process. Should be `NULL` or a numeric vector. If non-NULL, weighted least squares is used with weights `weights` (that is, minimizing `sum(w*e^2)`); otherwise ordinary least squares is used. See also ‘Details’,
`digits`	the number of significant digits to use when printing.
`...`	additional arguments to be passed to the low level regression fitting functions (see below).

The rlmb function produces iid samples for Bayesian generalized linear models. It has a more minimialistic interface than than the lmb function. Core required inputs for the function include the data vector, the design matrix and a prior specification. In addition, the dispersion parameter must currently be provided for the gaussian, Gamma, quasipoisson, and quasibinomial families (future implementations may incorporate a prior for these into the rlmb function).

Current implemented models include the gaussian family (identity link function), the poisson and quasipoisson families (log link function), the gamma family (log link function), as well as the binomial and quasibinomial families (logit, probit, and cloglog link functions). The function returns the simulated Bayesian coefficients and some associated outputs.

For the gaussian family, iid samples from the posterior density are genererated using standard simulation procedures for multivariate normal densities. For all other families, the samples are generated using accept-reject procedures by leveraging the likelihood subgradient approach of Nygren and Nygren (2006). This approach relies on tight enveloping functions that bound the posterior density from above. The Gridtype parameter is used to select the method used for determining the size of a Grid used to build the enveloping function. See EnvelopeBuild for details. Depending on the selection, the time to build the envelope and the acceptance rate during the simulation process may vary. The returned value iters contains the number of candidates generated before acceptance for each draw.

rlmb returns a object of class "rlmb". The function summary (i.e., summary.rglmb) can be used to obtain or print a summary of the results. The generic accessor functions coefficients, fitted.values, residuals, and extractAIC can be used to extract various useful features of the value returned by rglmb. An object of class "rlmb" is a list containing at least the following components:

`coefficients`	a matrix of dimension `n` by `length(mu)` with one sample in each row
`coef.mode`	a vector of `length(mu)` with the estimated posterior mode coefficients
`dispersion`	Either a constant provided as part of the call, or a vector of length `n` with one sample in each row.
`Prior`	A list with the priors specified for the model in question. Items in the list may vary based on the type of prior
`prior.weights`	a vector of weights specified or implied by the model
`y`	a vector with the dependent variable
`x`	a matrix with the implied design matrix for the model
`famfunc`	Family functions used during estimation process
`iters`	an `n` by `1` matrix giving the number of candidates generated before acceptance for each sample.
`Envelope`	the envelope that was used during sampling

Chambers, J.M.(1992) Linear models. Chapter 4 of Statistical Models in S eds J. M. Chambers and T. J. Hastie, Wadsworth & Brooks/Cole.

Wilkinson, G.N. and Rogers, C.E. (1973). Symbolic descriptions of factorial models for analysis of variance. Applied Statistics, 22, 392-399. doi: 10.2307/2346786.

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

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

The classical modeling functions lm and glm.

pfamily for documentation of pfamily functions used to specify priors.

Prior_Setup, Prior_Check for functions used to initialize and to check priors,

summary.glmb, predict.glmb, simulate.glmb, extractAIC.glmb, dummy.coef.glmb and methods(class="glmb") for methods inherited from class glmb and the methods and generic functions for classes glm and lm from which class lmb also inherits.

Other glmbayes modeling functions: glmb(), lmb(), rglmb()

## Annette Dobson (1990) "An Introduction to Generalized Linear Models".
## Page 9: Plant Weight Data.
ctl <- c(4.17,5.58,5.18,6.11,4.50,4.61,5.17,4.53,5.33,5.14)
trt <- c(4.81,4.17,4.41,3.59,5.87,3.83,6.03,4.89,4.32,4.69)
group <- gl(2, 10, 20, labels = c("Ctl","Trt"))
weight <- c(ctl, trt)

ps=Prior_Setup(weight ~ group)
x=ps$x
mu=ps$mu
V=ps$Sigma
mu[1,1]=mean(weight)

Prior_Check(weight ~ group,family =gaussian(),
            pfamily=dNormal(mu=mu,Sigma=V))
lm.D9 <- lm(weight ~ group,x=TRUE,y=TRUE)
y=lm.D9$y

## Dispersion for maximum likelihood estimate
disp_ML=sigma(lm.D9)^2
n_prior=2
shape=n_prior/2
rate= disp_ML*shape

# Two-Block Gibbs sampler
set.seed(180)
dispersion2=disp_ML
# Run 1000 burn-in iterations
for(i in 1:1000){
  out1=rlmb(n = 1, y=y, x=x, pfamily=dNormal(mu=mu,Sigma=V,dispersion=dispersion2))
  out2=rlmb(n = 1, y=y, x=x, pfamily=dGamma(shape=shape,rate=rate,beta=out1$coefficients[1,]))
  dispersions=out2$dispersion
}


# Run 1000 additional iterations and store output
beta_out<-matrix(0,nrow=1000, ncol=2)
disp_out=rep(0,1000)
for(i in 1:1000){
  out1=rlmb(n = 1, y=y, x=x, pfamily=dNormal(mu=mu,Sigma=V,dispersion=dispersion2))
  out2=rlmb(n = 1, y=y, x=x, pfamily=dGamma(shape=shape,rate=rate,beta=out1$coefficients[1,]))
  dispersions=out2$dispersion
  beta_out[i,1:2]=out1$coefficients[1,1:2]
  disp_out[i]=out2$dispersion
}

colMeans(beta_out)
mean(disp_out)

# Same model using Independent_Normal_Gamma_Prior
lmb.D9_v2=lmb(weight ~ group,dIndependent_Normal_Gamma(mu,V,shape=shape,rate=rate))
summary(lmb.D9_v2)