[R] case weights in coxph (survival)

[Robert Gentleman]

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
_._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._