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

Description Usage Arguments Details Value Examples

rnnorm_reg_std is used to generate iid samplers from Non-Gaussian Generalized Linear Models in standard form. The function should onlybe called after standardization of a Generalized Linear Model.

rnnorm_reg_std(
  n,
  y,
  x,
  mu,
  P,
  alpha,
  wt,
  f2,
  Envelope,
  family,
  link,
  progbar = 1L
)

`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`.
`mu`	a vector of length `p` giving the prior means of the variables in the design matrix.
`P`	a positive-definite symmetric matrix of dimension `p * p` specifying the prior precision matrix of the variable.
`alpha`	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`.
`wt`	an optional vector of ‘prior weights’ to be used in the fitting process. Should be NULL or a numeric vector.
`f2`	function used to calculate the negative of the log-posterior function
`Envelope`	an object of type `glmbenvelope`.
`family`	family used for simulation. Used that this is different from the family used in other functions.
`link`	link function used for simulation.
`progbar`	dummy for flagging if a progressbar should be produced during the call

This function uses the information contained in the constructed envelope list in order to sample from a model in standard form. The simulation proceeds as follows in order to generate each draw in the required number of samples.

1) A random number between 0 and 1 is generated and is used together with the information in the PLSD vector (from the envelope) in order to identify the part of the grid from which a candidate is to be generated.

2) For the part of the grid selected, the dimensions are looped through and a candidate component for each dimension is generated from a restricted normal using information from the Envelope (in particular, the values for logrt, loglt, and cbars corresponding to that the part of the grid selected and the dimension sampled)

3) The log-likelihood for the standardized model is evaluated for the generated candidate (note that the log-likelihood here includes the portion of the prior that was shifted to the log-likelihood)

4) An additional random number is generated and the log of this random number is compared to a log-acceptance rate that is calculated based on the candidate and the LLconst component from the Envelope component selected in order to determine if the candidate should be accepted or rejected

5) If the candidate was not accepted, the process above is repeated from step 1 until a candidate is accepted

A list consisting of the following:

`out`	A matrix with simulated draws from a model in standard form. Each row represents one draw from the density
`draws`	A vector with the number of candidates required before acceptance for each draw

data(menarche2)

summary(menarche2)
plot(Menarche/Total ~ Age, data=menarche2)

Age2=menarche2$Age-13

x<-matrix(as.numeric(1.0),nrow=length(Age2),ncol=2)
x[,2]=Age2

y=menarche2$Menarche/menarche2$Total
wt=menarche2$Total

mu<-matrix(as.numeric(0.0),nrow=2,ncol=1)
mu[2,1]=(log(0.9/0.1)-log(0.5/0.5))/3

V1<-1*diag(as.numeric(2.0))

# 2 standard deviations for prior estimate at age 13 between 0.1 and 0.9
## Specifies uncertainty around the point estimates

V1[1,1]<-((log(0.9/0.1)-log(0.5/0.5))/2)^2 
V1[2,2]=(3*mu[2,1]/2)^2  # Allows slope to be up to 1 times as large as point estimate 

#out<-rglmb(n = 1000, y=y, x=x, mu=mu, P=solve(V1), wt = wt, 
#           family = binomial(logit), Gridtype = 3) 
#summary(out)

famfunc<-glmbfamfunc(binomial(logit))

f1<-famfunc$f1
f2<-famfunc$f2  # Used in optim and glmbsim_cpp
f3<-famfunc$f3  # Used in optim
f5<-famfunc$f5
f6<-famfunc$f6

dispersion2<-as.numeric(1.0)
start <- mu
offset2=rep(as.numeric(0.0),length(y))
P=solve(V1)
n=1000


###### Adjust weight for dispersion

wt2=wt/dispersion2

######################### Shift mean vector to offset so that adjusted model has 0 mean

alpha=x%*%as.vector(mu)+offset2
mu2=0*as.vector(mu)
P2=P
x2=x


#####  Optimization step to find posterior mode and associated Precision

parin=start-mu

opt_out=optim(parin,f2,f3,y=as.vector(y),x=as.matrix(x),mu=as.vector(mu2),
              P=as.matrix(P),alpha=as.vector(alpha),wt=as.vector(wt2),
              method="BFGS",hessian=TRUE
)

bstar=opt_out$par  ## Posterior mode for adjusted model
bstar
bstar+as.vector(mu)  # mode for actual model
A1=opt_out$hessian # Approximate Precision at mode

## Standardize Model

Standard_Mod=glmb_Standardize_Model(y=as.vector(y), x=as.matrix(x),P=as.matrix(P),
                                    bstar=as.matrix(bstar,ncol=1), A1=as.matrix(A1))

bstar2=Standard_Mod$bstar2  
A=Standard_Mod$A
x2=Standard_Mod$x2
mu2=Standard_Mod$mu2
P2=Standard_Mod$P2
L2Inv=Standard_Mod$L2Inv
L3Inv=Standard_Mod$L3Inv

Env2=EnvelopeBuild(as.vector(bstar2), as.matrix(A),y, as.matrix(x2),
                   as.matrix(mu2,ncol=1),as.matrix(P2),as.vector(alpha),as.vector(wt2),
                   family="binomial",link="logit",Gridtype=as.integer(3), 
                   n=as.integer(n),sortgrid=TRUE)

## These now seem to match

Env2

#int n, NumericVector y, NumericMatrix x, NumericMatrix mu, NumericMatrix P, 
#NumericVector alpha, NumericVector wt, Function f2,
#Rcpp::List Envelope, Rcpp::CharacterVector family, 
#Rcpp::CharacterVector link, int progbar

### Note: getting the types correct here is important but potentially difficult for users
### May be better to call an R function wrapper that checks (and converts when possible) 
##  to correct types

sim=rnnorm_reg_std(n=as.integer(n),y=as.vector(y),x=as.matrix(x2),mu=as.matrix(mu2,ncol=1),
                   P=as.matrix(P2),alpha=as.vector(alpha),wt=as.vector(wt2),
                   f2=f2,Envelope=Env2,family="binomial",link="logit",as.integer(0))

out=L2Inv%*%L3Inv%*%t(sim$out)

for(i in 1:n){
  out[,i]=out[,i]+mu
}

summary(t(out))
mean(sim$draws)