Reviewing Model Predictions
We start with a look a model predictions, both in and out of sample. The predict function can be used for both of these purposes as we illustrate below.
Predictions for Actual Data
Generating in-sample predictions requires just a simple call to the predict function in order to extract the required information of the glmb class object (in this case we choose to extract the predictions for the “response). Once extracted, the randomized draws for the predictions can be analyzed further. Here we choose to so for my just looking at the mean predictions (we look at other outputs for the out of sample prediction below) and we then plot those for comparison to the underlying data.
As we can see from this plot, the predictions in general seems visually to match the underlying data relatively well although there may be some extent to which it overstates the percentage of girls who have had their first period for the lower age range. It is also hard to judge from this plot whether some of the points are outliers or whether they fall within a range where we would expect the observations to fall. We explore both of those issues further later in this chapter.
## Predicting for original data
pred1=predict(glmb.out1,type="response")
pred1_m=colMeans(pred1)
plot(Menarche/Total ~ Age, data=menarche,main="% of girls Who have had their first period")
lines(menarche$Age, pred1_m,lty=1)
Predictions for “New” Data
In order to understand the remaining uncertainty tied to our model and predictions for “new” observations, we use simulated value for the coefficients in order to simulate new data over a range of different ages so that we can provide “credible sets” around the model. We use the median sample size for the underlying observations in the the model. Let’s us briefly touch on each of the steps below.
In order to perform the simulation, we first set up a new dataset for the independent variables by creating a new age variable and setting the number of observations to simulate. We store the new independent variable in a a new data frame with the same variable names as in the original dataset.
## Pull Model Info
mod_Object=glmb.out1
n=nrow(mod_Object$coefficients)
obs_size=median(menarche$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=menarche$Age,
Menarche=menarche$Menarche,Total=menarche$Total,Age2=Age2)
# Should contain all variables used on RHS of model formula
newdata=data.frame(Age=Age_New,Age2=Age2_New)In order to perform the simulation, we first pass the output from the glmb, the newdata, the olddata, and the type to the predict function. We then set up a matrix into which we store the newly simulated values for the dependent variable and simulate from rbinom function using the predicted response values from the response function.
## 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)])
}The output is summarized by producing both means and credible set information for proportions based on the simulated coefficients and the simulated data (the intervals for the latter will be wider). The latter interval is more appropriate when examining if the actual data appears to correspond to outliers.
## 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),na.rm=TRUE)
quant2_m=apply(pred_menarche,2,FUN=quantile,probs=c(0.975),na.rm=TRUE)
quant1_m_y=apply(pred_y/obs_size,2,FUN=quantile,probs=c(0.025),na.rm=TRUE)
quant2_m_y=apply(pred_y/obs_size,2,FUN=quantile,probs=c(0.975),na.rm=TRUE)We now plot most of these summarized measures and the original data together in order to get a visualization of the various predictions and data. This plot suggests that the individual observations in the original data generally tends to fall between the lower and upper bounds credible sets based on the simulated data although it is hard to tell at lower end of the age scale based on this chart. We will explore the latter issue a bit more later in this chapter.
plot(Menarche/Total ~ Age, data=menarche,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)
We close this section with a comparison of our prior and posterior credible intervals. We store the related information for our posterior simulated credible intervals in a data frame and then look at these in comparison o the prior credible intervals. What we see is that the prior credible interval at age 13 is wider than the posterior credible interval and in fact contains the posterior interval. At age 15, the two posterior credible set is still contained within the prior credible set but skews heavily towards the upper end of the prior credible set.
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 credible set at age 13 was")
#> [1] "Prior credible set at age 13 was"
cbind(m1_lower,m1_upper)
#> m1_lower m1_upper
#> [1,] 0.3 0.7
print("Prior credible set at age 15 was")
#> [1] "Prior credible set at age 15 was"
cbind(m2_lower,m2_upper)
#> m2_lower m2_upper
#> [1,] 0.7 0.972
print("Posterior Predictions between ages 13 and 15 are")
#> [1] "Posterior Predictions between ages 13 and 15 are"
outdata[which(outdata$Age >= 13&outdata$Age <= 15 ),]
#> Age Age2 pred_m quant1_m quant2_m quant1_m_y quant2_m_y
#> 21 13.00 0.00 0.4950818 0.4631213 0.5269138 0.3997619 0.6000000
#> 22 13.25 0.25 0.5949412 0.5650843 0.6267817 0.5045238 0.6857143
#> 23 13.50 0.50 0.6875006 0.6593774 0.7176085 0.5904762 0.7809524
#> 24 13.75 0.75 0.7671624 0.7417636 0.7939361 0.6857143 0.8476190
#> 25 14.00 1.00 0.8314705 0.8096203 0.8547607 0.7523810 0.9047619
#> 26 14.25 1.25 0.8807580 0.8620487 0.9001824 0.8095238 0.9428571
#> 27 14.50 1.50 0.9170546 0.9016905 0.9321522 0.8571429 0.9714286
#> 28 14.75 1.75 0.9430077 0.9306520 0.9547214 0.8952381 0.9809524
#> 29 15.00 2.00 0.9611774 0.9516589 0.9700030 0.9142857 0.9904762
Residuals
Residuals
Residuals