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)Run the code above in your browser using DataLab