Zero-Altered Geometric Distribution
zageometric.RdFits a zero-altered geometric distribution based on a conditional model involving a Bernoulli distribution and a positive-geometric distribution.
Usage
zageometric(lpobs0 = "logitlink", lprob = "logitlink",
type.fitted = c("mean", "prob", "pobs0", "onempobs0"),
imethod = 1, ipobs0 = NULL, iprob = NULL, zero = NULL)
zageometricff(lprob = "logitlink", lonempobs0 = "logitlink",
type.fitted = c("mean", "prob", "pobs0", "onempobs0"),
imethod = 1, iprob = NULL, ionempobs0 = NULL, zero = "onempobs0")Arguments
- lpobs0
Link function for the parameter \(p_0\) or \(\phi\), called
pobs0orphihere. SeeLinksfor more choices.- lprob
Parameter link function applied to the probability of success, called
probor \(p\). SeeLinksfor more choices.- type.fitted
See
CommonVGAMffArgumentsandfittedvlmfor information.
Details
The response \(Y\) is zero with probability \(p_0\), or \(Y\) has a positive-geometric distribution with probability \(1-p_0\). Thus \(0 < p_0 < 1\), which is modelled as a function of the covariates. The zero-altered geometric distribution differs from the zero-inflated geometric distribution in that the former has zeros coming from one source, whereas the latter has zeros coming from the geometric distribution too. The zero-inflated geometric distribution is implemented in the VGAM package. Some people call the zero-altered geometric a hurdle model.
The input can be a matrix (multiple responses).
By default, the two linear/additive predictors
of zageometric
are \((logit(\phi), logit(p))^T\).
The VGAM family function zageometricff() has a few
changes compared to zageometric().
These are:
(i) the order of the linear/additive predictors is switched so the
geometric probability comes first;
(ii) argument onempobs0 is now 1 minus the probability of an observed 0,
i.e., the probability of the positive geometric distribution,
i.e., onempobs0 is 1-pobs0;
(iii) argument zero has a new default so that the pobs0
is intercept-only by default.
Now zageometricff() is generally recommended over
zageometric().
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,
and vgam.
The fitted.values slot of the fitted object,
which should be extracted by the generic function fitted, returns
the mean \(\mu\) (default) which is given by
$$\mu = (1-\phi) / p.$$
If type.fitted = "pobs0" then \(p_0\) is returned.
Warning
Convergence for this VGAM family function seems to depend quite strongly on providing good initial values.
Inference obtained from summary.vglm and summary.vgam
may or may not be correct. In particular, the p-values, standard errors
and degrees of freedom may need adjustment. Use simulation on artificial
data to check that these are reasonable.
Note
Note this family function allows \(p_0\) to be modelled as functions of the covariates. It is a conditional model, not a mixture model.
This family function effectively combines
binomialff and
posgeometric() and geometric into
one family function.
However, posgeometric() is not written because it
is trivially related to geometric.
Examples
zdata <- data.frame(x2 = runif(nn <- 1000))
zdata <- transform(zdata, pobs0 = logitlink(-1 + 2*x2, inverse = TRUE),
prob = logitlink(-2 + 3*x2, inverse = TRUE))
zdata <- transform(zdata, y1 = rzageom(nn, prob = prob, pobs0 = pobs0),
y2 = rzageom(nn, prob = prob, pobs0 = pobs0))
with(zdata, table(y1))
#> y1
#> 0 1 2 3 4 5 6 7 8 9 10 11 12 13 14 16 17 18 19 20
#> 549 154 91 64 32 19 26 6 12 11 3 7 1 6 5 2 1 3 1 2
#> 21 23 24
#> 1 2 2
fit <- vglm(cbind(y1, y2) ~ x2, zageometric, data = zdata, trace = TRUE)
#> Iteration 1: loglikelihood = -3369.0164
#> Iteration 2: loglikelihood = -3254.4975
#> Iteration 3: loglikelihood = -3230.5896
#> Iteration 4: loglikelihood = -3229.5065
#> Iteration 5: loglikelihood = -3229.5037
#> Iteration 6: loglikelihood = -3229.5037
coef(fit, matrix = TRUE)
#> logitlink(pobs01) logitlink(prob1) logitlink(pobs02)
#> (Intercept) -0.7900987 -1.943901 -1.073231
#> x2 2.0067554 2.894466 2.088027
#> logitlink(prob2)
#> (Intercept) -2.182907
#> x2 3.377365
head(fitted(fit))
#> y1 y2
#> 1 2.3191558 2.7993061
#> 2 1.4953437 1.7279799
#> 3 0.3387585 0.3702204
#> 4 4.4960752 5.8692034
#> 5 0.3276168 0.3582084
#> 6 1.8582491 2.1917790
head(predict(fit))
#> logitlink(pobs01) logitlink(prob1) logitlink(pobs02) logitlink(prob2)
#> [1,] -0.23746287 -1.1468000 -0.4982145 -1.2528216
#> [2,] 0.05268020 -0.7283088 -0.1963209 -0.7645114
#> [3,] 1.16003539 0.8688975 0.9558809 1.0991652
#> [4,] -0.66273067 -1.7601898 -0.9407052 -1.9685465
#> [5,] 1.18828091 0.9096377 0.9852704 1.1467023
#> [6,] -0.09216374 -0.9372261 -0.3470309 -1.0082835
summary(fit)
#>
#> Call:
#> vglm(formula = cbind(y1, y2) ~ x2, family = zageometric, data = zdata,
#> trace = TRUE)
#>
#> Coefficients:
#> Estimate Std. Error z value Pr(>|z|)
#> (Intercept):1 -0.79010 0.13428 -5.884 4.01e-09 ***
#> (Intercept):2 -1.94390 0.09952 -19.532 < 2e-16 ***
#> (Intercept):3 -1.07323 0.13809 -7.772 7.72e-15 ***
#> (Intercept):4 -2.18291 0.09429 -23.151 < 2e-16 ***
#> x2:1 2.00676 0.24131 8.316 < 2e-16 ***
#> x2:2 2.89447 0.23344 12.399 < 2e-16 ***
#> x2:3 2.08803 0.24169 8.639 < 2e-16 ***
#> x2:4 3.37737 0.22412 15.069 < 2e-16 ***
#> ---
#> Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
#>
#> Names of linear predictors: logitlink(pobs01), logitlink(prob1),
#> logitlink(pobs02), logitlink(prob2)
#>
#> Log-likelihood: -3229.504 on 3992 degrees of freedom
#>
#> Number of Fisher scoring iterations: 6
#>