I have data from a study in which male and female doctors each read 6 vignettes describing patient situations. The vignettes vary the age of the patient and also the level of pain the patient is in.

Each doctor reads the same 6 vignettes, and then makes a decision about whether to proceed with surgery or not.

Here is the table of counts

Doctor/Patient data

I've also made the raw data available here.

Ultimately what I'm interested in is the impact on the decision to go to surgery of

  1. Doctor's sex
  2. Patient pain level
  3. Patient age

It seems to me that there are some issues with the manipulation that was used, inasmuch as in many of the vignettes it's overwhelming that the doctor calls for surgery. For example, even in the mild pain vignette basically all doctors want surgery all the time, irrespective of whether the patient is old or young.

Under ideal circumstances I had hoped to say something like "There's an association between pain level and proceeding with surgery, but only if the doctor is female" or "there's an association between patient age and the doctor declining surgery, but only if the doctor is male".

It seems unlikely that such a conclusion can be reached from the observed data table, but I would like to do the analysis that could have justified such a conclusion if the results were different.


I think it would be perfectly reasonable to use a GLMM/mixed logistic regression model for this problem, i.e.

$$ \begin{split} y & \sim \textrm{Bernoulli}(\mu) \\ \mu & = 1/(1+\exp(-\eta)) \\ \eta & = \beta_0 + \beta_1 \cdot \textrm{docgender} + \beta_2 \cdot \textrm{patage} + \beta_3 \cdot \textrm{patpain}_1 + \beta_4 \cdot \textrm{patpain}_2 + \epsilon_{\textrm{doc}} \\ \epsilon_{\textrm{doc}} & \sim N(0,\sigma^2_{\textrm{doc}}) \end{split} $$ where the variables in the $\eta$ equation are dummy variables (we need two separate dummy variables for patient pain, e.g. $\textrm{patpain}_1$="mild vs. moderate pain" and $\textrm{patpain}_2$="mild vs. severe pain" since the predictor variable has three possible levels)

or in lme4 notation y ~ docgender + patage + patpain + (1|docID)

Note that the 'raw data' you gave isn't quite raw enough to solve this problem; aggregating across doctors to get the summary table loses the among-doctor variation, e.g. some doctors may be more likely to prescribe surgery across the board. In principle since every doctor reads all the vignettes you could estimate the among-doctor variance in response to the vignettes (not just in the average propensity to choose surgery), via

y ~ docgender + patage + patpain + (1+docgender+patage+patpain|docID)

but with 270 doctors and a Bernoulli response I think it would be a little optimistic to think you could estimate the 5x5 variance-covariance matrix implied by this formula (this is what's recommended by Barr et al in their "Keep it maximal" paper ...). Maybe worth a try.

Oops, now I notice that you're also interested in interactions. It would be reasonable to try

y ~ docgender*patage*patpain + (1|docID)

According to Frank Harrell's (Regression Modeling Strategies) rule of thumb, you should try not to have more than one parameter per 10-20 'data points', where in this case 'data points' corresponds to the minimum of the number of yeses and nos in your data set - there seem to be approx. 400 no values, so 20 parameters are OK, this interaction model has (2x2x3+1)=13 parameters, so seems reasonable.

  • $\begingroup$ Thanks for this great response. Just to clarify, is y ~ docgender*patage*patpain + (1|docID) what I should run in R if I have access to the full raw data, i.e. the data not aggregated across doctors? I'm pretty sure I can get access to this data. $\endgroup$ Sep 15 '16 at 23:54

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.