simulate.glmb: Simulate Responses

View source: R/simulate.glmb.R

simulate.glmbR Documentation

Simulate Responses

Description

Simulate responses from the posterior predictive distribution corresponding to a fitted glmb object \insertCiteGelman2013glmbayes.

Usage

## S3 method for class 'glmb'
simulate(object, nsim = 1, seed = NULL, ...)

Arguments

object

An object of class glmb, typically the result of a call to the function glmb.

nsim

Defunct (see below).

seed

an object specifying if and how the random number generator should be initialized (seeded).

...

Additional arguments passed to the function. Will frequently include a matrix pred of simulated predictions from the predict function, the family (e.g., binomial) and an optional vector of weights specifying prior.weights for the simulated values (default is 1)

Value

Simulated values for data corresponding to simulated model predictions that correspond either to the original data or to a newdata data frame provided to the predict function.

References

\insertAllCited

See Also

predict.glmb, glmb, glmbayes-package; rglmb, rlmb, lmb; simulate (e.g. simulate.glm, simulate.lm for classical fits); see \insertCiteglmbayesChapter06glmbayes for model statistics.

Examples

data(menarche,package="MASS")

## ----Analysis Setup-----------------------------------------------------------
## Number of variables in model
Age=menarche$Age
nvars=2
## Reference Ages for setting of priors and Age_Difference
ref_age1=13  # user can modify this
ref_age2=15  ## user can modify this
## Define variables used later in analysis
Age2=menarche$Age-ref_age1
Age_Diff=ref_age2-ref_age1
mu1=as.matrix(c(0,1.098612),ncol=1)
V1<-1*diag(nvars)
V1[1,1]=0.18687882
V1[2,2]=0.10576217
V1[1,2]=-0.03389182
V1[2,1]=-0.03389182
Menarche_Model_Data=data.frame(Age=menarche$Age,Total=menarche$Total,
                               Menarche=menarche$Menarche,Age2)
glmb.out1<-glmb(n=1000,cbind(Menarche, Total-Menarche) ~Age2,family=binomial(logit),
                pfamily=dNormal(mu=mu1,Sigma=V1),data=Menarche_Model_Data)

# Prediction from original model
pred1=predict(glmb.out1,type="response")

## Posterior predictive check (bayesplot): disabled in package examples.
## Former pp_check.glmb() lives in legacy_code/pp_check.glmb.R — source after
## install.packages("bayesplot") to restore, then:
##   bayesplot::pp_check(glmb.out1, ndraws = 100L, fun = "hist")

## Get Original Residuals, their means, and credible bounds
res_out=residuals(glmb.out1)
colMeans(res_out, na.rm=TRUE)

## Set up to simulate new data and residuals
res_mean=colMeans(res_out, na.rm=TRUE)
res_low1=apply(res_out,2,FUN=quantile,probs=c(0.025),na.rm=TRUE)
res_high1=apply(res_out,2,FUN=quantile,probs=c(0.975),na.rm=TRUE)

## Simulate new data and get residuals for simulated data

ysim1=simulate(glmb.out1,nsim=1,seed=10401L,pred=pred1,family="binomial",
               prior.weights=weights(glmb.out1))


res_ysim_out1=residuals(glmb.out1,ysim=ysim1)
res_low=apply(res_ysim_out1,2,FUN=quantile,probs=c(0.025),na.rm=TRUE)
res_high=apply(res_ysim_out1,2,FUN=quantile,probs=c(0.975),na.rm=TRUE)


oldpar <- par(no.readonly = TRUE)
par(mar = c(5, 4, 4, 2) + 0.1)  # Standard margin setup

# Plot Credible Interval bounds for Deviance Residuals

plot(res_mean~Age,ylim=c(-2.5,2.5),
main="Credible Interval Bound for Menarche - Logit Model Deviance Residuals",
xlab = "Age", ylab = "Avg. Dev. Res")
lines(Age, 0*res_mean,lty=1)
lines(Age, res_low,lty=1)
lines(Age, res_high,lty=1)
lines(Age, res_low1,lty=2)
lines(Age, res_high1,lty=2)

par(oldpar)


glmbayes documentation built on Aug. 5, 2026, 1:07 a.m.

Related to simulate.glmb in glmbayes...