Makeham Regression Family Function

Maximum likelihood estimation of the 3-parameter Makeham distribution.

Usage

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)

Arguments

nowarning

Logical. Suppress a warning? Ignored for VGAM 0.9-7 and higher.

lshape, lscale, lepsilon

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.

ishape, iscale, iepsilon

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!

gshape, gscale, gepsilon

See CommonVGAMffArguments.

nsimEIM, zero

See CommonVGAMffArguments. Argument probs.y is used only when imethod = 2.

oim.mean

To be currently ignored.

Details

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.

Value

An object of class "vglmff" (see vglmff-class). The object is used by modelling functions such as vglm, and vgam.

Author

T. W. Yee

Warning

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.

Examples

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 
#> 
 # \dontrun{}