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.
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 |
n |
number of draws to generate. If |
y |
a vector of observations of length |
x |
a design matrix of dimension |
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 |
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 |
digits |
Number of significant digits to use for printed outpute |
object |
an object of class |
... |
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 |
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()
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66 67 68 69 70 71 72 73 74 75 76 77 78 79 80 81 82 83 84 85 86 87 88 89 90 91 92 93 94 95 96 97 98 99 100 101 102 103 104 | ## 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
|
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.