combined visual presentation of logistic regression results

[Carrie Ngongo]

I am conducting logistic regressions with a dichotomous outcome (iatrogenic) and a single, continuous predictor (yeardevelop). I am running these regressions for different countries. I'm wondering about the best way to visually present the information. I could create formulas for each country and plot using Excel, but surely there must be a way to do this within Stata.

My input:

 logit iatrogenic yeardevelop if place==1

... place==2, place==3, etc.

A manual approach would plot:
log(p/1-p) = b0 + b1*yeardevelop
Where the years of interest are 1990-2020.

I will appreciate your guidance on (1) how to get Stata to visually present the curve and (2) how to save the information for multiple logistic regressions so that I can overlay the models for each of the countries, therefore demonstrating differences between the countries.

I will be most grateful for your help. Thank you.

[Nick Cox]

Consider this example:

webuse lbw, clear 

logit low age i.race 
label var low "birthweight < 2500 g"

set scheme s1color
regplot , sep(race) fitopts(lc(pink black blue)) mc(pink black blue) ms(Oh + Sh) jitter(3 ..) legend(order(2 3 1 5 6 4) col(1) pos(3)) name(G1, replace)

stripplot age, over(low) by(race, col(1) compact note(means and 95% confidence limits))  ms(sh) height(0.7) stack name(G2, replace) subtitle(, pos(9) fcolor(none) nobox nobexpand) bar

Here regplot is from the Stata Journal

. search regplot, sj

Search of official help files, FAQs, Examples, and Stata Journals

SJ-10-1 gr0009_1  . . . .  Software update for model diagnostic graph commands
        . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .  N. J. Cox
        (help anovaplot, indexplot, modeldiag, ofrtplot, ovfplot,
        qfrplot, racplot, rdplot, regplot, rhetplot, rvfplot2,
        rvlrplot, rvpplot2 if installed)
        Q1/10   SJ 10(1):164
        provides new command rbinplot for plotting means or medians
        of residuals by bins; provides new options for smoothing
        using restricted cubic splines; updates anova examples

SJ-4-4  gr0009  . . . . . . . . . . Speaking Stata: Graphing model diagnostics
        . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .  N. J. Cox
        (help anovaplot, indexplot, modeldiag, ofrtplot, ovfplot,
        qfrplot, racplot, rdplot, regplot, rhetplot, rvfplot2,
        rvlrplot, rvpplot2 if installed)
        Q4/04   SJ 4(4):449--475
        plotting diagnostic information calculated from residuals
        and fitted values from regression models with continuous
        responses

and stripplot is from SSC.

In your case I would start with one model, not several:

 
 logit iatrogenic yeardevelop i.place

[Carrie Ngongo]

Brilliant. Thanks so much for taking the time to direct me.

A couple of follow-up questions:
- Is it possible to suppress the data on regplot? As you'll see below, with my dichotomous outcome the data points along 0 and 1 are not particularly helpful. It would be ideal to keep the curves only, which would allow me to resize the y axis to, say, 0-0.4 (permitting a focus on the fitted curves).

- Curiously, the legend here includes only six countries, while my data set has nine. I see nine curves. Might you share a trick to have all nine names appear in the legend, please?

- All of this uses the code you recommended:

logit iatrogenic yeardevelop i.place

I assume here that the "i.place" approach is quicker but not statistically different from performing multiple logit regressions with the earlier approach of "if place==1", "if place==2," etc.

[Joseph Coveney]

I assume here that the "i.place" approach is quicker but not statistically different from performing multiple logit regressions with the earlier approach of "if place==1", "if place==2," etc.

No. Individual regression models will fit an individual (unique) temporal slope for each country. In order to obtain the same end with a single regression model, you'll need to include a country × year interaction term. See below. (Begin at the "Begin here" comment. The top part of the do-file just creates a fictitious dataset for illustration. I've made up random three-letter abbreviations to stand in for country names, because I don't know the names of the other three of your nine total countries.)

Your other questions—about plotting the fitted model predictions—are also illustrated below.

.ÿ
.ÿversionÿ16.1

.ÿ
.ÿclearÿ*

.ÿ
.ÿsetÿseedÿ`=strreverse("1567746")'

.ÿ
.ÿ//ÿNineÿcountries
.ÿquietlyÿsetÿobsÿ9

.ÿ
.ÿquietlyÿgenerateÿstrÿcnaÿ=ÿ""

.ÿforvaluesÿiÿ=ÿ1/3ÿ{
ÿÿ2.ÿÿÿÿÿquietlyÿreplaceÿcnaÿ=ÿcnaÿ+ÿchar(runiformint(65,ÿ90))
ÿÿ3.ÿ}

.ÿisidÿcna

.ÿencodeÿcna,ÿgenerate(cid)ÿlabel(Countries)

.ÿdrawnormÿcid_uÿcid×tim_u,ÿdoubleÿcorr(1ÿ-0.25ÿ\ÿ-0.25ÿ1)

.ÿ
.ÿ//ÿObservationÿintervalÿ1990ÿtoÿ2020,ÿinclusive
.ÿquietlyÿexpandÿ`=2020ÿ-ÿ1990ÿ+ÿ1'

.ÿbysortÿcid:ÿgenerateÿintÿtimÿ=ÿ1990ÿ+ÿ_nÿ-ÿ1

.ÿ
.ÿ//ÿIndividualÿobservationsÿ(20ÿpatientsÿperÿcountryÿperÿyear)
.ÿgenerateÿdoubleÿxbÿ=ÿcid_uÿ+ÿ(timÿ-ÿ15ÿ-ÿ1990)ÿ/ÿ31ÿ*ÿ(0.1ÿ+ÿcid×tim_u)

.ÿ
.ÿquietlyÿexpandÿ20

.ÿgenerateÿbyteÿiatÿ=ÿrbinomial(1,ÿinvlogit(xb))

.ÿ
.ÿ*
.ÿ*ÿBeginÿhere
.ÿ*
.ÿlogitÿiatÿi.cid##c.tim,ÿnolog

LogisticÿregressionÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿNumberÿofÿobsÿÿÿÿÿ=ÿÿÿÿÿÿ5,580
ÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿLRÿchi2(17)ÿÿÿÿÿÿÿ=ÿÿÿÿ1063.52
ÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿProbÿ>ÿchi2ÿÿÿÿÿÿÿ=ÿÿÿÿÿ0.0000
Logÿlikelihoodÿ=ÿ-3295.6788ÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿPseudoÿR2ÿÿÿÿÿÿÿÿÿ=ÿÿÿÿÿ0.1389

------------------------------------------------------------------------------
ÿÿÿÿÿÿÿÿÿiatÿ|ÿÿÿÿÿÿCoef.ÿÿÿStd.ÿErr.ÿÿÿÿÿÿzÿÿÿÿP>|z|ÿÿÿÿÿ[95%ÿConf.ÿInterval]
-------------+----------------------------------------------------------------
ÿÿÿÿÿÿÿÿÿcidÿ|
ÿÿÿÿÿÿÿÿDNNÿÿ|ÿÿÿ-9.38163ÿÿÿ33.69474ÿÿÿÿ-0.28ÿÿÿ0.781ÿÿÿÿ-75.42211ÿÿÿÿ56.65885
ÿÿÿÿÿÿÿÿFTCÿÿ|ÿÿ-184.6763ÿÿÿ29.95707ÿÿÿÿ-6.16ÿÿÿ0.000ÿÿÿÿÿ-243.391ÿÿÿ-125.9615
ÿÿÿÿÿÿÿÿNMWÿÿ|ÿÿÿ-138.246ÿÿÿ29.66232ÿÿÿÿ-4.66ÿÿÿ0.000ÿÿÿÿ-196.3831ÿÿÿ-80.10895
ÿÿÿÿÿÿÿÿOKWÿÿ|ÿÿ-169.3484ÿÿÿ29.77976ÿÿÿÿ-5.69ÿÿÿ0.000ÿÿÿÿ-227.7157ÿÿÿ-110.9812
ÿÿÿÿÿÿÿÿRDIÿÿ|ÿÿ-2.625427ÿÿÿ32.15121ÿÿÿÿ-0.08ÿÿÿ0.935ÿÿÿÿ-65.64064ÿÿÿÿ60.38978
ÿÿÿÿÿÿÿÿUQGÿÿ|ÿÿ-208.5959ÿÿÿ32.81811ÿÿÿÿ-6.36ÿÿÿ0.000ÿÿÿÿ-272.9183ÿÿÿ-144.2736
ÿÿÿÿÿÿÿÿYOEÿÿ|ÿÿ-88.24496ÿÿÿ29.37302ÿÿÿÿ-3.00ÿÿÿ0.003ÿÿÿÿÿ-145.815ÿÿÿ-30.67489
ÿÿÿÿÿÿÿÿZEIÿÿ|ÿÿÿ51.65138ÿÿÿ30.31722ÿÿÿÿÿ1.70ÿÿÿ0.088ÿÿÿÿÿ-7.76928ÿÿÿÿÿ111.072
ÿÿÿÿÿÿÿÿÿÿÿÿÿ|
ÿÿÿÿÿÿÿÿÿtimÿ|ÿÿ-.0389524ÿÿÿ.0115678ÿÿÿÿ-3.37ÿÿÿ0.001ÿÿÿÿ-.0616249ÿÿÿ-.0162799
ÿÿÿÿÿÿÿÿÿÿÿÿÿ|
ÿÿÿcid#c.timÿ|
ÿÿÿÿÿÿÿÿDNNÿÿ|ÿÿÿÿÿ.00619ÿÿÿÿ.016805ÿÿÿÿÿ0.37ÿÿÿ0.713ÿÿÿÿ-.0267472ÿÿÿÿ.0391271
ÿÿÿÿÿÿÿÿFTCÿÿ|ÿÿÿ.0927458ÿÿÿ.0149484ÿÿÿÿÿ6.20ÿÿÿ0.000ÿÿÿÿÿ.0634474ÿÿÿÿ.1220441
ÿÿÿÿÿÿÿÿNMWÿÿ|ÿÿÿÿ.069812ÿÿÿ.0148034ÿÿÿÿÿ4.72ÿÿÿ0.000ÿÿÿÿÿ.0407979ÿÿÿÿ.0988261
ÿÿÿÿÿÿÿÿOKWÿÿ|ÿÿÿ.0851389ÿÿÿ.0148605ÿÿÿÿÿ5.73ÿÿÿ0.000ÿÿÿÿÿ.0560128ÿÿÿÿÿ.114265
ÿÿÿÿÿÿÿÿRDIÿÿ|ÿÿÿ.0013748ÿÿÿ.0160494ÿÿÿÿÿ0.09ÿÿÿ0.932ÿÿÿÿ-.0300814ÿÿÿÿÿ.032831
ÿÿÿÿÿÿÿÿUQGÿÿ|ÿÿÿÿ.104097ÿÿÿ.0163644ÿÿÿÿÿ6.36ÿÿÿ0.000ÿÿÿÿÿ.0720233ÿÿÿÿ.1361707
ÿÿÿÿÿÿÿÿYOEÿÿ|ÿÿÿ.0447831ÿÿÿ.0146582ÿÿÿÿÿ3.06ÿÿÿ0.002ÿÿÿÿÿ.0160536ÿÿÿÿ.0735127
ÿÿÿÿÿÿÿÿZEIÿÿ|ÿÿÿ-.025213ÿÿÿ.0151311ÿÿÿÿ-1.67ÿÿÿ0.096ÿÿÿÿ-.0548694ÿÿÿÿ.0044434
ÿÿÿÿÿÿÿÿÿÿÿÿÿ|
ÿÿÿÿÿÿÿ_consÿ|ÿÿÿ76.67717ÿÿÿ23.17242ÿÿÿÿÿ3.31ÿÿÿ0.001ÿÿÿÿÿ31.26007ÿÿÿÿ122.0943
------------------------------------------------------------------------------

.ÿ
.ÿquietlyÿmarginsÿcid,ÿat(tim=(1990(2)2020))

.ÿmarginsplotÿ,ÿ///
>ÿÿÿÿÿÿÿÿÿtitle("")ÿnociÿ///
>ÿÿÿÿÿÿÿÿÿlegend(cols(4)ÿrows(2))ÿ///
>ÿÿÿÿÿÿÿÿÿplotopts(msymbol(none)ÿ///
>ÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿytitle(ProportionÿIatrogenicÿFistula)ÿ///
>ÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿylabel(ÿ,ÿformat(%03.1f)ÿangle(horizontal)ÿnogrid)ÿ///
>ÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿxtitle(YearÿofÿObservation)ÿxlabel(,ÿangle(45)))

ÿÿVariablesÿthatÿuniquelyÿidentifyÿmargins:ÿtimÿcid

.ÿquietlyÿgraphÿexportÿIatrogenicFistula.png,ÿreplace

.ÿ
.ÿexit

endÿofÿdo-file


.

The resulting graph is inserted below. Nine plots on a single graph is a lot. You might want to group the countries and create multiple graphs of fewer plots. And you might want to tweak the options for marginsplot in order to get something more to your liking in terms of appearance. (I, for one, would move the legend out of the graph completely and into a figure legend.)

[Carrie Ngongo]

I appreciate your guidance and code. I have a dichotomous outcome of interest. I want to look at it two ways: first, among all deliveries, and second, among operative deliveries (Cesarean section). And I'd like to do this for nine different countries.

Unfortunately, the approach that you propose does not result in smooth lines for each country. Part of the issue is that there are no observations for some countries in some years. Given that this is a 27-year time period (with observations between 1990 and 2017), interaction terms between country and year run into low sample sizes, even though the overall sample size is not small (~5,600 total deliveries, ~3,000 CS).

I recognize that my original approach was not elegant:

logit iatrogenic yeardevelop if place==1
logit iatrogenic yeardevelop if place==2
logit iatrogenic yeardevelop if place==3
...

Yet what I appreciate about this is that it gives me a clear story to tell for each country. I can see the comparative magnitude of the coefficients and which countries have statistically significant results. I am not able to glean a clear story from an approach that creates an interaction term between place and yeardevelop.

The overarching challenge for me is how best to visually convey the logit results. The best option that I've tried is Nick Cox's regplot, although I find busy and unhelpful the burst of color along the 0 and 1 lines. This is true even when I use regplot for a single country. Also, I'd appreciate a trick for having all nine countries show up in the legend for regplot rather than just the top six. Thanks.

[Joseph Coveney]

Unfortunately, the approach that you propose does not result in smooth lines for each country. Part of the issue is that there are no observations for some countries in some years. Given that this is a 27-year time period (with observations between 1990 and 2017), interaction terms between country and year run into low sample sizes, even though the overall sample size is not small (~5,600 total deliveries, ~3,000 CS).

. . . Yet what I appreciate about this is that it gives me a clear story to tell for each country. I can see the comparative magnitude of the coefficients and which countries have statistically significant results. I am not able to glean a clear story from an approach that creates an interaction term between place and yeardevelop.

. . . I'd appreciate a trick for having all nine countries show up in the legend for regplot rather than just the top six.

Study the following. You'll find that all of your requests and concerns are addressed using the approach that I gave above.

.ÿ
.ÿversionÿ16.1

.ÿ
.ÿclearÿ*

.ÿ
.ÿsetÿseedÿ`=strreverse("1567929")'

.ÿ
.ÿ//ÿNineÿcountries
.ÿquietlyÿsetÿobsÿ9

.ÿ
.ÿquietlyÿgenerateÿstrÿcnaÿ=ÿ""

.ÿforvaluesÿiÿ=ÿ1/3ÿ{
ÿÿ2.ÿÿÿÿÿquietlyÿreplaceÿcnaÿ=ÿcnaÿ+ÿchar(runiformint(65,ÿ90))
ÿÿ3.ÿ}

.ÿisidÿcna

.ÿencodeÿcna,ÿgenerate(cid)ÿlabel(Countries)

.ÿdrawnormÿcid_uÿcid×tim_u,ÿdoubleÿcorr(1ÿ-0.25ÿ\ÿ-0.25ÿ1)

.ÿ
.ÿ//ÿObservationÿintervalÿ1990ÿtoÿ2017,ÿinclusive
.ÿquietlyÿexpandÿ`=ÿ2017-ÿ1990ÿ+ÿ1'

.ÿbysortÿcid:ÿgenerateÿintÿtimÿ=ÿ1990ÿ+ÿ_nÿ-ÿ1

.ÿ
.ÿ//ÿIndividualÿobservationsÿ(20ÿpatientsÿperÿcountryÿperÿyear)
.ÿgenerateÿdoubleÿxbÿ=ÿcid_uÿ+ÿ(timÿ-ÿ15ÿ-ÿ1990)ÿ/ÿ31ÿ*ÿ(0.1ÿ+ÿcid×tim_u)

.ÿ
.ÿ//ÿBeginÿ"noÿobservationsÿforÿsomeÿcountriesÿinÿsomeÿyears"
.ÿcount
ÿÿ252

.ÿdropÿifÿruniform()ÿ>=ÿ0.95ÿ//ÿdropÿaboutÿ5%ÿofÿcountry-yearÿcombinations
(16ÿobservationsÿdeleted)

.ÿ//ÿEndÿ"noÿobservationsÿforÿsomeÿcountriesÿinÿsomeÿyears"
.ÿ
.ÿquietlyÿexpandÿ20

.ÿgenerateÿbyteÿiatÿ=ÿrbinomial(1,ÿinvlogit(xb))

.ÿ
.ÿ*
.ÿ*ÿBeginÿhere
.ÿ*
.ÿlogitÿiatÿibn.cidÿibn.cid#c.tim,ÿnoconstantÿnolog

LogisticÿregressionÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿNumberÿofÿobsÿÿÿÿÿ=ÿÿÿÿÿÿ4,720
ÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿWaldÿchi2(18)ÿÿÿÿÿ=ÿÿÿÿÿ648.07
Logÿlikelihoodÿ=ÿ-2884.9719ÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿProbÿ>ÿchi2ÿÿÿÿÿÿÿ=ÿÿÿÿÿ0.0000

------------------------------------------------------------------------------
ÿÿÿÿÿÿÿÿÿiatÿ|ÿÿÿÿÿÿCoef.ÿÿÿStd.ÿErr.ÿÿÿÿÿÿzÿÿÿÿP>|z|ÿÿÿÿÿ[95%ÿConf.ÿInterval]
-------------+----------------------------------------------------------------
ÿÿÿÿÿÿÿÿÿcidÿ|
ÿÿÿÿÿÿÿÿCHMÿÿ|ÿÿ-3.958147ÿÿÿ25.70204ÿÿÿÿ-0.15ÿÿÿ0.878ÿÿÿÿ-54.33321ÿÿÿÿ46.41692
ÿÿÿÿÿÿÿÿCMMÿÿ|ÿÿ-51.68557ÿÿÿÿ23.2675ÿÿÿÿ-2.22ÿÿÿ0.026ÿÿÿÿ-97.28902ÿÿÿ-6.082109
ÿÿÿÿÿÿÿÿEQYÿÿ|ÿÿÿ47.58423ÿÿÿ24.39646ÿÿÿÿÿ1.95ÿÿÿ0.051ÿÿÿÿ-.2319417ÿÿÿÿ95.40041
ÿÿÿÿÿÿÿÿKISÿÿ|ÿÿ-77.56061ÿÿÿÿ24.4504ÿÿÿÿ-3.17ÿÿÿ0.002ÿÿÿÿ-125.4825ÿÿÿ-29.63871
ÿÿÿÿÿÿÿÿUMHÿÿ|ÿÿÿ41.84074ÿÿÿ26.71574ÿÿÿÿÿ1.57ÿÿÿ0.117ÿÿÿÿ-10.52115ÿÿÿÿ94.20263
ÿÿÿÿÿÿÿÿURFÿÿ|ÿÿ-11.45371ÿÿÿ22.19142ÿÿÿÿ-0.52ÿÿÿ0.606ÿÿÿÿ-54.94809ÿÿÿÿ32.04068
ÿÿÿÿÿÿÿÿVEZÿÿ|ÿÿÿ41.63501ÿÿÿ21.41299ÿÿÿÿÿ1.94ÿÿÿ0.052ÿÿÿÿ-.3336696ÿÿÿÿÿ83.6037
ÿÿÿÿÿÿÿÿXNWÿÿ|ÿÿ-39.41971ÿÿÿ25.99171ÿÿÿÿ-1.52ÿÿÿ0.129ÿÿÿÿ-90.36252ÿÿÿÿ11.52311
ÿÿÿÿÿÿÿÿYLHÿÿ|ÿÿÿ-12.3842ÿÿÿ21.19096ÿÿÿÿ-0.58ÿÿÿ0.559ÿÿÿÿ-53.91772ÿÿÿÿ29.14932
ÿÿÿÿÿÿÿÿÿÿÿÿÿ|
ÿÿÿcid#c.timÿ|
ÿÿÿÿÿÿÿÿCHMÿÿ|ÿÿÿ.0013164ÿÿÿ.0128281ÿÿÿÿÿ0.10ÿÿÿ0.918ÿÿÿÿ-.0238263ÿÿÿÿ.0264591
ÿÿÿÿÿÿÿÿCMMÿÿ|ÿÿÿ.0254187ÿÿÿÿ.011612ÿÿÿÿÿ2.19ÿÿÿ0.029ÿÿÿÿÿ.0026597ÿÿÿÿ.0481778
ÿÿÿÿÿÿÿÿEQYÿÿ|ÿÿ-.0242233ÿÿÿ.0121772ÿÿÿÿ-1.99ÿÿÿ0.047ÿÿÿÿ-.0480901ÿÿÿ-.0003565
ÿÿÿÿÿÿÿÿKISÿÿ|ÿÿÿ.0386761ÿÿÿ.0122025ÿÿÿÿÿ3.17ÿÿÿ0.002ÿÿÿÿÿ.0147597ÿÿÿÿ.0625925
ÿÿÿÿÿÿÿÿUMHÿÿ|ÿÿ-.0215157ÿÿÿ.0133357ÿÿÿÿ-1.61ÿÿÿ0.107ÿÿÿÿ-.0476532ÿÿÿÿ.0046218
ÿÿÿÿÿÿÿÿURFÿÿ|ÿÿÿ.0056849ÿÿÿ.0110738ÿÿÿÿÿ0.51ÿÿÿ0.608ÿÿÿÿ-.0160194ÿÿÿÿ.0273892
ÿÿÿÿÿÿÿÿVEZÿÿ|ÿÿ-.0209552ÿÿÿÿ.010689ÿÿÿÿ-1.96ÿÿÿ0.050ÿÿÿÿ-.0419052ÿÿÿ-5.24e-06
ÿÿÿÿÿÿÿÿXNWÿÿ|ÿÿÿ.0190307ÿÿÿ.0129711ÿÿÿÿÿ1.47ÿÿÿ0.142ÿÿÿÿ-.0063921ÿÿÿÿ.0444535
ÿÿÿÿÿÿÿÿYLHÿÿ|ÿÿÿ.0061732ÿÿÿ.0105763ÿÿÿÿÿ0.58ÿÿÿ0.559ÿÿÿÿ-.0145558ÿÿÿÿ.0269023
------------------------------------------------------------------------------

.ÿ
.ÿquietlyÿmarginsÿcid,ÿat(tim=(1990(3)2017))

.ÿmarginsplotÿ,ÿ///
>ÿÿÿÿÿÿÿÿÿtitle("")ÿnociÿ///
>ÿÿÿÿÿÿÿÿÿlegend(cols(4)ÿrows(2))ÿ///
>ÿÿÿÿÿÿÿÿÿplotopts(msymbol(none)ÿ///
>ÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿytitle(ProportionÿIatrogenicÿFistula)ÿ///
>ÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿylabel(ÿ,ÿformat(%03.1f)ÿangle(horizontal)ÿnogrid)ÿ///
>ÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿÿxtitle(YearÿofÿObservation)ÿxlabel(,ÿangle(45)))

ÿÿVariablesÿthatÿuniquelyÿidentifyÿmargins:ÿtimÿcid

.ÿquietlyÿgraphÿexportÿIatrogenicFistula2.png,ÿreplace

.ÿ
.ÿexit

endÿofÿdo-file


.

[Carrie Ngongo]

I'm happy and grateful to report that I have managed to study your example and produce a lovely visual. Thanks so much, Joseph Coveney. Even though nine lines is indeed quite a lot, I really like how this visual representation shares a compelling story about how countries differ from one another over time.

I have two follow-up questions. First, might you have a trick to move the y-axis label further left so that it doesn't overlap with the axis tick labels? My country names also don't fit in the legend box, although I imagine that this can be solved if I use three-letter abbreviations as you do. I could also add labels at the right-hand end of each line (country name and perhaps p-value) using graphics software, which would make this easier for readers to interpret. I'll appreciate hearing if you would recommend an alternative solution.

Second, I would be most grateful if you could help me to describe your approach. I gather that you've kept all levels of the categorical place variable in the model, with no country treated as the base. You also created an interaction term between iatrogenic and each year, smoothing in three-year intervals (?). I would be most grateful if you could point me to an example of how someone has described this approach for publication. You know that I was previously relying on a bunch of separate logit regressions by country, so I want to be sure that I understand how this is different and how it should most appropriately be described.

Many thanks again for your help. I sincerely appreciate it.

[Carrie Ngongo]

Just as a quick postscript, I now catch that the three years just relate to the axis labels, not smoothing. And I would propose to put country names at the right-hand side of each line, but skip the p values. Thanks.

[Andrew Musau]

It is difficult to customize graph results from margins, so the best advice is to save these results and recreate the graph using a different command. Joseph Coveney has provided a reproducible example in #6, so here is a way to address the -ytitle- issue and to add labels within the graph.

version 16.1
clear
set seed `=strreverse("1567929")'
quietly set obs 9
quietly generate str cna=""
forvalues i=1/3{
    quietly replace cna=cna + char(runiformint(65, 90))
}
isid cna
encode cna, generate(cid) label(Countries)
drawnorm cid_u cid×tim_u, double corr(1 -0.25 \ -0.25 1)
quietly expand `=2017-1990+1'
bysort cid: generate int tim=1990+_n-1
generate double xb=cid_u+(tim-15-1990)/31*(0.1+cid×tim_u)
count
drop if runiform()>=0.95
quietly expand 20
generate byte iat=rbinomial(1,invlogit(xb))
logit iat ibn.cid ibn.cid#c.tim, noconstant nolog
quietly margins cid, at(tim=(1990(3)2017)) saving(myfile, replace)
use myfile, clear
xtset _m1 _at2
decode _m1, gen(country)
xtline _margin, overlay ytitle(Proportion Iatrogenic Fistula, margin(medium)) ///
ylabel(,format(%03.1f) angle(horizontal) nogrid) xtitle(Year of Observation) ///
xlabel(1990(3)2017, angle(45)) leg(off) scheme(s1color) ///
addplot(scatter _margin _at2 if _at2 ==2017, mcolor(none) mlab(country) xsc(r(. 2021)))

Res.:


With so many lines in one plot, readers will have difficulties following which line belongs to which country and labels may overlap and therefore be illegible. Therefore, I would recommend fabplot, from SSC by Nick Cox, which highlights data from one country with the rest as a backdrop. For the example above

*ssc install fabplot
fabplot line _margin _at2, by(country) xla(1990(3)2017, ///
format(%tyYY)) xtitle("") frontopts(lc(orange_red) lw(thick)) ///
ytitle(Proportion Iatrogenic Fistula) scheme(s1color)

Second, I would be most grateful if you could help me to describe your approach.

Finally, predictive margins and marginal effects are a pretty standard way of presenting results. See -help margins- and the manual entry of margins. Additionally, Richard Williams has a nice article on this in the Stata Journal (free to download).

https://journals.sagepub.com/doi/10....867X1201200209

[Andrew Musau]

">=" displays as ">" in the code in #9. Also, some issues with visibility of the graphs. So here is the code and graphs reposted.

version 16.1
clear
set seed `=strreverse("1567929")'
quietly set obs 9
quietly generate str cna=""
forvalues i=1/3{
    quietly replace cna=cna + char(runiformint(65, 90))
}
isid cna
encode cna, generate(cid) label(Countries)
drawnorm cid_u cid×tim_u, double corr(1 -0.25 \ -0.25 1)
quietly expand `=2017-1990+1'
bysort cid: generate int tim=1990+_n-1
generate double xb=cid_u+(tim-15-1990)/31*(0.1+cid×tim_u)
count
drop if runiform()>=0.95
quietly expand 20
generate byte iat=rbinomial(1,invlogit(xb))
logit iat ibn.cid ibn.cid#c.tim, noconstant nolog
quietly margins cid, at(tim=(1990(3)2017)) saving(myfile, replace)
use myfile, clear
xtset _m1 _at2
decode _m1, gen(country)
xtline _margin, overlay ytitle(Proportion Iatrogenic Fistula, margin(medium)) ///
ylabel(,format(%03.1f) angle(horizontal) nogrid) xtitle(Year of Observation) ///
xlabel(1990(3)2017, angle(45)) leg(off) scheme(s1color) ///
addplot(scatter _margin _at2 if _at2 ==2017, mcolor(none) mlab(country) xsc(r(. 2021)))
*ssc install fabplot
fabplot line _margin _at2, by(country) xla(1990(3)2017, ///
format(%tyYY)) xtitle("") frontopts(lc(orange_red) lw(thick)) ///
ytitle(Proportion Iatrogenic Fistula) scheme(s1color)

[Carrie Ngongo]

This is fabulous indeed. Thank you, Andrew Musau. I appreciate the good ideas and the helpful reference.

The fabplot graphs will be just what I need. I like your approach of having an overall graph above. When I produce a similar graph, however, I have many countries whose lines end close together. This results in country names being written over one another illegibly, like this:

Might you happen to have a trick for telling Stata to move the names a bit so that the order is clear and no names are written on top of others?

Many thanks.

[Andrew Musau]

You need to do this manually. Create another variable where you introduce a gap between observations for the last time period in the sample. Still, there remains the problem of identifying which line belongs to which label. Here is an even more extreme example. I use a condition to define the spacing, but you can insert the values by hand to get a better fit.

webuse grunfeld, clear
xtline invest, overlay addplot(scatter invest year if year==1954, mcolor(none) ///
mlabel(mvalue) leg(off) scheme(s1color) saving(gr1, replace) xtitle("") xsc(r(. 1957)))
sort invest
gen order=sum(year==1954)
gen invest2= cond(invest <300 & order>2, invest+ 22*order, invest)
xtline invest, overlay addplot(scatter invest2 year if year==1954, mcolor(none) ///
mlabel(mvalue) leg(off) scheme(s1color) saving(gr2, replace) xtitle("") xsc(r(. 1957)))
gr combine gr1.gph gr2.gph, scheme(s1color)

[Nick Cox]

I picked up on this thread only today and appreciate the mention of fabplot (SSC), which by a pleasant coincidence I have just revised, although I'll announce an update separately.

In #12 a logarithmic scale will help mightily. The same Grunfeld data appears extensively in many places, but most relevantly here in https://journals.sagepub.com/doi/abs...urnalCode=stja

In #11 log scale might help too although in some ways a logit scale is more natural. See also https://www.stata-journal.com/sjpdf....iclenum=gr0032 and mylabels (also from SSC);,

[Carrie Ngongo]

Andrew Musau, thank you for the code to try to fix labels manually. Unfortunately, I haven't managed to crack the overlapping text issue as you did. The first part of my code works fine:

logit iatrogenic ibn.place ibn.place#c.yeardevelop, noconstant nolog
quietly margins place, at(yeardevelop=(1990(3)2017)) saving(alldeliveries, replace)
use alldeliveries, clear
xtset _m1 _at2
decode _m1, gen(nation)
xtline _margin, overlay addplot(scatter _margin _at2 if _at2==2017, mcolor(none) mlabel(nation) leg(off) xtitle(Year of Observation) xlabel(1990(3)2020, angle(45)) leg(off) scheme(s1color) ytitle(Proportion Iatrogenic Fistula, margin(medium)) ylabel(, format(%03.1f) angle(horizontal) nogrid)  xsc(r(. 2017))) 
graph save gr1, replace

When you generate order, what does this mean? It doesn't seem to do the intuitive thing of listing countries by descending order of magnitude in 2017, as you'll see in the tabulation of order and nation below.

sort _margin
gen order=sum(_at2 ==2017)
tab order nation

gen _margin2= cond(_margin <0.3 & order >2, _margin+0.02*order, _margin)
xtline _margin, overlay addplot(scatter _margin2 _at2 if _at2==2017, mcolor(none) mlabel(nation) leg(off) xtitle(Year of Observation) xlabel(1990(3)2020, angle(45)) leg(off) scheme(s1color) ytitle(Proportion Iatrogenic Fistula, margin(medium)) ylabel(, format(%03.1f) angle(horizontal) nogrid) xsc(r(. 2017)))
graph save gr2, replace
gr combine gr1.gph gr2.gph, scheme(s1color)

I've managed to make it worse, as you can see below. South Sudan is nowhere near its line, and I haven't spaced out the other country names. I'm showing you my axis so that you can see what I'm working with. The labels for the three top lines are fine. Any chance that you could direct me to where I should tinker in the code to be able to avoid overlapping text?

Many thanks for your help.