Comparisons of means at timepoints (beginner question)

[Karen Gould]

Hi all,

I have fit the model below, and just wanted to make sure of the correct way to do a certain set of pairwise comparisons. I want to see whether the average on the outcome variable (adhdsev) was significantly higher than at the starting point (TimeWeight=0), for each of the four remaining timepoints (TimeWeight 1-4). I'd like to do this separately by i.adhdsubtype (3 possible subtypes, for a total of 12 pairwise comparisons). I would also like to do the same for overall means (not broken down by subtype, 4 pairwise comparisons). I was able to run margins with confidence intervals but I thought something like pwcompare might be necessary in order to do multiple comparisons. However I'm uncertain as to how to specify the levels for pwcompare for a continuous variable.

I also wanted to find out if there's any caveats in interpreting comparisons (or margins) given the presence of the additional interaction in the model of i.adhdsubtype#c.TimeWeight#c.TimeWeight. E.g. are the marginal means set for a certain value in the deceleration of the slope that might make interpretation difficult?

. mixed adhdsev c.TimeWeight c.TimeWeight#i.adhdsubtype c.TimeWeight#c.TimeWeight i.adhdsubtype#c.TimeWeight#c.TimeWeight, ///
>            || id: TimeWeight, variance mle covariance(unstructured) ///
>            residuals(independent,t(TimeWeight)), 
Note: t() not required for this residual structure; ignored

Performing EM optimization: 

Performing gradient-based optimization: 

Iteration 0:   log likelihood = -434.58713  
Iteration 1:   log likelihood = -432.68323  
Iteration 2:   log likelihood = -432.66667  
Iteration 3:   log likelihood = -432.66667  

Computing standard errors:

Mixed-effects ML regression                     Number of obs     =        244
Group variable: id                              Number of groups  =         93

                                                Obs per group:
                                                              min =          1
                                                              avg =        2.6
                                                              max =          4

                                                Wald chi2(6)      =      62.99
Log likelihood = -432.66667                     Prob > chi2       =     0.0000

-----------------------------------------------------------------------------------------------------------------
                                        adhdsev |      Coef.   Std. Err.      z    P>|z|     [95% Conf. Interval]
------------------------------------------------+----------------------------------------------------------------
                                     TimeWeight |  -.1340633    .360817    -0.37   0.710    -.8412517     .573125
                                                |
                       adhdsubtype#c.TimeWeight |
         ADHD, Predominantly Innattentive Type  |  -1.540448   .4016991    -3.83   0.000    -2.327764   -.7531322
ADHD, Predominantly Hyperactive-Impulsive Type  |  -2.915497   .9290948    -3.14   0.002    -4.736489   -1.094505
                                                |
                      c.TimeWeight#c.TimeWeight |  -.0286886   .1192213    -0.24   0.810    -.2623581    .2049809
                                                |
          adhdsubtype#c.TimeWeight#c.TimeWeight |
         ADHD, Predominantly Innattentive Type  |   .4004778   .1331424     3.01   0.003     .1395235    .6614322
ADHD, Predominantly Hyperactive-Impulsive Type  |   .7815377   .3276643     2.39   0.017     .1393275    1.423748
                                                |
                                          _cons |   4.752859   .1239529    38.34   0.000     4.509916    4.995802
-----------------------------------------------------------------------------------------------------------------

------------------------------------------------------------------------------
  Random-effects Parameters  |   Estimate   Std. Err.     [95% Conf. Interval]
-----------------------------+------------------------------------------------
id: Unstructured             |
               var(TimeWe~t) |   .2278224   .0975934      .0983927     .527509
                  var(_cons) |   .0480474   .0567667      .0047424    .4867911
         cov(TimeWe~t,_cons) |   .1046244   .0546602     -.0025076    .2117563
-----------------------------+------------------------------------------------
               var(Residual) |   1.463658   .1712223      1.163761    1.840836
------------------------------------------------------------------------------
LR test vs. linear model: chi2(3) = 21.75                 Prob > chi2 = 0.0001

Note: LR test is conservative and provided only for reference.

. 
.                    estimates store quadpredict ,

.                   lrtest final quadpredict

Likelihood-ratio test                                 LR chi2(2)  =     10.69
(Assumption: final nested in quadpredict)             Prob > chi2 =    0.0048

. margins i.adhdsubtype, at(TimeWeight=(0(1)4)) vsquish

Adjusted predictions                            Number of obs     =        244

Expression   : Linear prediction, fixed portion, predict()
1._at        : TimeWeight      =           0
2._at        : TimeWeight      =           1
3._at        : TimeWeight      =           2
4._at        : TimeWeight      =           3
5._at        : TimeWeight      =           4

-------------------------------------------------------------------------------------------------------------------
                                                  |            Delta-method
                                                  |     Margin   Std. Err.      z    P>|z|     [95% Conf. Interval]
--------------------------------------------------+----------------------------------------------------------------
                                  _at#adhdsubtype |
                           1#ADHD, Combined Type  |   4.752859   .1239529    38.34   0.000     4.509916    4.995802
         1#ADHD, Predominantly Innattentive Type  |   4.752859   .1239529    38.34   0.000     4.509916    4.995802
1#ADHD, Predominantly Hyperactive-Impulsive Type  |   4.752859   .1239529    38.34   0.000     4.509916    4.995802
                           2#ADHD, Combined Type  |   4.590107   .2435981    18.84   0.000     4.112664    5.067551
         2#ADHD, Predominantly Innattentive Type  |   3.450137   .1640803    21.03   0.000     3.128545    3.771728
2#ADHD, Predominantly Hyperactive-Impulsive Type  |   2.456148   .5955657     4.12   0.000      1.28886    3.623435
                           3#ADHD, Combined Type  |   4.369978   .3426055    12.76   0.000     3.698483    5.041473
         3#ADHD, Predominantly Innattentive Type  |   2.890993   .2343395    12.34   0.000     2.431696     3.35029
3#ADHD, Predominantly Hyperactive-Impulsive Type  |   1.665134   .7704188     2.16   0.031     .1551414    3.175128
                           4#ADHD, Combined Type  |   4.092471   .4908227     8.34   0.000     3.130477    5.054466
         4#ADHD, Predominantly Innattentive Type  |   3.075428   .3050299    10.08   0.000      2.47758    3.673275
4#ADHD, Predominantly Hyperactive-Impulsive Type  |   2.379819   1.083324     2.20   0.028     .2565441    4.503094
                           5#ADHD, Combined Type  |   3.757588   .9069076     4.14   0.000     1.980081    5.535094
         5#ADHD, Predominantly Innattentive Type  |   4.003441    .508255     7.88   0.000     3.007279    4.999602
5#ADHD, Predominantly Hyperactive-Impulsive Type  |   4.600202    2.19747     2.09   0.036     .2932402    8.907164
-------------------------------------------------------------------------------------------------------------------

. marginsplot, x(TimeWeight)

  Variables that uniquely identify margins: TimeWeight adhdsubtype

[Clyde Schechter]

It sounds like what you want requires nothing more than adding the -pwcompare- option (which is similar to but distinct from the -pwcompare- command) to your -margins- command. The specification of the levels for a continuous variable are just the same with -margins, pwcompare- as for -margins-. Thus:

 
 margins i.adhdsubtype, at(TimeWeight=(0(1)4)) pwcompare

will give you pairwise comparisons of all of the margins you show in your post.

[Karen Gould]

Thanks very much Clyde, seems like this will do just what I wanted.

[Karen Gould]

Is there any way to get p values on these comparisons (either with or without bonferroni correction)? The effects option does not seem to be allowed.

 
  margins i.adhdsubtype, at(TimeWeight=(0(1)4)) pwcompare mcompare(bonferroni)

[Karen Gould]

sorry, figured it out, was using effects as an option instead of a suboption. http://www.stata.com/manuals13/rmarginspwcompare.pdf

[Clyde Schechter]

I'm glad you figured out the solution to your problem.

That said, from the names of your variables, I'm guessing that this is some kind of clinical study of attention deficit hyperactivity disorder and that your outcome variable is some measure of severity. If that is the case, I would think that the actual differences between later times and baseline, by subtype, along with their confidence intervals would be far more useful and interesting than some p-value testing a straw-man null hypothesis. If I were reading or reviewing such a study, I'd much rather know how much the severity changed, and the associated level of precision with which that has been estimated, than see a statistic that fuses your sample size, the precision of your data, and the magnitude of the changes observed into a single jumbled statistic that cannot be disentangled to tell me any of those things separately.

You might want to check out Ronald L. Wasserstein & Nicole A. Lazar (2016): The ASA's statement on p-values: context, process, and purpose, The American Statistician, DOI:10.1080/00031305.2016.1154108 (which you can download from http://dx.doi.org/10.1080/00031305.2016.1154108) to see if (or where) p-values actually serve a real purpose in your study. (The supplementary materials accompanying that article also make for really interesting reading when you have a slow day to devote the time to them.)

[Karen Gould]

Thanks very much for your reply Clyde. Sorry for the late reply as I've only just noticed your post. I am very interested to read the study you've referenced and will dig in to the supplementary materials too when that slow day comes!

For the interaction, I have also plotted the margins with confidence intervals, which I hope illustrates the data more usefully as you've described, as well as reporting beta values and CIs for the subtypes in the interaction terms and overall significance of the interaction.

I am also reporting the slopes overall (not by subtype) at each timepoint, as I'm interested in the rate of change for the group as a whole. I am also interested in reporting the mean differences at each timepoint for the whole sample with confidence intervals (i.e. not broken down by subtype), and I would like to provide a pairwise comparison versus baseline for these to indicate whether overall, severity is significantly better than baseline at each of these times (I think p values will be requested here!). How do I do that for the sample as a whole?

(NB - after a group of participants had to be excluded, the subtype#time interaction is no longer significant in this analysis. I am interested to explore the change over time though. This data comes from a clinical study of another outcome, hence there is no control group data for this particular outcome in case you are wondering! No causal conclusions are to be drawn, but we are interested in describing these changes nonetheless as they may inspire experimental research)

[Clyde Schechter]

Sounds like a good plan.

As for contrasting time periods for the group as a whole, it sounds like you want:

margins, at(TimeWeight = (0(1)4)) pwcompare(pveffects)

[Karen Gould]

Thanks once again.