[Marina Soley]

Dear Joerg:

Sorry about bothering you again. I wonder if you could help me with one last question.

Is it possible to obtain expected costs under treatment and non-treatment for several values of one of the independent variables? In other words, in the predict mu statement, can I include an if statement specifying the values of one of the independent variables?

Thanks so much for all your help.
Marina

[Joerg Luedicke (StataCorp)]

Hi Marina,

You could do this with predict, cmean after stteffects or by fitting the gamma model with glm and then use margins. Using glm will give you the correct point estimates but if you want to estimate standard errors, they may only be valid for the simple regression adjustment model. Let's take a look at an example using our toy data:

clear
set seed 123
set obs 1000
gen double u0 = runiform()
gen double u1 = runiform()

* Gamma parameters:
local s0 = 1/exp(2*-.8)
local s1 = 1/exp(2*-.8)

* Covariates:
gen double x1 = runiform()
gen double x2 = rbeta(2,2)
gen double x3 = rnormal(1,0.5)
gen double x4 = runiform()
gen double x5 = rnormal(1,1)

* Treatment model:
local xb_t (0.5 + 0.8*x1 + 0.8*x2 + 0.4*x3)
gen double xb = 1/(1+exp(-(`xb_t')))
gen byte treat = rbinomial(1,xb)

* Outcome:
gen double xb0 = -1.5 + 0.50*(x1 + x2 + x3 + x4 + x5)
gen double xb1 = -1.0 + 0.90*(x1 + x2 + x3 + x4 + x5)
gen double y0 = invgammap(`s0', u0)*exp(xb0)/`s0'
gen double y1 = invgammap(`s1', u1)*exp(xb1)/`s1'
gen double y = treat*y1 + (1-treat)*y0

* True potential outcomes:
gen double pom0 = exp(xb0)
gen double pom1 = exp(xb1)

* stset outcome:
stset y

Say we wanted to compute hypothetical potential outcomes, evaluated at a number of points over the range of variable x4. Let's start by creating a copy of x4:

clonevar x4_copy = x4

We now fit the regression adjustment gamma model and then use predict, cmean while looping over six different values of x4 (which ranges from 0 to 1)

stteffects ra (x1 x2 x3 x4_copy x5, gamma) (treat)
forval i = 0(0.2)1 {
    replace x4_copy = `i'
    local u : disp %2.1f `i'
    local u = regexr("`u'","\.","_")
    predict cm0_`u' cm1_`u', cmean
}

We now have the observation level predictions for treated and controls, evaluated at some fixed values of x4. These are the averaged predictions:

summarize cm*

    Variable |        Obs        Mean    Std. Dev.       Min        Max
-------------+---------------------------------------------------------
     cm0_0_0 |      1,000    1.124215     .678133   .1635632   7.375252
     cm1_0_0 |      1,000    9.397557    13.58296   .2009165   247.8473
     cm0_0_2 |      1,000    1.274045    .7685117   .1853622   8.358193
     cm1_0_2 |      1,000    11.29478    16.32515   .2414784   297.8839
     cm0_0_4 |      1,000    1.443844    .8709356   .2100665   9.472137
-------------+---------------------------------------------------------
     cm1_0_4 |      1,000    13.57502    19.62095   .2902292    358.022
     cm0_0_6 |      1,000    1.636274    .9870101   .2380633   10.73454
     cm1_0_6 |      1,000    16.31561    23.58212   .3488221   430.3011
     cm0_0_8 |      1,000    1.854349    1.118555   .2697913    12.1652
     cm1_0_8 |      1,000    19.60948    28.34299   .4192439   517.1724
-------------+---------------------------------------------------------
     cm0_1_0 |      1,000    2.101489    1.267631    .305748   13.78652
     cm1_1_0 |      1,000    23.56834      34.065   .5038828   621.5815

If we were only interested in the averaged predictions, we could just fit the above model with glm and then use margins. Here is the regression adjustment model and an estimate of the average treatment effect:

glm y i.treat##c.(x1 x2 x3 x4 x5), family(gamma) link(log) vce(robust)
margins r.treat, vce(unconditional)

Contrasts of predictive margins

Expression   : Predicted mean y, predict()

------------------------------------------------
             |         df        chi2     P>chi2
-------------+----------------------------------
       treat |          1      356.26     0.0000
------------------------------------------------

--------------------------------------------------------------
             |            Unconditional
             |   Contrast   Std. Err.     [95% Conf. Interval]
-------------+------------------------------------------------
       treat |
   (1 vs 0)  |   14.07865   .7458975      12.61672    15.54058
--------------------------------------------------------------

Notice that the results are equivalent to our stteffects results. We can now use margins with the at() option to compute the same averaged predictions as above:

margins i.treat, vce(unconditional) at(x4=(0(0.2)1)) vsquish

Predictive margins                              Number of obs     =      1,000

Expression   : Predicted mean y, predict()
1._at        : x4              =           0
2._at        : x4              =          .2
3._at        : x4              =          .4
4._at        : x4              =          .6
5._at        : x4              =          .8
6._at        : x4              =           1

------------------------------------------------------------------------------
             |            Unconditional
             |     Margin   Std. Err.      z    P>|z|     [95% Conf. Interval]
-------------+----------------------------------------------------------------
   _at#treat |
        1 0  |   1.124215   .0791067    14.21   0.000     .9691683    1.279261
        1 1  |   9.397557   .5517552    17.03   0.000     8.316137    10.47898
        2 0  |   1.274045    .071325    17.86   0.000     1.134251     1.41384
        2 1  |   11.29478   .6073346    18.60   0.000     10.10443    12.48513
        3 0  |   1.443844   .0700351    20.62   0.000     1.306578    1.581111
        3 1  |   13.57502   .6896165    19.68   0.000      12.2234    14.92665
        4 0  |   1.636274   .0835553    19.58   0.000     1.472508    1.800039
        4 1  |   16.31561   .8182514    19.94   0.000     14.71187    17.91936
        5 0  |   1.854349   .1157792    16.02   0.000     1.627426    2.081272
        5 1  |   19.60948   1.017959    19.26   0.000     17.61432    21.60465
        6 0  |   2.101489    .165806    12.67   0.000     1.776515    2.426463
        6 1  |   23.56834   1.316372    17.90   0.000     20.98829    26.14838
------------------------------------------------------------------------------

Now, we could proceed similarly for the inverse-probability weighted regression adjustment model. Say we have the following model:

stteffects ipwra (x4, gamma) (treat x1 x2 x3, logit), vce(robust)

         failure _d:  1 (meaning all fail)
   analysis time _t:  y

Iteration 0:   EE criterion =  3.517e-24  
Iteration 1:   EE criterion =  1.865e-29  

Survival treatment-effects estimation           Number of obs     =      1,000
Estimator      : IPW regression adjustment
Outcome model  : gamma
Treatment model: logit
Censoring model: none
------------------------------------------------------------------------------
             |               Robust
          _t |      Coef.   Std. Err.      z    P>|z|     [95% Conf. Interval]
-------------+----------------------------------------------------------------
ATE          |
       treat |
   (1 vs 0)  |   14.05873   .7849306    17.91   0.000      12.5203    15.59717
-------------+----------------------------------------------------------------
POmean       |
       treat |
          0  |   1.510637   .0945599    15.98   0.000     1.325303    1.695971
------------------------------------------------------------------------------

We can replicate the point estimate of the average treatment effect using logit, glm, and margins:

logit treat x1 x2 x3
predict double ps
gen double ipw = 1.treat/ps + 0.treat/(1-ps)
glm y i.treat##c.(x4) [pw=ipw], family(gamma) link(log) vce(robust)
margins r.treat, vce(unconditional) noweights

Contrasts of predictive margins

Expression   : Predicted mean y, predict()

------------------------------------------------
             |         df        chi2     P>chi2
-------------+----------------------------------
       treat |          1      414.38     0.0000
------------------------------------------------

--------------------------------------------------------------
             |            Unconditional
             |   Contrast   Std. Err.     [95% Conf. Interval]
-------------+------------------------------------------------
       treat |
   (1 vs 0)  |   14.05873   .6906358      12.70511    15.41236
--------------------------------------------------------------

And here are the averaged predictions over the grid of evaluation points:

margins i.treat, vce(unconditional) at(x4=(0(0.2)1)) vsquish

Adjusted predictions                            Number of obs     =      1,000

Expression   : Predicted mean y, predict()
1._at        : x4              =           0
2._at        : x4              =          .2
3._at        : x4              =          .4
4._at        : x4              =          .6
5._at        : x4              =          .8
6._at        : x4              =           1

------------------------------------------------------------------------------
             |            Unconditional
             |     Margin   Std. Err.      z    P>|z|     [95% Conf. Interval]
-------------+----------------------------------------------------------------
   _at#treat |
        1 0  |    1.19795     .18933     6.33   0.000     .8268695     1.56903
        1 1  |    8.87711   1.095072     8.11   0.000     6.730807    11.02341
        2 0  |    1.30876   .1508259     8.68   0.000     1.013146    1.604373
        2 1  |   10.89789   .9718971    11.21   0.000     8.993008    12.80278
        3 0  |    1.42982   .1133837    12.61   0.000     1.207592    1.652048
        3 1  |   13.37868   .8154319    16.41   0.000     11.78047     14.9769
        4 0  |   1.562079   .0987974    15.81   0.000     1.368439    1.755718
        4 1  |    16.4242   .8215155    19.99   0.000     14.81406    18.03434
        5 0  |   1.706571   .1367028    12.48   0.000     1.438638    1.974503
        5 1  |     20.163   1.324372    15.22   0.000     17.56728    22.75872
        6 0  |   1.864429    .216928     8.59   0.000     1.439257      2.2896
        6 1  |    24.7529    2.36804    10.45   0.000     20.11163    29.39417
------------------------------------------------------------------------------

We see that the standard error of the average treatment effect is quite a bit smaller than the one from stteffects. This is generally to be expected because here we ignore the fact that the weights are estimated. If you want to estimate consistent standard errors for the averaged predictions, you may want to bootstrap the whole thing.

I hope this helps,

Joerg

[ye zhang]

Dear Joerg,

I have read your response to a topic regarding "treatment effect estimate with binary treatment and continuous outcome" and it helps me a lot.

Now my study is about survival treatment effects estimate which means "with binary treatment and categorical outcome". When I run command "stteffects ipwra", lots of iterations comes. Then I add "iteration (5)" option at the end of my "stteffects ipwra" and I got results. However, I am not sure if these results are correct.

After I read your comments for "continuous outcome", I also want to check my analysis step by step. That means 1) estimate propensity score using logit model; 2) then calculate weights. 3) estimate outcome for control group and treatment group. However, I am not sure how to estimate outcome after calculating weights. Could you please give me some help about this and how to deal with the "iterations"?

Thanks in advance!

Looking forward to your reply.

Best,

Ye