# Example 1: simulated ZIP data
zdata <- data.frame(x2 = runif(nn <- 2000))
zdata <- transform(zdata, phi1 = logit(-0.5 + 1*x2, inverse = TRUE),
phi2 = logit( 0.5 - 1*x2, inverse = TRUE),
Phi1 = logit(-0.5 , inverse = TRUE),
Phi2 = logit( 0.5 , inverse = TRUE),
lambda1 = loge( 0.5 + 2*x2, inverse = TRUE),
lambda2 = loge( 0.5 + 2*x2, inverse = TRUE))
zdata <- transform(zdata, y1 = rzipois(nn, lambda1, Phi1),
y2 = rzipois(nn, lambda2, Phi2))
with(zdata, table(y1)) # Eyeball the data
with(zdata, table(y2))
fit1 <- vglm(y1 ~ x2, zipoisson(zero = 1), zdata, crit = "coef")
fit2 <- vglm(y2 ~ x2, zipoisson(zero = 1), zdata, crit = "coef")
coef(fit1, matrix = TRUE) # These should agree with the above values
coef(fit2, matrix = TRUE) # These should agree with the above values
# Fit all two simultaneously, using a different parameterization:
fit12 <- vglm(cbind(y1, y2) ~ x2, zipoissonff, zdata, crit = "coef")
coef(fit12, matrix = TRUE) # These should agree with the above values
# Example 2: McKendrick (1926). Data from 223 Indian village households
cholera <- data.frame(ncases = 0:4, # Number of cholera cases,
wfreq = c(168, 32, 16, 6, 1)) # Frequencies
fit <- vglm(ncases ~ 1, zipoisson, wei = wfreq, cholera, trace = TRUE)
coef(fit, matrix = TRUE)
with(cholera, cbind(actual = wfreq,
fitted = round(dzipois(ncases, lambda = Coef(fit)[2],
phi = Coef(fit)[1]) *
sum(wfreq), dig = 2)))
# Example 3: data from Angers and Biswas (2003)
abdata <- data.frame(y = 0:7, w = c(182, 41, 12, 2, 2, 0, 0, 1))
abdata <- subset(abdata, w > 0)
fit <- vglm(y ~ 1, zipoisson(lphi = probit, iphi = 0.3),
abdata, weight = w, trace = TRUE)
fit@misc$prob0 # Estimate of P(Y = 0)
coef(fit, matrix = TRUE)
Coef(fit) # Estimate of phi and lambda
fitted(fit)
with(abdata, weighted.mean(y, w)) # Compare this with fitted(fit)
summary(fit)
# Example 4: This RR-ZIP is known as a COZIGAM or COZIVGLM-ZIP
rrzip <- rrvglm(Alopacce ~ bs(WaterCon), zipoissonff(zero = NULL),
hspider, trace = TRUE)
coef(rrzip, matrix = TRUE)
Coef(rrzip)
summary(rrzip)
plotvgam(rrzip, lcol = "blue")
Run the code above in your browser using DataLab