st: Gray's test. R results in Stata ouput.

[Allen Buxton <ABUXTON@childrensoncologygroup.org>]

Since I have not seen in Statalist or elsewhere that Stata will calculate 
Gray's test, I offer this too long posting. Simply to have Stata keep
data transfer and the running of R organized helps me, though nothing is
returned in r().
This is to accomplish Gray's test comparing cumulative incidence curves. 
The following are some notes and the contents of two files for the PERSONAL: 
directory, in Windows 'c:\ado\personal/'.
 stgrays.ado and stgrays.R. 
Last, is an example of the output listed in the Stata results on typing:
'stgrays treatment status _t'.
re: http://www.stata.com/statalist/archive/2008-05/msg01260.html
re: cmprsk: Subdistribution Analysis of Competing Risks
http://cran.r-project.org/web/packages/cmprsk/index.html
thank you,
Allen
*******************************************************************************
* note                                                                        *
*******************************************************************************
. di `"`=c(sysdir_personal)'"'
c:\ado\personal/
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.
This is written with file names for Windows. C:\... , paths must be revised 
per requirements of the OS, type in Stata: display `"`=c(sysdir_personal)'"'
The Stata lines may need revision:
saveold C:\data\stgraydt.dta,replace
local using    `"C:\ado\personal\stgrays.R"'
local rpath    `"C:\Program Files\R\R-2.12.1\bin\r"'
The R line may also need revision for the filename:
dta1<-read.dta("C:\\data\\stgraydt.dta",convert.dates=TRUE,
  convert.factors=TRUE,missing.type=FALSE,convert.underscore=TRUE,
  warn.missing.labels=FALSE)
*******************************************************************************
* end note                                                                    *
*******************************************************************************
*******************************************************************************
*  c:\ado\personal\stgrays.ado                                                *
*  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.   *
*******************************************************************************
*! stgrays abuxton jan2011 attn lines: saveold, local using, local rpath
*! stgrays varlist [if] [in]
*! varlist(exactly 3 variables) for in order 'varx status time' 
program stgrays,sortpreserve
version 11.1
syntax varlist(min=3 max=3 numeric) [if] [in]
   tempvar touse
   mark `touse' `if' `in'
   qui count if `touse'
   if r(N)==0 {
      di as error `"[if] [in] yield"'
      error 2000 /* no observations */
   } 
*   cap qui count if _st&`touse'
*   if r(N)==0 | _rc!=0 {
*      di as error `"_st is not found please see stset, or [if] [in] yield"'
*      error 2000 /* no observations */
*   } 
*
preserve
qui keep if `touse'
qui keep `varlist'
tempvar group status time
foreach i of numlist 1/3 {
local vname `=word(`"`varlist'"',`i')'
if `i'==1 {
clonevar `group' =`vname'
drop `vname'
}
else if `i'==2 {
clonevar `status'=`vname'
drop `vname'
}
else if `i'==3 {
clonevar `time'=`vname'
drop `vname'
}
}
rename `group' group
rename `status' status
if `"`:value label status'"'==`""' {
label define status 0 "censored"
label val status status
}
else {
label define `:value label status' 0 "censored",modify
}
rename `time' t
saveold C:\data\stgraydt.dta,replace
*rsource using `"C:\ado\personal\stgrays.R"',rpath(C:\Program Files\R\R-2.12.1\bin\r) roptions(--slave)
*Execute Rterm
*di `"`=c(sysdir_personal)'"'
local using    `"C:\ado\personal\stgrays.R"'
local rpath    `"C:\Program Files\R\R-2.12.1\bin\r"'
local roptions `"--slave"'
tempfile tempsource templis
copy `"`using'"' `"`tempsource'"'
local Rcommand `""`rpath'" `roptions' < "`tempsource'" > "`templis'""'
shell `Rcommand'
*List R output
*
  disp _n as text "Beginning of R output from source file: " as result `"`using'"'
  type `"`templis'"'
  disp _n as text "End of R output from source file: " as result `"`using'"'
end
*******************************************************************************
* end c:\ado\personal\stgrays.ado                                             *
*******************************************************************************
#******************************************************************************
#C:\ado\personal\stgrays.R                                                    *
#******************************************************************************
#ref C:\ado\personal\stgrays.R
#ref 23jan2011 use print() and formatC()
rm(list=ls(all=TRUE))
autoload("cuminc", "cmprsk")    #load cuminc
autoload("read.dta","foreign")  #load - to read Stata v8 dataset
#options(digits=5)
                
dta1<-read.dta("C:\\data\\stgraydt.dta",convert.dates=TRUE,convert.factors=TRUE,missing.type=FALSE,convert.underscore=TRUE,warn.missing.labels=FALSE)
#dta1
summary(dta1) #check dataset to see that things look okay
attach(dta1)  #prepare vectors for cuminc function
ftime<-t
fstatus<-status
grp<-group
detach()
status.grp<-table(fstatus,grp) #check crosstabs to verify the data
print(status.grp)
print(margin.table(status.grp,1) )
print(margin.table(status.grp,2))
print(margin.table(table(grp)))
result1<-cuminc(ftime,fstatus,grp,rho=0,cencode="censored",na.action=na.omit)
print(formatC(result1$Tests[,],digits=3,format="f",drop0trailing=TRUE),quote=FALSE)
print(result1,format="fg",digits=4)
#summary(result1)
#******************************************************************************
#end C:\ado\personal\stgrays.R                                                *
#******************************************************************************
*******************************************************************************
*Example Output Re:
*Coviello, V. and Boggess, M. (2004). 
*Cumulative incidence estimation in the presence of competing risks.
*Stata Journal, 4(2):103–112(10). 
*******************************************************************************
. use "http://www.stata-journal.com/software/sj4-2/st0059/prostatecancer.dta";
. label define treatment 0 "0:not" 1 "1:trt"
. label val treatment treatment
. *keep if touse==1
. stset time,f(status==1)
(omitted output)
. stgrays treatment status _t
file C:\data\stgraydt.dta saved
Beginning of R output from source file: C:\ado\personal\stgrays.R
   group          status          t        
 0:not:253   censored:150   Min.   : 0.10  
 1:trt:253   Cancer  :155   1st Qu.:14.25  
             CVD     :141   Median :34.50  
             Other   : 60   Mean   :36.17  
                            3rd Qu.:57.00  
                            Max.   :76.00  
          grp
fstatus    0:not 1:trt
  censored    63    87
  Cancer      91    64
  CVD         63    78
  Other       36    24
fstatus
censored   Cancer      CVD    Other 
     150      155      141       60 
grp
0:not 1:trt 
  253   253 
[1] 506
       stat  pv    df
Cancer 6.216 0.013 1 
CVD    2.655 0.103 1 
Other  2.274 0.132 1 
Tests:
           stat         pv df
Cancer 6.215563 0.01266321  1
CVD    2.655360 0.10320136  1
Other  2.273801 0.13157686  1
Estimates and Variances:
$est
                  20     40      60
0:not Cancer 0.15810 0.3004 0.34214
1:trt Cancer 0.12253 0.1937 0.23874
0:not CVD    0.11462 0.2055 0.24828
1:trt CVD    0.14625 0.2609 0.31066
0:not Other  0.04743 0.0909 0.13729
1:trt Other  0.04348 0.0830 0.09195
$var
                    20        40        60
0:not Cancer 0.0005278 0.0008326 0.0009011
1:trt Cancer 0.0004263 0.0006192 0.0007478
0:not CVD    0.0004025 0.0006471 0.0007543
1:trt CVD    0.0004949 0.0007639 0.0008605
0:not Other  0.0001793 0.0003280 0.0004843
1:trt Other  0.0001650 0.0003020 0.0003367
End of R output from source file: C:\ado\personal\stgrays.R
*******************************************************************************

*
*   For searches and help try:
*   http://www.stata.com/help.cgi?search
*   http://www.stata.com/support/statalist/faq
*   http://www.ats.ucla.edu/stat/stata/