ppplot and its if statement

[Chen Samulsion]

Dear Stata users,

I tried using Nick Cox's -ppplot- command (SSC) to compare the cumulative distributions of variable x across multiple groups. The command runs perfectly when I do not use an if qualifier. However, when I try to include an if qualifier in the command to select specific groups, ppplot outputs a blank graph. Does ppplot not support the if qualifier, or is there something I am misunderstanding about how to use this command? Thank you.

use2 lbw
tabulate race
ppplot line bwt, by(race) plot(function y=x, clp(dash)) name(g1)
ppplot line bwt if race<=2, by(race) plot(function y=x, clp(dash)) name(g2)
keep if race<=2
ppplot line bwt, by(race) plot(function y=x, clp(dash)) name(g3)

[Nick Cox]

Thanks for providing an example, but I don't get past a reference to use2 -- what is that please? Also where is lbw.dta?

[Chen Samulsion]

Sorry Nick, use2 is my own small program to load dataset in my personal directory, in order to obtain lbw, you can just type

webuse lbw

[Nick Cox]

Thanks. You found a bug! The separate command should be

separate `varlist' if `touse', by(`by') gen(`stub') `missing' short

The command goes back to 2004. What the help most needs is, I think, an ack't to Chen and a plug for qplot (Stata Journal).

[Chen Samulsion]

Sorry, the problem seems aroused by the data structure itself. Maybe the second category (race==black, N=26) has too less observations?

webuse lbw
tabulate race
ppplot line bwt, by(race) plot(function y=x, clp(dash)) name(gr123)
ppplot line bwt if race==1 | race==2, by(race) plot(function y=x, clp(dash)) name(gr12)
ppplot line bwt if race==1 | race==3, by(race) plot(function y=x, clp(dash)) name(gr13)
ppplot line bwt if race==2 | race==3, by(race) plot(function y=x, clp(dash)) name(gr23)
graph combine gr123 gr12 gr13 gr23

. tab race

       race |      Freq.     Percent        Cum.
------------+-----------------------------------
      white |         96       50.79       50.79
      black |         26       13.76       64.55
      other |         67       35.45      100.00
------------+-----------------------------------
      Total |        189      100.00

[Chen Samulsion]

The reason I have returned to this (ancient) command developed in 2004 is that I have seen references in the literature to the use of P-P (probability-probability) plots for testing the proportional odds (PO) assumption. I believe that P-P plots can be very useful. Statistical analysis of Likert-based ordinal scales: a guide for clinical trialists

[Chen Samulsion]

Fig. 2
The probability-probability (P-P) plot displays the cumulative probabilities of ordinal scores (0–10) for the control (x-axis) and intervention (y-axis) groups. All data points lie below the reference line (y = x), indicating that the control group consistently has higher cumulative probabilities across most scores. The varying vertical deviations (0.08 to 0.21) demonstrate inconsistent differences between groups, suggesting that the PO assumption is likely violated.

[Chen Samulsion]

Thank you so much Nick Cox. I correct the code that you pointed out, now it works! My speculation in #5 was completely wrong.

Thanks. You found a bug! The separate command should be

separate `varlist' if `touse', by(`by') gen(`stub') `missing' short

The command goes back to 2004. What the help most needs is, I think, an ack't to Chen and a plug for qplot (Stata Journal).

[Nick Cox]

Thanks for that closure.

I readily admit that you have to use any particular graph type repeatedly before you have learned how to use them easily and effectively! That applies even to simple kinds such as bar charts, histograms and scatter plots, although you've probably forgotten whatever small struggles were needed for them to become familiar.

That said, I have never warmed to these P-P plots and have rarely used this command since writing it.

There is, I think, a general reason and a specific reason for my disinclination.

Generally, I find them more indirect than the alternatives. If the issue is comparison of probability distributions. there are several more direct possibilities.

Specifically, I find that such plots are better at comparing the middles of distributions than the tails, just as Kolmogorov-Smirnov tests and the like are not so good at comparing tails. I drew this pithy quotation to the attention of StataCorp and it can be found at [R] diagnostic plots.

The idea behind qnorm is recommended strongly by Miller (1997): he calls it probit plotting. His
recommendations from much practical experience should interest many users. “My recommendation for
detecting nonnormality is probit plotting” (Miller 1997, 10). “If a deviation from normality cannot be
spotted by eye on probit paper, it is not worth worrying about. I never use the Kolmogorov–Smirnov
test (or one of its cousins) or the [chi-square] test as a preliminary test of normality. They do not tell you how the
sample is differing from normality, and I have a feeling they are more likely to detect irregularities in the
middle of the distribution than in the tails” (Miller 1997, 13–14).

Let's use the dataset of this thread and look at one alternative, quantile plots. (There are good cases for dot plots, histograms with bin widths not too small, and some others.)

A quantile plot is just the ordered data plotted as such. There are many other terms for the same idea, and quantile plots were used avant la lettre from at least the early 19th century onwards. A recent book https://www.routledge.com/Explorator.../9781032814803 uses the term sorted data plot, which may appeal.

Quantile plots may and often should be complemented by other displays. A simple but subtle enough philosophy is to aim at a display that tries to do two things at once, show detail that might be interesting or important, and show a summary of the data too. Here is one example.


If this seems over-complicated, it is not compulsory either. There isn't anything particularly exotic or unusual about any of the choices.

1. The quantile plots as said show the ordered data (and all of them). We get to see a fair amount of detail, including a quiet signal about sample sizes (rather few black babies) and a direct display of every data point, including moderate outliers in the case of this dataset. There is no need for arbitrary choices about binning or to resort to devices like jittering the data to shake points apart.

2. There is a box plot with median and quartiles. I tend not to use the Tukey convention of plotting data points individually if and only if they lie more than 1.5 IQR from the nearer quartile. All too often that is used but not explained, or explained quite incorrectly, and my wicked guess is that even statistical people who have seen it often used would find it hard to explain the exact rule without checking a reference. Here I just use 5 and 95% points as a simpler criterion for where spikes stop, but by far the larger point is to choose your own rule according to taste and whatever suits the data and the problem. Showing medians as markers or point symbols is in fact a convention that Tukey used in later papers. I don't know if he explained why, but my motive is to be able to refer to the medians directly. Otherwise you have some rigmarole to explain: the box plot shows three lines for the quartiles and the median, and so on. But as in #1 you are showing all the data points any way, so exactly what you show to summarize what is going on in the tails is not so important, and up to you always.

3. I show means as horizontal lines. If you like, there is an echo here of the fact that the mean is the area under the quantile function implied by the array of quantiles. Graphically the larger imperative is that showing means as extra points or markers added to the box plots can just complicate or confuse interpretation of those box plots, even if it turns out that mean and median are close.

4. I try to remember to explain what I am doing, using note() for that. Look carefully at books and papers and you'll find that to be uncommon, surprisingly or not.

Here is some code. As flagged within I use qplot 2.5.0 which is in press at the Stata Journal, but in any previous version of qplot you could get the same x axis labels 0 0.25 0.5 0.75 1 (or any others you liked) through option choices.

Using pctilesets to get a data set of percentiles and then merging with the main datasets may seem more complicated than you want, but the attitude is to use a toolset to get what you want. Other people will have different preferences, which is what it should be. Once you've worked your way to a serviceable script you can recycle code in a do-file or ado-file.

webuse lbw, clear

* qplot 2.5.0 Stata Journal 26(3)
qplot bwt , by(race, row(1) compact)
egen mean = mean(bwt), by(race)
bysort race (bwt) : gen x = cond(_n == 1, 0, cond(_n == _N, 1, 0))
qplot bwt , by(race, row(1) compact note("lines show means") legend(off)) addplot(line mean x) xtitle(Fraction of data) name(G1, replace)

* pctilesets from Stata Journal 26(2)
pctilesets bwt, over(race) pctile(5 25 50 75 95) saving(pctiles, replace)
clonevar origgvar=race
merge m:1 origgvar using pctiles
gen where = 1.1

qplot bwt , by(race, row(1) compact note("lines show means; diamonds show medians; boxes show quartiles; spikes to 5% and 95% points") legend(off)) xtitle(Fraction of data) addplot(line mean x || scatter p50 where, ms(Dh) msize(medlarge) pstyle(p2) || rbar p25 p75 where, barw(0.1) pstyle(p2) fcol(none) || rspike p75 p95 where, pstyle(p2) || rspike p25 p5 where, pstyle(p2))  name(G2, replace)

[Nick Cox]

Thanks as ever to KitBaum , ppplot has been updated on SSC. I fixed the bug reported by Chen Samulsion and took the opportunity to make some other small changes. The most visible change is that the probability labels on each axis default to 0 0.25 0.5 0.75 1. If you dislike this choice, or at least prefer something else, you could override that with a call to xlabel() ylabel().

The .hlp is now .sthlp. There is a new do-file that runs the examples in the help. That includes the example in #9 in this thread. In short, the help for ppplot now suggests that you use something else instead.

Thanks again to Chen for the problem report!