# NOT RUN {
avec <- c(10, 20, 30) # Alter these values
tvec <- 0 # Truncate this value
pobs.a <- logitlink(-(2:4), inverse = TRUE) # Between 0.02 and 0.12
size1 <- exp(1)
gdata <- data.frame(x2 = runif(nn <- 1000))
gdata <- transform(gdata, lambda1 = exp(2 + 0.5 * x2))
gdata <- transform(gdata,
y1 = rgaitnbinom.mlm(nn, size1, mu = lambda1, pobs.a = pobs.a,
truncate = tvec, byrow = TRUE, alter = avec))
gatnbinomial.mlm(alter = avec)
(ty1 <- with(gdata, table(y1)))
propn1 <- c(ty1) / sum(ty1)
plot(as.numeric(names(ty1)), propn1, las = 1, xlab = "y",
yaxs = "i", ylim = c(0, max(propn1) * 1.1), main = "Heaped data",
ylab = "Proportion", lwd = 3, type = "h", col = "blue")
fit1 <- vglm(y1 ~ x2, trace = TRUE, data = gdata, crit = "coef",
gatnbinomial.mlm(alter = avec, truncate = tvec,
zero = c("size", "pobs")))
head(fitted(fit1))
head(predict(fit1))
coef(fit1, matrix = TRUE)
summary(fit1)
# }
Run the code above in your browser using DataLab