Notice: On April 23, 2014, Statalist moved from an email list to a forum, based at statalist.org.
| From | Allen Buxton <ABUXTON@childrensoncologygroup.org> |
| To | "statalist@hsphsun2.harvard.edu" <statalist@hsphsun2.harvard.edu> |
| Subject | st: Gray's test. R results in Stata ouput. |
| Date | Mon, 31 Jan 2011 08:16:44 -0800 |
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/