Description Usage Arguments Details Value Examples
View source: R/Envrnnorm_reg_std.R
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.
1 2 3 4 5 6 7 8 9 10 11 12 13 14 | rnnorm_reg_std(
n,
y,
x,
mu,
P,
alpha,
wt,
f2,
Envelope,
family,
link,
progbar = 1L
)
|
n |
number of draws to generate. If |
y |
a vector of observations of length |
x |
a design matrix of dimension |
mu |
a vector of length |
P |
a positive-definite symmetric matrix of dimension |
alpha |
this can be used to specify an a priori known component
to be included in the linear predictor during fitting. This should be
|
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 |
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 |
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 105 106 107 108 109 110 111 112 | 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)
|
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.