[R] survreg.distributions() error

Nick Reich nick.reich at gmail.com
Mon Jun 1 05:47:37 CEST 2009


Hi there.

I am receiving an unexpected error message when creating a new
distribution for the survreg() function in the survival package.  I
understand the survival.distributions() function and have been
following the Cauchy example provided in the help file.

My goal is to use survreg to fit a gamma distribution to interval
censored data.

Here is a simple example of what I'm trying to do.  This code
recreates the error on my machine.

require(survival)

## create the gamma distribution for survreg()
## paramter vector is parms = c(shape, scale)
mygamma <- list(name='gamma',
		init = function(x, weights, ...) {
			c(median(x), mad(x)) },
		density = function(x, parms) {
			shape <- parms[1]  ## "k" in wikipedia definition
			scale <- parms[2]  ## "theta" in wikipedia
			cbind(pgamma(x, shape=shape, scale=scale),
			      1-pgamma(x, shape=shape, scale=scale),
			      dgamma(x, shape=shape, scale=scale),
			      (shape-1)/x - 1/scale,
			      (shape-1)*(shape-2)/x^2
			         - 2*(shape-1)/(x*scale)
			         + 1/scale^2)},
		quantile = function(p, parms) {
			qgamma(p, shape=parms[1], scale=parms2)},
		deviance= function(...) stop('deviance residuals not defined'),
		parms=c(2,2)
		)

## create some interval censored data
xx <- Surv(time=rep(c(.5, 1, 1.5), each=5),
	   time2=rep(2, 15),
	   event=rep(3, 15),
	   type="interval")

## lognormal model runs fine
survreg(xx~1, dist="lognormal")

## mygamma model gives an error
survreg(xx~1, dist=mygamma)

This is the error message that I get:
Error in density(z2) : argument "parms" is missing, with no default

I've tracked this down to a single line in the survreg.fit() function.
 This function defines another function derfun() whose first few lines
look like this:

derfun <- function(y, eta, sigma, density, parms) {
        ny <- ncol(y)
        status <- y[, ny]
        z <- (y[, 1] - eta)/sigma
        dmat <- density(z, parms)  ## line A: THIS COMMAND DOES NOT
GIVE AN ERROR
        dtemp <- dmat[, 3] * dmat[, 4]
        if (any(status == 3)) {
            z2 <- (y[, 2] - eta)/sigma
            dmat2 <- density(z2)    ## line B: THIS IS THE COMMAND
THAT GIVES THE ERROR
        }
        else {
            dmat2 <- matrix(0, 1, 5)
            z2 <- 0
        }

>From playing with debug() a little bit, it looks to me like if line B
(see labels above) called "density(z2, parms)", similarly to line A,
then the error message might be avoided.  Instead, it calls
"density(z2)".  Maybe there's a fundamental error in my setting up of
the gamma distribution in the first place, but this seems like
possibly just as simple as a typo mistake.

So...  two questions:
1.  Can anyone verify that this is indeed a bug that's as simple as a
typo in the package code?
2.  I tried cutting and pasting the code for survreg.fit() into a new
file, making the edit myself and then re-sourcing the function code to
overwrite the "bad" function.  But this isn't working.  Is there a
simple way to test out a new version of this function?

Thanks,
Nick




More information about the R-help mailing list