Grey's test

[Andrzej Niemierko]

Has Grey's test (to compare cumulative incidence curves) been implemented in Stata? Or, is there an alternative to Grey's test implemented in Stata?
Andrzej

[Roman Mostazir]

May be 'stcrreg' is the option you are looking for. See

help stcrreg

[Andrzej Niemierko]

Thanks, Roman. But I don't see an alternative to `lrtest' after using `stcrreg'.

[Roman Mostazir]

I don't do survival analysis much therefore, can't help much in this regard. Hope someone who is experienced in time to event data analysis will come up to suggest you something better. But a quick look at the manual at pag-9 reveals the following:

The log pseudolikelihood is used as a maximization criterion to obtain parameter estimates, but is not representative of the distribution of the data. For this reason, likelihood-ratio (LR) tests (the lrtest command) are not valid after stcrreg. Use Wald tests (the test command) instead.

[Allen Buxton]

I have not think that Stata will calculate Gray's test. To have Stata keep data transfer and the running of R organized helps me, though nothing is returned in r().

Two files go in the personal directory, PERSONAL. (di `"`=c(sysdir_personal)'"')
1) stgrays.ado and
2) stgrays.R. (attached as stgrays_Rprogram.txt)

************************************************** *****************************
The ado file listed, stgrays.ado, takes as arguments of any name, however the meaning must be ordered as
1) group - describe the groups to be compared,
2) status - 0=censored,
3) time.
stgrays <varlist>

A three variable Stata dataset is created C:\data\stgraydt.dta.
This dataset is read into R where Gray's test is done and the results are typed in the Stata log / output.

The Stata lines (stgrays.ado) may need revision:
local using `"C:\ado\personal\stgrays.R"'
local rpath `"C:\Program Files\R\R-3.1.0\bin\x64\r"'

The R line may also need revision for the filename:
dta1<-read.dta("C:\\data\\stgraydt.dta",convert.dates=TR UE,convert.factors=TRUE,missing.type=FALSE,convert .underscore=TRUE,warn.missing.labels=FALSE)

The foreign package in R reads a Stata version 12 dta file.

please note that this is inspired by rsource * Newson, R. (2007). *
Rsource: Stata module to run r from inside stata using an r source file. *
************************************************** *****************************
* example

use "http://www.stata-journal.com/software/sj4-2/st0059/prostatecancer.dta",clear
label define treatment 0 "0:not" 1 "1:trt"
label val treatment treatment
stset time,f(status==1)
stgrays treatment status _t

* with stcrreg -similar test
stcrreg i.treatment, compete(status==2/3)
stset time,f(status==2)
stcrreg i.treatment, compete(status==1 3)
stset time,f(status==3)
stcrreg i.treatment, compete(status==1/2)

[Allen Buxton]

I do not think that Stata will calculate Gray's test.

[Allen Buxton]

Here is an update for the posting #5 (18 May 2015, 12:39 )

To have Stata keep data transfer and the running of R organized helps me, though nothing is returned in r().

Two files go in the personal directory, PERSONAL. (di `"`=c(sysdir_personal)'"')
1) stgrays.ado and
2) stgrays.R. (attached as stgrays_Rprogram.txt)

************************************************** *****************************
The ado file listed, stgrays.ado, takes as arguments of any name, however the meaning must be ordered as
1) group - describe the groups to be compared,
2) status - 0=censored,
3) time.
stgrays <varlist>

A three variable Stata dataset is created C:\data\stgraydt.dta.
This dataset is read into R where Gray's test is done and the results are typed in the Stata log / output.

The Stata lines (stgrays.ado) may need revision:
local using `"C:\ado\personal\stgrays.R"'
local rpath `"C:\Program Files\R\R-3.5.3\bin\x64\r"'

The R line may also need revision for the filename:
dta1<-read.dta("C:\\data\\stgraydt.dta",convert.dates=TR UE,convert.factors=TRUE,missing.type=FALSE,convert .underscore=TRUE,warn.missing.labels=FALSE)

The foreign package in R reads a Stata version 12 dta file.

please note that this is inspired by rsource * Newson, R. (2007). *
Rsource: Stata module to run r from inside stata using an r source file. *
************************************************** *****************************
* example

use "http://www.stata-journal.com/software/sj4-2/st0059/prostatecancer.dta",clear
label define treatment 0 "0:not" 1 "1:trt"
label val treatment treatment
stset time,f(status==1)
stgrays treatment status _t

* with stcrreg -similar test
stcrreg i.treatment, compete(status==2/3)
stset time,f(status==2)
stcrreg i.treatment, compete(status==1 3)
stset time,f(status==3)
stcrreg i.treatment, compete(status==1/2)

[Allen Buxton]

Here is another update with an alternative way to read Stata .dta version 14 into the R program with better handling of unlabeled factors (if any).

To have Stata keep data transfer and the running of R organized helps me, though nothing is returned in r().

Two files go in the personal directory, PERSONAL. (di `"`=c(sysdir_personal)'"')
1) stgrays.ado and
2) stgrays.R. (attached as stgrays_Rprogram.txt)

************************************************** *****************************
The ado file listed, stgrays.ado, takes as arguments of any name, however the meaning must be ordered as
1) group - describe the groups to be compared,
2) status - 0=censored,
3) time.
stgrays <varlist>

A three variable Stata dataset is created C:\data\stgraydt.dta.
This dataset is read into R where Gray's test is done and the results are typed in the Stata log / output.

The Stata lines (stgrays.ado) may need revision:
local using `"C:\ado\personal\stgrays.R"'
local rpath `"C:\Program Files\R\R-3.5.3\bin\x64\r"'

The R program line 9 may also need revision for the filename:
dta1<-as_factor(as_tibble(read_dta("C:\\data\\stgraydt.d ta"), "labels"))

This reads a Stata version 14 dataset (or proir Stata version) and requires R functions from the haven and tibble packages (R program lines 5 to 7):
autoload("read_dta","haven") #load - to read Stata
autoload("as_factor","haven") #load - to read Stata
autoload("as_tibble","tibble") #load - to read Stata

please note that this is inspired by rsource * Newson, R. (2007). *
Rsource: Stata module to run r from inside stata using an r source file. *
************************************************** *****************************
* example

use "http://www.stata-journal.com/software/sj4-2/st0059/prostatecancer.dta",clear
label define treatment 0 "0:not" 1 "1:trt"
label val treatment treatment
stset time,f(status==1)
stgrays treatment status _t
* Even if the event status is not fully-factored (all levels labeled) the R results, as coded, ///
  will yield correct results:
clonevar status0 = status
label val status0 .
stgrays treatment status0 _t

* with stcrreg -similar test
stcrreg i.treatment, compete(status==2/3)
stset time,f(status==2)
stcrreg i.treatment, compete(status==1 3)
stset time,f(status==3)
stcrreg i.treatment, compete(status==1/2)

[Allen Buxton]

I attached the .ado and R program.

[Allen Buxton]

Here is another update with an option for providing a variable for a stratified Gray's test. As previously noted, this reads Stata dta into the R program with correct handling of unlabeled factors levels (if any).

To have Stata keep data transfer and the running of R organized helps me, though nothing is returned in r().

Two files go in the personal directory, PERSONAL. (di `"`=c(sysdir_personal)'"')
1) stgrays.ado and
2) stgrays.R. (attached as stgraysRprogram.txt)

************************************************** *****************************
The ado file listed, stgrays.ado, takes as arguments of any name, however the meaning must be ordered as
1) group - describe the groups to be compared,
2) status - 0=censored,
3) time,
4) optionally, strata.

stgrays <varlist> [ , strata( <varname> ) ]

A four variable Stata dataset is created C:\data\stgraydt.dta.
This dataset is read into R where Gray's test is done and the results are typed in the Stata log / output.

The Stata lines (stgrays.ado) may need revision:
local using `"C:\ado\personal\stgrays.R"'
local rpath `"C:\Program Files\R\R-3.5.3\bin\x64\r"'

The R program line 9 may also need revision for the filename:
dta1<-as_factor(as_tibble(read_dta("C:\\data\\stgraydt.d ta"), "labels"))

This reads a Stata version 14 dataset (or prior Stata version) and requires R functions from the haven and tibble packages (R program lines 5 to 7):
autoload("read_dta","haven") #load - to read Stata
autoload("as_factor","haven") #load - to read Stata
autoload("as_tibble","tibble") #load - to read Stata

please note that this is inspired by rsource * Newson, R. (2007). *
Rsource: Stata module to run r from inside stata using an r source file. *
************************************************** *****************************
* example

use "http://www.stata-journal.com/software/sj4-2/st0059/prostatecancer.dta",clear
label define treatment 0 "0:not" 1 "1:trt"
label val treatment treatment
stset time,f(status==1)
stgrays treatment status _t

* with strata - as example
* stgrays treatment status _t if !missing(AG),strata(AG)

* Even if the event status is not fully-factored (all levels labeled) the R results, as coded, ///
  will yield correct results:
clonevar status0 = status
label val status0 .
stgrays treatment status0 _t

* with stcrreg -similar test
stcrreg i.treatment, compete(status==2/3)
stset time,f(status==2)
stcrreg i.treatment, compete(status==1 3)
stset time,f(status==3)
stcrreg i.treatment, compete(status==1/2)

[Allen Buxton]

Here is an update to stgrays.ado to calculate Gray's test. This also saves in Stata datasets CI estimates from the R cmprsk package.
stgraystest.dta - Gray's test results as applicable
stgrayslist.dta - Estimated CI, se, and confidence intervals at specific timepoints.
stgraysplot.dta - Estimated CI, se, and confidence intervals at all timepoints for plotting CI curves in Stata.

Two files go in the personal directory, PERSONAL.

di `"`=c(sysdir_personal)'"'

1) stgrays.ado and
2) stgrays.R. (attached as stgraysRprogram.txt)

. which stgrays

C:\Users\abuxton\ado\personal\stgrays.ado
*! stgrays abuxton sep2023 (R-4.2.3) attn lines 12, 13 and line 14 rpath
*! stgrays varlist [if] [in] [ , strata(varname) at(numlist) level(#) ]
*! varlist(exactly 3 variables) for in order 'varx status time'
*! censored status must be equal to 0.

Please edit line 14 of stgrays.ado to define the R path rpath & lines 12 and 13 as preferred.

scalar `sc_stgraydt'="C:\data\stgraydt.dta"                 // line 12 for working dataset
scalar `sc_stgrayR' ="`=c(sysdir_personal)'stgrays.R"       // line 13 for path to R program
scalar `sc_rpath'   ="C:\Users\abuxton\R\R-4.2.3\bin\x64\r" // line 14 for path to R exe

The stgrays.R code sets the working directory in R call to:
setwd("C:/data") ## line 9 for working dataset location and the three saved datasets noted above.


************************************************** *****************************
The ado file listed, stgrays.ado, takes as arguments of any name, however the meaning must be ordered as
1) group - describe the groups to be compared,
2) status - 0=censored,
3) time,
4) optionally, strata.

stgrays <varlist> [if] [in] [ , strata(varname) at(numlist) level(#) ]

A four variable Stata dataset is created: C:\data\stgraydt.dta.
This dataset is read into R where Gray's test is done and the results are typed in the Stata log / output.

The Stata lines (stgrays.ado) may need revision:

scalar `sc_stgraydt'="C:\data\stgraydt.dta"                 // line 12 for working dataset
scalar `sc_stgrayR' ="`=c(sysdir_personal)'stgrays.R"       // line 13 for path to R program
scalar `sc_rpath'   ="C:\Users\abuxton\R\R-4.2.3\bin\x64\r" // line 14 for path to R exe

The R program line 9 may also need revision for the filename:
setwd("C:/data") ## line 9 for working dataset location and the three saved datasets noted above.

This reads a Stata version 14 dataset (or prior Stata version) and requires R functions from the haven, tibble, and tidyverse packages (R program lines 5 to 9):
lines 5 to 9:
autoload("read_dta","haven") #load - to read Stata
autoload("as_factor","haven") #load - to read Stata
autoload("as_tibble","tibble") #load - to read Stata
library(tidyverse)

The R program creates or overwrites three datasets
"C:\data\stgraystest.dta"
"C:\data\stgrayslist.dta"
"C:\data\stgraysplot.dta"

please note that this is inspired by rsource * Newson, R. (2007). *
Rsource: Stata module to run r from inside stata using an r source file. *
************************************************** *****************************
* example

frame reset
use "http://www.stata-journal.com/software/sj4-2/st0059/prostatecancer.dta",clear
label define treatment 0 "0:not" 1 "1:trt"
label val treatment treatment
stset time,f(status==1)

frame copy default dtaA, replace
stgrays treatment status _t , at(12(12)60)

* with strata - as example
*frame copy dtaA default , replace
*stgrays treatment status _t if !missing(AG),strata(AG) at(12(12)60)

cwf default
frame default : use "C:\data\stgraystest.dta" , clear
format status %40.0g
list , noobs

frame default : use "C:\data\stgrayslist.dta" , clear
format group_status %40.0g
bysort group_status: list time est se lci* uci* if time<=60, noobs

frame default : use "C:\data\stgraysplot.dta" , clear
drop if time>60
#delimit;
twoway
(connected est time if inlist(group_status,1), sort connect(stairstep) msymbol(none) lpattern("dash") lcolor(black)  lwidth(thin))
(connected est time if inlist(group_status,2), sort connect(stairstep) msymbol(none) lpattern("solid") lcolor(black) lwidth(thin))
(connected est time if inlist(group_status,3), sort connect(stairstep) msymbol(none) lpattern("dash") lcolor(red) lwidth(thin))
(connected est time if inlist(group_status,4), sort connect(stairstep) msymbol(none) lpattern("solid") lcolor(red) lwidth(thin))
(connected est time if inlist(group_status,5), sort connect(stairstep) msymbol(none) lpattern("dash") lcolor(blue) lwidth(thin))
(connected est time if inlist(group_status,6), sort connect(stairstep) msymbol(none) lpattern("solid") lcolor(blue) lwidth(thin))
,
xlabel(0(6)60, nogrid labsize(*1.2))
ylabel(0(0.05)0.35, nogrid labsize(*1.2) angle(0))
ytitle("cumulative incidence")
xtitle("time to event")
graphregion(style(none) margin(medium) fcolor(white) ifcolor(white))
plotregion(style(none) margin(small) ifcolor(white))
title(`""', color(black) size(medium))
legend(label(1 "not Cancer") label(2 "trt Cancer") label(3 "not CVD") label(4 "trt CVD") label(5 "not Other") label(6 "trt Other")
       rows(2) bm(zero) region(lw(none)))
note(`"R cmprsk"')
scheme(tufte)
;
#delimit cr
cwf dtaA

[Allen Buxton]

Here is a log output listing (smcl) of the example code and an image of the cumulative incidence plot.

[Allen Buxton]

Here is an example on the use of the cuminc function of Robert Gray's cmprsk R package from within in Stata (stgrays.ado re: recent posting-replies #11 and #12). It is featured in section 5.3 of Competing Risks A Practical Perspective by Melania Pintilie, John Wiley & Sons, 2006.