This function fit a Zero Inflated Bivariate Poisson model under Geoffroy et. al (2021).

ZIBP_Geoffroy(
  l1.fo,
  l2.fo,
  psi.fo,
  data,
  initial.values = NULL,
  optimizer = "nlminb"
)

Arguments

l1.fo

a formula object to explain the parameter \(\lambda_1\), the first response \(X_1\) on the left of an ~ operator, and the terms, separated by operators, on the right.

l2.fo

a formula object to explain the parameter \(\lambda_2\), the second response \(X_2\) on the left of an ~ operator, and the terms, separated by operators, on the right.

psi.fo

a formula object to explain the \(\psi\) parameter. An example could be: ~ age + salary.

data

a data frame containing the variables occurring in the formula.

initial.values

a vector with possible initial values for the parameters.

optimizer

the optimizer to be used: nlminb, optim.

Value

Returns a object with class "ZIBPGeoffroy".

References

Kouakou, K. J. G., Hili, O., & Dupuy, J. F. (2021). Estimation in the zero-inflated bivariate Poisson model with an application to health-care utilization data. Afrika Statistika, 16(2), 2767-2788.

Author

Freddy Hernandez-Barajas, fhernanb@unal.edu.co

Examples

# Example 1 ---------------------------------------------------------------
# Estimating parameters as a regression model

l1 <- 5
l2 <- 3
psi <- 0.20
l0 <- 1.20

set.seed(12345)
data1 <- rZIBP_Geoffroy(n=500, l1=l1, l2=l2, l0=l0, psi=psi)
data1 <- as.data.frame(data1)
colnames(data1) <- c("y1", "y2")

mod <- ZIBP_Geoffroy(l1.fo=y1~1,
                     l2.fo=y2~1,
                     psi.fo=~1,
                     data=data1)

summary(mod)
#> ---------------------------------------------------------------
#> Fixed effects for log(l1) 
#> ---------------------------------------------------------------
#>             Estimate Std. Error z value  Pr(>|z|)    
#> (Intercept) 1.552910   0.054075  28.718 < 2.2e-16 ***
#> ---
#> Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
#> ---------------------------------------------------------------
#> Fixed effects for log(l2) 
#> ---------------------------------------------------------------
#>             Estimate Std. Error z value  Pr(>|z|)    
#> (Intercept) 1.043192   0.086603  12.046 < 2.2e-16 ***
#> ---
#> Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
#> ---------------------------------------------------------------
#> Fixed effects for logit(psi) 
#> ---------------------------------------------------------------
#>             Estimate Std. Error z value  Pr(>|z|)    
#> (Intercept) -1.26625    0.10801 -11.724 < 2.2e-16 ***
#> ---
#> Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
#> ---------------------------------------------------------------
#> Estimation for l0 
#> ---------------------------------------------------------------
#>             Estimate Std. Error z value  Pr(>|z|)    
#> (Intercept)  1.39196    0.23802  5.8482 4.969e-09 ***
#> ---
#> Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
#> ---------------------------------------------------------------

c(exp(mod$par[1:3]), mod$par[4])
#>    l1_hat    l2_hat    l0_hat   psi_hat 
#> 4.7251993 2.8382636 0.2818859 1.3919595 

moments_estim_ZIBP_Geoffroy(data1)
#>    l1_hat    l2_hat    l0_hat   psi_hat 
#> 4.8858732 3.3560129 1.1656262 0.2114351 


# Example 2 ---------------------------------------------------------------
# Replicating application from section 5.1 from Geoffroy et. al (2021).
if (FALSE) { # \dontrun{
  library(AER)
  data("NMES1988")

  # To obtain some statistics as in the paper
  y1 <- NMES1988$nvisits
  y2 <- NMES1988$novisits

  cor(y1, y2)
  cor.test(y1, y2)
  cov(y1, y2)
  mean(y1==0 & y2==0)

  # Model as in section 5.1 from Geoffroy et. al (2021).
  mod3 <- NULL
  mod3 <- ZIBP_Geoffroy(l1.fo=nvisits ~health+chronic+age+gender+
                          married+school+income+medicaid,
                        l2.fo=novisits~health+chronic+age+gender+
                          married+school+income+medicaid,
                        psi.fo=~chronic+gender+school+medicaid,
                        data=NMES1988)

  summary(mod3)
} # }