# Out-of-sample predictions for mixed model are the same as naive model (ignoring the random effects)

I have a dataset that consists of subjects coming into the clinic (for treatment of another disease) and they are screened for Tuberclosis (as they are a high risk population). Every time they are screened, they are either screen positive or negative (1 or 0) and then for each visit, there is an assessment of TB disease. My colloborator wants to assess the utility of this screening test.

The problem I’m running into goes like this - Since each subject is screened multiple times, I incluide a random intercept for each subject and a fixed effect for the screen result and run a mixed model logistic regression

library(rstanarm)
m1<-stan_glmer(disease_binary ~ screen_binary + (1|chart_number_2),data=df,family="binomial")


There’s about 11000 subjects with a total observations of 55000. This took me about 30 minutes running all 4 chains in parrallel.

To assess the screening test, I would like to calculate an AUC which is done using below -

screen_preds <- posterior_linpred(m1,newdata=df,transform=TRUE) %>% tbl_df

prob_data <- screen_preds %>%
map_dfr(~median(.x))%>%
gather(.,key=id,value="prob")%>%
mutate(id=as.numeric(id))%>%
inner_join(.,dta_bayes,by="id") .  ### I take the median of the predicted probabilities and join them with the original dataframe which has the 'truth' (dta_bayes is the original df)

## Using the probabiilites and the true responses, I calculate the AUROC
library(pROC)
roc_random<-roc(prob_data$$disease_binary,prob_data$$prob,ci=TRUE,boot.n=100,print.auc=TRUE,plot=TRUE)
auc(roc_random)


I get an AUC of 0.998. This is quite surprising and bizarre to me.

Instead, I took a different approach. I split the data into training and testing - randomly sampled 80% of the subjects and took all of their data (so 80% of the groups were sampled), and left 20% as a hold out test.

I repeated the above, instead fitting the model on the 80% data.

m_train <- stan_glmer(disease_binary~screen_binary+ (1|chart_number_2),data=train,family="binomial",cores=4)

## Using the m_train model to now predict on the hold out test set
test_preds <- posterior_linpred(m_train,newdata = test,transform = TRUE)%>%tbl_df  ## test set probabilities

## Using these test set probabilities -
## take the median of the posterior predictions
median_preds<-test_preds%>%
map_dfr(.,~median(.x))

### Map the probabilities to the test df
test_prob<-median_preds%>%
gather(.,key=id,value="prob")%>%
mutate(id=as.numeric(id))%>%
inner_join(.,test,by="id")       ### The test set probabilites are now mapped to the test data (with the true responses)


AUC

test_roc<-roc(test_prob$$disease_binary,test_prob$$prob,ci=TRUE,boot.n=100,print.auc=TRUE,plot=TRUE)
auc(roc_random)


The result is 0.603 (0.544-0.662) - quite a drastic drop from 0.998.

Alright, so we maybe overfit our model with so many random effects - how about we do the same train-test split but this time, ignore the correlation and treat each individual visit as an independent response?

train_naive <- glm(disease_binary~screen_binary,data=train,family="binomial")
prob_naive <- predict(train_naive,newdata = test,type="response")
prob_naive <- as.data.frame(prob_naive)

test_prob <- prob_naive%>%
tbl_df %>%
mutate(id=seq(1,nrow(.),by=1)) %>%
inner_join(.,test_prob,by="id")

test_roc_naive<-roc(test_prob$$disease_binary,test_prob$$prob_naive,ci=TRUE,boot.n=100,print.auc=TRUE,plot=TRUE)

auc(test_roc_naive)


I get the same exact AUC - 0.603 (0.544-0.662).

The naive and the random effect model produce exactly the same estimate, which makes me wonder that the random effect model is essentially using the same coefficient(s) as the naive model.

What am I missing here?

screen_binary<-rbinom(1000,1,0.5)
diseas_binary<-rbinom(1000,1,0.05)
chart_number_2<-seq(1,150)

df<-as.data.frame(cbind(screen_binary,diseas_binary,chart_number_2))%>%
tbl_df%>%
arrange(.,chart_number_2)

Reproducible example - Added a separate slope for screen_binary too here
screen_binary<-rbinom(1000,1,0.5)
disease_binary<-rbinom(1000,1,0.05)
chart_number_2<-seq(1,150)

df<-as.data.frame(cbind(screen_binary,disease_binary,chart_number_2))%>%
tbl_df%>%
arrange(.,chart_number_2)%>%
mutate(id=seq(1,nrow(.),by=1))

library(rstanarm)
m1<-stan_glmer(disease_binary ~ screen_binary + (1 + screen_binary | chart_number_2),data=df,family="binomial")
screen_preds <- posterior_linpred(m1,newdata=df,transform=TRUE) %>% tbl_df

prob_data <- screen_preds %>%
map_dfr(~median(.x))%>%
gather(.,key=id,value="prob")%>%
mutate(id=as.numeric(id))%>%
inner_join(.,df,by="id")

## Using the probabilites and the true responses, I calculate the AUROC
library(pROC)
roc_random<-roc(prob_data$$disease_binary,prob_data$$prob,ci=TRUE,boot.n=100,print.auc=TRUE,plot=TRUE)
auc(roc_random)


I get an AUC of 0.86


## train and test
# split the observations into training and test
subjs<-df%>%
distinct(.,chart_number_2)%>%
sample_frac(.,0.8)

train<-df%>%
inner_join(.,subjs,by="chart_number_2")
test<-df%>%
anti_join(.,train,by="chart_number_2")%>%
mutate(id=seq(1,nrow(.),by=1))

## Repeat trainign on train and instead of predicting on train, predict on test
library(rstanarm)
m1<-stan_glmer(disease_binary ~ screen_binary + (1 + screen_binary | chart_number_2),data=train,family="binomial")
screen_preds <- posterior_linpred(m1,newdata=test,transform=TRUE) %>% tbl_df

prob_data <- screen_preds %>%
map_dfr(~median(.x))%>%
gather(.,key=id,value="prob")%>%
mutate(id=as.numeric(id))%>%
inner_join(.,test,by="id")

## Using the probabilites and the true responses, I calculate the AUROC
library(pROC)
roc_random<-roc(prob_data$$disease_binary,prob_data$$prob,ci=TRUE,print.auc=TRUE,plot=TRUE)
auc(roc_random)


I get an AUC of 0.503

The screen_binary is a slope — there’s an implied 1 + in your formula, which is the population intercept and your screen_binary is a slope. So perhaps you should also have a slope (1 + screen_binary | chart_number_2) in your random effect?
Also have you made sure your screen_binary and disease_binary are factors rather than an integer? Not sure about rstanarm, but integers instead of factors can wreak havoc with some packages/functions.