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