---
title: "Comparison to other software"
output: rmarkdown::html_vignette
vignette: >
%\VignetteIndexEntry{Comparison to other software}
%\VignetteEngine{knitr::rmarkdown}
%\VignetteEncoding{UTF-8}
---
```{r, include = FALSE}
knitr::opts_chunk$set(
collapse = TRUE,
comment = "#>"
)
```
## R code and results
```{r setup}
library(eventglm)
colon.cifit <- cumincglm(Surv(time, status) ~ rx, time = 2500, data = colon)
summary(colon.cifit)
confint(colon.cifit)
```
## Stata code and results
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;
```
<div align="center">
<table class="table" cellspacing="0" cellpadding="5" rules="all" frame="box" bordercolor="#C1C1C1" summary="Procedure Reg: Parameter Estimates">
<colgroup>
<col>
<col>
</colgroup>
<colgroup>
<col>
<col>
<col>
<col>
<col>
<col>
<col>
</colgroup>
<thead>
<tr>
<th class="c b header" colspan="9" scope="colgroup">Parameter Estimates</th>
</tr>
<tr>
<th class="l b header" rowspan="2" scope="col">Variable</th>
<th class="r b header" rowspan="2" scope="col">DF</th>
<th class="r b header" rowspan="2" scope="col">Parameter<br/>Estimate</th>
<th class="r b header" rowspan="2" scope="col">Standard<br/>Error</th>
<th class="r b header" rowspan="2" scope="col">t Value</th>
<th class="r b header" rowspan="2" scope="col">Pr > |t|</th>
<th class="c b header" colspan="3" scope="colgroup">Heteroscedasticity Consistent</th>
</tr>
<tr>
<th class="r b header" scope="col">Standard<br/>Error</th>
<th class="r b header" scope="col">t Value</th>
<th class="r b header" scope="col">Pr > |t|</th>
</tr>
</thead>
<tbody>
<tr>
<th class="l rowheader" scope="row">Intercept</th>
<th class="r data">1</th>
<td class="r data">0.54345</td>
<td class="r data">0.02928</td>
<td class="r data">18.56</td>
<td class="r data"><.0001</td>
<td class="r data">0.02936</td>
<td class="r data">18.51</td>
<td class="r data"><.0001</td>
</tr>
<tr>
<th class="l rowheader" scope="row">rxlev</th>
<th class="r data">1</th>
<td class="r data" nowrap>-0.02907</td>
<td class="r data">0.04158</td>
<td class="r data" nowrap>-0.70</td>
<td class="r data">0.4846</td>
<td class="r data">0.04160</td>
<td class="r data" nowrap>-0.70</td>
<td class="r data">0.4847</td>
</tr>
<tr>
<th class="l rowheader" scope="row">rxlevplus</th>
<th class="r data">1</th>
<td class="r data" nowrap>-0.13176</td>
<td class="r data">0.04179</td>
<td class="r data" nowrap>-3.15</td>
<td class="r data">0.0017</td>
<td class="r data">0.04172</td>
<td class="r data" nowrap>-3.16</td>
<td class="r data">0.0016</td>
</tr>
</tbody>
</table>
</div>