3
$\begingroup$

I'm trying to convey some findings in which a score from 1-10 seems to predict disease status (binary).

I predict yhat and plot yhat with a quadratic fit against my predictor (score).

It looks accurate but when I add a confidence interval to the quadratic fit, it's VERY narrow. Too narrow for me to believe it.

Does a predicted plot deflate the confidence interval of the original data? If so, does anyone have a suggestion on how I convey my data in a similar format (i.e. for every increase in score, percentage of success increases by y) with a confidence interval?

First plot:

This plot is obtained by entering twoway qfitci effect score where effect is a binary variable denoting whether or not the patient had the desired effect, 1 being effect, and score being a nominal/continuous variable where 1 is the lowest and 10 is the highest score, with the hypothesis that a higher score increases probability of effect. qfitci is a quadratic fitting plot with CI in gray.

CI of original data plot twoway qfitci effect score:

enter image description here

2nd plot:

This plot is obtained by running a logistic regression model: logistic effect score age gender in stata. This command returns OR's as opposed to the logit command which returns coefficients in e. This model is then predicted using predict yhat in stata which creates a new variable yhat with the predicted probabilities of the model.

DATA (CSV): https://gofile.io/?c=sxNnuM or: https://easyupload.io/fnp6r8

CI of prediction plot twoway qfitci yhat score:

enter image description here

Improve this question
$\endgroup$
6
  • $\begingroup$ There's not enough information in what you've provided to say much at all. Exactly how were the CI bands for the original plot and for the prediction plot determined? Note also that the y axis for neither plot is well-defined; are they probabilities, log-odds, or what? It also might help to have both y axes on the same magnification scale; presently one covers about 3 times the range of the other. $\endgroup$
    – EdM
    Commented Jan 24, 2020 at 20:07
  • $\begingroup$ I had hoped I had defined them by saying one is a binary outcome and raw data, so it's a continuous predictor (score) plotted against a binary variable. The other is the prediction of the model (predict "variableName" in Stata) so it is in percentages. The model used is a multiple logistic regression using OR's and age+gender as control variables. $\endgroup$
    – Paze
    Commented Jan 24, 2020 at 20:09
  • $\begingroup$ The top graph might be based on a model of a binary outcome, but it has a continuous y-axis scale suggesting that the y axis is either log-odds (consistent with the plot) or a probability (except that the CI extend beyond 0 and 1). Is the second graph perhaps a probability scale instead? (Not likely to be percentages with values below 1 unless this is a very rare event.) Please provide details on how the plots and confidence intervals were obtained. $\endgroup$
    – EdM
    Commented Jan 24, 2020 at 20:16
  • 1
    $\begingroup$ Data added. I have included the yhat prediction variable in the data, but it could also be generated in stata by first regressing: "regress effect score age gender" and then "predict yhat", I don't know the syntaxes for R or SPSS, sorry. $\endgroup$
    – Paze
    Commented Jan 24, 2020 at 20:49
  • 1
    $\begingroup$ Hm I get the same problem with the excel file. The information is as follows: A = age B = effect C = score D = predicted model E = gender I'm sorry for the confusion, the CSV file seems to work better, I'll reupload it another place: easyupload.io/fnp6r8 Does this work? There is some missing data, yes, it's normal in this case. $\endgroup$
    – Paze
    Commented Jan 24, 2020 at 21:24

1 Answer 1

Reset to default
5
$\begingroup$

I'm afraid I don't know Stata (at all...), so I'm making some guesses here.

Your real data differ somehow, or Stata is doing something I can't divine, because you have yhats for two patients with missing responses. Using complete case analysis, I replicated the logistic regression model in R. Your yhats are predicted probabilities from a standard logistic regression model with additive (on the linear scale) effects of score, age, and gender. Your top plot seems to treat the 0/1 effect data as a response and fits a linear (OLS) regression model with a quadratic on score, and uses normal theory to add a confidence band. Your second plot seems to treat those predicted probabilities as though they were the raw response data and does the same thing. Neither of these is correct.

You presumably want to marginalize over the other variables somehow and plot the predicted values with their SE's from the logistic regression model. There are various ways to do this, e.g., you could use 'least squares means'. A simple way to get something is to solve the model equation at specified values of your variables. For example, you could plot two lines for the two sexes at the mean of age for the different values of score and plot that. Here's an example, coded in R (I don't have Stata):

m = glm(effect~score+age+gender, d, family=binomial)
round(coef(summary(m)), digits=3)
#             Estimate Std. Error z value Pr(>|z|)
# (Intercept)   -2.322      3.372  -0.689    0.491
# score          0.485      0.154   3.144    0.002
# age            0.007      0.045   0.157    0.875
# gender         0.249      0.782   0.319    0.750

pred.mat.0 = predict(m, newdata=data.frame(score=0:10, age=75, gender=0), 
                     type="link", se.fit=T)
pred.mat.1 = predict(m, newdata=data.frame(score=0:10, age=75, gender=1),
                     type="link", se.fit=T)
lo2p = function(lo){ exp(lo)/(1+exp(lo)) }
pred.mat.0.probs = data.frame(yhat=lo2p(pred.mat.0$fit),
                               low=lo2p(pred.mat.0$fit - 2*pred.mat.0$se.fit),
                                up=lo2p(pred.mat.0$fit + 2*pred.mat.0$se.fit) )
pred.mat.1.probs = data.frame(yhat=lo2p(pred.mat.1$fit),
                               low=lo2p(pred.mat.1$fit - 2*pred.mat.1$se.fit),
                                up=lo2p(pred.mat.1$fit + 2*pred.mat.1$se.fit) )

set.seed(1)
plot(effect~jitter(score, amount=.25), d, xlim=c(0,10), 
     col=rgb(.5,.5,.5,alpha=.5), pch=16)
lines(0:10, pred.mat.0.probs$yhat, col="blue3", lwd=2)
lines(0:10, pred.mat.0.probs$low,  col="blue1")
lines(0:10, pred.mat.0.probs$up,   col="blue1")
lines(0:10, pred.mat.1.probs$yhat, col="red3", lwd=2)
lines(0:10, pred.mat.1.probs$low,  col="red1")
lines(0:10, pred.mat.1.probs$up,   col="red1")

enter image description here

Improve this answer
$\endgroup$
6
  • $\begingroup$ Thank you Gung that explains it. I should have known that there is no easy way to graphically interpret a multiple model (one of my pet peeves with statistics is how important multiple models are in my field and how difficult it is to portray this in a graphic). As you suggested I can graph a few different scenarios and display those. Perhaps what I am looking for is simply a ROC curve or confusion matrix? It would portray the model fit in a way and the implications/coefficients of the variable of interest? $\endgroup$
    – Paze
    Commented Jan 25, 2020 at 0:57
  • $\begingroup$ @Paze, I don't know Stata. It may well have an easy method for you. It's just that the methods you used don't do that. Gender isn't doing much in your model, you could probably just pick one, if you wanted a single line. What is the 'right' way to present your model depends on what you want to show about it. I doubt a ROC curve or a confusion matrix would be best for you, though. $\endgroup$
    – gung - Reinstate Monica
    Commented Jan 25, 2020 at 1:46
  • $\begingroup$ @gung-ReinstateMonica Thank you for your answer here. I know this is old, but how would you expand this if score had a quadratic component (e.g., B1 * score + B2 * score2) and you wanted to plot the combined effects with a confidence band (I'm also using R)? I can make a new post as well. $\endgroup$
    – aarsmith
    Commented Apr 25, 2023 at 17:04
  • $\begingroup$ @aarsmith, you should probably ask a new question, but it's going to be the same thing. Once you've fit a model with whatever variables, the rest is identical. $\endgroup$
    – gung - Reinstate Monica
    Commented Apr 25, 2023 at 18:15
  • $\begingroup$ @gung-ReinstateMonica the same as in you could just set score = 0:10 and scoresquared = 0:10^2? $\endgroup$
    – aarsmith
    Commented Apr 25, 2023 at 19:59

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

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