source("04 BP Laksh 1999.R")Bivariate Poisson distribution
Bivariate Poisson
Lakshminarayana, Pandit, and Srinivasa Rao (1999) proposed a version of the bivariate Poisson that can deal with negative or positive values for the correlation of the two variables \(X_1\) and \(X_2\). This distribution has the joint probability mass function:
\[ f_{BP}(x_1, x_2 | \lambda_1, \lambda_2, \alpha) = \frac{e^{-(\lambda_1+\lambda_2)} \lambda_1^{x_1} \lambda_2^{x_2} }{x_1! \, x_2!} \left[ 1+\alpha (e^{-x_1}-e^{-\lambda_1 c}) (e^{-x_2}-e^{-\lambda_2 c}) \right], \]
where the distribution parameters are \(\lambda_1>0\), \(\lambda_2>0\) and \(|\alpha| \leq \frac{1}{(1-A)(1-B)}\) with \(A=\exp(-\lambda_1c)\), \(B=\exp(-\lambda_2 c)\) and \(c=1-\exp(-1)\). The parameter \(\alpha\) must fall within this interval to ensure that \(f(x_1, x_2 | \lambda_1, \lambda_2, \alpha)\) constitutes a proper probability mass function.
Some characteristics of this distribution include \(E(X_1)=\lambda_1\), \(E(X_2)=\lambda_2\), \(Cov(X_1,X_2)=\alpha \lambda_1 \lambda_2 c^2 A B\), and \(Cor(X_1,X_2)=\alpha \sqrt{\lambda_1 \lambda_2} c^2 e^{-(\lambda_1 + \lambda_2)c}\). The covariance and correlation can take on negative or positive values depending on the sign of the \(\alpha\) parameter.
The dBP function can be used to calculate probabilities, the structure function is given below.
dBP(x, l1, l2, alpha, log=FALSE)Next we present an example using the dBP function to obtain the density for three pairs \((x_1, x_2)\) given the parameters.
X <- matrix(c(2, 1,
1, 3,
3, 4), ncol=2, byrow=TRUE)
dBP(x=X, l1=3, l2=4, alpha=1)[1] 0.01634399 0.02898963 0.04404024Another useful function is probability_grid_BP that can be used to obtain a grid or matrix with the probabilities for each combination of the parameters. Next, we present an example for the \(BP(\lambda_1=3, \lambda_2=4, \alpha=-1.23)\) obtaining the probabilities for the combinations of \(X_1=0, 1, 2, 3\) and \(X_2=0, 1, 2, 3, 4, 5\).
probability_grid_BP(l1=3, l2=4, alpha=-1.23,
max_val_x1=3,
max_val_x2=5) 0 1 2 3 4 5
0 0.0000346869 0.002549011 0.006871393 0.01003171 0.01035171 0.008375543
1 0.0020613570 0.010098167 0.021559503 0.02941465 0.02966063 0.023800894
2 0.0041721090 0.016499834 0.032860902 0.04374647 0.04372143 0.034969775
3 0.0045694459 0.016997422 0.033052805 0.04360833 0.04343834 0.034700647The probability_grid_BP function can be used to create heat maps with the probabilities for any \(BP(\lambda_1, \lambda_2, \alpha)\). Above you can find an example.
freq_table <- probability_grid_BP(l1=3, l2=4, alpha=-1.23,
max_val_x1 = 10,
max_val_x2 = 10)
filled.contour(x=as.numeric(rownames(freq_table)),
y=as.numeric(colnames(freq_table)),
z=freq_table,
nlevels=20,
xlab="X1", ylab="X2",
color = topo.colors)
As the parameter \(\alpha\) depends on the other parameters \(\lambda_1\) and \(\lambda_2\), it is important to obtain the interval of admissible values for \(\alpha\). The correct_alpha_BP function can be used to obtain automatically the interval, next there is an example assuming \(\lambda_1=3\) and \(\lambda_2=4\).
correct_alpha_BP(l1=3, l2=4)$min_alpha
[1] -1.278638
$max_alpha
[1] 1.278638Another useful function is rBP which can be used to generate random samples from the BP distribution given the parameter \(\lambda_1\), \(\lambda_2\) and \(\alpha\). Next we present an example using the rBP function to simulate some pairs \((x_1, x_2)\) given the parameters.
rBP(n=5, l1=3, l2=4, alpha=-1.23) X1 X2
1 1 2
2 3 5
3 3 2
4 4 3
5 4 5The moments_estim_BP function calculates the moment estimators for the BP. Next, we present an example simulating 5000 values from a \(BP(\lambda_1=3, \lambda_2=4, \alpha=-1.23)\) and then we use the moments_estim_BP function to obtain moment estimators.
x <- rBP(n=5000, l1=3, l2=4, alpha=-1.23)
moments_estim_BP(x) l1_hat l2_hat alpha_hat
3.0344 4.0066 -0.8860 From the last output we observe that estimated values are close to \((\lambda_1=3, \lambda_2=4, \alpha=-1.23)^\top\).
The llBP is the log-likelihood function that can be used with the optim function to estimate parameters using the maximum likelihood method.
In the next example we used the last random sample x. We use the moments estimators as initial values for the optim function.
# To obtain reasonable values for alpha
est <- as.numeric(moments_estim_BP(x))
# To obtain reasonable limits for alpha
min_alpha <- correct_alpha_BP(l1=est[1], l2=est[2])$min_alpha
max_alpha <- correct_alpha_BP(l1=est[1], l2=est[2])$max_alpha
# Estimating parameters
optim(fn = llBP,
par = moments_estim_BP(x),
lower = c(0.001, 0.001, min_alpha),
upper = c( Inf, Inf, max_alpha),
method = "L-BFGS-B",
control = list(maxit=100000, fnscale=-1),
x=x)$par
l1_hat l2_hat alpha_hat
3.036020 4.006553 -1.273339
$value
[1] -20038.48
$counts
function gradient
12 12
$convergence
[1] 0
$message
[1] "CONVERGENCE: REL_REDUCTION_OF_F <= FACTR*EPSMCH"From the last output we observe that estimated values are close to \((\lambda_1=3, \lambda_2=4, \alpha=-1.23)^\top\).
Zero Inflated Bivariate Poisson
In this section we propose a version of the bivariate Poisson inflated with zeros, we denote this distribution as ZIBP. This distribution has the joint probability mass function:
\[ f_{ZIBP}(x_1, x_2 | \lambda_1, \lambda_2, \alpha, \psi) = \begin{cases} \psi + (1-\psi) e^{-(\lambda_1+\lambda_2)} \left[ 1+\alpha (1-e^{-\lambda_1 c}) (1-e^{-\lambda_2 c}) \right] & \text{for } (x_1, x_2)=(0, 0) \\ (1-\psi) \frac{e^{-(\lambda_1+\lambda_2)} \lambda_1^{x_1} \lambda_2^{x_2} }{x_1! \, x_2!} \left[ 1+\alpha (e^{-x_1}-e^{-\lambda_1 c}) (e^{-x_2}-e^{-\lambda_2 c}) \right] & \text{for} (x_1, x_2) \neq (0, 0) \end{cases} \]
The new parameter for this distribution is \(\psi\) and corresponds to the probability of excess of zeros for an usual BP distribution. The parameter \(\psi\) has assume values in \(0 \leq \psi \leq 1\).
Some characteristics of this distribution include \(E(X_1)=\lambda_1 - 2\lambda_1 \psi + \lambda_1 \psi^2\) and \(E(X_2)=\lambda_2 - 2\lambda_2 \psi + \lambda_2 \psi^2\).
source("04 ZIBP Laksh 1999.R")The dZIBP function can be used to calculate probabilities, the structure function is given below.
dZIBP(x, l1, l2, alpha, psi, log=FALSE)Next we present an example using the dZIBP function to obtain the density for three pairs \((x_1, x_2)\) given the parameters.
X <- matrix(c(2, 1,
1, 3,
3, 4), ncol=2, byrow=TRUE)
dZIBP(x=X, l1=3, l2=4, alpha=1, psi=0.5)[1] 0.008171995 0.014494817 0.022020122Another useful function is probability_grid_ZIBP that can be used to obtain a grid or matrix with the probabilities for each combination of the parameters. Next, we present an example for the \(ZIBP(\lambda_1=3, \lambda_2=4, \alpha=-1.23, \psi=0.07)\) obtaining the probabilities for the combinations of \(X_1=0, 1, 2, 3\) and \(X_2=0, 1, 2, 3, 4, 5\).
probability_grid_ZIBP(l1=3, l2=4, alpha=-1.23, psi=0.07,
max_val_x1=3,
max_val_x2=5) 0 1 2 3 4 5
0 0.070032259 0.002370580 0.006390396 0.009329488 0.009627088 0.007789255
1 0.001917062 0.009391295 0.020050338 0.027355622 0.027584383 0.022134832
2 0.003880061 0.015344846 0.030560639 0.040684218 0.040660931 0.032521891
3 0.004249585 0.015807602 0.030739109 0.040555749 0.040397660 0.032271602The probability_grid_ZIBP function can be used to create heat maps with the probabilities for any \(ZIBP(\lambda_1, \lambda_2, \alpha, \psi)\). Above you can find an example.
freq_table <- probability_grid_ZIBP(l1=3, l2=4,
alpha=-1.23, psi=0.07,
max_val_x1 = 10,
max_val_x2 = 10)
# Para ver la densidad conjunta
filled.contour(x=as.numeric(rownames(freq_table)),
y=as.numeric(colnames(freq_table)),
z=freq_table,
nlevels=20,
xlab="X1", ylab="X2",
color = topo.colors)
Another useful function is rZIBP which can be used to generate random samples from the ZIBP distribution given the parameter \(\lambda_1\), \(\lambda_2\), \(\alpha\) and \(\psi\). Next we present an example using the rZIBP function to simulate some pairs \((x_1, x_2)\) given the parameters.
rZIBP(n=10, l1=3, l2=4, alpha=-1.23, psi=0.20) X1 X2
1 1 2
2 2 4
3 0 0
4 5 6
5 0 0
6 3 2
7 4 3
8 6 2
9 5 7
10 3 4The moments_estim_ZIBP function calculates the moment estimators for the BP. Next, we present an example simulating 5000 values from a \(ZIBP(\lambda_1=3, \lambda_2=4, \alpha=-1.23, \psi=0.15)\) and then we use the moments_estim_ZIBP function to obtain moment estimators.
x <- rZIBP(n=5000, l1=3, l2=4, alpha=-1.23, psi=0.15)
moments_estim_ZIBP(x) l1_hat l2_hat alpha_hat psi_hat
3.0117 4.0030 0.4537 0.0751 From the last output we observe that estimated values are close to \((\lambda_1=3, \lambda_2=4, \alpha=-1.23, \psi=0.15)^\top\).
The llZIBP is the log-likelihood function that can be used with the optim function to estimate parameters using the maximum likelihood method.
In the next example we used the last random sample x. We use the moments estimators as initial values for the optim function.
# To obtain reasonable initial values for theta
theta <- moments_estim_ZIBP(x)
# To replace the third moment estimator by zero
start_param <- c(theta[1:2], 0, theta[4])
# To obtain reasonable limits for alpha
min_alpha <- correct_alpha_BP(l1=theta[1], l2=theta[2])$min_alpha
max_alpha <- correct_alpha_BP(l1=theta[1], l2=theta[2])$max_alpha
# Estimating parameters
res <- optim(fn = llZIBP,
par = start_param,
lower = c(0.001, 0.001, min_alpha, 0.0001),
upper = c( Inf, Inf, max_alpha, 0.9999),
method = "L-BFGS-B",
control = list(maxit=100000, fnscale=-1),
x=x)
res$par
l1_hat l2_hat psi_hat
3.0115340 4.0027334 -1.2767662 0.1446092
$value
[1] -19259.52
$counts
function gradient
26 26
$convergence
[1] 0
$message
[1] "CONVERGENCE: REL_REDUCTION_OF_F <= FACTR*EPSMCH"From the last output we observe that estimated values are close to \((\lambda_1=3, \lambda_2=4, \alpha=-1.23, \psi=0.15)^\top\).
Marginal ZIBP model
In this section we present the Marginal Bivariate Poisson regression model to model the parameters \(\mu_1\), \(\mu_2\) and \(\psi\) using covariates for the ZIBP model. The model defined as follows:
\[\begin{align*} (Y_1, Y_2) &\sim~ ZIBP(\lambda_1, \lambda_2, \alpha, \psi) \\ log(\mu_1) &= \beta_0 + \beta_1 X_1 + \ldots + \beta_p X_p \\ log(\mu_2) &= \gamma_0 + \gamma_1 X_1 + \ldots + \gamma_q X_q \\ logit(\psi) &= \delta_0 + \delta_1 X_1 + \ldots + \delta_r X_r \\ \alpha &= \text{constant}, \end{align*}\]
where \(\mu_1=E(Y_1)\) and \(\mu_2=E(Y_2)\).
The main function is marZIBP, this function has the next structure:
marZIBP(mu1.fo, m2.fo, psi.fo,
data, initial.values=NULL)The arguments of this function are:
m1.fo: formula to indicate the variables to model the parameter \(\mu_1\) that corresponds to the \(E(X_1)\).mu2.fo: formula to indicate the variables to model the parameter \(\mu_2\) that corresponds to the \(E(X_2)\).psi.fo: formula to indicate the variables to model the parameter \(\psi\).data: an optional data frame, list or environment containing the variables in the model.initial.values: a vector with initial values for \(\boldsymbol{\Theta}\).
To load the main function and auxiliar functions we can use the next code.
source("04 marZIBP model.R")Example 1
# Example 1 - toy example with a random sample
datos1 <- rZIBP(n=120, l1=3, l2=2, alpha=1, psi=0.15)
datos1 <- as.data.frame(datos1)
head(datos1) X1 X2
1 5 0
2 5 1
3 3 2
4 4 0
5 2 1
6 0 1# To fit the model
mod1 <- marZIBP(mu1.fo=X1~1,
mu2.fo=X2~1,
psi.fo=~1,
data=datos1)
# To obtain the usual summary table
summary(mod1)---------------------------------------------------------------
Fixed effects for log(mu1)
---------------------------------------------------------------
Estimate Std. Error z value Pr(>|z|)
(Intercept) 0.7418288 0.0097813 75.841 < 2.2e-16 ***
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
---------------------------------------------------------------
Fixed effects for log(mu2)
---------------------------------------------------------------
Estimate Std. Error z value Pr(>|z|)
(Intercept) 0.249954 0.010609 23.561 < 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.402122 0.021411 -65.487 < 2.2e-16 ***
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
---------------------------------------------------------------
Estimation for alpha
---------------------------------------------------------------
Estimate Std. Error z value Pr(>|z|)
(Intercept) -0.66227 0.21757 -3.0439 0.002335 **
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
---------------------------------------------------------------# To explore the estimations of mu1, mu2, mu and p
# To obtain E(Y1)=mu1 and E(Y2)=mu2
mod1$fitted.mu1[1][1] 2.099772mod1$fitted.mu2[1][1] 1.283966# To compare sample means with mu1 and mu2
colMeans(datos1) X1 X2
2.616667 1.600000 # To obtain psi and alpha
mod1$fitted.psi[1][1] 0.1974796mod1$fitted.alpha[1] -0.6622689# To obtain l1 and l2
mod1$fitted.l1[1][1] 3.260318mod1$fitted.l2[1][1] 1.993616Example 2
# Example 2 - toy example using covariates
datos2 <- data.frame(y1=c(0, 5, 3, 2, 1),
y2=c(3, 2, 4, 1, 3),
x1=c(3, 4, 6, 9, 8),
x2=c(0, 3, 2, 0, 5))
datos2 y1 y2 x1 x2
1 0 3 3 0
2 5 2 4 3
3 3 4 6 2
4 2 1 9 0
5 1 3 8 5mod2 <- marZIBP(mu1.fo=y1~x1+x2,
mu2.fo=y2~x1+x2,
psi.fo=~x1+x2,
data=datos2)
summary(mod2)---------------------------------------------------------------
Fixed effects for log(mu1)
---------------------------------------------------------------
Estimate Std. Error z value Pr(>|z|)
(Intercept) 0.81934305 0.39003019 2.1007 0.03567 *
x1 -0.00042732 0.05980424 -0.0071 0.99430
x2 -0.00650096 0.07198976 -0.0903 0.92805
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
---------------------------------------------------------------
Fixed effects for log(mu2)
---------------------------------------------------------------
Estimate Std. Error z value Pr(>|z|)
(Intercept) 0.9652925 0.3635994 2.6548 0.007935 **
x1 -0.0028542 0.0559839 -0.0510 0.959339
x2 0.0028183 0.0669909 0.0421 0.966443
---
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) 0.60230 97.17597 0.0062 0.9951
x1 -1.73517 32.34637 -0.0536 0.9572
x2 -0.27982 12.59478 -0.0222 0.9823
---------------------------------------------------------------
Estimation for alpha
---------------------------------------------------------------
Estimate Std. Error z value Pr(>|z|)
(Intercept) -0.9881 4.4472 -0.2222 0.8242
---------------------------------------------------------------References
Lakshminarayana, J, SN Narahari Pandit, and K Srinivasa Rao. 1999. βOn a Bivariate Poisson Distribution.β Communications in Statistics-Theory and Methods 28 (2): 267β76.