knitr::opts_chunk$set( collapse = TRUE, comment = "#>" )
library(eventglm) colon.cifit <- cumincglm(Surv(time, status) ~ rx, time = 2500, data = colon) summary(colon.cifit) confint(colon.cifit)
This uses the st0202_1
package, available from here: https://www.stata-journal.com/article.html?article=st0202_1
```{stata, include = TRUE, eval = FALSE} . 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. ```{sas, eval = FALSE} 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 |
Any scripts or data that you put into this service are public.
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.