[R] Error in step() (or stepAIC) for Cox model

[Prof Brian Ripley]

On Mon, 24 Oct 2005, Li, Jia wrote:

> Hello all,
>
> I am trying to use stepwise procedure to select covariates in Cox model
> and use bootstrap to repeat stepwise selection,  then record how many
> times  variables  are chosen by step() in bootstrap replications. When I
> use step() (or stepAIC) to do model selection, I got errors. Here is the
> part of my code
>
> for (j in 1:mm){    #<--mm=10
>
> for (b in 1:nrow(reg.bs)){ #<--bootstrap 10 times
>
> mi<-data.frame(tlfup,cen,complete(imp,j)) #<--completed data sets after
> MI
>
> in.star<-sample(1:n,n,replace=T) #<--to sample id number 1-1851.
>
> data.star<-mi[in.star,]
>
> M<-coxph(Surv(tlfup,cen)~mi$trt+mi$nodes+mi$htypeed1+mi$htypeed2+mi$ngra
> ded2+mi$agem40
>
> +mi$agem40sq+mi$er+mi$pr,data=data.star)
>
> reg.model<-step(M) #<--do stepwise selection
>
> reg.step[[b]]<-c((reg.model$coef)) #<-store selected variables
>
> chosen.vb[[b]]<-names(reg.step[[b]]) #<--store names of selected
> variables
>
> }
>
> tot.vb[[j]]<-c(unlist(chosen.vb))
>
> }
>
> Error in reg.step[[b]] : subscript out of bounds
>
>
> varibles<-unlist(tot.vb) #<--change to a vector
>
> table(varibles)           #<--how many times the names have been
> selected
>
> Error in sort(unique.default(x), na.last = TRUE) :
>        'x' must be atomic
>
>
> I figure the reason may be that when stepwise procedure selects the
> chosen model with no varibles, that is
>
> Start:  AIC= 12436.85
> Surv(tlfup, cen) ~ mi$trt + mi$nodes + mi$htypeed1 + mi$htypeed2 +
>    mi$ngraded2 + mi$agem40 + mi$agem40sq + mi$er + mi$pr
>
>              Df   AIC
> - mi$pr        1 12435
> - mi$trt       1 12435
> - mi$agem40sq  1 12435
> - mi$agem40    1 12435
> - mi$htypeed2  1 12435
> - mi$nodes     1 12435
> - mi$er        1 12436
> - mi$ngraded2  1 12436
> - mi$htypeed1  1 12436
> <none>           12437
>
> Step:  AIC= 12425.91
> Surv(tlfup, cen) ~ mi$htypeed1 + mi$ngraded2 + mi$er
>
>              Df   AIC
> - mi$er        1 12425
> - mi$ngraded2  1 12425
> - mi$htypeed1  1 12426
> <none>           12426
>
> Step:  AIC= 12424.77
> Surv(tlfup, cen) ~ 1
>
> then reg.model$coef dosen't exist anymore (am I right?), but I need to
> count the variables that are selected by step() in bootstrap
> replications, what should I do?

Specify the lower scope suitably.

-- 
Brian D. Ripley,                  ripley at stats.ox.ac.uk
Professor of Applied Statistics,  http://www.stats.ox.ac.uk/~ripley/
University of Oxford,             Tel:  +44 1865 272861 (self)
1 South Parks Road,                     +44 1865 272866 (PA)
Oxford OX1 3TG, UK                Fax:  +44 1865 272595