Comparing the difference in case/control medians across three groups

[Nick Vigala]

Hi,

Thank you all for your time. I'm trying to show that there's increasing stratification of an outcome variable between cases and controls as the number of certain genetic variants increases. Basically, I'm looking to test the delta medians across three groups. I can run a median or kruskal-wallis test to compare the medians across the groups, but I'm really trying to test the difference in the medians. If I simply calculate the delta medians, I only end up with three values (case minus control for each group), which is an insufficient number of observations for a statistical test. I would be grateful for any advice!

Thank you very much!

[Carlo Lazzaro]

Nick:
welcome to the list.
I would take a look at http://www.stata-journal.com/sjpdf.h...iclenum=st0007, especially pages 59-61.
An extensive search of SSC (please, see -help SSC-) user-written packages by Roger Newson is another (and related) research option.

[Nick Vigala]

Hi Carlo,

Thanks for the reply! I'm looking for something a little bit different. Newsom's cendif and somersd packages only let me compare two groups at a time and so give me multiple p-values with my data, but I'm looking for a single value concerning the following null hypothesis.

Ho: The difference in medians between cases and controls does not change across groups, where the groups are defined by number of genetic variants possessed.

So I have three groups (by number of variants) with two subgroups each (cases and controls).

Group 1: Cases and controls with no variants
Group 2: Cases and controls with 1 variant
Group 3: Cases and controls with 2 variants

Between cases and controls in each group, there is a difference in the median value of the outcome variable. That difference in the medians increases from Group 1 to Group 2 and again from Group 2 to Group 3. I'd like to know whether that increase is statistically significant. What test is appropriate for this question?

Thanks again for your help!!!!

[Steve Samuels]

You can get a CI for the difference in differences (DID) by bootstrapping, below. As an aside, cendif and somersd will not estimate a difference in two medians. Rather, cendif will compute a Hodges-Lehman median difference: if X is the variable in group 1 and Y is the same variable in group 2, the H-L median difference is the median of X-Y differences in the population of possible pairs of observations, one from each group.

sysuse auto, clear
gen group = rep78
recode group 1/3=1 4=2 5=3

/* Define a program to return medians */
program define qq, rclass
    centile mpg if group==1
    scalar  m1 = r(c_1)
    centile mpg if group==2
    scalar  m2 = r(c_1)
    centile mpg if group==3
    scalar m3 = r(c_1)
    return scalar m1  = m1
    return scalar m2  = m2
    return scalar m3  = m3
end

set seed 701427
#delim;
bootstrap
    m1 =  r(m1)
    m2 =  r(m2)
    m3 =  r(m3)
    d21 = (r(m2)-r(m1))
    d32 = (r(m3)-r(m2))
    DID =  (r(m3)-2*r(m2)+r(m1)),
    ties strata(group) rep(1000) nodots saving(bs01, replace): qq;
#delim cr
estat bootstrap

Excerpted Results from bootstrap

------------------------------------------------------------------------------
------------------------------------------------------------------------------
             |   Observed   Bootstrap                         Normal-based
             |      Coef.   Std. Err.      z    P>|z|     [95% Conf. Interval]
-------------+----------------------------------------------------------------
          m1 |         19   .5939931    31.99   0.000     17.83579    20.16421
          m2 |       22.5   1.627773    13.82   0.000     19.30962    25.69038
          m3 |         30   5.532102     5.42   0.000     19.15728    40.84272
         d21 |        3.5   1.716505     2.04   0.041     .1357116    6.864288
         d32 |        7.5   5.812611     1.29   0.197    -3.892508    18.89251
         DID |          4   6.490542     0.62   0.538    -8.721229    16.72123
------------------------------------------------------------------------------

Results from estat bootstrap:

-----------------------------------------------------------------------------
             |    Observed               Bootstrap
             |       Coef.       Bias    Std. Err.  [95% Conf. Interval]
-------------+----------------------------------------------------------------
          m1 |          19       .095   .59399312          18         20  (BC)
          m2 |        22.5     -.3475   1.6277725          18         25  (BC)
          m3 |          30     -2.357   5.5321022          18         35  (BC)
         d21 |         3.5     -.4425   1.7165052         -.5          6  (BC)
         d32 |         7.5    -2.0095    5.812611          -6         14  (BC)
         DID |           4     -1.567   6.4905422         -10       15.5  (BC)
------------------------------------------------------------------------------
(BC)   bias-corrected confidence interval (adjusted for ties)

These BC confidence intervals are the ones that I would quote. I note that you asked for a hypothesis test . By definition, a p-value for a test is a probability calculated under the null hypothesis \(H_0:\). The z statistics for d21, d32, and DID in the bootstrap are not computed under the null hypothesis, so the P>|z| column does not contain p-values for a hypothesis test per se. You may not go far wrong if you quote them, but I suggest that you check. See Example 3 in the Manual Entry for bootstrap, which shows how to estimate the Achieved Significance Level (ASL) in a two-sample problem,
:

Efron and Tibshirani (1993, 224) describe an alternative to Satterthwaite’s approximation that estimates the ASL by bootstrapping the statistic from the test of equal means. Their idea is to recenter the two samples to the combined sample mean so that the data now conform to the null hypothesis but that the variances within the samples remain unchanged.

You'd need to do an entirely new bootstrap run on the recentered data (d32 = d21).

[Steve Samuels]

Centering the data to get the Achieved Significance Level for the bootstrap test of differences in differences turns out to be easy. I assume that the original medians are naturally ordered, as you say they are: \(m1 \leq m2 \leq m3\). Before bootstrap

scalar midpt = (m1 + m3)/2  //  midpoint
scalar list midpt
replace mpg = mpg  + midpt -m2  if group ==2
table group, c(p50 mpg)  // group 2 median = midpt
----------------------
    group |   med(mpg)
----------+-----------
        1 |         19
        2 |       24.5
        3 |         30
----------------------

[Nick Vigala]

Steve,

Thanks so much for your help! Unfortunately, I'm still having trouble. I have six medians (case and control for each of three groups) and tried running the following version of your program with no success. I received the error "m1 command not found," which doesn't appear when I run your exact program with the sysauto dataset. Any ideas? Thanks!

. /* Define a program to return medians */
. program define aa, rclass
1. centile outcome if group==1
2. scalar m1 = r(c_1)
3. centile outcome if group==2
4. scalar m2 = r(c_1)
5. centile outcome if group==3
6. scalar m3 = r(c_1)
7. centile outcome if group==4
8. scalar m4 = r(c_1)
9. centile outcome if group==5
10. scalar m5 = r(c_1)
11. centile outcome if group==6
12. scalar m6 = r(c_1)
13. return scalar m1 = m1
14. return scalar m2 = m2
15. return scalar m3 = m3
16. return scalar m4 = m4
17. return scalar m5 = m5
18. return scalar m6 = m6
19. end

.
. set seed 701427

. #delim;
delimiter now ;
. bootstrap
> m1 = r(m1)
> m2 = r(m2)
> m3 = r(m3)
> m4 = r(m4)
> m5 = r(m5)
> m6 = r(m6)
> d41 = (r(m4)-r(m1))
> d52 = (r(m5)-r(m2))
> d63 = (r(m6)-r(m3))
> DID = ((r(m6)-r(m3))-(r(m5)-r(m2))-(r(m4)-r(m1)),
> ties strata(group) rep(1000) nodots saving(bs01, replace): aa;
m1 command not found

[Steve Samuels]

Sorry I missed this post, Nick. You must write programs in do files, text files that are fed to Stata, not at the Stata command line. do files can be written and saved in Stata's built in do file editor, which is very capable.. See the help for do and doedit. You can also use certain text editors instead of the built in editor. See, for example https://www.maxmasnick.com/2015/08/12/real-text-editor/ and http://huebler.blogspot.com/2008/04/stata.html

,I don't know how to define a single DID for more than two groups. If there is one treatment group and two controls groups, I think you need two DIDs, one to compare the treatment with each of the two controls. You might save some typing if you define all the differences and DIDs inside the program to save the work of defining them in the bootstrap statement as I did.
[
Please read FAQ 12 which describes how to put code and results inside CODE delimiters. Jiust before the program define definition, I always put the line

capture program drop _all

Otherwise, revisions of the program won't be recognized.

[Roman Mostazir]

I did not understand your null hypothesis, but you may have a better picture of it than I. However, if you are trying to test the median between groups, 'qreg' is perhaps an option. 'qreg' stands for quantile regression and default is 50th percentile/median or you can test any percentile.

Example:

gen double y = (50-10)*runiform()+10 //generate outcome variable

gen group = floor((1-0+1)*runiform()+0) ////generate group variable

data:

list in 501/520,clean

              y   group  
501.    29.8959       1  
502.   29.98235       0  
503.   30.01979       1  
504.   30.13703       0  
505.   30.17403       0  
506.   30.18587       0  
507.   30.27591       1  
508.   30.27891       0  
509.   30.28846       0  
510.   30.31692       0  
511.   30.33744       0  
512.   30.43347       0  
513.   30.44278       1  
514.   30.47454       0  
515.   30.47853       1  
516.   30.56844       0  
517.   30.56936       1  
518.   30.57958       0  
519.   30.62828       0  
520.   30.65757       1  

==================================


H0: There is no difference in median between the two groups

Regression:

qreg y group

Results:


Median regression                                    Number of obs =      1000
  Raw sum of deviations 5013.394 (about 29.879847)
  Min sum of deviations 4983.536                     Pseudo R2     =    0.0060

------------------------------------------------------------------------------
           y |      Coef.   Std. Err.      t    P>|t|     [95% Conf. Interval]
-------------+----------------------------------------------------------------
       group |    3.12377   1.218894     2.56   0.011     .7318799     5.51566
       _cons |   28.24497   .8567016    32.97   0.000     26.56383    29.92611
------------------------------------------------------------------------------

The overall median for outcome is presented in the parenthesis above the result (29.879847) .
Median for Group-0 = 28.25 and the between group median difference is 3.12 of
which significance and 95% CIs are presented.

Search help for 'qreg' for details.

[Steve Samuels]

I agree that qreg would fit the standard regression estimate for a DID: the coefficient of the interaction of Treatment and Period variables. In fact bsqreg would save the trouble of bootstrapping a program.

I did not recommend qreg for two reasons:

1. To get standard errors for the medians and their pairwise differences one would need several qreg models or several lincom statements after the full model. This extra work seems to negate the original advantage.

Perhaps more important:

2. qreg estimates the 50th percentile, not the median. These need not be the same when \(n\) is even. In the following example, the median is 5.5 (middle of two middle observations when n is even); the 50th percentile is 2.

. input y
             y
  1. 1
  2. 2
  3. 9
  4. 10
  5. end
. expand 2 // qreg needs n>4   to work
(4 observations created)
. centile y //

                                                       -- Binom. Interp. --
    Variable |       Obs  Percentile    Centile        [95% Conf. Interval]
-------------+-------------------------------------------------------------
           y |         8         50         5.5               1          10

. qui qreg y
. di "p50 ="  _b[_cons]
p50 =2

The potential discrepancy between median and p50 can also affect the DID calculation with qreg.

sysuse auto, clear
gen group = rep78
recode group 1/2= 1 3=2 4 =3 5 =4
sort group
/* Get DID from centile */
qui forvalues i = 1/4{
  centile mpg if group==`i'
  scalar m`i' = r(c_1)
  }
scalar did_med =(m4 - m3) -(m2 - m1)
/* Set up qreg regression with 0-1 treatment & period variables */
gen treat =  inlist(group,3,4)
gen period = inlist(group,2,4)
egen tag = tag(group)
list group treat period if tag  // check

qui qreg  mpg i.treat##i.period
scalar  did_qreg = _b[1.treat#1.period]
scalar list did_med did_qreg

The results are:

 list group treat period if tag  // check

     +------------------------+
     | group   treat   period |
     |------------------------|
  1. |     1       0        0 |
 11. |     2       0        1 |
 41. |     3       1        0 |
 59. |     4       1        1 |
     +------------------------+
.
. qui qreg  mpg i.treat##i.period

. scalar  did_qreg = _b[1.treat#1.period]

. scalar list did_med did_qreg
   did_med =        6.5
  did_qreg =          8

[Steve Samuels]

I didn't mean to leave the impression that dislike qreg. p50, as computed by qreg is a reasonable target parameter in itself and will be close to the median in most data sets. I should have noted that with qreg can control for exogenous covariates. For Nick's problem (3 groups: 1 treatment, 2 control), qreg has an additional dvantage: one can get individual tests, corrected for multiplicity, that each DID is zero and the joint test that both are zero.

sysuse auto, clear
gen id = _n
set seed 557141
gen group = rep78
/* Create 6 groups */
recode group 1/2= 1 3=2 4 =3 5 =4
/* Create two new groups */
expand 2 if inlist(group,3,4)
bys id: gen ct =_n
replace group = 5 if group==3 & ct==2
replace group = 6 if group==4 & ct==2
replace mpg = rnormal(mpg) if inlist(group,5,6)
replace mpg = mpg +5 if group ==6  


/* Set up qreg regression */

gen     treat = 0 if inlist(group,1,2)
replace treat = 1 if inlist(group,3,4)
replace treat = 2 if inlist(group,5,6)
gen period = inlist(group,2,4,6)


bsqreg  mpg i.treat##i.period, reps(500) // interaction coefficients are DIDs

/* Individual and joint tests that both DIDs are zero */
test 1.treat#1.period 2.treat#1.period, mtest(sidak)