In knygren/glmbayes: Bayesian Generalized Linear Models (iid Samples)

references: - id: nygren2006 title: Likelihood Subgradient Densities author: - family: Nygren given: Kjell - family: Nygren given: Lan Ma container-title: Journal of the American Statistical Association volume: 101 URL: 'https://doi.org/10.1198/016214506000000357' DOI: 10.1198/016214506000000357 issue: 475 publisher: Taylor and Francis, Ltd page: 1144-1156 type: article-journal issued: year: 2006 month: 9

id: Venables2002 title: Modern Applied Statistics with S. author:
family: Venables given: W. N.
family: Ripley given: B. D. URL: 'https://www.springer.com/gp/book/9780387954578' DOI: 10.1007/978-0-387-21706-2 publisher: Springer-Verlag New York page: 1-498 type: book issued: year: 2002
id: Hamilton2004 title: Lecture notes obstetrics and gynaecology author:
family: Hamilton-Fairley given: Diana URL: 'https://www.bookdepository.com/Lecture-Notes-Obstetrics-Gynaecology-Diana-Hamilton-Fairley/9781405120661' publisher: Blackwell Pub. type: book issued: year: 2004 month: 11
id: Evans2006 title: Checking for Prior-Data Conflict author:
family: Evans given: Michael
family: Moshonov given: Hadas container-title: Bayesian Analysis volume: 1 URL: 'http://dx.doi.org/10.1214/06-BA129' DOI: 10.1214/06-BA129 issue: 4 page: 893-914 type: article-journal issued: year: 2006 month:
id: Spiegelhalter2002 title: Bayesian measures of model complexity and fit author:
family: Spiegelhalter given: David J.
family: Best given: Nicola G.
family: Carlin given: Bradley P.
family: Van Der Linde given: Angelika container-title: Journal of the Royal Statistical Society Series B (Statistical Methodology) volume: 64 URL: 'https://doi.org/10.1111/1467-9868.00353' DOI: 10.1111/1467-9868.00353 issue: 4 page: 583-639 type: article-journal issued: year: 2002 month: 9

nocite: | @nygren2006 ...

knitr::opts_chunk$set(
  collapse = TRUE,
  comment = "#>"
)

library(glmbayes)

## Load menarche data
data(menarche2)
head(menarche2, 5)

## Number of variables in model
Age=menarche2$Age
nvars=2
set.seed(333)

## 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=Age-ref_age1
Age_Diff=ref_age2-ref_age1

## Point estimates at reference ages
m1=0.5  
m2=0.9

## Lower bound of prior credible intervals for point estimates
m1_lower=0.3
m2_lower=0.7

## Assumed correlation between the two (on link scale)
m_corr=0.4

## Set up link function and initialize prior mean and Variance-Covariance matrices
bi_logit <- binomial(link="logit")
mu1<-matrix(0,nrow=nvars,ncol=1)
rownames(mu1)=c("Intercept","Age2")
colnames(mu1)=c("Prior Mean")
V1<-1*diag(nvars)
rownames(V1)=c("Intercept","Age2")
colnames(V1)=c("Intercept","Age2")

## Prior mean for intercept is set to point estimate 
## at reference age1 (on logit scale)
mu1[1,1]=bi_logit$linkfun(m1)

## Prior mean for slope is set to difference in point estimates
## on logit scale divided by Age_Diff

mu1[2,1]=(bi_logit$linkfun(m2) -bi_logit$linkfun(m1))/Age_Diff 
print(mu1)

## Implied standard deviations for point estimates on logit scale

sd_m1= (bi_logit$linkfun(m1) -bi_logit$linkfun(m1_lower))/1.96
sd_m2= (bi_logit$linkfun(m2) -bi_logit$linkfun(m2_lower))/1.96

## Also compute implied estimate for upper bound of confidence intervals

m1_upper=bi_logit$linkinv(bi_logit$linkfun(m1)+sd_m1*1.96)
m2_upper=bi_logit$linkinv(bi_logit$linkfun(m2)+sd_m2*1.96)
print("m1_upper is:")
m1_upper
print("m2_upper is:")
m2_upper


## Implied Standard deviation for slope (using variance formula for difference between two variables)
a=(1/Age_Diff)
sd_slope=sqrt((a*sd_m1)^2+(a*sd_m2)^2-2*a*a*(sd_m1*sd_m2*m_corr))

#Cov(m1,slope)=cov(m1, a*(m2-m1)) =a*E[(m1-E[m1])((m2-m1)-E[m2-m1])]
#   =a*E[(m1-E[m1])(m2-E[m2])]- a* E[(m1-E[m1])(m1-E[m1])]
##   =a*Cov[m1,m2] - a*Var[m1]
##  =a*sd_m1*sd_m2*m_corr-a* sd_m1*sd_m1
cov_V1=a*sd_m1*sd_m2*m_corr-a* sd_m1*sd_m1

# Set covariance matrix
V1[1,1]=sd_m1^2
V1[2,2]=sd_slope^2
V1[1,2]=cov_V1
V1[2,1]=V1[1,2]
print("V1 is:")
print(V1)

Menarche_Model_Data=data.frame(Age=menarche2$Age,Total=menarche2$Total,Menarche=menarche2$Menarche,Age2)
prior1=list(mu=mu1,Sigma=V1)
#glmb.out1<-glmb(n=1000,cbind(Menarche, Total-Menarche) ~ #Age2,family=binomial(logit),mu=mu1,Sigma=V1,data=Menarche_Model_Data)
glmb.out1<-glmb(n=1000,cbind(Menarche, Total-Menarche) ~ Age2,family=binomial(logit),pfamily=dNormal(mu=mu1,Sigma=V1),data=Menarche_Model_Data)

Introduction

In this section, we will be continuing our discussion of the Bayesian model using the Menarche data that we covered in chapter 4. Recall that our main model (with the carefully specified prior) had the summary output below. The focus of this chapter will be to dig deeper and to examine model fit, residuals, outliers, and the DIC statistic (for model comparisons).

Reviewing Model Predictions

Fits and Predictions for Actual Data

## Predicting for original data

pred1=predict(glmb.out1,type="response")
pred1_m=colMeans(pred1)

plot(Menarche/Total ~ Age, data=menarche2,main="% of girls Who have had their first period")
lines(menarche2$Age, pred1_m,lty=1)

Predictions for Simulated Data

## Pull Model Info
mod_Object=glmb.out1
n=nrow(mod_Object$coefficients)
obs_size=median(menarche2$Total) ## Median of Total counts

# Generate New Data for The Independent Variable
Age_New <- seq(8, 20, 0.25)
Age2_New=Age_New-13

# olddata - Should contain all variables used in original model formula
olddata=data.frame(Age=menarche2$Age,
Menarche=menarche2$Menarche,Total=menarche2$Total,Age2=Age2)

# Should contain all variables used on RHS of model formula
newdata=data.frame(Age=Age_New,Age2=Age2_New)

## Predict and Simulate for newdata

pred_menarche=predict(mod_Object,newdata=newdata,olddata=olddata,type="response")
pred_y=matrix(0,nrow=n,ncol=length(Age_New))
for(i in 1:n){
  pred_y[i,1:length(Age_New)]=rbinom(length(Age_New),size=obs_size,
  prob=pred_menarche[i,1:length(Age_New)])
}

## Produce mean and quantile information

pred_m=colMeans(pred_menarche)
pred_y_m=colMeans(pred_y/obs_size)
quant1_m=apply(pred_menarche,2,FUN=quantile,probs=c(0.025))
quant2_m=apply(pred_menarche,2,FUN=quantile,probs=c(0.975))
quant1_m_y=apply(pred_y/obs_size,2,FUN=quantile,probs=c(0.025))
quant2_m_y=apply(pred_y/obs_size,2,FUN=quantile,probs=c(0.975))

plot(Menarche/Total ~ Age, data=menarche2,main="Percentage of girls Who have had their first period")
lines(Age_New, pred_m,lty=1)
lines(Age_New, quant1_m,lty=2)
lines(Age_New, quant2_m,lty=2)
lines(Age_New, quant1_m_y,lty=2)
lines(Age_New, quant2_m_y,lty=2)

outdata=data.frame(Age=Age_New,Age2=Age2_New,pred_m=pred_m, quant1_m, quant2_m,quant1_m_y,quant2_m_y)

print("Prior Confidence Interval at age 13 was")
cbind(m1_lower,m1_upper)

print("Prior Confidence Interval at age 15 was")
cbind(m2_lower,m2_upper)

print("Posterior Predictions between ages 13 and 15 are")
outdata[which(outdata$Age >= 13&outdata$Age <= 15 ),]

Reviewing Model Residuals

## Get Residuals for model, their means, and credible bounds
res_out=residuals(glmb.out1)
res_mean=colMeans(res_out)
sqrt(mean(res_mean^2))

res_low1=apply(res_out,2,FUN=quantile,probs=c(0.025))
res_high1=apply(res_out,2,FUN=quantile,probs=c(0.975))

## Simulate new data and get residuals for simulated data

ysim1=simulate(glmb.out1,nsim=1,seed=NULL,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))
res_high=apply(res_ysim_out1,2,FUN=quantile,probs=c(0.975))

# Plot Credible Interval bounds for Deviance Residuals

plot(res_mean~Age,ylim=c(-2.5,2.5),
  main="Credible Interval Bounds for Logit Model Deviance Residuals",
  xlab = "Age", ylab = "Dev. Residuals")
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)

# Look a the raw data
menarche2

Using the DIC Statistic

The DIC Statistic is mainly useful when comparing the predictions from several models with the same likelihood function but with different specifications (say different link-functions or different model variables). Here, we will explore alternative link functions for the Menarche data by using the Probit and Clog-Log Models and will then compare the models using the DIC Statistic and related outputs.

Running The Probit Model

## ----Probit: set up link function info and initialize prior matrices-----------

## Set up link function and initialize prior mean and Variance-Covariance matrices
bi_probit <- binomial(link="probit")
mu2<-matrix(0,nrow=nvars,ncol=1)
rownames(mu2)=c("Intercept","Age2")
colnames(mu2)=c("Prior Mean")
V2<-1*diag(nvars)
rownames(V2)=c("Intercept","Age2")
colnames(V2)=c("Intercept","Age2")


## ----Probit:set prior means----------------------------------------------------
## Prior mean for intercept is set to point estimate 
## at reference age1 (on logit scale)
mu2[1,1]=bi_probit$linkfun(m1)

## Prior mean for slope is set to difference in point estimates
## on logit scale divided by Age_Diff

mu2[2,1]=(bi_probit$linkfun(m2) -bi_probit$linkfun(m1))/Age_Diff 
print(mu2)


## ----Probit:set prior Variance Covariance matrix-------------------------------
## Implied standard deviations for point estimates on logit scale

sd_m1= (bi_probit$linkfun(m1) -bi_probit$linkfun(m1_lower))/1.96
sd_m2= (bi_probit$linkfun(m2) -bi_probit$linkfun(m2_lower))/1.96

## Implied Standard deviation for slope (using variance formula for difference between two variables)
a=(1/Age_Diff)
sd_slope=sqrt((a*sd_m1)^2+(a*sd_m2)^2-2*a*a*(sd_m1*sd_m2*m_corr))

#Cov(m1,slope)=cov(m1, a*(m2-m1)) =a*E[(m1-E[m1])((m2-m1)-E[m2-m1])]
#   =a*E[(m1-E[m1])(m2-E[m2])]- a* E[(m1-E[m1])(m1-E[m1])]
##   =a*Cov[m1,m2] - a*Var[m1]
##  =a*sd_m1*sd_m2*m_corr-a* sd_m1*sd_m1
cov_V2=a*sd_m1*sd_m2*m_corr-a* sd_m1*sd_m1

# Set covariance matrix
V2[1,1]=sd_m1^2
V2[2,2]=sd_slope^2
V2[1,2]=cov_V2
V2[2,1]=V2[1,2]
print(V2)


## ----Run Probit,results = "hide"-----------------------------------------------
Menarche_Model_Data=data.frame(Age=menarche2$Age,Total=menarche2$Total,
                               Menarche=menarche2$Menarche,Age2)

glmb.out2<-glmb(n=1000,cbind(Menarche, Total-Menarche) ~Age2,family=binomial(probit),
                pfamily=dNormal(mu=mu2,Sigma=V2),data=Menarche_Model_Data)

print(glmb.out2)
summary(glmb.out2)

Running The Clog-Log Model

## ----cloglog: set up link function info and initialize prior matrices-----------

## Set up link function and initialize prior mean and Variance-Covariance matrices
bi_cloglog <- binomial(link="cloglog")
mu3<-matrix(0,nrow=nvars,ncol=1)
rownames(mu3)=c("Intercept","Age2")
colnames(mu3)=c("Prior Mean")
V3<-1*diag(nvars)
rownames(V3)=c("Intercept","Age2")
colnames(V3)=c("Intercept","Age2")


## ----cloglog:set prior means----------------------------------------------------
## Prior mean for intercept is set to point estimate 
## at reference age1 (on logit scale)
mu3[1,1]=bi_cloglog$linkfun(m1)

## Prior mean for slope is set to difference in point estimates
## on logit scale divided by Age_Diff

mu3[2,1]=(bi_cloglog$linkfun(m2) -bi_cloglog$linkfun(m1))/Age_Diff 
print(mu3)


## ----Probit:set prior Variance Covariance matrix-------------------------------
## Implied standard deviations for point estimates on logit scale

sd_m1= (bi_cloglog$linkfun(m1) -bi_cloglog$linkfun(m1_lower))/1.96
sd_m2= (bi_cloglog$linkfun(m2) -bi_cloglog$linkfun(m2_lower))/1.96

## Implied Standard deviation for slope (using variance formula for difference between two variables)
a=(1/Age_Diff)
sd_slope=sqrt((a*sd_m1)^2+(a*sd_m2)^2-2*a*a*(sd_m1*sd_m2*m_corr))

#Cov(m1,slope)=cov(m1, a*(m2-m1)) =a*E[(m1-E[m1])((m2-m1)-E[m2-m1])]
#   =a*E[(m1-E[m1])(m2-E[m2])]- a* E[(m1-E[m1])(m1-E[m1])]
##   =a*Cov[m1,m2] - a*Var[m1]
##  =a*sd_m1*sd_m2*m_corr-a* sd_m1*sd_m1
cov_V3=a*sd_m1*sd_m2*m_corr-a* sd_m1*sd_m1

# Set covariance matrix
V3[1,1]=sd_m1^2
V3[2,2]=sd_slope^2
V3[1,2]=cov_V3
V3[2,1]=V3[1,2]
print(V3)


## ----Run cloglog,results = "hide"-----------------------------------------------
Menarche_Model_Data=data.frame(Age=menarche2$Age,Total=menarche2$Total,
                               Menarche=menarche2$Menarche,Age2)

glmb.out3<-glmb(n=1000,cbind(Menarche, Total-Menarche) ~Age2,family=binomial(cloglog),
                pfamily=dNormal(mu=mu3,Sigma=V3),data=Menarche_Model_Data)

print(glmb.out3)
summary(glmb.out3)

Running The Reverse Clog-Log Model

## ----cloglog: set up link function info and initialize prior matrices-----------

## Set up link function and initialize prior mean and Variance-Covariance matrices
bi_cloglog <- binomial(link="cloglog")
mu4<-matrix(0,nrow=nvars,ncol=1)
rownames(mu4)=c("Intercept","Age2")
colnames(mu4)=c("Prior Mean")
V4<-1*diag(nvars)
rownames(V4)=c("Intercept","Age2")
colnames(V4)=c("Intercept","Age2")


## ----cloglog:set prior means----------------------------------------------------
## Prior mean for intercept is set to point estimate 
## at reference age1 (on logit scale)
mu4[1,1]=bi_cloglog$linkfun((1-m1))

## Prior mean for slope is set to difference in point estimates
## on logit scale divided by Age_Diff

mu4[2,1]=(bi_cloglog$linkfun((1-m2))-bi_cloglog$linkfun((1-m1)))/Age_Diff 
print(mu4)


## ----Probit:set prior Variance Covariance matrix-------------------------------
## Implied standard deviations for point estimates on logit scale

sd_m1= (bi_cloglog$linkfun((1-m1_lower))-bi_cloglog$linkfun((1-m1)))/1.96
sd_m2= (bi_cloglog$linkfun((1-m2_lower))-bi_cloglog$linkfun((1-m2)))/1.96

## Implied Standard deviation for slope (using variance formula for difference between two variables)
a=(1/Age_Diff)
sd_slope=sqrt((a*sd_m1)^2+(a*sd_m2)^2-2*a*a*(sd_m1*sd_m2*m_corr))

#Cov(m1,slope)=cov(m1, a*(m2-m1)) =a*E[(m1-E[m1])((m2-m1)-E[m2-m1])]
#   =a*E[(m1-E[m1])(m2-E[m2])]- a* E[(m1-E[m1])(m1-E[m1])]
##   =a*Cov[m1,m2] - a*Var[m1]
##  =a*sd_m1*sd_m2*m_corr-a* sd_m1*sd_m1
cov_V4=a*sd_m1*sd_m2*m_corr-a* sd_m1*sd_m1

# Set covariance matrix
V4[1,1]=sd_m1^2
V4[2,2]=sd_slope^2
V4[1,2]=cov_V4
V4[2,1]=V4[1,2]
print(V4)


## ----Run cloglog,results = "hide"-----------------------------------------------
Menarche_Model_Data=data.frame(Age=menarche2$Age,Total=menarche2$Total,
                               Menarche=menarche2$Menarche,Age2)

glmb.out4<-glmb(n=1000,cbind(Total-Menarche,Menarche) ~Age2,family=binomial(cloglog),
                pfamily=dNormal(mu=mu4,Sigma=V4),data=Menarche_Model_Data)

print(glmb.out4)
summary(glmb.out4)

Comparing the DIC Statistics

DIC_Out=rbind(
  extractAIC(glmb.out1),
  extractAIC(glmb.out2),
  extractAIC(glmb.out3),
  extractAIC(glmb.out4))

rownames(DIC_Out)=c("Logit","Probit","Clog-Log","Reverse Clog-Log")
DIC_Out

glm.out1=glm(cbind(Menarche, Total-Menarche) ~Age2,family=binomial(logit),data=Menarche_Model_Data)
glm.out2=glm(cbind(Menarche, Total-Menarche) ~Age2,family=binomial(probit),data=Menarche_Model_Data)
glm.out3=glm(cbind(Menarche, Total-Menarche) ~Age2,family=binomial(cloglog),data=Menarche_Model_Data)
glm.out4=glm(cbind(Total-Menarche,Menarche ) ~Age2,family=binomial(cloglog),data=Menarche_Model_Data)

AIC_Out=rbind(extractAIC(glm.out1),
extractAIC(glm.out2),
extractAIC(glm.out3),
extractAIC(glm.out4))

rownames(AIC_Out)=c("Logit","Probit","Clog-Log","Reverse Clog-Log")
colnames(AIC_Out)=c("edf","AIC")
AIC_Out

Plotting Predictions and Residuals For the Additional Models

Probit Model

Predictions

## Predict and Simulate for newdata

mod_Object=glmb.out2
pred1=predict(glmb.out2,type="response")
pred1_m=colMeans(pred1)
n=nrow(mod_Object$coefficients)
pred_menarche=predict(mod_Object,newdata=newdata,olddata=olddata,type="response")
pred_y=matrix(0,nrow=n,ncol=length(Age_New))
for(i in 1:n){
  pred_y[i,1:length(Age_New)]=rbinom(length(Age_New),size=obs_size,
  prob=pred_menarche[i,1:length(Age_New)])
}

## Produce mean and quantile information

pred_m=colMeans(pred_menarche)
pred_y_m=colMeans(pred_y/obs_size)
quant1_m=apply(pred_menarche,2,FUN=quantile,probs=c(0.025))
quant2_m=apply(pred_menarche,2,FUN=quantile,probs=c(0.975))
quant1_m_y=apply(pred_y/obs_size,2,FUN=quantile,probs=c(0.025))
quant2_m_y=apply(pred_y/obs_size,2,FUN=quantile,probs=c(0.975))

plot(Menarche/Total ~ Age, data=menarche2,main="Percentage of girls Who have had their first period")
lines(Age_New, pred_m,lty=1)
lines(Age_New, quant1_m,lty=2)
lines(Age_New, quant2_m,lty=2)
lines(Age_New, quant1_m_y,lty=2)
lines(Age_New, quant2_m_y,lty=2)

Residuals

## Get Residuals for model, their means, and credible bounds
res_out=residuals(glmb.out2)
res_mean=colMeans(res_out)
sqrt(mean(res_mean^2))
res_low1=apply(res_out,2,FUN=quantile,probs=c(0.025))
res_high1=apply(res_out,2,FUN=quantile,probs=c(0.975))

## Simulate new data and get residuals for simulated data

ysim1=simulate(glmb.out2,nsim=1,seed=NULL,pred=pred1,family="binomial",
                        prior.weights=weights(glmb.out2))
res_ysim_out1=residuals(glmb.out2,ysim=ysim1)
res_low=apply(res_ysim_out1,2,FUN=quantile,probs=c(0.025))
res_high=apply(res_ysim_out1,2,FUN=quantile,probs=c(0.975))

# Plot Credible Interval bounds for Deviance Residuals

plot(res_mean~Age,ylim=c(-2.5,2.5),
  main="Credible Interval Bounds for Probit Model Deviance Residuals",
  xlab = "Age", ylab = "Dev. Residuals")
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)

clog-log Model

Predictions

## Predict and Simulate for newdata
mod_Object=glmb.out3
pred1=predict(glmb.out3,type="response")
pred1_m=colMeans(pred1)
n=nrow(mod_Object$coefficients)
pred_menarche=predict(mod_Object,newdata=newdata,olddata=olddata,type="response")
pred_y=matrix(0,nrow=n,ncol=length(Age_New))
for(i in 1:n){
  pred_y[i,1:length(Age_New)]=rbinom(length(Age_New),size=obs_size,
  prob=pred_menarche[i,1:length(Age_New)])
}

## Produce mean and quantile information

pred_m=colMeans(pred_menarche)
pred_y_m=colMeans(pred_y/obs_size)
quant1_m=apply(pred_menarche,2,FUN=quantile,probs=c(0.025))
quant2_m=apply(pred_menarche,2,FUN=quantile,probs=c(0.975))
quant1_m_y=apply(pred_y/obs_size,2,FUN=quantile,probs=c(0.025))
quant2_m_y=apply(pred_y/obs_size,2,FUN=quantile,probs=c(0.975))

plot(Menarche/Total ~ Age, data=menarche2,main="Percentage of girls Who have had their first period")
lines(Age_New, pred_m,lty=1)
lines(Age_New, quant1_m,lty=2)
lines(Age_New, quant2_m,lty=2)
lines(Age_New, quant1_m_y,lty=2)
lines(Age_New, quant2_m_y,lty=2)

Residuals

## Get Residuals for model, their means, and credible bounds
res_out=residuals(glmb.out3)
res_mean=colMeans(res_out)
sqrt(mean(res_mean^2))
res_low1=apply(res_out,2,FUN=quantile,probs=c(0.025))
res_high1=apply(res_out,2,FUN=quantile,probs=c(0.975))

## Simulate new data and get residuals for simulated data

ysim1=simulate(glmb.out3,nsim=1,seed=NULL,pred=pred1,family="binomial",
                        prior.weights=weights(glmb.out3))
res_ysim_out1=residuals(glmb.out3,ysim=ysim1)
res_low=apply(res_ysim_out1,2,FUN=quantile,probs=c(0.025))
res_high=apply(res_ysim_out1,2,FUN=quantile,probs=c(0.975))

# Plot Credible Interval bounds for Deviance Residuals

plot(res_mean~Age,ylim=c(-4.0,4.0),
  main="Credible Interval Bounds for cloglog Model Deviance Residuals",
  xlab = "Age", ylab = "Dev. Residuals")
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)

Reverse clog-log Model

Predictions

## Predict and Simulate for newdata
mod_Object=glmb.out4
pred1=predict(glmb.out4,type="response")
pred1_m=colMeans(pred1)
n=nrow(mod_Object$coefficients)
pred_menarche=predict(mod_Object,newdata=newdata,olddata=olddata,type="response")
pred_y=matrix(0,nrow=n,ncol=length(Age_New))
for(i in 1:n){
  pred_y[i,1:length(Age_New)]=rbinom(length(Age_New),size=obs_size,
  prob=pred_menarche[i,1:length(Age_New)])
}

## Produce mean and quantile information

pred_m=colMeans(pred_menarche)
pred_y_m=colMeans(pred_y/obs_size)
quant1_m=apply(pred_menarche,2,FUN=quantile,probs=c(0.025))
quant2_m=apply(pred_menarche,2,FUN=quantile,probs=c(0.975))
quant1_m_y=apply(pred_y/obs_size,2,FUN=quantile,probs=c(0.025))
quant2_m_y=apply(pred_y/obs_size,2,FUN=quantile,probs=c(0.975))

plot(Menarche/Total ~ Age, data=menarche2,main="Percentage of girls Who have had their first period")
lines(Age_New, 1-pred_m,lty=1)
lines(Age_New, 1-quant1_m,lty=2)
lines(Age_New, 1-quant2_m,lty=2)
lines(Age_New, 1-quant1_m_y,lty=2)
lines(Age_New, 1-quant2_m_y,lty=2)

Residuals

## Get Residuals for model, their means, and credible bounds
res_out=residuals(glmb.out4)
res_mean=colMeans(res_out)
sqrt(mean(res_mean^2))
res_low1=apply(res_out,2,FUN=quantile,probs=c(0.025))
res_high1=apply(res_out,2,FUN=quantile,probs=c(0.975))

## Simulate new data and get residuals for simulated data

ysim1=simulate(glmb.out4,nsim=1,seed=NULL,pred=pred1,family="binomial",
                        prior.weights=weights(glmb.out4))
res_ysim_out1=residuals(glmb.out4,ysim=ysim1)
res_low=apply(res_ysim_out1,2,FUN=quantile,probs=c(0.025))
res_high=apply(res_ysim_out1,2,FUN=quantile,probs=c(0.975))

# Plot Credible Interval bounds for Deviance Residuals

plot(-res_mean~Age,ylim=c(-3.0,3.0),
  main="Credible Interval Bounds - Rev. clog-log Model Deviance Residuals",
  xlab = "Age", ylab = "Dev. Residuals")
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)

References

knygren/glmbayes documentation built on Sept. 4, 2020, 4:39 p.m.

rdrr.io home R language documentation Run R code online

CRAN packages Bioconductor packages R-Forge packages GitHub packages

Note that we can't provide technical support on individual packages. You should contact the package authors for that.

knygren/glmbayes
Bayesian Generalized Linear Models (iid Samples)

In knygren/glmbayes: Bayesian Generalized Linear Models (iid Samples)

Introduction

Reviewing Model Predictions

Fits and Predictions for Actual Data

Predictions for Simulated Data

Reviewing Model Residuals

Using the DIC Statistic

Running The Probit Model

Running The Clog-Log Model

Running The Reverse Clog-Log Model

Comparing the DIC Statistics

Plotting Predictions and Residuals For the Additional Models

Probit Model

clog-log Model

Reverse clog-log Model

References

R Package Documentation

Browse R Packages

We want your feedback!

knygren/glmbayes Bayesian Generalized Linear Models (iid Samples)

In knygren/glmbayes: Bayesian Generalized Linear Models (iid Samples)

Introduction

Reviewing Model Predictions

Fits and Predictions for Actual Data

Predictions for Simulated Data

Reviewing Model Residuals

Using the DIC Statistic

Running The Probit Model

Running The Clog-Log Model

Running The Reverse Clog-Log Model

Comparing the DIC Statistics

Plotting Predictions and Residuals For the Additional Models

Probit Model

clog-log Model

Reverse clog-log Model

References

R Package Documentation

Browse R Packages

We want your feedback!

knygren/glmbayes
Bayesian Generalized Linear Models (iid Samples)