Makeham Regression Family Function
Maximum likelihood estimation of the 3-parameter Makeham distribution.
makeham(lscale = "loglink", lshape = "loglink", lepsilon = "loglink",
iscale = NULL, ishape = NULL, iepsilon = NULL,
gscale = exp(-5:5),gshape = exp(-5:5), gepsilon = exp(-4:1),
nsimEIM = 500, oim.mean = TRUE, zero = NULL, nowarning = FALSE)nowarning |
Logical. Suppress a warning? Ignored for VGAM 0.9-7 and higher. |
lshape, lscale, lepsilon |
Parameter link functions applied to the
shape parameter |
ishape, iscale, iepsilon |
Optional initial values.
A |
gshape, gscale, gepsilon |
|
nsimEIM, zero |
See |
oim.mean |
To be currently ignored. |
The Makeham distribution, which adds another parameter to the Gompertz distribution, has cumulative distribution function
F(y; alpha, beta, epsilon) = 1 - exp(-y * epsilon + (alpha / beta) * [1 - e^(beta * y)])
which leads to a probability density function
f(y; alpha, beta, epsilon) = (epsilon + alpha * e^(beta y) ) * exp(-y * epsilon + (alpha / beta) * [1 - e^(beta * y)])
for alpha > 0,
beta > 0,
epsilon >= 0,
y > 0.
Here, beta is called the scale parameter scale,
and alpha is called a shape parameter.
The moments for this distribution do
not appear to be available in closed form.
Simulated Fisher scoring is used and multiple responses are handled.
An object of class "vglmff" (see vglmff-class).
The object is used by modelling functions such as vglm,
and vgam.
A lot of care is needed because
this is a rather difficult distribution for parameter estimation,
especially when the shape parameter is large relative to the
scale parameter.
If the self-starting initial values fail then try experimenting
with the initial value arguments, especially iepsilon.
Successful convergence depends on having very good initial values.
More improvements could be made here.
Also, monitor convergence by setting trace = TRUE.
A trick is to fit a gompertz distribution and use
it for initial values; see below.
However, this family function is currently numerically fraught.
T. W. Yee
## Not run: set.seed(123)
mdata <- data.frame(x2 = runif(nn <- 1000))
mdata <- transform(mdata, eta1 = -1,
ceta1 = 1,
eeta1 = -2)
mdata <- transform(mdata, shape1 = exp(eta1),
scale1 = exp(ceta1),
epsil1 = exp(eeta1))
mdata <- transform(mdata,
y1 = rmakeham(nn, shape = shape1, scale = scale1, eps = epsil1))
# A trick is to fit a Gompertz distribution first
fit0 <- vglm(y1 ~ 1, gompertz, data = mdata, trace = TRUE)
fit1 <- vglm(y1 ~ 1, makeham, data = mdata,
etastart = cbind(predict(fit0), log(0.1)), trace = TRUE)
coef(fit1, matrix = TRUE)
summary(fit1)
## End(Not run)Please choose more modern alternatives, such as Google Chrome or Mozilla Firefox.