zipoisson: Zero-Inflated Poisson Distribution Family Function

Description Usage Arguments Details Value Warning Note Author(s) References See Also Examples

View source: R/family.zeroinf.R

Description

Fits a zero-inflated or zero-deflated Poisson distribution by full maximum likelihood estimation.

Usage

1
2
3
4
5
6
7
8
9
zipoisson(lpstr0 = "logitlink", llambda = "loglink", type.fitted =
          c("mean", "lambda", "pobs0", "pstr0", "onempstr0"), ipstr0 =
          NULL, ilambda = NULL, gpstr0 = NULL, imethod = 1,
          ishrinkage = 0.95, probs.y = 0.35, parallel = FALSE, zero = NULL)
zipoissonff(llambda = "loglink", lonempstr0 = "logitlink", type.fitted =
            c("mean", "lambda", "pobs0", "pstr0", "onempstr0"),
            ilambda = NULL, ionempstr0 = NULL, gonempstr0 = NULL,
            imethod = 1, ishrinkage = 0.95, probs.y = 0.35, zero =
            "onempstr0")

Arguments

lpstr0, llambda

Link function for the parameter phi and the usual lambda parameter. See Links for more choices; see CommonVGAMffArguments for more information. For the zero-deflated model see below.

ipstr0, ilambda

Optional initial values for phi, whose values must lie between 0 and 1. Optional initial values for lambda, whose values must be positive. The defaults are to compute an initial value internally for each. If a vector then recycling is used.

lonempstr0, ionempstr0

Corresponding arguments for the other parameterization. See details below.

type.fitted

Character. The type of fitted value to be returned. The first choice (the expected value) is the default. The estimated probability of an observed 0 is an alternative, else the estimated probability of a structural 0, or one minus the estimated probability of a structural 0. See CommonVGAMffArguments and fittedvlm for more information.

imethod

An integer with value 1 or 2 which specifies the initialization method for lambda. If failure to converge occurs try another value and/or else specify a value for ishrinkage and/or else specify a value for ipstr0. See CommonVGAMffArguments for more information.

ishrinkage

How much shrinkage is used when initializing lambda. The value must be between 0 and 1 inclusive, and a value of 0 means the individual response values are used, and a value of 1 means the median or mean is used. This argument is used in conjunction with imethod. See CommonVGAMffArguments for more information.

zero

Specifies which linear/additive predictors are to be modelled as intercept-only. If given, the value can be either 1 or 2, and the default is none of them. Setting zero = 1 makes phi a single parameter. See CommonVGAMffArguments for more information.

gpstr0, gonempstr0, probs.y

Details at CommonVGAMffArguments.

parallel

Details at CommonVGAMffArguments, but unlikely to be practically used actually.

Details

These models are a mixture of a Poisson distribution and the value 0; it has value 0 with probability phi else is Poisson(lambda) distributed. Thus there are two sources for zero values, and phi is the probability of a structural zero. The model for zipoisson() can be written

P(Y = 0) = phi + (1-phi) * exp(-lambda),

and for y=1,2,…,

P(Y = y) = (1-phi) * exp(-lambda) * lambda^y / y!.

Here, the parameter phi satisfies 0 < phi < 1. The mean of Y is (1-phi)*lambda and these are returned as the fitted values, by default. The variance of Y is (1-phi)*lambda*(1 + phi lambda). By default, the two linear/additive predictors of zipoisson() are (logit(phi), log(lambda))^T.

The VGAM family function zipoissonff() has a few changes compared to zipoisson(). These are: (i) the order of the linear/additive predictors is switched so the Poisson mean comes first; (ii) onempstr0 is now 1 minus the probability of a structural 0, i.e., the probability of the parent (Poisson) component, i.e., onempstr0 is 1-pstr0; (iii) argument zero has a new default so that the onempstr0 is intercept-only by default. Now zipoissonff() is generally recommended over zipoisson() (and definitely recommended over yip88). Both functions implement Fisher scoring and can handle multiple responses.

Value

An object of class "vglmff" (see vglmff-class). The object is used by modelling functions such as vglm, rrvglm and vgam.

Warning

Numerical problems can occur, e.g., when the probability of zero is actually less than, not more than, the nominal probability of zero. For example, in the Angers and Biswas (2003) data below, replacing 182 by 1 results in nonconvergence. Half-stepping is not uncommon. If failure to converge occurs, try using combinations of imethod, ishrinkage, ipstr0, and/or zipoisson(zero = 1) if there are explanatory variables. The default for zipoissonff() is to model the structural zero probability as an intercept-only.

Note

This family function can be used to estimate the 0-deflated model, hence pstr0 is not to be interpreted as a probability. One should set, e.g., lpstr0 = "identitylink". Likewise, the functions in Zipois can handle the zero-deflated Poisson distribution too. Although the iterations might fall outside the parameter space, the validparams slot should keep them inside. A (somewhat) similar alternative for zero-deflation is to try the zero-altered Poisson model (see zapoisson).

The use of this VGAM family function with rrvglm can result in a so-called COZIGAM or COZIGLM. That is, a reduced-rank zero-inflated Poisson model (RR-ZIP) is a constrained zero-inflated generalized linear model. See COZIGAM. A RR-ZINB model can also be fitted easily; see zinegbinomial. Jargon-wise, a COZIGLM might be better described as a COZIVGLM-ZIP.

Author(s)

T. W. Yee

References

Thas, O. and Rayner, J. C. W. (2005). Smooth tests for the zero-inflated Poisson distribution. Biometrics, 61, 808–815.

Data: Angers, J-F. and Biswas, A. (2003). A Bayesian analysis of zero-inflated generalized Poisson model. Computational Statistics & Data Analysis, 42, 37–46.

Cameron, A. C. and Trivedi, P. K. (1998). Regression Analysis of Count Data. Cambridge University Press: Cambridge.

Yee, T. W. (2014). Reduced-rank vector generalized linear models with two linear predictors. Computational Statistics and Data Analysis, 71, 889–902.

See Also

gaitpoisson, zapoisson, Zipois, yip88, spikeplot, rrvglm, zipebcom, rpois, simulate.vlm, hdeff.vglm.

Examples

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
# Example 1: simulated ZIP data
zdata <- data.frame(x2 = runif(nn <- 1000))
zdata <- transform(zdata, pstr01  = logitlink(-0.5 + 1*x2, inverse = TRUE),
                          pstr02  = logitlink( 0.5 - 1*x2, inverse = TRUE),
                          Ps01    = logitlink(-0.5       , inverse = TRUE),
                          Ps02    = logitlink( 0.5       , inverse = TRUE),
                          lambda1 =  loglink(-0.5 + 2*x2, inverse = TRUE),
                          lambda2 =  loglink( 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), data = zdata, crit = "coef")
fit2 <- vglm(y2 ~ x2, zipoisson(zero = 1), data = 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, data = zdata, crit = "coef")
coef(fit12, matrix = TRUE)  # These should agree with the above values

# For the first observation compute the probability that y1 is
# due to a structural zero.
(fitted(fit1, type = "pstr0") / fitted(fit1, type = "pobs0"))[1]


# 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), digits = 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 = probitlink, ipstr0 = 0.8),
            data = abdata, weight = w, trace = TRUE)
fitted(fit, type = "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 intercept-only data
zdata <- transform(zdata, lambda3 = loglink(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 / 2)  # Not too near the boundary
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 = "identitylink"),
             data = zdata, trace = TRUE, crit = "coef")
coef(fit3, matrix = TRUE)
# Check how accurate it was:
zdata[1, "usepstr0"]  # Answer
coef(fit3)[1]         # Estimate
Coef(fit3)
vcov(fit3)  # Is positive-definite

# Example 5: This RR-ZIP is known as a COZIGAM or COZIVGLM-ZIP
set.seed(123)
rrzip <- rrvglm(Alopacce ~ sm.bs(WaterCon, df = 3), zipoisson(zero = NULL),
                data = hspider, trace = TRUE, Index.corner = 2)
coef(rrzip, matrix = TRUE)
Coef(rrzip)
summary(rrzip)
## Not run: plotvgam(rrzip, lcol = "blue")

Example output

Loading required package: stats4
Loading required package: splines
y1
  0   1   2   3   4   5   6   7   8  10 
542 154 106 102  42  29  14   4   5   2 
y2
  0   1   2   3   4   5   6   7   8   9  10  11  12  13  14  15  16  17 
660  30  49  48  43  42  30  15  19  11  16  11   6   6   5   6   2   1 
            logitlink(pstr0) loglink(lambda)
(Intercept)       -0.2922953      -0.3177837
x2                 0.0000000       1.7820549
            logitlink(pstr0) loglink(lambda)
(Intercept)        0.6058409       0.5763848
x2                 0.0000000       1.8579954
            loglink(lambda1) logitlink(onempstr01) loglink(lambda2)
(Intercept)       -0.3177837             0.2922953        0.5763848
x2                 1.7820549             0.0000000        1.8579954
            logitlink(onempstr02)
(Intercept)            -0.6058409
x2                      0.0000000
[1] 0.7655351
VGLM    linear loop  1 :  loglikelihood = -180.0534
VGLM    linear loop  2 :  loglikelihood = -179.3605
VGLM    linear loop  3 :  loglikelihood = -179.3477
VGLM    linear loop  4 :  loglikelihood = -179.3477
            logitlink(pstr0) loglink(lambda)
(Intercept)         0.419289     -0.02821652
     actual fitted
[1,]    168 168.00
[2,]     32  32.53
[3,]     16  15.81
[4,]      6   5.12
[5,]      1   1.25
VGLM    linear loop  1 :  loglikelihood = -191.2889
VGLM    linear loop  2 :  loglikelihood = -190.4574
VGLM    linear loop  3 :  loglikelihood = -190.437
VGLM    linear loop  4 :  loglikelihood = -190.437
VGLM    linear loop  5 :  loglikelihood = -190.437
       [,1]
1 0.7583333
2 0.7583333
3 0.7583333
4 0.7583333
5 0.7583333
8 0.7583333
            probitlink(pstr0) loglink(lambda)
(Intercept)         0.1944214      -0.1657264
    pstr0    lambda 
0.5770771 0.8472781 
       [,1]
1 0.3583333
2 0.3583333
3 0.3583333
4 0.3583333
5 0.3583333
8 0.3583333
[1] 0.3583333

Call:
vglm(formula = y ~ 1, family = zipoisson(lpstr0 = probitlink, 
    ipstr0 = 0.8), data = abdata, weights = w, trace = TRUE)

Pearson residuals:
  probitlink(pstr0) loglink(lambda)
1            7.0568          -2.864
2          -13.7264          -3.660
3           -3.8216           6.901
4           -0.0887           6.443
5            1.3828          10.069
8            4.0993          14.811

Coefficients: 
              Estimate Std. Error z value Pr(>|z|)
(Intercept):1   0.1944     0.1741   1.117    0.264
(Intercept):2  -0.1657     0.1786  -0.928    0.353

Names of linear predictors: probitlink(pstr0), loglink(lambda)

Log-likelihood: -190.437 on 10 degrees of freedom

Number of Fisher scoring iterations: 5 

No Hauck-Donner effect found in any of the estimates

         x2    pstr01    pstr02      Ps01      Ps02   lambda1   lambda2 y1 y2
1 0.3966344 0.4741816 0.5258184 0.3775407 0.6224593 1.3408031  3.644681  1  5
2 0.3092640 0.4524600 0.5475400 0.3775407 0.6224593 1.1258384  3.060346  0  2
3 0.2475031 0.4372090 0.5627910 0.3775407 0.6224593 0.9950186  2.704741  3  0
4 0.2826161 0.4458670 0.5541330 0.3775407 0.6224593 1.0674068  2.901513  0  0
5 0.9746158 0.6164757 0.3835243 0.3775407 0.6224593 4.2598396 11.579445  0  0
6 0.9034912 0.5995262 0.4004738 0.3775407 0.6224593 3.6950066 10.044069  4  6
  lambda3 deflat.limit   usepstr0 y3
1       1   -0.5819767 -0.2909884  0
2       1   -0.5819767 -0.2909884  1
3       1   -0.5819767 -0.2909884  1
4       1   -0.5819767 -0.2909884  0
5       1   -0.5819767 -0.2909884  1
6       1   -0.5819767 -0.2909884  1
y3
  0   1   2   3   4   5   6 
173 480 244  76  22   4   1 
VGLM    linear loop  1 :  coefficients = -0.197569862,  0.088324347
VGLM    linear loop  2 :  coefficients = -0.2999035477,  0.0042826213
VGLM    linear loop  3 :  coefficients = -0.3059282301,  0.0031083281
VGLM    linear loop  4 :  coefficients = -0.3059294725,  0.0031121105
VGLM    linear loop  5 :  coefficients = -0.3059294725,  0.0031121105
                 pstr0 loglink(lambda)
(Intercept) -0.3059295      0.00311211
[1] -0.2909884
(Intercept):1 
   -0.3059295 
     pstr0     lambda 
-0.3059295  1.0031170 
              (Intercept):1 (Intercept):2
(Intercept):1   0.001405118   0.001381882
(Intercept):2   0.001381882   0.001821518
RR-VGLM    linear loop  1 :  loglikelihood = -102.86636
   Alternating iteration 1 ,   Convergence criterion  =  1 
    ResSS  =  138.31 
   Alternating iteration 2 ,   Convergence criterion  =  0.0001542638 
    ResSS  =  138.2886 
   Alternating iteration 3 ,   Convergence criterion  =  7.621223e-06 
    ResSS  =  138.2876 
   Alternating iteration 4 ,   Convergence criterion  =  3.7527e-07 
    ResSS  =  138.2875 
   Alternating iteration 5 ,   Convergence criterion  =  1.846195e-08 
    ResSS  =  138.2875 
RR-VGLM    linear loop  2 :  loglikelihood = -87.831774
   Alternating iteration 1 ,   Convergence criterion  =  1 
    ResSS  =  181.07 
   Alternating iteration 2 ,   Convergence criterion  =  0.001859694 
    ResSS  =  180.7339 
   Alternating iteration 3 ,   Convergence criterion  =  0.0006638747 
    ResSS  =  180.614 
   Alternating iteration 4 ,   Convergence criterion  =  0.0001545886 
    ResSS  =  180.5861 
   Alternating iteration 5 ,   Convergence criterion  =  2.792251e-05 
    ResSS  =  180.581 
   Alternating iteration 6 ,   Convergence criterion  =  4.466665e-06 
    ResSS  =  180.5802 
   Alternating iteration 7 ,   Convergence criterion  =  6.787567e-07 
    ResSS  =  180.5801 
RR-VGLM    linear loop  3 :  loglikelihood = -87.467439
   Alternating iteration 1 ,   Convergence criterion  =  1 
    ResSS  =  218.2467 
   Alternating iteration 2 ,   Convergence criterion  =  0.0001374247 
    ResSS  =  218.2168 
   Alternating iteration 3 ,   Convergence criterion  =  3.636422e-05 
    ResSS  =  218.2088 
   Alternating iteration 4 ,   Convergence criterion  =  1.028386e-05 
    ResSS  =  218.2066 
   Alternating iteration 5 ,   Convergence criterion  =  3.00539e-06 
    ResSS  =  218.2059 
   Alternating iteration 6 ,   Convergence criterion  =  8.934244e-07 
    ResSS  =  218.2057 
   Alternating iteration 7 ,   Convergence criterion  =  2.680223e-07 
    ResSS  =  218.2057 
RR-VGLM    linear loop  4 :  loglikelihood = -87.287751
   Alternating iteration 1 ,   Convergence criterion  =  1 
    ResSS  =  242.7126 
   Alternating iteration 2 ,   Convergence criterion  =  1.259938e-10 
    ResSS  =  242.7126 
RR-VGLM    linear loop  5 :  loglikelihood = -87.285865
   Alternating iteration 1 ,   Convergence criterion  =  1 
    ResSS  =  233.5506 
   Alternating iteration 2 ,   Convergence criterion  =  9.663082e-09 
    ResSS  =  233.5506 
RR-VGLM    linear loop  6 :  loglikelihood = -87.285862
   Alternating iteration 1 ,   Convergence criterion  =  1 
    ResSS  =  233.4989 
   Alternating iteration 2 ,   Convergence criterion  =  1.689959e-10 
    ResSS  =  233.4989 
RR-VGLM    linear loop  7 :  loglikelihood = -87.285862
                         logitlink(pstr0) loglink(lambda)
(Intercept)                    -2.5901076       2.7264838
sm.bs(WaterCon, df = 3)1        0.8885816      -0.4154863
sm.bs(WaterCon, df = 3)2       -0.7413091       0.3466240
sm.bs(WaterCon, df = 3)3        5.6709152      -2.6516274
A matrix:
                    latvar
logitlink(pstr0) -2.138655
loglink(lambda)   1.000000

C matrix:
                             latvar
sm.bs(WaterCon, df = 3)1 -0.4154863
sm.bs(WaterCon, df = 3)2  0.3466240
sm.bs(WaterCon, df = 3)3 -2.6516274

B1 matrix:
            logitlink(pstr0) loglink(lambda)
(Intercept)        -2.590108        2.726484

Call:
rrvglm(formula = Alopacce ~ sm.bs(WaterCon, df = 3), family = zipoisson(zero = NULL), 
    data = hspider, trace = TRUE, Index.corner = 2)

Pearson residuals:
                     Min       1Q    Median      3Q     Max
logitlink(pstr0) -2.7421 -0.41788 -0.290225 0.70877  1.7665
loglink(lambda)  -3.2187 -1.06371 -0.039912 0.20871 10.7146

Coefficients: 
                         Estimate Std. Error z value  Pr(>|z|)    
I(latvar.mat)            -2.13876    0.98189 -2.1782 0.0146951 *  
(Intercept):1            -2.59011    1.04668 -2.4746 0.0066694 ** 
(Intercept):2             2.72648    0.15195 17.9432 < 2.2e-16 ***
sm.bs(WaterCon, df = 3)1 -0.41549    0.75021 -0.5538 0.2898476    
sm.bs(WaterCon, df = 3)2  0.34662    0.86246  0.4019 0.3438791    
sm.bs(WaterCon, df = 3)3 -2.65163    0.80844 -3.2799 0.0005192 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Names of linear predictors: logitlink(pstr0), loglink(lambda)

Log-likelihood: -87.28586 on 49 degrees of freedom

Number of Fisher scoring iterations: 7 

VGAM documentation built on Jan. 16, 2021, 5:21 p.m.