Generalized linear mixed model for assessing diagnostic accuracy of doctors

Question

I am trying to determine how a doctor's expertise (either expert:1 or not an expert:0) affects their diagnostic accuracy of medical images. Additionally, I'd like to examine how the hospital at which the images were taken affects the diagnostic accuracy, both for expert and non-expert doctors.

I have thousands of medical images coming from 30 hospitals around the world, each of which has been diagnosed by 10-20 doctors, some of which are experts and some of which are not. The doctors assessed random images not coming from their own hospital. Here is an example of how the data looks:

Image	Doctor	Expertise	Hospital	Correct Diagnosis
1	A	1	Italy	1
1	B	0	Italy	0
1	C	1	Italy	1
2	D	0	USA	0
2	A	1	USA	1
3	E	0	France	0

I am wondering if it suitable to use a generalized linear mixed model (GLMM) to understand the diagnostic accuracy of doctors based on their expertise and where the images were taken? Also, is it appropriate to model the hospital and expertise as fixed effects and the Doctor as a random affect, to account for the correlation caused by the same doctor assessing different images?

Here is an example of code I have used in R:

model <- glmer(
Correct Diagnosis ~ Expertise + Hospital + (1 | Doctor),
data = df, 
family = binomial, 
control = glmerControl(optimizer="bobyqa")
)

Thanks for the help!

Answer 1

Why not model hospitals as random effects? Your doctors are clearly nested within hospitals, and you have enough hospitals (30) to use them as random effects. Referring to Page 7 of the original paper from lme4's creator, you could instead use a model like the one below:

model <- glmer(
Correct Diagnosis ~ Expertise + (1 | Hospital / Doctor),
data = df, 
family = binomial, 
control = glmerControl(optimizer="bobyqa")
)

Here you are now estimating the random intercepts of doctors as well as the random intercepts for hospitals. This allows you to see how correct diagnoses vary among doctors and hospitals. Keep in mind though that the more complex your random effects become, the more likely it will not converge if there isn't actually any variation among the random effects. If this happens, you can keep the model you already have to simplify things and allow the model to converge.