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

Description Usage Arguments Details Value References See Also Examples

rGamma_reg is used to generate iid samples for the dispersion parameter in Bayesian Generalized Linear models. The model is specified by providing the model structure, regression coefficient, and a Gamma prior.

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

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

## S3 method for class 'rGamma_reg'
summary(object, ...)

`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 prior parameters (shape and rate) for the prior distribution and a vector with the assumed value (beta) for the regression coefficients.
`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. This can be a character string naming a family function, a family function or the result of a call to a family function. (See `family` for details of family functions.)
`digits`	Number of significant digits to use for printed outpute
`object`	an object of class `"glmb_dispersion"` that is to be summarized
`...`	further arguments passed to or from other methods

The rGamma_reg function produces iid samples for the dispersion parameter in Bayesian generalized linear models (gaussian and Gamma families only). Core required inputs for the function include the data vector, the design matrix, an estimate for the regression coefficients, and a prior gamma distribution specification. The function returns the simulated Bayesian conditional dispersion and the prior specification. The most practical use of this function is likely to be in conjunction with the rglmb function as part of block Gibbs sampling procedues.

For the gaussian family, iid samples from the posterior density are genererated using standard simulation procedures for gamma distributions. For the Gamma family, , the samples are generated using accept-reject procedures by leveraging the likelihood subgradient approach of Nygren and Nygren (2006).

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

`dispersion`	an `n` by `1` matrix with simulated values for the dispersion
`Prior`	A list with two components. The first being the prior mean vector and the second the prior precision matrix

A reference

Other family simulation functions: rNormal_Gamma_reg(), rNormal_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)
lm.D9 <- lm(weight ~ group,x=TRUE,y=TRUE)

lm_summary=summary(lm.D9)

lm_summary
lm.D9$coefficients

n<-10000
y<-lm.D9$y
x<-as.matrix(lm.D9$x)

### Set old regression and dispersion coefficients

b_old=lm.D9$coefficients
v_old=lm_summary$sigma^2

#### Set up for rglmb_dispersion

n0=0.1
shape=n0/2
rate=shape*v_old

rate/(shape-1)
rate/shape

#a1<-shape+n1/2
#b1<-rate+sum(SS)/2

#out<-1/rgamma(n,shape=a1,rate=b1)
## v0=sum(SS)/n1 ~ (2*rate+sum(SS))/(2*shape+n1)=[(2*rate+sum(SS))/2]/((2*shape+n1)/2) 
n1=length(y)
SS=v_old*n1

n1 # This is equal to 20


n0=2 # Prior observations
v_prior=v_old  # Prior point estimate for variance (the mean of (1/dispersion=1/v_prior))
wt0=(n0/n1)  

## set shape=0.01*(n1/2)
## set rate= 0.01*SS/2

shape=wt0*(n1/2)   ###  Shape is prior observations /2
rate=shape*v_prior  ### rate is essentiall prior SS - V in rmultireg should be this
rate/shape ## Should match v_prior (currently also v_old)

## We see that this currently matches 
### (test different v_prior with various prior observations below) 

prior_list=list(beta=b_old,shape=shape,rate=rate)

dispout<-rGamma_reg(n=n,y,x,prior_list=prior_list,
offset= rep(0, length(y)),family=gaussian())

mean(dispout$dispersion) 
v_prior
v_old
#summary(dispout)  ## Summary function not working - need to write

### Set up test rglmb regressions without and with prior for dispersion 

mu<-c(0,0)
mu=b_old  ### For testing purposes, set prior=b_old
P<-0.1*diag(2)
wt2=rep(1,length(y))

### Check
prior=list(mu=mu,Sigma=solve(P),dispersion=v_old)
outtemp1<-glmb(n = 1000, weight ~ group, family = gaussian(),
pfamily=dNormal(mu=mu,Sigma=solve(P),dispersion=v_old))
## Could use a residuals function here -- For now, maybe run the glmb function

summary(outtemp1)
mean(colMeans(residuals(outtemp1)^2))
v_old
colMeans((outtemp1$coefficients))
b_old
prior=list(mu=mu,P=P,dispersion=v_old)
outtemp2<-rglmb(n = 1000, y, x, family = gaussian(),
pfamily=dNormal(mu=mu,Sigma=solve(P),dispersion=v_old),
offset = rep(0, length(y)), weights = wt2)
summary(outtemp2)
colMeans((outtemp2$coefficients))
b_old

prior=list(mu=mu,Sigma=solve(P),shape=shape,rate=rate)
outtemp3<-glmb(n = 10000, weight ~ group,family = gaussian(),  
               dNormal_Gamma(mu=mu,Sigma=solve(P),shape=shape,rate=rate))
## Could use a residuals function here -- For now, maybe run the glmb function

summary(outtemp3)
mean(colMeans(residuals(outtemp3)^2))
v_old
colMeans((outtemp3$coefficients))
b_old
mean(outtemp3$dispersion)  ## Seems slightly smaller --> Needs qc
v_old