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