Can I ask how to get variability around a prediction after fitting a random intercept model in stan_glmer please? I have fit the following model where my event is 0/1 and x1 and x2 are continuous covariates:
model <- stan_glmer(binary event ~ x1 + x2 +(1 | group), family="binomial"
I attach an example – but essentially my goal is to plot the predicted probability for a new individual (who doesn’t belong to an observed group) and therefore I marginalise over the random effects distribution, for different values of x2 whilst holding x1 at its mean value, say. I can plot this but without variability. I can bootstrap my mcmc sample of 4000 draws but I’m not sure this is adding the right amount of variability. So my question is how/can I get variability on my prediction plot or is it just not possible/needed?
Any help would be greatly appreciated.
#generate some data as an example: 4 groups of people, 20 people in each group each with covariates x1 and x2
groupn <- 15 # sample size for new study
group1 <- data.frame(group=“1”,
id=1:15, # id’s
x1=seq(40, 70, length.out = 15),
x2=seq(0, 10, length.out = 15),
Event=c(0,1,0,0,1, 0,0,1,0,0, 0,0,0,0,0))
group2 <- data.frame(group=“2”,
id=16:30, # id’s
x1=seq(30, 60, length.out = 15),
x2=seq(5, 15, length.out = 15),
Event=c(1,0,0,0,0, 0,0,0,0,0, 0,0,0,0,0))
group3 <- data.frame(group=“3”,
id=31:45, # id’s
x1=seq(20, 40, length.out = 15),
x2=seq(0, 5, length.out = 15),
Event=c(0,0,0,0,1, 0,0,0,0,0, 0,0,0,0,0))
group4 <- data.frame(group=“4”,
id=46:60, # id’s
x1=seq(30, 50, length.out = 15),
x2=seq(2, 9, length.out = 15),
Event=c(0,1,0,0,1, 0,0,0,0,0, 0,0,0,0,0))
dat<-rbind(group1, group2, group3, group4)
set.seed(123)
model <- stan_glmer(Event ~ x1 + x2 +(1 | group), data=dat, family=“binomial”,
prior_intercept=cauchy(0, 10),
prior=cauchy(0, 2.5),
prior_covariance=decov(regularization=1, concentration=1, shape=1, scale=1),
cores=4)
loomod <- loo(model)
print(model, digits=3)
print(loomod)
#7/60=0.116
#Estimates:
#Median MAD_SD
#(Intercept) -0.874 2.067
#x1 0.013 0.062
#x2 -0.403 0.232
#Error terms:
#Groups Name Std.Dev.
#group (Intercept) 1.111
#Num. levels: group 4
#Sample avg. posterior predictive
#distribution of y (X = xbar):
#Median MAD_SD
#mean_PPD 0.117 0.049
summary(dat$x1)
summary(dat$x2)
new_n<-1
out<-prop<-NULL
x<-seq(0,15,0.01)
new_n<-1
for(i in (x)){
nd <- data.frame(group=“NEW group”, # to identify prediction over new study
id=1:new_n, # patient id’s
Event=rep(0, new_n)) # column is needed for posterior_predict,
#actual values are ignored
#vary x2 while keeping x1 fixed at the mean
nd$x2 <- c(i)
nd$x1<- c(42.5)
#Number of simulations from the predictive distribution
draws<-4000
ppred1 <- posterior_predict(model, newdata=nd, re.form=~(1|group), draws=draws , seed=123)
##Look at the individual probabilities - columns
sumcol <- data.frame(iter=1:1, count=apply(ppred1, 2, sum) )
totalcol<- filter(sumcol)%>% mutate(prop=count/draws)
prop <- cbind(prop, totalcol$prop)
out <- cbind(out, as.vector(ppred1))
}
probc<-t(prop)
x2var<-cbind(x, probc)
x2var<-data.frame(x2var)
#plot
ggplot(x2var, aes(x=x, y=V2)) + geom_point(size=1, color=‘red’)+
scale_y_continuous(breaks = round(seq(0, 1.2, by = 0.2),1)) +
scale_x_continuous(breaks = round(seq(0, 15, by = 1),1))+
labs(x=“X2”, y=“Probability of an event”)+
coord_cartesian(ylim=c(0, 1), xlim=c(0, 15))