Dear R-help folks:
Here is my platform:
> version
platform sparc-sun-solaris2.7
arch sparc
os solaris2.7
system sparc, solaris2.7
status
major 1
minor 5.0
year 2002
month 04
day 29
language R
I have a function called FitModels(), which simply takes in the
names of a data.frame and two variable names within that data.frame, and
fits and returns a list of objects from 2 coxph() fits, one main effects and
one interaction model. Sometimes the two variables are such that there is a
warning message: convergence has not been reached, or the X matrix is
singular, or some such. When a warning occurs, then, I would like to play
it safe and have FitModels() return a list at least one of whose elements
L[[i]] inherits(L[[i]],"try-error"); so, I guess, I want to set
options(warn=2), so that, "all warnings are turned into errors," as
the help
says. Then, I hope, try() will return "try-error"'s. So I start
an R
session, which has warn=0 and show.error.messages=TRUE by default, and here
is what happens (question after R session):
################################# Start R session:
> library(survival)
> ## Here are 2 calls with different pairs of variables, both of which
yield warnings:> qq <-
FitModels(vnames=c("acd11apcd3p","drpcd8p"),dframe="ALLmth12",survname="fail
tm12", censname="failcens")
Warning message:
Loglik converged before variable 2 ; beta may be infinite. in: fitter(X,
Y, strats, offset, init, control, weights = weights, > ## warnings are just warnings when options()$warn=0:
> lapply(qq,class)
$cph.x1x2
[1] "coxph"
$cph.x1x2.int
[1] "coxph"> qq <-
FitModels(vnames=c("cd3ncd56p","adrpcd8p"),dframe="ALLmth12",survname="failt
m12", censname="failcens")
Warning message:
Loglik converged before variable 1 ; beta may be infinite. in: fitter(X,
Y, strats, offset, init, control, weights = weights, > lapply(qq,class)
$cph.x1x2
[1] "coxph"
$cph.x1x2.int
[1] "coxph"> ## Now attempt to turn all warnings into errors:
> options(warn=2)
> ## same model fitting done:
> qq <-
FitModels(vnames=c("acd11apcd3p","drpcd8p"),dframe="ALLmth12",survname="fail
tm12", censname="failcens")
Error in fitter(X, Y, strats, offset, init, control, weights = weights, :
(converted from warning) Loglik converged before variable 2 ; beta
may be infinite. > ## The interaction model had a warning, which now yields a try-error
(what I wanted!!!):> lapply(qq,class)
$cph.x1x2
[1] "coxph"
$cph.x1x2.int
[1] "try-error"> ## Here is the same second model:
> qq <-
FitModels(vnames=c("cd3ncd56p","adrpcd8p"),dframe="ALLmth12",survname="failt
m12", censname="failcens")> ## But now the warning has vanished and no try-error is produced. Why
not???:> lapply(qq,class)
$cph.x1x2
[1] "coxph"
$cph.x1x2.int
[1] "coxph"> ## Check the options:
> options()$warn
[1] 2> options()$show.err
[1] TRUE
################################# End R session:
Shouldn't the second call to FitModels() also yield a try-error for one of
the models? If I do the calls in the reverse order, the same thing happens:
the first call always yields a try-error, the second does not. What am I
missing about the behavior of the warn option? Is there a way to do what I
want, given that I really want to do this for lots of pairs of variables and
easily flag the ones with such fitting problems? Thanks in advance.
Tom Richards
P.S., And here is the simple function FitModels(), if you care to see it.
It was extracted from a much larger function that computes several other
covariates from x1 and x2:
FitModels <-
function(vnames,dframe,survname,censname,censval="Progressed"){
## Get covariate x1 from data frame:
x1x1 <- eval(parse(text=paste(dframe,vnames[1],sep="$")))
## Get covariate x2 from data frame:
x2x2 <- eval(parse(text=paste(dframe,vnames[2],sep="$")))
## Get survival time:
response <- eval(parse(text=paste(dframe,survname,sep="$")))
## Get censoring variable:
censred <- eval(parse(text=paste(dframe,censname,sep="$")))
censred <- as.numeric(censred == censval)
## Take out observations with missing x1 or x2.
SUBSET.to.use <- (!is.na(x1x1))&(!is.na(x2x2))
## Take only x1/x2 values specified in SUBSET.to.use:
x1 <- x1x1[(SUBSET.to.use)]
x2 <- x2x2[(SUBSET.to.use)]
## Properly subset the binary response variable:
response <- response[(SUBSET.to.use)]
censred <- censred[(SUBSET.to.use)]
## Make new data frame to use for modelling:
DD <- data.frame(srvtime=response,censored=censred,x1=x1,x2=x2)
#######################################################################
## Fit MAIN EFFECTS model with both covariates x1 and x2:
cph.x1x2 <- try(coxph(Surv(srvtime,censored) ~ x1+x2,
data=DD,na.action=na.exclude,x=TRUE))
## Fit MAIN EFFECTS + INTERACTION LR model with covariates
x1 and x2:
cph.x1x2.int <- try(coxph(Surv(srvtime,censored) ~ x1+x2+x1:x2,
data=DD,na.action=na.exclude,x=TRUE))
## Return models fitted to these data:
retval <- list(
cph.x1x2 = cph.x1x2,
cph.x1x2.int = cph.x1x2.int)
invisible(retval)
}
-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-
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
_._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._
