Discrepancy in number of subjects after using stset with inverse probability weights. Does anyone know where I'm going wrong?

[Gianfranco Di Gennaro]

Dear all, I have a study with 3 treatments and a time-to-event outcome measured multiple times over time for each patient.
I am in an observational setting and I calculated the IPW (inverse probability weights) after calculating the propensity scores via a multinomial model.

When I do

stset stset followup [pweight=group_w], id(id) failure(OutYN==0)

STATA in the output correctly tells me that I have 838 subjects.

. stset followup [pweight=group_w], id(id) failure(OutYN==0)  

Survival-time data settings

           ID variable: id
         Failure event: OutYN==0
Observed time interval: (followup[_n-1], followup]
     Exit on or before: failure
                Weight: [pweight=group_w]

--------------------------------------------------------------------------
      2,554  total observations
         20  weights invalid                                PROBABLE ERROR
--------------------------------------------------------------------------
      2,534  observations remaining, representing
        838  subjects
        205  failures in single-failure-per-subject data
   20,208.4  total analysis time at risk and under observation
                                                At risk from t =         0
                                     Earliest observed entry t =         0
                                          Last observed exit t =        36

However, when I go to do my cox model STATA sees each observation as a single patient and says I've got 4,899 patients!

. xi: stcox ib2.AED_4_trattamenti

        Failure _d: OutYN==0
  Analysis time _t: followup
       ID variable: id
            Weight: [pweight=group_w]

(sum of wgt is 14,692.6996538639)
Iteration 0:   log pseudolikelihood = -9614.0512
Iteration 1:   log pseudolikelihood = -9585.6003
Iteration 2:   log pseudolikelihood = -9585.3606
Iteration 3:   log pseudolikelihood = -9585.3606
Refining estimates:
Iteration 0:   log pseudolikelihood = -9585.3606

Cox regression with Breslow method for ties

No. of subjects =       4,899                           Number of obs =  2,534
No. of failures =       1,159
Time at risk    = 116,995.778
                                                        Wald chi2(2)  =   5.05
Log pseudolikelihood = -9585.3606                       Prob > chi2   = 0.0800

                                        (Std. err. adjusted for 838 clusters in id)
-----------------------------------------------------------------------------------
                  |               Robust
               _t | Haz. ratio   std. err.      z    P>|z|     [95% conf. interval]
------------------+----------------------------------------------------------------
AED_4_trattamenti |
    Brivaracetam  |    1.65829   .4509027     1.86   0.063     .9732213     2.82559
      Perampanel  |   1.658218   .3839784     2.18   0.029     1.053263    2.610638
-----------------------------------------------------------------------------------

Do you know where I am wrong?

Thank you so much!
Gianfranco

[Carlo Lazzaro]

Gianfranco:
what if you include the -strata()- option in your -stcox- code?
In addition, please note that the -xi.- prefix looks redundant with -fvvarlist- notation.

[Gianfranco Di Gennaro]

Dear Carlo Lazzaro , thank you.
Anyway, no: -strata()- option does NOT solve the problem.
Thank you in advance. The solution is probably obvious and right under my nose, but I just can't see it!

Please see my data:

* Example generated by -dataex-. For more info, type help dataex
clear
input float(id followup) int Età byte Sesso float(aedprebin AED_4_trattamenti group_w OutYN)
 1         6 39 0 1 1 1.7021646 1
 1        12 39 0 1 1 1.7021646 1
 1 12.866667 39 0 1 1 1.7021646 1
 1        36 39 0 1 1 1.7021646 1
 2         6 28 1 1 3  7.868358 1
 2        12 28 1 1 3  7.868358 1
 2 20.766666 28 1 1 3  7.868358 1
 2        36 28 1 1 3  7.868358 1
 3         6 31 1 1 3  5.920408 1
 3        12 31 1 1 3  5.920408 1
 3      21.3 31 1 1 3  5.920408 1
 3        36 31 1 1 3  5.920408 1
 4         6 31 1 1 2  4.582242 1
 4        12 31 1 1 2  4.582242 1
 4      14.2 31 1 1 2  4.582242 1
 4        36 31 1 1 2  4.582242 1
 5         6 21 0 1 3  5.496475 1
 5        12 21 0 1 3  5.496475 1
 5        24 21 0 1 3  5.496475 1
 5 33.466667 21 0 1 3  5.496475 1
 6         6 28 1 1 3  7.932226 1
 6        12 28 1 1 3  7.932226 1
 6        20 28 1 1 3  7.932226 0
 7         6 27 1 1 3  6.886532 1
 7        12 27 1 1 3  6.886532 1
 7        24 27 1 1 3  6.886532 1
 7 33.466667 27 1 1 3  6.886532 1
 8         6 30 1 1 2  4.560753 1
 8        12 30 1 1 2  4.560753 1
 8 33.466667 30 1 1 2  4.560753 1
 8        36 30 1 1 2  4.560753 1
 9         6 31 0 1 3  6.210071 1
 9        12 31 0 1 3  6.210071 1
 9        24 31 0 1 3  6.210071 1
 9 33.466667 31 0 1 3  6.210071 1
15         6 41 1 1 3  4.759955 1
15        12 41 1 1 3  4.759955 1
15        24 41 1 1 3  4.759955 1
15        36 41 1 1 3  4.759955 1
16         6 22 1 1 3  7.087768 1
16        12 22 1 1 3  7.087768 1
16        24 22 1 1 3  7.087768 1
16        36 22 1 1 3  7.087768 1
17         6 58 1 1 3 11.030135 1
17        12 58 1 1 3 11.030135 1
17        24 58 1 1 3 11.030135 1
17        36 58 1 1 3 11.030135 1
18         6 38 1 1 3  6.868028 1
18        12 38 1 1 3  6.868028 1
18        24 38 1 1 3  6.868028 1
18        36 38 1 1 3  6.868028 1
19         6 22 0 0 2  4.820758 1
19        12 22 0 0 2  4.820758 1
19        24 22 0 0 2  4.820758 1
19 33.466667 22 0 0 2  4.820758 1
20         6 18 1 1 2  5.176184 1
20        12 18 1 1 2  5.176184 1
20        24 18 1 1 2  5.176184 1
20        36 18 1 1 2  5.176184 1
21         6 20 1 0 2 2.8915505 1
21        12 20 1 0 2 2.8915505 1
21        24 20 1 0 2 2.8915505 1
21 33.266666 20 1 0 2 2.8915505 1
22         6 29 0 1 2  7.763869 1
22        12 29 0 1 2  7.763869 1
22        24 29 0 1 2  7.763869 1
22 33.466667 29 0 1 2  7.763869 1
23         6 59 1 0 2  8.511524 1
23        12 59 1 0 2  8.511524 1
23        24 59 1 0 2  8.511524 1
23        36 59 1 0 2  8.511524 1
24         6 19 0 0 2  2.985111 1
24        12 19 0 0 2  2.985111 1
24        24 19 0 0 2  2.985111 1
24 33.466667 19 0 0 2  2.985111 1
25         6 43 1 1 2  7.950534 1
25        12 43 1 1 2  7.950534 1
25        13 43 1 1 2  7.950534 0
26         6 40 0 0 2  2.687587 1
26        12 40 0 0 2  2.687587 1
26        13 40 0 0 2  2.687587 0
27         6 28 1 0 2  5.554121 1
27        12 28 1 0 2  5.554121 1
27        24 28 1 0 2  5.554121 1
27        36 28 1 0 2  5.554121 1
28         6 29 1 0 2  3.227032 1
28        12 29 1 0 2  3.227032 1
28        24 29 1 0 2  3.227032 1
28 33.466667 29 1 0 2  3.227032 1
29         6 27 0 1 1  2.518501 1
29        12 27 0 1 1  2.518501 1
29        24 27 0 1 1  2.518501 1
29 33.466667 27 0 1 1  2.518501 1
30         6 40 1 1 1    2.7129 1
30        12 40 1 1 1    2.7129 1
30        19 40 1 1 1    2.7129 0
31         6 33 1 1 1  2.089566 1
31        12 33 1 1 1  2.089566 1
31        16 33 1 1 1  2.089566 0
32         6 55 1 1 1 2.2411718 1
end
label values AED_4_trattamenti varl2
label def varl2 1 "Brivaracetam", modify
label def varl2 2 "Lacosamide", modify
label def varl2 3 "Perampanel", modify

And then my script:

. stset followup [pweight=group_w], id(id) failure(OutYN==0)  

Survival-time data settings

           ID variable: id
         Failure event: OutYN==0
Observed time interval: (followup[_n-1], followup]
     Exit on or before: failure
                Weight: [pweight=group_w]

--------------------------------------------------------------------------
      2,554  total observations
         20  weights invalid                                PROBABLE ERROR
--------------------------------------------------------------------------
      2,534  observations remaining, representing
        838  subjects
        205  failures in single-failure-per-subject data
   20,208.4  total analysis time at risk and under observation
                                                At risk from t =         0
                                     Earliest observed entry t =         0
                                          Last observed exit t =        36

The simplest Cox-Model withou STRATA. The number of subjects, as you see, gives me 4899.
However I see that the output gives me "Std. err. adjusted for 838 clusters in id".
Can I trust this model (obviously to be developed with the appropriate covariates)?

. xi: stcox ib1.AED_4_trattamenti

        Failure _d: OutYN==0
  Analysis time _t: followup
       ID variable: id
            Weight: [pweight=group_w]

(sum of wgt is 14,692.6996538639)
Iteration 0:   log pseudolikelihood = -9614.0512
Iteration 1:   log pseudolikelihood = -9585.6003
Iteration 2:   log pseudolikelihood = -9585.3606
Iteration 3:   log pseudolikelihood = -9585.3606
Refining estimates:
Iteration 0:   log pseudolikelihood = -9585.3606

Cox regression with Breslow method for ties

No. of subjects =       4,899                           Number of obs =  2,534
No. of failures =       1,159
Time at risk    = 116,995.778
                                                        Wald chi2(2)  =   5.05
Log pseudolikelihood = -9585.3606                       Prob > chi2   = 0.0800

                                        (Std. err. adjusted for 838 clusters in id)
-----------------------------------------------------------------------------------
                  |               Robust
               _t | Haz. ratio   std. err.      z    P>|z|     [95% conf. interval]
------------------+----------------------------------------------------------------
AED_4_trattamenti |
      Lacosamide  |    .603031   .1639691    -1.86   0.063     .3539084    1.027516
      Perampanel  |   .9999571   .2072112    -0.00   1.000     .6661857    1.500954
-----------------------------------------------------------------------------------

And the same model with -strata(id)- option, that collapse coefficients to 1

. xi: stcox ib1.AED_4_trattamenti  , strata(id)

        Failure _d: OutYN==0
  Analysis time _t: followup
       ID variable: id
            Weight: [pweight=group_w]

(sum of wgt is 14,692.6996538639)
Iteration 0:   log pseudolikelihood = -2239.4306
Refining estimates:
Iteration 0:   log pseudolikelihood = -2239.4306

Stratified Cox regression with no ties
Strata variable: id

No. of subjects =       4,899                           Number of obs =  2,534
No. of failures =       1,159
Time at risk    = 116,995.778
                                                        Wald chi2(0)  =      .
Log pseudolikelihood = -2239.4306                       Prob > chi2   =      .

                                        (Std. err. adjusted for 838 clusters in id)
-----------------------------------------------------------------------------------
                  |               Robust
               _t | Haz. ratio   std. err.      z    P>|z|     [95% conf. interval]
------------------+----------------------------------------------------------------
AED_4_trattamenti |
      Lacosamide  |          1  (omitted)
      Perampanel  |          1  (omitted)
-----------------------------------------------------------------------------------

[Carlo Lazzaro]

Gianfranco:
as you write, -strata()- does not explain anything.
The discrepancy you noticed seems to rest on -pweight-, as in the following toy-example:

. stset followup [pweight=group_w], id(id) failure(OutYN==0)

Survival-time data settings

           ID variable: id
         Failure event: OutYN==0
Observed time interval: (followup[_n-1], followup]
     Exit on or before: failure
                Weight: [pweight=group_w]

--------------------------------------------------------------------------
        100  total observations
          0  exclusions
--------------------------------------------------------------------------
        100  observations remaining, representing
         27  subjects
          5  failures in single-failure-per-subject data
        820  total analysis time at risk and under observation
                                                At risk from t =         0
                                     Earliest observed entry t =         0
                                          Last observed exit t =        36

. stcox ib1.AED_4_trattamenti

        Failure _d: OutYN==0
  Analysis time _t: followup
       ID variable: id
            Weight: [pweight=group_w]

(sum of wgt is 538.0457727909088)
Iteration 0:  Log pseudolikelihood = -114.11296
Iteration 1:  Log pseudolikelihood =  -110.4191
Iteration 2:  Log pseudolikelihood = -110.19043
Iteration 3:  Log pseudolikelihood = -110.19014
Refining estimates:
Iteration 0:  Log pseudolikelihood = -110.19014

Cox regression with Breslow method for ties

No. of subjects =       142                             Number of obs =    100
No. of failures =        23
Time at risk    = 4,477.536
                                                        Wald chi2(2)  =   3.40
Log pseudolikelihood = -110.19014                       Prob > chi2   = 0.1823

                                         (Std. err. adjusted for 27 clusters in id)
-----------------------------------------------------------------------------------
                  |               Robust
               _t | Haz. ratio   std. err.      z    P>|z|     [95% conf. interval]
------------------+----------------------------------------------------------------
AED_4_trattamenti |
      Lacosamide  |   .2819603   .3011478    -1.19   0.236     .0347586    2.287252
      Perampanel  |   .1632966   .1704014    -1.74   0.082     .0211223    1.262447
-----------------------------------------------------------------------------------


. bysort id : gen wanted= group_w if _n==1


. total wanted

Total estimation                            Number of obs = 27

--------------------------------------------------------------
             |      Total   Std. err.     [95% conf. interval]
-------------+------------------------------------------------
    wanted   |   142.0355   12.48628      116.3696    167.7014
--------------------------------------------------------------

[Gianfranco Di Gennaro]

Thank you so much, Carlo. This puts me in a tough spot because, frankly, I can't understand how the total number of subjects is being recognized as the sum of the weights. I wonder if, at this point, it would make more sense to simply include the weight as a covariate in the model:

stcox ib1.AED_4_trattamenti group_w

[Carlo Lazzaro]

Gianfranco:
I think you have two options (I assume you want to disseminate your research via a research report and/or an article):
1) using a "simpler" semi-parametric survival regression specification:

stcox ib1.AED_4_trattamenti

2) using (and comparing) the coefficients of both specifications, explaining how, when -pweight- is used, the number of subjects changes.

[Rich Goldstein]

I'm a little confused here - since you have pweights, don't you want to use the "svy" prefix; see

help svy_estimation