Custom distributions using ActuDistns

[RossHamil]

Hi Chris

I am trying to perform analysis using the custom distributions functionality and the package ActuDistns. In this instance I am using the Perks survival model. Are you able to show me how to use these two packages together / where I am going wrong with the code below?

Also separately is there an easy way to generate the p-values with significance level stars as per the standard glm output?

#Using eha example for ActuDistns if (require("ActuDistns")) { custom.perks <- list(name="perks", pars=c("alpha","beta"), location="alpha", transforms=c(log, log), inv.transforms=c(exp, exp), inits=function(t){ c(0.001, 1/median(t)) } ) } fitperks <- flexsurvreg(formula = Surv(modelentryage, modelexitage,d) ~ modelentryage + Sex, data = rdf_data, dist=custom.perks) fitperks

And I get the following error: Forming integrated hazard function... Forming integrated rmst function... Forming integrated mean function... Error in (function (x, alpha = 1, beta = 1) : unused argument (log = TRUE) Some(<code style = 'font-size:10pt'> Error in (function (x, alpha = 1, beta = 1) : unused argument (log = TRUE) </code>) Error in (function (x, alpha = 1, beta = 1) : unused argument (log = TRUE)

[chjackson]

This is hard for me to check, because the ActuDistns package doesn't seem to be available for current R. It was removed from CRAN in 2017 . I don't know anything about the Perks distribution, but there appears to be one called that in the VGAM package which is available.

I'm not a fan of p values, but you can derive them from the confidence intervals.

[chjackson]

Though I wonder if this is the problem - see this bit in help(flexsurvreg), "Density functions should also have an argument log, after the parameters, which when TRUE, computes the log density, using a numerically stable additive formula if possible.". Do the dperks function not have a log argument? If so, just write your own wrapper around it that does. See the flexsurv vignette Section 4 for an example of a similar wrapper.

[RossHamil]

Thanks Chris, I'll take a look at the VGAM package. I managed to get hold of ActuDistns through an old snapshot:

install.packages("ActuDistns", repos = "https://cran.microsoft.com/snapshot/2017-02-04/")
library(ActuDistns)

Alternatively could you show me an example of how to define the distribution without the package, eg. looking at the code for the perks model:

#Perks probability density function

dperks=function(x,alpha=1,beta=1)
{
    ret = ifelse(x<=0 | alpha<=0 | beta<=0,NaN,
           alpha*exp(beta*x)*(1+alpha)/
           (1+alpha*exp(beta*x))**2)
    return(ret)
}

#Perks hazard rate function

hperks=function(x,alpha=1,beta=1)
{
    ret = ifelse(x<=0 | alpha<=0 | beta<=0,NaN,
           alpha*exp(beta*x) /
           (1+alpha*exp(beta*x)))
    return(ret)
}

#Perks integrated hazard rate function

iperks=function(x,t=1,alpha=1,beta=1)
{
    ret = ifelse(x<=0 | t<=0 | length(x)!=length(t) | alpha<=0 | beta<=0,
           NaN, (1/beta)*log(
           (1+alpha*exp(beta*(x+t))) /
           (1+alpha*exp(beta*x))))
    return(ret)
}

I had an attempt at doing this using the examples given in the instructions, but it runs very slowly (is still going)

## Custom distribution: supply the hazard function only
hperks2 <- function (x, alpha = 1, beta = 1) 
{
    ret = ifelse(x <= 0 | alpha <= 0 | beta <= 0, NaN, alpha * 
        exp(beta * x)/(1 + alpha * exp(beta * x)))
    return(ret)
}
Hperks2 <- function (x, t = 1, alpha = 1, beta = 1) 
{
    ret = ifelse(x <= 0 | t <= 0 | length(x) != length(t) | alpha <= 
        0 | beta <= 0, NaN, (1/beta) * log((1 + alpha * exp(beta * 
        (x + t)))/(1 + alpha * exp(beta * x))))
    return(ret)
}
hperks2 <- Vectorize(hperks2)
Hperks2 <- Vectorize(Hperks2)
custom.perks2 <- list(name="perks2", pars=c("alpha","beta"), location="alpha",
transforms=c(log,log), inv.transforms=c(exp,exp),
inits= function (t, mf, mml, aux) 
{
    c(0.001, 1/mean(t))
})
fitperks2 <- flexsurvreg(Surv(modelentryage, modelexitage, d) ~ modelentryage + Sex, data = rdf_data, dist=custom.perks2)
fitperks2