*Built using Zelig version 5.1.2.9000*

Bivariate Logistic Regression for Two Dichotomous Dependent Variables with `blogit`

from ZeligChoice.

Use the bivariate logistic regression model if you have two binary dependent variables \((Y_1, Y_2)\), and wish to model them jointly as a function of some explanatory variables. Each pair of dependent variables \((Y_{i1}, Y_{i2})\) has four potential outcomes, \((Y_{i1}=1, Y_{i2}=1)\), \((Y_{i1}=1, Y_{i2}=0)\), \((Y_{i1}=0, Y_{i2}=1)\), and \((Y_{i1}=0, Y_{i2}=0)\). The joint probability for each of these four outcomes is modeled with three systematic components: the marginal Pr \((Y_{i1} = 1)\) and Pr \((Y_{i2} = 1)\), and the odds ratio \(\psi\), which describes the dependence of one marginal on the other. Each of these systematic components may be modeled as functions of (possibly different) sets of explanatory variables.

In every bivariate logit specification, there are three equations which correspond to each dependent variable (\(Y_1\), \(Y_2\)), and \(\psi\), the odds ratio. You should provide a list of formulas for each equation or, you may use `cbind()`

if the right hand side is the same for both equations

`formulae <- list(cbind(Y1, Y2), X1 + X2)`

which means that all the explanatory variables in equations 1 and 2 (corresponding to \(Y_1\) and \(Y_2\)) are included, but only an intercept is estimated (all explanatory variables are omitted) for equation 3 (\(\psi\)).

**Anticipated feature, not currently enabled:**

You may use the function `tag()`

to constrain variables across equations:

`formulae <- list(mu1 = y1 x1 + tag(x3, "x3"), mu2 = y2 ~ x2 + tag(x3, "x3"))`

where `tag()`

is a special function that constrains variables to have the same effect across equations. Thus, the coefficient for x3 in equation mu1 is constrained to be equal to the coefficient for x3 in equation mu2.

Load the data and estimate the model:

`data(sanction)`

```
z.out1 <- zelig(cbind(import, export) ~ coop + cost + target,
model = "blogit", data = sanction)
```

```
## How to cite this model in Zelig:
## Thomas W. Yee. 2007.
## blogit: Bivariate Logit Regression for Dichotomous Dependent Variables
## in Christine Choirat, Christopher Gandrud, James Honaker, Kosuke Imai, Gary King, and Olivia Lau,
## "Zelig: Everyone's Statistical Software," http://zeligproject.org/
```

`summary(z.out1)`

```
## Model:
##
## Call:
## z5$zelig(formula = cbind(import, export) ~ coop + cost + target,
## data = sanction)
##
##
## Pearson residuals:
## Min 1Q Median 3Q Max
## logit(mu1) -2.184 -0.5255 -0.27618 0.4669 2.112
## logit(mu2) -4.924 -0.4599 0.08869 0.6242 3.079
## loge(oratio) -2.122 -0.4145 0.09651 0.2451 4.750
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept):1 -2.88048 1.09504 -2.630 0.008527
## (Intercept):2 -3.28639 1.05862 -3.104 0.001907
## (Intercept):3 -0.63660 0.73944 -0.861 0.389283
## coop:1 0.28705 0.33854 0.848 0.396501
## coop:2 -0.06214 0.33798 -0.184 0.854117
## cost:1 2.43481 0.66010 3.689 0.000226
## cost:2 2.72058 0.70009 3.886 0.000102
## target:1 -1.23994 0.48258 -2.569 0.010188
## target:2 -0.49957 0.44868 -1.113 0.265533
##
## Number of linear predictors: 3
##
## Names of linear predictors: logit(mu1), logit(mu2), loge(oratio)
##
## Residual deviance: 143.9813 on 225 degrees of freedom
##
## Log-likelihood: -71.9906 on 225 degrees of freedom
##
## Number of iterations: 8
##
## Odds ratio: 0.5291
## Next step: Use 'setx' method
```

By default, zelig() estimates two effect parameters for each explanatory variable in addition to the odds ratio parameter; this formulation is parametrically independent (estimating unconstrained effects for each explanatory variable), but stochastically dependent because the models share an odds ratio. Generate baseline values for the explanatory variables (with cost set to 1, net gain to sender) and alternative values (with cost set to 4, major loss to sender):

Simulate fitted values and first differences:

`s.out1 <- sim(z.out1, x = x.low, x1 = x.high)`

`summary(s.out1)`

```
##
## sim x :
## -----
## ev
## mean sd 50% 2.5% 97.5%
## Pr(Y1=0, Y2=0) 0.766093739 0.075317101 0.77358310 0.602784319 0.88645025
## Pr(Y1=0, Y2=1) 0.150966892 0.067414468 0.13886736 0.048405075 0.30918751
## Pr(Y1=1, Y2=0) 0.074055394 0.043283243 0.06410452 0.018120382 0.18385203
## Pr(Y1=1, Y2=1) 0.008883975 0.008944314 0.00607899 0.000926747 0.03495784
## pv
## 0 1
## (Y1=0, Y2=0) 0.229 0.771
## (Y1=0, Y2=1) 0.863 0.137
## (Y1=1, Y2=0) 0.918 0.082
## (Y1=1, Y2=1) 0.990 0.010
##
## sim x1 :
## -----
## ev
## mean sd 50% 2.5%
## Pr(Y1=0, Y2=0) 6.856656e-05 0.000289293 8.438768e-06 1.130919e-07
## Pr(Y1=0, Y2=1) 2.190235e-02 0.037090302 9.189970e-03 6.227473e-04
## Pr(Y1=1, Y2=0) 6.068223e-03 0.013777284 1.676166e-03 6.898759e-05
## Pr(Y1=1, Y2=1) 9.719609e-01 0.038579035 9.849495e-01 8.793918e-01
## 97.5%
## Pr(Y1=0, Y2=0) 0.0005408712
## Pr(Y1=0, Y2=1) 0.1141765919
## Pr(Y1=1, Y2=0) 0.0465000275
## Pr(Y1=1, Y2=1) 0.9983380073
## pv
## 0 1
## (Y1=0, Y2=0) 1.000 0.000
## (Y1=0, Y2=1) 0.989 0.011
## (Y1=1, Y2=0) 0.989 0.011
## (Y1=1, Y2=1) 0.022 0.978
## fd
## mean sd 50% 2.5% 97.5%
## Pr(Y1=0, Y2=0) -0.76602517 0.07539928 -0.77345473 -0.8864501 -0.60273253
## Pr(Y1=0, Y2=1) -0.12906454 0.08036197 -0.12107038 -0.3010280 0.01301308
## Pr(Y1=1, Y2=0) -0.06798717 0.04695905 -0.05992501 -0.1749069 -0.00240991
## Pr(Y1=1, Y2=1) 0.96307689 0.04264955 0.97673796 0.8552438 0.99647840
```

`plot(s.out1)`

For each observation, define two binary dependent variables, \(Y_1\) and \(Y_2\), each of which take the value of either 0 or 1 (in the following, we suppress the observation index). We model the joint outcome \((Y_1\), \(Y_2)\) using a marginal probability for each dependent variable, and the odds ratio, which parameterizes the relationship between the two dependent variables. Define \(Y_{rs}\) such that it is equal to 1 when \(Y_1=r\) and \(Y_2=s\) and is 0 otherwise, where \(r\) and \(s\) take a value of either 0 or 1. Then, the model is defined as follows,

- The
*stochastic component*is

\[ \begin{aligned} Y_{11} &\sim& \textrm{Bernoulli}(y_{11} \mid \pi_{11}) \\ Y_{10} &\sim& \textrm{Bernoulli}(y_{10} \mid \pi_{10}) \\ Y_{01} &\sim& \textrm{Bernoulli}(y_{01} \mid \pi_{01}) \end{aligned} \]

where \(\pi_{rs}=\Pr(Y_1=r, Y_2=s)\) is the joint probability, and \(\pi_{00}=1-\pi_{11}-\pi_{10}-\pi_{01}\).

- The
*systematic components*model the marginal probabilities, \(\pi_j=\Pr(Y_j=1)\), as well as the odds ratio. The odds ratio is defined as \(\psi = \pi_{00} \pi_{01}/\pi_{10}\pi_{11}\) and describes the relationship between the two outcomes. Thus, for each observation we have

\[ \begin{aligned} \pi_j & = & \frac{1}{1 + \exp(-x_j \beta_j)} \quad \textrm{ for} \quad j=1,2, \\ \psi &= & \exp(x_3 \beta_3). \end{aligned} \]

- The expected values (qi$ev) for the bivariate logit model are the predicted joint probabilities. Simulations of \(\beta_1\), \(\beta_2\), and \(\beta_3\) (drawn from their sampling distributions) are substituted into the systematic components \((\pi_1, \pi_2, \psi)\) to find simulations of the predicted joint probabilities:

\[ \begin{aligned} \pi_{11} & = & \left\{ \begin{array}{ll} \frac{1}{2}(\psi - 1)^{-1} - {a - \sqrt{a^2 + b}} & \textrm{for} \; \psi \ne 1 \\ \pi_1 \pi_2 & \textrm{for} \; \psi = 1 \end{array} \right., \\ \pi_{10} &=& \pi_1 - \pi_{11}, \\ \pi_{01} &=& \pi_2 - \pi_{11}, \\ \pi_{00} &=& 1 - \pi_{10} - \pi_{01} - \pi_{11}, \end{aligned} \]

where \(a = 1 + (\pi_1 + \pi_2)(\psi - 1)\), \(b = -4 \psi(\psi - 1) \pi_1 \pi_2\), and the joint probabilities for each observation must sum to one. For \(n\) simulations, the expected values form an \(n \times 4\) matrix for each observation in x.

The predicted values (

`qi$pr`

) are draws from the multinomial distribution given the expected joint probabilities.The first differences (

`qi$fd`

) for each of the predicted joint probabilities are given by

\[ \textrm{FD}_{rs} = \Pr(Y_1=r, Y_2=s \mid x_1)-\Pr(Y_1=r, Y_2=s \mid x). \]

- The risk ratio (
`qi$rr`

) for each of the predicted joint probabilities are given by

\[ \textrm{RR}_{rs} = \frac{\Pr(Y_1=r, Y_2=s \mid x_1)}{\Pr(Y_1=r, Y_2=s \mid x)} \]

- In conditional prediction models, the average expected treatment effect (
`att.ev`

) for the treatment group is

\[ \frac{1}{\sum_{i=1}^n t_i}\sum_{i:t_i=1}^n \left\{ Y_{ij}(t_i=1) - E[Y_{ij}(t_i=0)] \right\} \textrm{ for } j = 1,2, \]

where \(t_i\) is a binary explanatory variable defining the treatment (\(t_i=1\)) and control (\(t_i=0\)) groups. Variation in the simulations are due to uncertainty in simulating \(E[Y_{ij}(t_i=0)]\), the counterfactual expected value of \(Y_{ij}\) for observations in the treatment group, under the assumption that everything stays the same except that the treatment indicator is switched to \(t_i=0\).

- In conditional prediction models, the average predicted treatment effect (att.pr) for the treatment group is

\[ \frac{1}{\sum_{i=1}^n t_i}\sum_{i:t_i=1}^n \left\{ Y_{ij}(t_i=1) - \widehat{Y_{ij}(t_i=0)} \right\} \textrm{ for } j = 1,2, \]

where \(t_i\) is a binary explanatory variable defining the treatment (\(t_i=1\)) and control (\(t_i=0\)) groups. Variation in the simulations are due to uncertainty in simulating \(\widehat{Y_{ij}(t_i=0)}\), the counterfactual predicted value of \(Y_{ij}\) for observations in the treatment group, under the assumption that everything stays the same except that the treatment indicator is switched to \(t_i=0\).

The output of each Zelig command contains useful information which you may view. For example, if you run `z.out <- zelig(y ~ x, model = "blogit", data)`

, then you may examine the available information in `z.out`

by using `names(z.out)`

, see the coefficients by using z.out$coefficients, and obtain a default summary of information through summary(z.out). Other elements available through the $ operator are listed below.

From the

`zelig()`

output object`z.out`

, you may extract:coefficients: the named vector of coefficients.

fitted.values: an \(n \times 4\) matrix of the in-sample fitted values.

predictors: an \(n \times 3\) matrix of the linear predictors \(x_j \beta_j\).

residuals: an \(n \times 3\) matrix of the residuals.

df.residual: the residual degrees of freedom.

df.total: the total degrees of freedom.

rss: the residual sum of squares.

y: an \(n \times 2\) matrix of the dependent variables.

zelig.data: the input data frame if

`save.data = TRUE`

.From summary(z.out), you may extract:

coef3: a table of the coefficients with their associated standard errors and \(t\)-statistics.

cov.unscaled: the variance-covariance matrix.

pearson.resid: an \(n \times 3\) matrix of the Pearson residuals.

The bivariate logit function is part of the VGAM package by Thomas Yee . In addition, advanced users may wish to refer to `help(vglm)`

in the VGAM library.

Yee TW (2015). *Vector Generalized Linear and Additive Models: With an Implementation in R*. Springer, New York, USA.

Yee TW and Wild CJ (1996). “Vector Generalized Additive Models.” *Journal of Royal Statistical Society, Series B*, **58** (3), pp. 481-493.

Yee TW (2010). “The VGAM Package for Categorical Data Analysis.” *Journal of Statistical Software*, **32** (10), pp. 1-34. <URL: http://www.jstatsoft.org/v32/i10/>.

Yee TW and Hadi AF (2014). “Row-column interaction models, with an R implementation.” *Computational Statistics*, **29** (6), pp. 1427-1445.

Yee TW (2017). *VGAM: Vector Generalized Linear and Additive Models*. R package version 1.0-3, <URL: https://CRAN.R-project.org/package=VGAM>.

Yee TW (2013). “Two-parameter reduced-rank vector generalized linear models.” *Computational Statistics and Data Analysis*. <URL: http://ees.elsevier.com/csda>.

Yee TW, Stoklosa J and Huggins RM (2015). “The VGAM Package for Capture-Recapture Data Using the Conditional Likelihood.” *Journal of Statistical Software*, **65** (5), pp. 1-33. <URL: http://www.jstatsoft.org/v65/i05/>.