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

Description Usage Arguments Details Value Author(s) References See Also Examples

lmb is used to fit Bayesian linear models, specified by giving a symbolic descriptions of the linear predictor and the prior distribution.

lmb(
  formula,
  pfamily,
  n = 1000,
  data,
  subset,
  weights,
  na.action,
  method = "qr",
  model = TRUE,
  x = TRUE,
  y = TRUE,
  qr = TRUE,
  singular.ok = TRUE,
  contrasts = NULL,
  offset,
  ...
)

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

`formula`	an object of class `"formula"` (or one that can be coerced to that class): a symbolic description of the model to be fitted. The details of model specification are given under ‘Details’.
`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).
`n`	number of draws to generate. If `length(n) > 1`, the length is taken to be the number required.
`data`	an optional data frame, list or environment (or object coercible by `as.data.frame` to a data frame) containing the variables in the model. If not found in `data`, the variables are taken from `environment(formula)`, typically the environment from which `lm` is called.
`subset`	an optional vector specifying a subset of observations to be used in the fitting process.
`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’,
`na.action`	a function which indicates what should happen when the data contain `NA`s. The default is set by the `na.action` setting of `options`, and is `na.fail` if that is unset. The ‘factory-fresh’ default is `stats{na.omit}`. Another possible value is `NULL`, no action. Value `stats{na.exclude}` can be useful.
`method`	the method to be used in fitting the classical model during a call to `glm`. The default method `glm.fit` uses iteratively reweighted least squares (IWLS): the alternative `"model.frame"` returns the model frame and does no fitting. User-supplied fitting functions can be supplied either as a function or a character string naming a function, with a function which takes the same arguments as `glm.fit`. If specified as a character string it is looked up from within the stats namespace.
`model`	logicals. If `TRUE` the corresponding components of the fit (the model frame, the model matrix, the response, the QR decomposition) are returned.
`x`	logicals. If `TRUE` the corresponding components of the fit (the model frame, the model matrix, the response, the QR decomposition) are returned.
`y`	logicals. If `TRUE` the corresponding components of the fit (the model frame, the model matrix, the response, the QR decomposition) are returned.
`qr`	logicals. If `TRUE` the corresponding components of the fit (the model frame, the model matrix, the response, the QR decomposition) are returned.
`singular.ok`	logical. If `FALSE` (the default in S but not in R) a singular fit is an error.
`contrasts`	an optional list. See the `contrasts.arg` of `model.matrix.default`.
`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`.
`...`	additional arguments to be passed to the low level regression fitting functions (see below).
`digits`	the number of significant digits to use when printing.

The function lmb is a Bayesian version of the classical lm function. Setup of the models mirrors those for the lm function but add the required argument pfamily with a prior formulation. The function generates random iid samples from the posterior density associated with the model.

The function returns the output from a call to the function lm as well as the simulated Bayesian coefficients and associated outputs. In addition, the function returns model diagnostic information related to the DIC, a Bayesian Information Criterion similar to the AIC for classical models.

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.

glmb returns an object of class "glmb". The function summary (i.e., summary.glmb) 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 codeglmb.

An object of class "lmb" is a list containing at least the following components:

`lm`	an object of class `"lm"` containing the output from a call to the function `lm`
`coefficients`	a matrix of dimension `n` by `length(mu)` with one sample in each row
`coef.means`	a vector of `length(mu)` with the estimated posterior mean coefficients
`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 list may vary based on the type of prior
`residuals`	a matrix of dimension `n` by `length(y)` with one sample for the deviance residuals in each row
`fitted.values`	a matrix of dimension `n` by `length(y)` with one sample for the fitted values in each row
`linear.predictors`	an `n` by `length(y)` matrix with one sample for the linear fit on the link scale in each row
`deviance`	an `n` by `1` matrix with one sample for the deviance in each row
`pD`	An Estimate for the effective number of parameters
`Dbar`	Expected value for minus twice the log-likelihood function
`Dthetabar`	Value of minus twice the log-likelihood function evaluated at the mean value for the coefficients
`DIC`	Estimated Deviance Information criterion
`weights`	a vector of weights specified or implied by the model
`prior.weights`	a vector of weights specified or implied by the model
`y`	a vector of observations of length `m`.
`x`	a design matrix of dimension `m * p`
`model`	if requested (the default),the model frame
`call`	the matched call
`formula`	the formula supplied
`terms`	the `terms` object used
`data`	the `data argument`
`famfunc`	family functions used during estimation and post processing
`iters`	an `n` by `1` matrix giving the number of candidates generated before acceptance for each sample.
`contrasts`	(where relevant) the contrasts used.
`xlevels`	(where relevant) a record of the levels of the factors used in fitting
`pfamily`	the prior family specified
`digits`	the number of significant digits to use when printing.

In addition, non-empty fits will have (yet to be implemented) components qr, R and effects relating to the final weighted linear fit for the posterior mode. Objects of class "glmb" are normall of class c("glmb","glm","lm"), that is inherit from classes glm and lm and well-designed methods from those classed will be applied when appropriate.

If a binomial glmb model was specified by giving a two-column response, the weights returned by prior.weights are the total number of cases (factored by the supplied case weights) and the component of y of the result is the proportion of successes.

The R implementation of lmb has been written by Kjell Nygren and was built to be a Bayesian version of the lm function and hence tries to mirror the features of the lm function to the greatest extent possible while also taking advantage of some of the method functions developed for the glmb function. For details on the author(s) for the lm function see the documentation for lm.

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(), rglmb(), rlmb()

## 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)
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)
disp_ML=sigma(lm.D9)^2
n_prior=2
shape=n_prior/2
rate= disp_ML*shape

# Conjugate Normal_Gamma Prior 
lmb.D9=lmb(weight ~ group,dNormal_Gamma(mu,V/disp_ML,shape=shape,rate=rate))
summary(lmb.D9)

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