# Example 1: simulated ZIP data
zdata <- data.frame(x2 = runif(nn <- 2000))
zdata <- transform(zdata, pstr01 = logit(-0.5 + 1*x2, inverse = TRUE),
pstr02 = logit( 0.5 - 1*x2, inverse = TRUE),
Ps01 = logit(-0.5 , inverse = TRUE),
Ps02 = 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, lambda = lambda1, pstr0 = Ps01),
y2 = rzipois(nn, lambda = lambda2, pstr0 = Ps02))
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],
pstr0 = 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(lpstr0 = probit, ipstr0 = 0.8),
abdata, weight = w, trace = TRUE)
fit@misc$pobs0 # Estimate of P(Y = 0)
coef(fit, matrix = TRUE)
Coef(fit) # Estimate of pstr0 and lambda
fitted(fit)
with(abdata, weighted.mean(y, w)) # Compare this with fitted(fit)
summary(fit)
# Example 4: zero-deflated model for an intercept-only data
zdata <- transform(zdata, lambda3 = loge( 0.0 , inverse = TRUE))
zdata <- transform(zdata, deflat_limit = -1 / expm1(lambda3)) # Boundary
# The 'pstr0' parameter is negative and in parameter space:
zdata <- transform(zdata, usepstr0 = deflat_limit / 1.5)
zdata <- transform(zdata, y3 = rzipois(nn, lambda3, pstr0 = usepstr0))
head(zdata)
with(zdata, table(y3)) # A lot of deflation
fit3 <- vglm(y3 ~ 1, zipoisson(zero = -1, lpstr0 = identity),
zdata, trace = TRUE, crit = "coef")
coef(fit3, matrix = TRUE)
# Check how accurate it was:
zdata[1, 'usepstr0'] # Answer
coef(fit3)[1] # Estimate
Coef(fit3)
# Example 5: This RR-ZIP is known as a COZIGAM or COZIVGLM-ZIP
set.seed(123)
rrzip <- rrvglm(Alopacce ~ bs(WaterCon, df = 3), zipoisson(zero = NULL),
hspider, trace = TRUE, Index.corner = 2)
coef(rrzip, matrix = TRUE)
Coef(rrzip)
summary(rrzip)
plotvgam(rrzip, lcol = "blue")
Run the code above in your browser using DataLab