rNormal_reg: The Bayesian Generalized Linear Model with Normal Prior...
In knygren/glmbayes: Bayesian Generalized Linear Models (iid Samples)

Description Usage Arguments Details Value References See Also Examples

rNormal_reg is used to generate iid samples from Bayesian Generalized linear models with a normal prior. The model is specified by providing a data vector, a design matrix, and 2 prior constants (mu and Sigma) for the normal prior.

rNormal_reg(
  n,
  y,
  x,
  prior_list,
  offset = NULL,
  weights = 1,
  family = gaussian()
)

`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`	a design matrix of dimension `m * p`.
`prior_list`	a list with the two parameters (mu and Sigma) for the prior distribution.
`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 of length equal to the number of cases. One or more offset terms can be included in the formula instead or as well, and if more than one is specified their sum is used. See `model.offset`.
`weights`	an optional vector of ‘prior weights’ to be used in the fitting process. Should be NULL or a numeric vector.
`family`	a description of the error distribution and link function to be used in the model. For `glm` this can be a character string naming a family function, a family function or the result of a call to a family function. For `glm.fit` only the third option is supported. (See `family` for details of family functions.)

The rNormal_reg function produces iid samples for Bayesian generalized linear models. with multivariate-normal priors. Core required inputs for the function include the data vector, the design matrix and a prior specification (provided in a list when this function is called directly). Specifying the pfamily as dNormal in the lmb or glmb is equivalent to calling this function directly.

The function returns the simulated Bayesian coefficients and some associated outputs. While it is possible to call this function directly, it is generally recommended that the lmb, rlmb, glmb or rglmb functions be used instead as those functions have more documentation and the resulting objects come with more methods.

When the specified family is gaussian, the multivariate normal is the conjugate prior distribution for the likelihood function and the posterior distribution is therefore also a multivariate normal density. Fairly standard procedures are applied in order to generate samples in that case. Currently, this includes using a Cholesky decomposition. Later implementations may switch to a QR decomposition to increase consistency with the lm function and to increase associated numerical accuracy.

When the specified family is any log-concave non-gaussian family, then the estimation uses the Likelihood subgradient density approach of Nygren and Nygren (2006). This approaches uses tangencies to the log-likelihood function in order to construct an enveloping function from which candidate draws are generated and then either accepted/rejected using accept/reject methods. The core C function performing this simulation essentially goes through the following steps:

1) The model is standardized to have prior mean vector equal to 0 (i.e., offsets and any prior mean are combined into a constant term).

2) The posterior mode for this transformed model is found. Currently this uses the optim function. Later implementations may replace this with iteratively reweighted least squares (IWLS) to increase consistency with the glm function and to enhance numerical accuracy.

3) The model is further standardized so that (a) the precision matrix at the posterior mode is diagonal and (b) the prior variance-covariance matrix is the identity matrix (see glmb_Standardize_Model for details).

4) An enveloping function is built for the the standardized model, containing constants needed during simulation (see EnvelopeBuild for details).

5) Samples for the standardized model are generated using accept-reject methods (see rnnorm_reg_std for details).

6) The output from the standardized model are transformed back to the original scale by reversing the two eigenvalue decompositions and by adding back the prior mean.

rNormal_reg returns a object of class "rglmb". 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 rNormal_reg. An object of class "rglmb" is a list containing at least the following components:

`coefficients`	a `n` by `length(mu)` matrix with one sample in each row
`PostMode`	a vector of `length(mu)` with the estimated posterior mode coefficients
`Prior`	A list with two components. The first being the prior mean vector and the second the prior precision matrix
`iters`	an `n` by `1` matrix giving the number of candidates generated before acceptance for each sample.
`famfunc`	an object of class `"famfunc"`
`Envelope`	an object of class `"envelope"`
`dispersion`	an `n` by `1` matrix with simulated values for the dispersion
`loglike`	a `n` by `1` matrix containing the negative loglikelihood for each sample.

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.

Other family simulation functions: rGamma_reg(), rNormal_Gamma_reg(), rindependent_norm_gamma_reg()

## 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))

## Will move this step inside the Prior_Check
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

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

n<-10000
y<-lmb.D9$y
x<-as.matrix(lmb.D9$x)
prior_list=list(mu=mu,Sigma=V/disp_ML,shape=shape,rate=rate)
ngamma.D9=rNormal_Gamma_reg(n=1000,y=y,x=x,
  prior_list=prior_list)
summary(ngamma.D9$coefficients)