makeham.RdMaximum likelihood estimation of the 3-parameter Makeham distribution.
Logical. Suppress a warning? Ignored for VGAM 0.9-7 and higher.
Parameter link functions applied to the
shape parameter shape,
scale parameter scale, and
other parameter epsilon.
All parameters are treated as positive here
(cf. dmakeham allows epsilon = 0, etc.).
See Links for more choices.
Optional initial values.
A NULL means a value is computed internally.
A value must be given for iepsilon currently, and this
is a sensitive parameter!
See CommonVGAMffArguments.
Argument probs.y is used only when imethod = 2.
To be currently ignored.
The Makeham distribution, which adds another parameter
to the Gompertz distribution,
has cumulative distribution function
$$F(y; \alpha, \beta, \varepsilon) =
1 - \exp
\left\{
-y \varepsilon + \frac {\alpha}{\beta}
\left[ 1 - e^{\beta y} \right]
\right\}
$$
which leads to a probability density function
$$f(y; \alpha, \beta, \varepsilon) =
\left[
\varepsilon + \alpha e^{\beta y} \right]
\;
\exp
\left\{
-y \varepsilon + \frac {\alpha}{\beta}
\left[ 1 - e^{\beta y} \right]
\right\},
$$
for \(\alpha > 0\),
\(\beta > 0\),
\(\varepsilon \geq 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.
if (FALSE) 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)
#> Iteration 1: loglikelihood = -247.82362
#> Iteration 2: loglikelihood = -246.16387
#> Iteration 3: loglikelihood = -246.15013
#> Iteration 4: loglikelihood = -246.15009
#> Iteration 5: loglikelihood = -246.15009
fit1 <- vglm(y1 ~ 1, makeham, data = mdata,
etastart = cbind(predict(fit0), log(0.1)), trace = TRUE)
#> Iteration 1: loglikelihood = -16380.129
#> Iteration 2: loglikelihood = NaN
#> Taking a modified step....................
coef(fit1, matrix = TRUE)
#> loglink(scale) loglink(shape) loglink(epsilon)
#> (Intercept) -142.7856 177.0155 2.27214
summary(fit1)
#>
#> Call:
#> vglm(formula = y1 ~ 1, family = makeham, data = mdata, etastart = cbind(predict(fit0),
#> log(0.1)), trace = TRUE)
#>
#> Coefficients:
#> Estimate Std. Error z value Pr(>|z|)
#> (Intercept):1 -142.786 292.008 -0.489 0.625
#> (Intercept):2 177.016 366.996 0.482 0.630
#> (Intercept):3 2.272 2.430 0.935 0.350
#>
#> Names of linear predictors: loglink(scale), loglink(shape), loglink(epsilon)
#>
#> Log-likelihood: -4.56584e+79 on 2997 degrees of freedom
#>
#> Number of Fisher scoring iterations: 2
#>
#> No Hauck-Donner effect found in any of the estimates
#>
# \dontrun{}