The bivariate response comprises \(Y_1\)
  from the linear model having parameters
  mean and sd for
  \(\mu_1\) and \(\sigma_1\),
  and the binary
  \(Y_2\) having parameter
  prob for its mean \(\mu_2\).
  The
  joint probability density/mass function is
  \(P(y_1, Y_2 = 0) = \phi_1(y_1; \mu_1, \sigma_1)
       (1 - \Delta)\)
     and
  \(P(y_1, Y_2 = 1) = \phi_1(y_1; \mu_1, \sigma_1)
       \Delta\)
     where \(\Delta\) adjusts \(\mu_2\)
     according to the association parameter
     \(\alpha\).
     The quantity \(\Delta\) is
     \(\Phi((\Phi^{-1}(\mu_2) - \alpha Z_1)/
       \sqrt{1 - \alpha^2})\).
     The quantity \(Z_1\) is \((Y_1-\mu_1) / \sigma_1\).
     Thus there is an underlying bivariate normal
     distribution, and a copula is used to bring the
     two marginal distributions together.
  Here,
  \(-1 < \alpha < 1\), and
  \(\Phi\) is the
  cumulative distribution function
  pnorm
  of a standard univariate normal.
The first marginal
  distribution is a normal distribution
  for the linear model.
  The second column of the response must
  have values 0 or 1,
  e.g.,
  Bernoulli random variables.
  When \(\alpha = 0\)
  then \(\Delta=\mu_2\).
  Together, this family function combines
  uninormal and
  binomialff.
  If the response are correlated then
  a more efficient joint analysis
  should result.
This VGAM family function cannot handle
  multiple responses. Only a two-column
  matrix is allowed.
  The two-column fitted
  value matrix has columns \(\mu_1\)
  and \(\mu_2\).