# Creating a custom distribution with flexsurvreg

I'm interested in fitting a parametric survival model but would like to explore the use of a Beta distribution for this purpose rather than a Weibull, Exponential, etc. model.

The Beta distribution is useful in that it can be parameterized with a min and a max parameter which means that, for my purposes, I can constrain my resulting survival curve to only take values within those limits. My hope is that the Beta is flexible enough to fit a survival curve that is comparable to other, more common distributions.

According to the flexsurvreg reference manual on the use of custom distributions, I need only create an input list that has a few features: a name, names parameters, etc. My understanding is that I also need to define a probability density function and a cumulative density function for my distribution.

A Beta density parameterized with $$\theta_{1}$$/$$\theta_{2}$$/min/max parameters can be written as: $$Pr(X=x)=\frac{1}{Beta(\theta_{1},\theta_{2})}\cdot \frac{(x-Min)^{\theta_{1}-1} \cdot(Max-x)^{\theta_{2}-1}}{(Max-Min)^{\theta_{1}+\theta_{2}-1}}$$

The next thing flexsurveg( ) needs is a list with some metadata about my distribution (initializing values, etc). I'm trying to get a simple example working before I start adding predictors, etc., but I can't get a basic model to fit without hearing about "non-finite" function values. I currently suspect that, when $$\theta_1$$ and/or $$\theta_{2}$$ are possibly large values, the exponents are "blowing up" and causing this issue.

Here's a workable example:

#Simulate some right-censored data
set.seed(8765309)
myData <- data.frame(i = 1:100,
time = ceiling(runif(100, min = 1, max = 60)),
cens = sample(0:1, 100, replace = TRUE),
age = ceiling(rnorm(100, mean = 45, sd = 15)))

#Define the beta density
dbetamax <- Vectorize(function(x, theta1, theta2, log = FALSE){
d_y <- ifelse(x <= 0, 0,
ifelse(x >= 120, 0,
(1 / beta(theta1, theta2)) * ( (x - 0)^(theta1- 1) ) * ( (120 - x)^(theta2-1) ) * (1 / ((120 - 0)^(theta1+theta2-1)))))
ifelse(log, log(d_y), d_y)
})

#Define the custom distribution
custom.betamax <- list(name = "betamax",
pars = c("theta1", "theta2"),
location = "theta1",
transforms = c(log, log),
inv.transforms = c(exp, exp),
inits = function(t){
c(1, 1)
})

#Pass to flexsurv
flexsurv::flexsurvreg(survival::Surv(myData$$time, myData$$cens) ~ age, data = myData, dist = custom.betamax)


Which results in the following:

Forming cumulative distribution function...

Forming integrated rmst function...

Forming integrated mean function...

Error in integrate(f = function (u, ...) : non-finite function value

Error: $operator is invalid for atomic vectors Does anybody have much experience with custom distributions in flexsurvreg( ). If so? If so, do you know what the nature of this error is? And, if so, how does one avoid integration issues like these? • Your uses of "shape" and "rate" to describe the Beta parameters are strange, because they correspond to neither property. This leads one to suspect that specifying location = "shape" is going to cause problems. – whuber Commented Feb 17, 2020 at 21:37 • I edited the question to use the same terminology in R that I used in the written density just above it. The flexsurvreg documentation indicates that that particular element of the list, "location", should identify the "main parameter governing the mean of the distribution." This does not have a clear analog for the 4-parameter Beta (whose mean is a function of all 4 parameters). Do you have a sense of what the nature of these "problems" are? Commented Feb 17, 2020 at 21:58 • Knowing nothing about this package, my first guess would be that it's feeding an invalid parameter to the density function. That could easily happen if it thinks either of the Beta parameters is a location. A location parameter has a strict meaning: it must modify the PDF only by appearing in the form$f(x;\ldots,\mu) = f(x-\mu;\ldots).\$
– whuber
Commented Feb 17, 2020 at 22:04

I know this is not flexsurvreg, or even R, but if you are interested in Python, you can do this in Python's lifelines package¹ also using custom model creation.

Instead of defining a PDF, we define a cumulative hazard (often easier to express and reason about vs a PDF). It's all wrapped into a class, subclassed from ParametricUnivariateFitter.

from lifelines.fitters import ParametricUnivariateFitter

from autograd import numpy as np

# create some fake data
n = 100
T_actual = 10 * np.random.random(n)**2
T_censor = 10 * np.random.random(n)**2
E = T_actual < T_censor
T_obs = np.minimum(T_actual, T_censor)

MAX = 120
MIN = 0

class BetaFitter(ParametricUnivariateFitter):
_fitted_parameter_names = ['alpha_', 'beta_']
_bounds = [(0, None), (0, None)]

def _cumulative_density(self, params, times):
alpha_, beta_ = params
return betainc(alpha_, beta_, (times - MIN) / (MAX-MIN))

def _cumulative_hazard(self, params, times):
return -np.log(1-self._cumulative_density(params, times))


Inference is done:

bf = BetaFitter()
bf.fit(T_obs, E)
bf.print_summary()


Plotting:

bf.plot_survival_function()


More documentation on custom models here.

¹I'm the author of lifelines.