Comparison to other software • eventglm

R code and results

library(eventglm)

colon.cifit <- cumincglm(Surv(time, status) ~ rx, time = 2500, data = colon)
summary(colon.cifit)
#> 
#> Call:
#> cumincglm(formula = Surv(time, status) ~ rx, time = 2500, data = colon)
#> 
#> Coefficients:
#>             Estimate Std. Error z value Pr(>|z|)    
#> (Intercept)  0.54345    0.02946  18.449  < 2e-16 ***
#> rxLev       -0.02907    0.04173  -0.697  0.48596    
#> rxLev+5FU   -0.13176    0.04186  -3.148  0.00165 ** 
#> ---
#> Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
#> 
#> (Dispersion parameter for quasi family taken to be 1)
#> 
#>     Null deviance: 253.10  on 928  degrees of freedom
#> Residual deviance: 250.15  on 926  degrees of freedom
#> AIC: NA
#> 
#> Number of Fisher Scoring iterations: 2
confint(colon.cifit)
#>                  2.5 %      97.5 %
#> (Intercept)  0.4857153  0.60118745
#> rxLev       -0.1108637  0.05271372
#> rxLev+5FU   -0.2137964 -0.04971921

Stata code and results

This uses the st0202_1 package, available from here: https://www.stata-journal.com/article.html?article=st0202_1

. import delimited "colon.csv", clear
(18 vars, 929 obs)

. 
. stset time, failure(status==1)

     failure event:  status == 1
obs. time interval:  (0, time]
 exit on or before:  failure

------------------------------------------------------------------------------
        929  total observations
          0  exclusions
------------------------------------------------------------------------------
        929  observations remaining, representing
        452  failures in single-record/single-failure data
  1,551,389  total analysis time at risk and under observation
                                                at risk from t =         0
                                     earliest observed entry t =         0
                                          last observed exit t =     3,329

. 
. // requires st0202_1 install (search stpci)
. stpci, at(2500)
Pseudo-observations for the cumulative incidence function.
Competing risks: (none).
Computing pseudo-observations (progress dots indicate percent completed).
----+--- 1 ---+--- 2 ---+--- 3 ---+--- 4 ---+--- 5 
..................................................    50
..................................................   100
Generated variable: pseudo.

. tabulate rx, gen(rxdum)

         rx |      Freq.     Percent        Cum.
------------+-----------------------------------
        Lev |        310       33.37       33.37
    Lev+5FU |        304       32.72       66.09
        Obs |        315       33.91      100.00
------------+-----------------------------------
      Total |        929      100.00

. glm pseudo rxdum1 rxdum2, vce(robust)

Iteration 0:   log pseudolikelihood =  -708.7556  

Generalized linear models                         Number of obs   =        929
Optimization     : ML                             Residual df     =        926
                                                  Scale parameter =    .270145
Deviance         =   250.154308                   (1/df) Deviance =    .270145
Pearson          =   250.154308                   (1/df) Pearson  =    .270145

Variance function: V(u) = 1                       [Gaussian]
Link function    : g(u) = u                       [Identity]

                                                  AIC             =   1.532305
Log pseudolikelihood = -708.7556003               BIC             =   -6078.23

------------------------------------------------------------------------------
             |               Robust
      pseudo |      Coef.   Std. Err.      z    P>|z|     [95% Conf. Interval]
-------------+----------------------------------------------------------------
      rxdum1 |   -.029075   .0416186    -0.70   0.485    -.1106459    .0524959
      rxdum2 |  -.1317578   .0417443    -3.16   0.002    -.2135752   -.0499404
       _cons |   .5434514     .02938    18.50   0.000     .4858676    .6010352
------------------------------------------------------------------------------

SAS code and results

Assuming you have loaded the colon data in your workspace.

proc lifetest data=colon noprint plots=none timelist=2500 reduceout outsurv=sall;
time time*status(0);
run;

data sall;
set sall;
theta = survival;
keep theta;
run;

data sout;
set colon;
keep id;
run;

%macro pseudosurv;

%do ip=1 %to 929;
data coloni;
set colon;
where id ^= &ip;
run;

proc lifetest data=coloni noprint plots=none timelist=2500 reduceout outsurv=salli;
time time*status(0);
run;

data salli;
set salli;
thetamini = survival;
id = &ip;
keep id thetamini;
run;

data souti;
merge salli sall;
run;

data sout; 
merge sout souti;
by id;
run;
%end;
%mend pseudosurv;

%pseudosurv;


data sout2;
set sout;
pseudoci = 1 - (929 * theta - (929 - 1) * thetamini);
run;

data colon2;
merge colon sout2;
by id;
if rx='Lev' then rxlev = 1;
else rxlev = 0;
if rx='Lev+5FU' then rxlevplus = 1;
else rxlevplus = 0;
run;

proc reg data = colon2;
model pseudoci = rxlev rxlevplus / white;
run;

Parameter Estimates
Variable	DF	Parameter Estimate	Standard Error	t Value	Pr > \|t\|	Heteroscedasticity Consistent
Variable	DF	Parameter Estimate	Standard Error	t Value	Pr > \|t\|	Standard Error	t Value	Pr > \|t\|
Intercept	1	0.54345	0.02928	18.56	<.0001	0.02936	18.51	<.0001
rxlev	1	-0.02907	0.04158	-0.70	0.4846	0.04160	-0.70	0.4847
rxlevplus	1	-0.13176	0.04179	-3.15	0.0017	0.04172	-3.16	0.0016