3
$\begingroup$

Stata has a 'contrast' function after using a regression on a categorical variable. What this essentially does is compare each level of the categorical variable to the mean of the dependent variable.

So say I am interested in tumor size in 10 different locations in the brain. My research question is:

Does any location in the brain harbor statistically significantly smaller tumors than the rest?

Many would suggest running a binomial regression with locations as the dependents, or logistic regression with locations as a categorical independent. But this does not answer the research question. This answers the question "does any location in the brain harbor small tumors compared to the base in question."

Myself I would probably run 10 different regressions by dichotomizing location, but I also know many would shake their fists at this approach.

So I've looked at using contrast (this is what it's called in stata, it may have a different name in whatever environment you work in or in statistics generally, and is described above). However I'm posting this question because this function is new to me.

My question is essentially: Does this function answer my research question sufficiently? If 3 locations pop out with p values lower than my alpha, can I report these as locations that harbor significantly smaller (or larger) aneurysms than the mean?

Improve this question
$\endgroup$
3
  • $\begingroup$ Are locations ordered in any way? $\endgroup$ – Dimitriy V. Masterov Aug 11 '20 at 20:34
  • $\begingroup$ How it tumor size measured? From your question, it sometimes sounds like it is binary (small versus big, hence logit), but it can often be continuous variable. $\endgroup$ – Dimitriy V. Masterov Aug 11 '20 at 20:46
  • $\begingroup$ Locations are NOT ordered. Tumor size can either be binary (dichotomized) or gradual, it doesn't really matter for my original question. Only difference would be whether the results would be reported as logarithmic means or true means. I am aware gradual data would yield more information. $\endgroup$ – Paze Aug 11 '20 at 23:22
1
$\begingroup$

At first, it seems that "global" contrast with the grand mean does not make sense here, since the questions asks you about "the rest of them." On the other hand, I've also heard that it does not matter (like in the Statalist thread that @JTS365 linked to). I think the intuition for this is that the hypothesis that, say, the first mean is the same as the grand mean is

$$\mu_1 = \frac{1}{K}\sum_{i=1}^{K} \mu_i,$$ which is algebraically equivalent to $$\mu_1 = \frac{1}{K-1}\sum_{i=2}^{K} \mu_i$$ since $$\mu_1 = \frac{1}{K}\sum_{i=1}^{K} \mu_i=\frac{1}{K}\mu_1 + \frac{1}{K}\sum_{i=2}^{K}\mu_i \implies \mu_1 - \frac{1}{K-1}\sum_{i=2}^{K}\mu_i =0.$$

In any case, I did both contrasts below and they turn out the same judging by the F statistic and its denominator. The leave-one-out, user-defined contrasts are a pain to define, so since the grand mean contrasts are functionally the same, this is probably why the nice folks in College Station did not bother to give us a separate canned option for it.

We will use a dataset of pig weights. This will be our surrogate for tumor size, with week standing in for location. We will treat week as an unordered categorical variable in what follows to keep it similar to your problem. The data show a clear trend of growing weight (so we will be likely to reject):

enter image description here

We start with a longitudinal panel from which we sample each pig at some point in his lifecycle (N=48) to get a cross-section:

. /* Date Step */
. webuse pig, clear
(Longitudinal analysis of pig weights)

. xtset id week
       panel variable:  id (strongly balanced)
        time variable:  week, 1 to 9
                delta:  1 unit

. xtdescribe

      id:  1, 2, ..., 48                                     n =         48
    week:  1, 2, ..., 9                                      T =          9
           Delta(week) = 1 unit
           Span(week)  = 9 periods
           (id*week uniquely identifies each observation)

Distribution of T_i:   min      5%     25%       50%       75%     95%     max
                         9       9       9         9         9       9       9

     Freq.  Percent    Cum. |  Pattern
 ---------------------------+-----------
       48    100.00  100.00 |  111111111
 ---------------------------+-----------
       48    100.00         |  XXXXXXXXX

. set seed 08112020

. sample 1, by(id) count // sample a pig in a random week
(384 observations deleted)

. isid id

Now we move on to the model, which will be a simple het-robust regression:

. /* Estimate Model */
. regress weight i.week, robust

Linear regression                               Number of obs     =         48
                                                F(8, 39)          =     132.50
                                                Prob > F          =     0.0000
                                                R-squared         =     0.9381
                                                Root MSE          =     4.6974

------------------------------------------------------------------------------
             |               Robust
      weight |      Coef.   Std. Err.      t    P>|t|     [95% Conf. Interval]
-------------+----------------------------------------------------------------
        week |
          2  |   8.833333   1.467599     6.02   0.000     5.864835    11.80183
          3  |   14.16667    1.71303     8.27   0.000     10.70174     17.6316
          4  |   19.36667   2.430893     7.97   0.000     14.44972    24.28361
          5  |   27.66667   2.165253    12.78   0.000     23.28703     32.0463
          6  |         34   2.093267    16.24   0.000     29.76597    38.23403
          7  |    40.7381   1.935656    21.05   0.000     36.82286    44.65333
          8  |   45.59524   2.466316    18.49   0.000     40.60664    50.58383
          9  |   55.16667   6.708416     8.22   0.000     41.59761    68.73572
             |
       _cons |   23.33333    1.36292    17.12   0.000     20.57657     26.0901
------------------------------------------------------------------------------

Now for some predictions:

. /* Calculate and compare expected weights */
. margins // global mean

Predictive margins                              Number of obs     =         48
Model VCE    : Robust

Expression   : Linear prediction, predict()

------------------------------------------------------------------------------
             |            Delta-method
             |     Margin   Std. Err.      t    P>|t|     [95% Conf. Interval]
-------------+----------------------------------------------------------------
       _cons |   50.84375   .6780132    74.99   0.000     49.47234    52.21516
------------------------------------------------------------------------------

. margins week // expected weight in each week 

Adjusted predictions                            Number of obs     =         48
Model VCE    : Robust

Expression   : Linear prediction, predict()

------------------------------------------------------------------------------
             |            Delta-method
             |     Margin   Std. Err.      t    P>|t|     [95% Conf. Interval]
-------------+----------------------------------------------------------------
        week |
          1  |   23.33333    1.36292    17.12   0.000     20.57657     26.0901
          2  |   32.16667   .5443311    59.09   0.000     31.06565    33.26768
          3  |       37.5   1.037749    36.14   0.000     35.40095    39.59905
          4  |       42.7   2.012882    21.21   0.000     38.62856    46.77144
          5  |         51   1.682489    30.31   0.000     47.59684    54.40316
          6  |   57.33333   1.588778    36.09   0.000     54.11973    60.54694
          7  |   64.07143   1.374486    46.61   0.000     61.29127    66.85159
          8  |   68.92857    2.05552    33.53   0.000     64.77089    73.08625
          9  |       78.5   6.568508    11.95   0.000     65.21394    91.78606
------------------------------------------------------------------------------

First we will do the global mean comparison, where the first four weeks will be below the mean (negative) and then will be above for the next five (positive):

. /* compare expected weight in each week to the global mean */
. margins g.week, mcompare(sidak)

Contrasts of adjusted predictions               Number of obs     =         48
Model VCE    : Robust

Expression   : Linear prediction, predict()

-----------------------------------------------------------
             |                                        Sidak
             |         df           F        P>F        P>F
-------------+---------------------------------------------
        week |
(1 vs mean)  |          1      336.94     0.0000     0.0000
(2 vs mean)  |          1      342.18     0.0000     0.0000
(3 vs mean)  |          1      107.38     0.0000     0.0000
(4 vs mean)  |          1       16.00     0.0003     0.0025
(5 vs mean)  |          1        0.05     0.8242     1.0000
(6 vs mean)  |          1       16.55     0.0002     0.0020
(7 vs mean)  |          1       81.07     0.0000     0.0000
(8 vs mean)  |          1       82.81     0.0000     0.0000
(9 vs mean)  |          1       22.66     0.0000     0.0002
      Joint  |          8      132.50     0.0000
             |
 Denominator |         39
-----------------------------------------------------------
Note: Sidak-adjusted p-values are reported for tests on
      individual contrasts only.

---------------------------
             |    Number of
             |  Comparisons
-------------+-------------
        week |            9
---------------------------

--------------------------------------------------------------
             |            Delta-method           Sidak
             |   Contrast   Std. Err.     [95% Conf. Interval]
-------------+------------------------------------------------
        week |
(1 vs mean)  |  -27.28148   1.486245     -31.63183   -22.93113
(2 vs mean)  |  -18.44815   .9973037     -21.36733   -15.52896
(3 vs mean)  |  -13.11481    1.26561     -16.81935   -9.410277
(4 vs mean)  |  -7.914815   1.978757     -13.70679   -2.122842
(5 vs mean)  |   .3851852   1.722171      -4.65574    5.426111
(6 vs mean)  |   6.718519   1.651497      1.884461    11.55258
(7 vs mean)  |   13.45661   1.494507      9.082078    17.83115
(8 vs mean)  |   18.31376    2.01256      12.42284    24.20468
(9 vs mean)  |   27.88519   5.858465      10.73701    45.03336
--------------------------------------------------------------

Here all but the fifth week are significantly different, and we reject the joint null that they are all the same. This matches what we see in the graph, where the grand mean is the dashed line.

Now we do the manual version where we use the mean of all other weeks instead of the global mean:

. // Compare the weights in each week to the average of all OTHER weeks
. // Here 1/8 = .125, b/c H0: mu_1 - (1/8)(mu_2 + mu_3 + ... + m_9) == 0
. contrast ///
> {week +1.00 -.125 -.125 -.125 -.125 -.125 -.125 -.125 -.125} ///
> {week -.125 +1.00 -.125 -.125 -.125 -.125 -.125 -.125 -.125} ///
> {week -.125 -.125 +1.00 -.125 -.125 -.125 -.125 -.125 -.125} ///
> {week -.125 -.125 -.125 +1.00 -.125 -.125 -.125 -.125 -.125} ///
> {week -.125 -.125 -.125 -.125 +1.00 -.125 -.125 -.125 -.125} ///
> {week -.125 -.125 -.125 -.125 -.125 +1.00 -.125 -.125 -.125} ///
> {week -.125 -.125 -.125 -.125 -.125 -.125 +1.00 -.125 -.125} ///
> {week -.125 -.125 -.125 -.125 -.125 -.125 -.125 +1.00 -.125} ///
> {week -.125 -.125 -.125 -.125 -.125 -.125 -.125 -.125 +1.00} ///
> , effects mcompare(sidak)

Contrasts of marginal linear predictions

Margins      : asbalanced

-----------------------------------------------------------
             |                                        Sidak
             |         df           F        P>F        P>F
-------------+---------------------------------------------
        week |
        (1)  |          1      336.94     0.0000     0.0000
        (2)  |          1      342.18     0.0000     0.0000
        (3)  |          1      107.38     0.0000     0.0000
        (4)  |          1       16.00     0.0003     0.0025
        (5)  |          1        0.05     0.8242     1.0000
        (6)  |          1       16.55     0.0002     0.0020
        (7)  |          1       81.07     0.0000     0.0000
        (8)  |          1       82.81     0.0000     0.0000
        (9)  |          1       22.66     0.0000     0.0002
      Joint  |          8      132.50     0.0000
             |
 Denominator |         39
-----------------------------------------------------------
Note: Sidak-adjusted p-values are reported for tests on
      individual contrasts only.

---------------------------
             |    Number of
             |  Comparisons
-------------+-------------
        week |            9
---------------------------

------------------------------------------------------------------------------
             |                              Sidak                Sidak
             |   Contrast   Std. Err.      t    P>|t|     [95% Conf. Interval]
-------------+----------------------------------------------------------------
        week |
        (1)  |  -30.69167   1.672025   -18.36   0.000    -35.58581   -25.79752
        (2)  |  -20.75417   1.121967   -18.50   0.000    -24.03825   -17.47008
        (3)  |  -14.75417   1.423811   -10.36   0.000    -18.92177   -10.58656
        (4)  |  -8.904167   2.226101    -4.00   0.002    -15.42014   -2.388197
        (5)  |   .4333333   1.937442     0.22   1.000    -5.237708    6.104374
        (6)  |   7.558333   1.857934     4.07   0.002     2.120018    12.99665
        (7)  |   15.13869    1.68132     9.00   0.000     10.21734    20.06004
        (8)  |   20.60298   2.264131     9.10   0.000     13.97569    27.23026
        (9)  |   31.37083   6.590774     4.76   0.000     12.07913    50.66253
------------------------------------------------------------------------------

We no longer get the same differences, but the individual p-values and F stats are all the same as before.

Three other issues are worth mentioning. We are making 9 comparisons, so I adjusted for this with the Sidak correction to the individual contrasts only. I believe this correction has a conservative FWER when contrasts are positively dependent, as they are here.

Second, the pig dataset is pretty balanced, so it's not a big deal that each contrast assumes an equal number of observations in each level of each factor. If, however, our data were not balanced, we might prefer that contrast use the actual cell frequencies from our data in computing the marginal means. You will need to change the multipliers above to be unequal. With margins g.week, you can do that with margins gw.week. I have no idea if this is the case with tumors.

Third, the question arguably calls for a one-sided hypothesis since it asks about locations that are smaller (rather than unequal). This means your null is $H_0: \delta \ge 0$ against $H_a: \delta < 0$. Everything above was two-sided, so we need to divide the p-values by 2 (and/or calculate 90% CIs). This does not really make a difference here.

Code:

cls
/* Date Step */
webuse pig, clear
xtset id week
xtdescribe
set seed 08112020
sample 1, by(id) count // sample a pig in a random week
isid id

/* Estimate Model */
regress weight i.week, robust

/* Calculate and compare expected weights */
margins // global mean
margins week // expected weight in each week 

/* compare expected weight in each week to the global mean */
margins g.week, mcompare(sidak)


// Compare the weights in each week to the average of all OTHER weeks
// Here 1/8 = .125, b/c H0: mu_1 - (1/8)(mu_2 + mu_3 + ... + m_9) == 0
contrast ///
{week +1.00 -.125 -.125 -.125 -.125 -.125 -.125 -.125 -.125} ///
{week -.125 +1.00 -.125 -.125 -.125 -.125 -.125 -.125 -.125} ///
{week -.125 -.125 +1.00 -.125 -.125 -.125 -.125 -.125 -.125} ///
{week -.125 -.125 -.125 +1.00 -.125 -.125 -.125 -.125 -.125} ///
{week -.125 -.125 -.125 -.125 +1.00 -.125 -.125 -.125 -.125} ///
{week -.125 -.125 -.125 -.125 -.125 +1.00 -.125 -.125 -.125} ///
{week -.125 -.125 -.125 -.125 -.125 -.125 +1.00 -.125 -.125} ///
{week -.125 -.125 -.125 -.125 -.125 -.125 -.125 +1.00 -.125} ///
{week -.125 -.125 -.125 -.125 -.125 -.125 -.125 -.125 +1.00} ///
, effects mcompare(sidak)
Improve this answer
$\endgroup$
0
$\begingroup$

See if this link helps (it talks about getting the contrast for each category vs. the other categories):

https://www.statalist.org/forums/forum/general-stata-discussion/general/1461598-getting-the-contrast-for-each-category-vs-the-other-categories

Improve this answer
$\endgroup$

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

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