Binary logistic regression with two dependent variables

I have a continuous independent variable, X, and two binary dependent variables, Y and Z. I'm trying to run a binary logistic regression that models the correlation between Y and Z. X = age, while Y and Z are two pass/fail tasks. I want to determine the probability of an individual passing a given task at a given age, but I also want to test whether individuals who pass one task tend to also pass the other.

How should I go about doing this? I've tried the following in R:

  glm(cbind(Z,Y)~X, family="binomial", data= my_data)


This gives me a coefficient of -0.02335 for X. However, when I reverse the order of (Z,Y) in the formula above, the X coefficient becomes positive. Why is that?

Does anyone know how to run a binary logistic regression in R with correlated dependent variables?

Logistic regression is an univariate model, if you have two target variables, then you need to use a different model (there are a number of choices, including the bivariate logistic regression model given by Palmgren (1989), "Regression Models for Bivariate Binary Responses", UW Biostatistics Working Paper Series, 101, and implemented by binom2.or in the VGAM package).