2
$\begingroup$

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?

Adding a made up df 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

$\endgroup$
0
$\begingroup$

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.

$\endgroup$
  • $\begingroup$ Thanks for the repsonse. I re-ran a random made up df by also adding a separate slope for each subject - still have the same issue. I agree with the factor/integer but rstanarm has no problem dealing with them as far as I know. $\endgroup$ – Dilsher Singh Dhillon Apr 4 at 22:06

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.