# Prediction plot and confidence intervals problems. Is this the best way to plot the relationship here?

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:

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:

• 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.
– EdM
Jan 24 '20 at 20:07
• 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.
– Paze
Jan 24 '20 at 20:09
• 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.
– EdM
Jan 24 '20 at 20:16
• 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.
– Paze
Jan 24 '20 at 20:49
• 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.
– Paze
Jan 24 '20 at 21:24

## 1 Answer

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")


• 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?
– Paze
Jan 25 '20 at 0:57
• @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. Jan 25 '20 at 1:46
{"extraction_info": {"found_math": true, "script_math_tex": 0, "script_math_asciimath": 0, "math_annotations": 0, "math_alttext": 0, "mathml": 0, "mathjax_tag": 0, "mathjax_inline_tex": 0, "mathjax_display_tex": 0, "mathjax_asciimath": 1, "img_math": 0, "codecogs_latex": 0, "wp_latex": 0, "mimetex.cgi": 0, "/images/math/codecogs": 0, "mathtex.cgi": 0, "katex": 0, "math-container": 8, "wp-katex-eq": 0, "align": 0, "equation": 0, "x-ck12": 0, "texerror": 0, "math_score": 0.541089653968811, "perplexity": 1397.627280555644}, "config": {"markdown_headings": true, "markdown_code": true, "boilerplate_config": {"ratio_threshold": 0.18, "absolute_threshold": 10, "end_threshold": 15, "enable": false}, "remove_buttons": true, "remove_image_figures": true, "remove_link_clusters": true, "table_config": {"min_rows": 2, "min_cols": 3, "format": "plain"}, "remove_chinese": true, "remove_edit_buttons": true, "extract_latex": true}, "warc_path": "s3://commoncrawl/crawl-data/CC-MAIN-2021-43/segments/1634323588284.71/warc/CC-MAIN-20211028100619-20211028130619-00362.warc.gz"}