Postestimation using test command following bootstrap

[Vicky Taxiarchi]

Dear all,

Please allow me to ask some advice from you as I have looked everywhere else I have thought of and could not find an answer to my problem.

I use Stata 16 and I have trying to implement constraint longitudinal analysis. Although I have seen that Stata 17 has a constraint option incorporated in the mixed command, this is not available at e earlier option and therefore I follow the methodology introduced at a previous post here: https://www.statalist.org/forums/for...-data-analysis

My dataset is comprised of the outcome variable, time as 0 at baseline and 3 follow up points at 3, 6 and 12 months, and the treatment arm with 2 categories.
My code therefore is:

mixed Outcome i.Treatmentarm##i.time || ParticipantIDA: , covariance(unstructured) level(95)
test 1.Treatmentarm=0, coef


( 1)  [Outcome]1.Treatmentarm = 0

chi2(  1) =    2.37
Prob > chi2 =    0.1234


Constrained coefficients

                
Coef.   Std. Err.    z    P>z    [95% Conf.    Interval]
                
Outcome      
Treatmentarm 
Intervention            0  (omitted)

time 
3    -3.872614   2.191517    -1.77    0.077    -8.167908    .4226801
6    -4.298432   2.108544    -2.04    0.041    -8.431101    -.1657617
12     1.164081   2.083127    0.56    0.576    -2.918774    5.246935

Treatmentarm#time 
Intervention# 3    -.1910133   2.967006    -0.06    0.949    -6.006239    5.624213
Intervention# 6     1.806498   2.809967    0.64    0.520    -3.700937    7.313932
Intervention#12    -4.000512   2.704179    -1.48    0.139    -9.300606    1.299581

_cons    23.07878       1.74    13.26    0.000    19.66845    26.48912
                
lns1_1_1          
_cons    2.414024   .0991465    24.35    0.000    2.2197    2.608347
                
lnsig_e           
_cons    2.141509   .0539483    39.70    0.000    2.035772    2.247246

Due to failure of the linearity assumption, I considered bootstrapping as:

bootstrap, reps(1000): mixed Outcome Treatmentarm##i.time || ParticipantIDA: , covariance(unstructured) 
test 1.Treatmentarm=0, coef

Constrained coefficients

-----------------------------------------------------------------------------------
                  |              Bootstrap
                  |      Coef.   Std. Err.      z    P>|z|     [95% Conf. Interval]
------------------+----------------------------------------------------------------
Outcome |
     Treatmentarm |
    Intervention  |          0  (omitted)
                  |
             time |
               3  |  -4.949825   1.866378    -2.65   0.008    -8.607858   -1.291791
               6  |  -5.096868   1.986886    -2.57   0.010    -8.991094   -1.202642
              12  |   .0389688   2.021361     0.02   0.985    -3.922826    4.000764
                  |
Treatmentarm#time |
 Intervention# 3  |   2.824151   2.720397     1.04   0.299     -2.50773    8.156031
 Intervention# 6  |   4.213531   2.355054     1.79   0.074    -.4022901    8.829353
 Intervention#12  |  -.7735578    2.49915    -0.31   0.757    -5.671801    4.124685
                  |
            _cons |     22.597   1.180223    19.15   0.000     20.28381    24.91019
------------------+----------------------------------------------------------------
lns1_1_1          |
            _cons |   2.421146   .0414222    58.45   0.000      2.33996    2.502332
------------------+----------------------------------------------------------------
lnsig_e           |
            _cons |   2.147335   .0934787    22.97   0.000      1.96412     2.33055

My question here is: Why do the effects obtained following the postestimation test command differ before and after bootstrap? (As expected, and as per Stata's manual, the effects following the initial mixed model before and after bootstrap are the same -not presented here).

Please consider that I tested this to another dataset after a simple regression model as:

 
 regress Outcome i.Country   
 test 2.Country==0, coef

and

 
 bootstrap, reps(1000) seed(45645): regress Outcome i.Country   
 test 2.Country==0, coef

where the variable country has 3 categories. The effects obtained before and after bstrap after the test command, were as before different.

Any advice or any thoughts would be very helpful. Many thanks in advance.

[Andrew Musau]

My question here is: Why do the effects obtained following the postestimation test command differ before and after bootstrap? (As expected, and as per Stata's manual, the effects following the initial mixed model before and after bootstrap are the same -not presented here). Please consider that I tested this to another dataset after a simple regression model as: Code:

regress Outcome i.Country
test 2.Country==0, coef

and

Code:
bootstrap, reps(1000) seed(45645): regress Outcome i.Country
test 2.Country==0, coef

where the variable country has 3 categories. The effects obtained before and after bstrap after the test command, were as before different. Any advice or any thoughts would be very helpful. Many thanks in advance.

Your p values from the test commands are the same as those on your specified coefficients. In other words, you are testing the same hypothesis. So your question is really why bootstrap standard errors differ from conventional standard errors. For this, you have to read up on what bootstrap standard errors are and why you would want to compute them. A good place to start is the PDF manual entry of bootstrap in the Stata manual.

help bootstrap

[Vicky Taxiarchi]

Dear Andrew,

Thank you very much for your time on it and apologies for my late reply.

I am sorry, I don't see where the p-values are the same. For example, the p-value for the interaction of intervention with time=3 without bootstrap is 0.949 and the one after bootstrap is 0.299. Please if you could have a second look at my result and correct me if I am wrong, but I see all the p-values to differ. In addition to that all the effect estimates also differ, which according to the bootstrap manual should not, although this or the bootstrap postestimation manual does not refer to estimation of effects following the test postestimation, rather only for direct estimations within the bootstrap command. But yes, I understand why the standard errors differ and how these are extracted.

Any further advice is more than welcome.

Best wishes
Vicky

[Andrew Musau]

You do not provide a reproducible example, but it is not difficult to show what I mean in #2.

webuse grunfeld, clear
bootstrap, reps(1000) seed(45645) nodots: regress invest i.company 
test 2.company==0

See the highlighted:

Res.:

. bootstrap, reps(1000) seed(45645) nodots: regress invest i.company 

Linear regression                               Number of obs     =        200
                                                Replications      =      1,000
                                                Wald chi2(9)      =    1009.64
                                                Prob > chi2       =     0.0000
                                                R-squared         =     0.7602
                                                Adj R-squared     =     0.7489
                                                Root MSE          =   108.6848

------------------------------------------------------------------------------
             |   Observed   Bootstrap                         Normal-based
      invest |      Coef.   Std. Err.      z    P>|z|     [95% Conf. Interval]
-------------+----------------------------------------------------------------
     company |
          2  |   -197.545   75.75836    -2.61   0.009    -346.0287   -49.06134
          3  |    -505.73   70.48405    -7.18   0.000    -643.8762   -367.5838
          4  |  -521.8965     70.438    -7.41   0.000    -659.9524   -383.8406
          5  |  -546.2175   70.01702    -7.80   0.000    -683.4483   -408.9867
          6  |   -552.609   70.07906    -7.89   0.000    -689.9614   -415.2566
          7  |  -560.4245   70.15277    -7.99   0.000    -697.9214   -422.9276
          8  |  -565.1285   70.07039    -8.07   0.000    -702.4639   -427.7931
          9  |   -566.131   69.87813    -8.10   0.000    -703.0896   -429.1724
         10  |  -604.9355   69.79858    -8.67   0.000    -741.7382   -468.1328
             |
       _cons |     608.02    69.7879     8.71   0.000     471.2382    744.8018
------------------------------------------------------------------------------

. 
. test 2.company==0

 ( 1)  2.company = 0

           chi2(  1) =    6.80
         Prob > chi2 =    0.0091

[Andrew Musau]

I will expand a bit on #2 and #4

1. Bootstrapping does not affect the estimated coefficients, just the standard errors. Compare coefficients in

webuse grunfeld, clear
regress invest i.company
bootstrap, reps(1000) seed(45645) nodots: regress invest i.company

2. The t/z statistic is computed as \(\frac{coefficient}{std. error}.\)

Given that you accept 1 and 2, then differences in t/z statistics and the ensuing levels of significance before and after bootstrap are due to differences in the estimated standard errors.

3. Finally, for a coefficient in a regression, the null hypothesis is that the slope parameter is equal to zero. This is exactly the same hypothesis that you are testing when you run your Wald test. The difference here is that the test statistic that the commands output differ: t/z under the null hypothesis for regress and Chi-square under the null hypothesis for the Wald test. You obtain the same p value as the hypotheses are the same.

4. Taking 1-3, it is apparent that the test results before and after bootstrap differ because the bootstrap standard error differs from the conventional standard error.

[Vicky Taxiarchi]

Dear Andrew,

Thank you so much for keeping on this conversation with me.

Yes, you are absolutely right, the p-value there is the same which of course would be expected.

My interest lies to the coefficients of the interaction terms of intervention and time, which I am sorry that I have not been clear about. I use the test command to constraint for the baseline outcome to be the same between the intervention and the control and examining whether the outcome differs between the two groups at the follow up time points. I'd certainly expect the confidence interval limits to differ before and after bootstrapping given the bootstrap standard error, which is the reason for bootstrapping anyway. But still, what about the coefficients? Do you mean to say that these are affected by the bootstrapped standard errors, or the bootstrap in general, and if so is there any documentation to show this?

Best wishes
Vicky

[Andrew Musau]

My interest lies to the coefficients of the interaction terms of intervention and time, which I am sorry that I have not been clear about. I use the test command to constraint for the baseline outcome to be the same between the intervention and the control and examining whether the outcome differs between the two groups at the follow up time points. I'd certainly expect the confidence interval limits to differ before and after bootstrapping given the bootstrap standard error, which is the reason for bootstrapping anyway. But still, what about the coefficients?

My first point in #5 was

1. Bootstrapping does not affect the estimated coefficients, just the standard errors.

The example using the Grunfeld dataset illustrates this. Therefore, either using your dataset or one of Stata's datasets, you need to provide a reproducible example that replicates your issue. If using your dataset, utilize dataex to provide the data example as recommended in FAQ Advice #12. Otherwise, see

help mixed

to access one of Stata's datasets. If you need to, within the Stata dataset, you can create a 4 level categorical variable and an indicator as follows:

set seed 12152021
gen indicator1 = runiformint(0, 1)
gen indicator2= rnormal(20, 30)>0
gen cat4 = runiformint(1, 4)

[Vicky Taxiarchi]

Dear Andrew,

I totally agree with your point that bootstrap does not affect the estimated coefficients, but please spare another minute to look at the results I have sent already at my initial post.

Please see there that the coef estimates of the interaction of treatment and time after the postestimation are different when bootstrapping is not utilized to when it is. For example, you can see above that the estimate of intervention#3 is -0.19 at the first case and 2.82 at the second. Exactly because I agree with you is the reason I am so confused with the results I got. This change occurs after the postestimation test command, whilst after the initial mixed command the coef estimates are the same as they should be and which is the reason I have not included them here -on top of that I am not interested in them, as my interest is a constraint longitudinal analysis where the baseline score is constrained to be the same between the treatment groups. I have noticed that Stata v17 has the option to constraint within the mixed command, but Stata 16 does not. I wonder whether the results would be as such at Stata 17 but unfortunately I have no means to check it.

Best wishes
Vicky

[Vicky Taxiarchi]

Dear Andrew,

Please consider this question resolved -the solution to this lies on how the test command is run which is outlined at the manual for the respective command. Many thanks for all the time you spent on this.

[Andrew Musau]

You do not provide a reproducible example as I asked in #7, but with some work I am able to replicate your problem and I am at a loss explaining it. Your issue is that the constrained coefficients reported by test after bootstrap differ from the constrained coefficients also reported by test after estimation using the standard maximum-likelihood variance estimator. I do not know how test computes these coefficients or whether it is a bug but if no one else comes up with an answer, send an email to Tech Support and ask. However, I am able to provide a way to get the constrained coefficients and in this way, bootstrap does not appear to behave inconsistently.

webuse pig, clear
set seed 01052022
gen Treat= rnormal(20, 30)>0
gen time = runiformint(1, 4)
mixed weight Treat##i.time ||  id: , covariance(unstructured) level(95)
test 1.Treat=0, coef
bootstrap, reps(1000) nodots: mixed weight Treat##i.time || id: , covariance(unstructured)
test 1.Treat=0, coef

Res.:

. mixed weight Treat##i.time ||  id: , covariance(unstructured) level(95)
Note: single-variable random-effects specification in id equation; covariance structure set to identity

Performing EM optimization:

Performing gradient-based optimization:

Iteration 0:   log likelihood = -1824.8004  
Iteration 1:   log likelihood = -1823.7127  
Iteration 2:   log likelihood = -1823.7111  
Iteration 3:   log likelihood = -1823.7111  

Computing standard errors:

Mixed-effects ML regression                     Number of obs     =        432
Group variable: id                              Number of groups  =         48

                                                Obs per group:
                                                              min =          9
                                                              avg =        9.0
                                                              max =          9

                                                Wald chi2(7)      =       7.06
Log likelihood = -1823.7111                     Prob > chi2       =     0.4228

------------------------------------------------------------------------------
      weight |      Coef.   Std. Err.      z    P>|z|     [95% Conf. Interval]
-------------+----------------------------------------------------------------
     1.Treat |  -2.502247   3.289592    -0.76   0.447    -8.949729    3.945234
             |
        time |
          2  |  -3.697059    3.97019    -0.93   0.352    -11.47849     4.08437
          3  |   -7.92069    4.14016    -1.91   0.056    -16.03525    .1938754
          4  |      -4.55   4.102251    -1.11   0.267    -12.59026    3.490264
             |
  Treat#time |
        1 2  |   5.284363   4.684174     1.13   0.259    -3.896449    14.46518
        1 3  |   8.468525   4.918645     1.72   0.085    -1.171842    18.10889
        1 4  |   1.660581   4.941015     0.34   0.737    -8.023632    11.34479
             |
       _cons |       53.3   2.786931    19.12   0.000     47.83772    58.76228
------------------------------------------------------------------------------

------------------------------------------------------------------------------
  Random-effects Parameters  |   Estimate   Std. Err.     [95% Conf. Interval]
-----------------------------+------------------------------------------------
id: Identity                 |
                  var(_cons) |   4.22e-22   2.37e-21      7.07e-27    2.52e-17
-----------------------------+------------------------------------------------
               var(Residual) |   271.8444   18.49676      237.9048    310.6259
------------------------------------------------------------------------------
LR test vs. linear model: chibar2(01) = 0.00          Prob >= chibar2 = 1.0000

. test 1.Treat=0, coef

 ( 1)  [weight]1.Treat = 0

           chi2(  1) =    0.58
         Prob > chi2 =    0.4469


Constrained coefficients

------------------------------------------------------------------------------
             |      Coef.   Std. Err.      z    P>|z|     [95% Conf. Interval]
-------------+----------------------------------------------------------------
weight       |
     1.Treat |          0  (omitted)
             |
        time |
          2  |  -1.901091    3.19182    -0.60   0.551    -8.156944    4.354762
          3  |  -6.124722   3.400917    -1.80   0.072     -12.7904    .5409534
          4  |  -2.754032   3.354664    -0.82   0.412    -9.329054    3.820989
             |
  Treat#time |
        1 2  |   2.782116   3.334677     0.83   0.404    -3.753731    9.317963
        1 3  |   5.966278   3.656728     1.63   0.103    -1.200776    13.13333
        1 4  |  -.8416667   3.686763    -0.23   0.819    -8.067589    6.384256
             |
       _cons |   51.50403    1.48064    34.78   0.000     48.60203    54.40603
-------------+----------------------------------------------------------------
lns1_1_1     |
       _cons |  -24.60847   2.805316    -8.77   0.000    -30.10679   -19.11015
-------------+----------------------------------------------------------------
lnsig_e      |
       _cons |   2.802615   .0340209    82.38   0.000     2.735935    2.869295
------------------------------------------------------------------------------

.
. bootstrap, reps(1000) nodots: mixed weight Treat##i.time || id: , covariance(unstructured)

Mixed-effects ML regression                     Number of obs     =        432
Group variable: id                              Number of groups  =         48

                                                Obs per group:
                                                              min =          9
                                                              avg =        9.0
                                                              max =          9

                                                Wald chi2(7)      =       7.46
Log likelihood = -1823.7111                     Prob > chi2       =     0.3822

------------------------------------------------------------------------------
             |   Observed   Bootstrap                         Normal-based
      weight |      Coef.   Std. Err.      z    P>|z|     [95% Conf. Interval]
-------------+----------------------------------------------------------------
     1.Treat |  -2.502247   2.985471    -0.84   0.402    -8.353663    3.349169
             |
        time |
          2  |  -3.697059   3.529956    -1.05   0.295    -10.61565    3.221528
          3  |   -7.92069   3.888919    -2.04   0.042    -15.54283   -.2985494
          4  |      -4.55   4.314856    -1.05   0.292    -13.00696    3.906962
             |
  Treat#time |
        1 2  |   5.284363   4.434182     1.19   0.233    -3.406473     13.9752
        1 3  |   8.468525   4.679923     1.81   0.070    -.7039555    17.64101
        1 4  |   1.660581   5.195752     0.32   0.749    -8.522907    11.84407
             |
       _cons |       53.3   2.363008    22.56   0.000     48.66859    57.93141
------------------------------------------------------------------------------

------------------------------------------------------------------------------
                             |   Observed   Bootstrap         Normal-based
  Random-effects Parameters  |   Estimate   Std. Err.     [95% Conf. Interval]
-----------------------------+------------------------------------------------
id: Identity                 |
                  var(_cons) |   4.22e-22   2.56e-21      2.92e-27    6.09e-17
-----------------------------+------------------------------------------------
               var(Residual) |   271.8444   14.12827      245.5171    300.9948
------------------------------------------------------------------------------
LR test vs. linear model: chibar2(01) = 0.00          Prob >= chibar2 = 1.0000

.
. test 1.Treat=0, coef

 ( 1)  [weight]1.Treat = 0

           chi2(  1) =    0.70
         Prob > chi2 =    0.4020


Constrained coefficients

------------------------------------------------------------------------------
             |              Bootstrap
             |      Coef.   Std. Err.      z    P>|z|     [95% Conf. Interval]
-------------+----------------------------------------------------------------
weight       |
     1.Treat |          0  (omitted)
             |
        time |
          2  |  -2.144115   3.004591    -0.71   0.475    -8.033005    3.744775
          3  |  -6.352354   3.409146    -1.86   0.062    -13.03416    .3294488
          4  |  -3.051587   3.927063    -0.78   0.437    -10.74849    4.645315
             |
  Treat#time |
        1 2  |   2.790067   3.287169     0.85   0.396    -3.652666    9.232801
        1 3  |   6.002177   3.639032     1.65   0.099    -1.130194    13.13455
        1 4  |  -.7897427   4.295215    -0.18   0.854     -9.20821    7.628724
             |
       _cons |   51.71689   1.419893    36.42   0.000     48.93395    54.49982
-------------+----------------------------------------------------------------
lns1_1_1     |
       _cons |  -24.69169   3.029039    -8.15   0.000     -30.6285   -18.75489
-------------+----------------------------------------------------------------
lnsig_e      |
       _cons |   2.802405   .0259847   107.85   0.000     2.751476    2.853334
------------------------------------------------------------------------------

I have noticed that Stata v17 has the option to constraint within the mixed command, but Stata 16 does not. I wonder whether the results would be as such at Stata 17 but unfortunately I have no means to check it.

You can estimate the model using meglm which allows constraints, and doing so reveals that you have a valid point as bootstrap does not affect the estimated coefficients.

webuse pig, clear
set seed 01052022
gen Treat= rnormal(20, 30)>0
gen time = runiformint(1, 4)
constraint define 1 _b[1.Treat]=0
meglm weight Treat##i.time ||  id: , constraint(1) nolog
bootstrap, reps(1000) nodots: meglm weight Treat##i.time || id: , constraint(1) nolog

Res.:

. meglm weight Treat##i.time ||  id: , constraint(1) nolog

Mixed-effects GLM                               Number of obs     =        432
Family:                Gaussian
Link:                  identity
Group variable:              id                 Number of groups  =         48

                                                Obs per group:
                                                              min =          9
                                                              avg =        9.0
                                                              max =          9

Integration method: mvaghermite                 Integration pts.  =          7

                                                Wald chi2(6)      =       6.47
Log likelihood = -1824.0002                     Prob > chi2       =     0.3725
 ( 1)  [weight]1.Treat = 0
-------------------------------------------------------------------------------
       weight |      Coef.   Std. Err.      z    P>|z|     [95% Conf. Interval]
--------------+----------------------------------------------------------------
      1.Treat |          0  (omitted)
              |
         time |
           2  |  -1.901091   3.193957    -0.60   0.552    -8.161132     4.35895
           3  |  -6.124722   3.403194    -1.80   0.072    -12.79486    .5454157
           4  |  -2.754032    3.35691    -0.82   0.412    -9.333455    3.825391
              |
   Treat#time |
         1 2  |   2.782116    3.33691     0.83   0.404    -3.758106    9.322339
         1 3  |   5.966278   3.659176     1.63   0.103    -1.205574    13.13813
         1 4  |  -.8416667   3.689231    -0.23   0.820    -8.072426    6.389093
              |
        _cons |   51.50403   1.481631    34.76   0.000     48.60009    54.40798
--------------+----------------------------------------------------------------
id            |
    var(_cons)|   1.46e-30   2.76e-15                             .           .
--------------+----------------------------------------------------------------
 var(e.weight)|   272.2085   18.52144                      238.2236    311.0417
-------------------------------------------------------------------------------

. bootstrap, reps(1000) nodots: meglm weight Treat##i.time || id: , constraint(1) nolog

Mixed-effects GLM                               Number of obs     =        432
                                                Replications      =      1,000
                                                Wald chi2(6)      =       6.43
Log likelihood = -1824.0002                     Prob > chi2       =     0.3771

 ( 1)  [weight]1.Treat = 0
-------------------------------------------------------------------------------
              |   Observed   Bootstrap                         Normal-based
       weight |      Coef.   Std. Err.      z    P>|z|     [95% Conf. Interval]
--------------+----------------------------------------------------------------
      1.Treat |          0  (omitted)
              |
         time |
           2  |  -1.901091   3.013442    -0.63   0.528    -7.807329    4.005147
           3  |  -6.124722    3.41928    -1.79   0.073    -12.82639    .5769433
           4  |  -2.754032   3.943378    -0.70   0.485    -10.48291    4.974847
              |
   Treat#time |
         1 2  |   2.782116   3.286608     0.85   0.397    -3.659518     9.22375
         1 3  |   5.966278   3.638369     1.64   0.101    -1.164793    13.09735
         1 4  |  -.8416667   4.296283    -0.20   0.845    -9.262227    7.578894
              |
        _cons |   51.50403     1.4423    35.71   0.000     48.67718    54.33089
--------------+----------------------------------------------------------------
id            |
    var(_cons)|   1.46e-30   8.727489                             .           .
--------------+----------------------------------------------------------------
 var(e.weight)|   272.2085   13.10288                      247.7015    299.1401
-------------------------------------------------------------------------------

ADDED IN EDIT: These estimations took a while to run, so I missed the response in #9. I am glad that you figured it out. It would be great if you can point out the specific reference in the manual for the benefit of others who will encounter the same problem in the future.