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)
#>
#> Deviance Residuals:
#> Min 1Q Median 3Q Max
#> -0.5875 -0.4902 -0.3467 0.4863 2.1103
#>
#> 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.04971921This 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
------------------------------------------------------------------------------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 | ||
|
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 |