On Tue, Sep 18, 2001 at 04:20:22PM +0800, Nicholas Lewin-Koh
wrote:> Hi,
> I am having trouble with the survival library, particualrily the coxph
> function.
>
> the following works
> coxph(jtree9$cph.call,z,rep(1,dim(z)[1]))
> Call:
> coxph(formula = jtree9$cph.call, data = z, weights = rep(1, dim(z)[1]))
>
>
> coef exp(coef) se(coef) z p
> SM 0.2574 1.294 0.0786 3.274 1.1e-03
> Sex -0.1283 0.880 0.1809 -0.709 4.8e-01
> RACE -0.3226 0.724 0.0817 -3.950 7.8e-05
> MHHT:MHDM 0.0524 1.054 0.0574 0.913 3.6e-01
> FHHT:FHCAD:FHDM 0.0361 1.037 0.0412 0.877 3.8e-01
>
> Likelihood ratio test=43.8 on 5 df, p=2.55e-08 n= 646
>
> but
> > coxph(jtree9$cph.call,z,rep(1,dim(z)[1]))$weights
> NULL
>
> while the documentation says that a coxph object should have a weights
> element if supplied. Further (what I am really trying to do) if I put this
> in a function
>
> optimize.W<-function(W,G,Groups,cph.call,z){
> n<-length(Groups)
> grp.wt<-rep(0,n)
> for(i in 1:length(G)){
> ind<-Groups == G[i]
> grp.wt[ind]<-W[i]
> }
> mod<-coxph(cph.call,z,grp.wt,na.action=na.omit,singular.ok=T)
> sum(mod$residuals^2)
> }
>
> >
optimize.W(rep(1,13),unique(jtree9$members[,jtree9$depth]),jtree9$members[,jtree9$depth],
> jtree9$cph.call,z)
>
> Error in eval(expr, envir, enclos) : Object "grp.wt" not found
>
> I'm puzzled, is this some sort of environment problem?
> Here's my R particulars
Nope, and if you search the R-help mail archive you will find a
number of instances of this (and different solutions).
The problem is that at some point in time it was decided that
modelling functions would evaluate some of their arguments in
special ways. Usually this is a good thing but sometimes it isn't.
When you make a formula it captures the local environment and hence
bindings for any variables specified (you make yours outside of the
call to your optimize.W function so it captures something out
there).
When evaluating the arguments (the weights one as well as some
others) the first place that is searched is the data argument (if
it's a dataframe then its column names are used), so you can get
what you want by adding a column to z with names grp.wt (I think).
You could also get what you want by doing something like
assign("grp.wt", grp.wt, env=environment(cph.call))
to explicitly put this variable in a place that is searched for a
binding.
At this point there is no clean solution. We probably have to decide
at some point that this rather peculiar evaluation model needs to
disappear but that is hard because it is not fairly ingrained.
There are alternate views and solutions...
A simple example:
data(aml)
cph1 <- Surv(time,status)~x
foo <- function(cphf, z) {
wts <- rexp(nrow(z))
mod <-coxph(cphf,z, wts)
}
foo(cph1,aml)
#Error in eval(expr, envir, enclos) : Object "wts" not found
foo2 <- function(cphf, z) {
wts <- rexp(nrow(z))
assign("wts",wts,env=environment(cphf))
mod<-coxph(cphf,z,wts)
}
foo2(cph1, aml)
>
> > R.Version()
> $platform
> [1] "i686-pc-linux-gnu"
> $arch
> [1] "i686"
> $os
> [1] "linux-gnu"
> $system
> [1] "i686, linux-gnu"
> $status
> [1] ""
> $major
> [1] "1"
> $minor
> [1] "3.1"
> $year
> [1] "2001"
> $month
> [1] "08"
> $day
> [1] "31"
> $language
> [1] "R"
>
> I am also using the most recent version of survival from cran.
>
> Thanks
> Nicholas
>
>
>
> CH3
> |
> N Nicholas Lewin-Koh
> / \ Dept of Statistics
> N----C C==O Program in Ecology and Evolutionary Biology
> || || | Iowa State University
> || || | Ames, IA 50011
> CH C N--CH3 http://www.public.iastate.edu/~nlewin
> \ / \ / nlewin at iastate.edu
> N C
> | || Currently
> CH3 O Graphics Lab
> School of Computing
> National University of Singapore
> The Real Part of Coffee kohnicho at comp.nus.edu.sg
>
>
-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-
> r-help mailing list -- Read http://www.ci.tuwien.ac.at/~hornik/R/R-FAQ.html
> Send "info", "help", or "[un]subscribe"
> (in the "body", not the subject !) To: r-help-request at
stat.math.ethz.ch
>
_._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._
--
+---------------------------------------------------------------------------+
| Robert Gentleman phone : (617) 632-5250 |
| Associate Professor fax: (617) 632-2444 |
| Department of Biostatistics office: M1B28
| Harvard School of Public Health email: rgentlem at jimmy.dfci.harvard.edu |
+---------------------------------------------------------------------------+
-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-
r-help mailing list -- Read http://www.ci.tuwien.ac.at/~hornik/R/R-FAQ.html
Send "info", "help", or "[un]subscribe"
(in the "body", not the subject !) To: r-help-request at
stat.math.ethz.ch
_._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._