Marginal Zero Inflated Bivariate Poisson Distribution - Geoffroy

This function fit a marginal Zero Inflated Bivariate Poisson distribution under the of Geoffroy et. al (2021).

marZIBP_Geoffroy(
  mu1.fo,
  mu2.fo,
  psi.fo,
  data,
  initial.values = NULL,
  optimizer = "nlminb"
)

Arguments

mu1.fo

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

mu2.fo

a formula object to explain the marginal mean \(\mu_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 formula2.

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 "marZIBPGeoffroy".

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 ---------------------------------------------------------------
l1 <- 3
l2 <- 4
l0 <- 1.5
psi <- 0.20

set.seed(12345678)
data1 <- rZIBP_Geoffroy(n=500, l1=l1, l2=l2, l0=l0, psi=psi)
#> Error in rZIBP_Geoffroy(n = 500, l1 = l1, l2 = l2, l0 = l0, psi = psi): could not find function "rZIBP_Geoffroy"
data1 <- as.data.frame(data1)
#> Error in eval(expr, envir, enclos): object 'data1' not found

# To fit the model
mod1 <- NULL
mod1 <- marZIBP_Geoffroy(mu1.fo=X1~1,
                         mu2.fo=X2~1,
                         psi.fo=~1,
                         data=data1)
#> Error in marZIBP_Geoffroy(mu1.fo = X1 ~ 1, mu2.fo = X2 ~ 1, psi.fo = ~1,     data = data1): could not find function "marZIBP_Geoffroy"

# To obtain the usual summary table
summary(mod1)
#> Length  Class   Mode 
#>      0   NULL   NULL 

# To explore the estimations of l1, l2, mu and p

# To obtain E(Y1)=v1 and E(Y2)=v2
mod1$fitted.mu1[1]
#> NULL
mod1$fitted.mu2[1]
#> NULL

# To compare sample means with v1 and v2
colMeans(data1)
#> Error in eval(expr, envir, enclos): object 'data1' not found

# To obtain l0 and psi
mod1$fitted.l0
#> NULL
mod1$fitted.psi[1]
#> NULL

# To obtain l1 and l2
mod1$fitted.l1[1]
#> NULL
mod1$fitted.l2[1]
#> NULL

# Example 2 ---------------------------------------------------------------

gen_data_ZIBP_Geoffroy <- function(n=100) {
  # To generate the covariates
  x1 <- runif(n=n)
  x2 <- runif(n=n)

  # To generate the means
  mu1 <- exp(1 + 1.3 * x1)
  mu2 <- exp(1 + 2.1 * x2)

  # To generate the psi
  logit_inv <- function(x) exp(x) / (1+exp(x))
  psi <- logit_inv(-2.4 + 2.1 * x2)

  # The third lambda
  l0 <- 1

  # To obtain lambdas
  l1 <- mu1 / (1-psi) - l0
  l2 <- mu2 / (1-psi) - l0

  # To generate Y1 and Y2
  y <- NULL
  for (i in 1:n)
    y <- rbind(y, rZIBP_Geoffroy(n=1, l1=l1[i], l2=l2[i],
                                 l0=l0, psi=psi[i]))

  # To create the dataset
  dataset <- data.frame(y1=y[,1], y2=y[,2],
                        x1=x1, x2=x2,
                        mu1=mu1, mu2=mu2,
                        l0=l0, psi=psi,
                        l1=l1, l2=l2)

  return(dataset)
}

set.seed(12)
data2 <- gen_data_ZIBP_Geoffroy(n=300)
#> Error in rZIBP_Geoffroy(n = 1, l1 = l1[i], l2 = l2[i], l0 = l0, psi = psi[i]): could not find function "rZIBP_Geoffroy"
head(data2, n=8)
#> Error in eval(expr, envir, enclos): object 'data2' not found

mod2 <- NULL
mod2 <- marZIBP_Geoffroy(mu1.fo=y1~x1,
                         mu2.fo=y2~x2,
                         psi.fo=~x2,
                         data=data2)
#> Error in marZIBP_Geoffroy(mu1.fo = y1 ~ x1, mu2.fo = y2 ~ x2, psi.fo = ~x2,     data = data2): could not find function "marZIBP_Geoffroy"

summary(mod2)
#> Length  Class   Mode 
#>      0   NULL   NULL