Bivariate Probit Model

Density and random generation for a bivariate probit model. The correlation parameter rho is the measure of dependency.

rbinom2.rho(n, mu1, mu2 = if (exchangeable) mu1 else
  stop("'mu2' not specified"), rho = 0,
  exchangeable = FALSE, twoCols = TRUE,
  colnames = if (twoCols) c("y1","y2") else
  c("00", "01", "10", "11"), ErrorCheck = TRUE)
dbinom2.rho(mu1, mu2 = if (exchangeable) mu1 else
  stop("'mu2' not specified"), rho = 0, exchangeable = FALSE,
  colnames = c("00", "01", "10", "11"), ErrorCheck = TRUE)

Arguments

n

number of observations. Same as in runif. The arguments mu1, mu2, rho are recycled to this value.

mu1, mu2

The marginal probabilities. Only mu1 is needed if exchangeable = TRUE. Values should be between 0 and 1.

rho

The correlation parameter. Must be numeric and lie between \(-1\) and \(1\). The default value of zero means the responses are uncorrelated.

exchangeable

Logical. If TRUE, the two marginal probabilities are constrained to be equal.

twoCols

Logical. If TRUE, then a \(n\) \(\times\) \(2\) matrix of 1s and 0s is returned. If FALSE, then a \(n\) \(\times\) \(4\) matrix of 1s and 0s is returned.

colnames

The dimnames argument of matrix is assigned list(NULL, colnames).

ErrorCheck

Logical. Do some error checking of the input parameters?

Details

The function rbinom2.rho generates data coming from a bivariate probit model. The data might be fitted with the VGAM family function binom2.rho.

The function dbinom2.rho does not really compute the density (because that does not make sense here) but rather returns the four joint probabilities.

Value

The function rbinom2.rho returns either a 2 or 4 column matrix of 1s and 0s, depending on the argument twoCols.

The function dbinom2.rho returns a 4 column matrix of joint probabilities; each row adds up to unity.

Author

T. W. Yee

See also

binom2.rho.

Examples

(myrho <- rhobitlink(2, inverse = TRUE))  # Example 1
#> [1] 0.7615942
nn <- 2000
ymat <- rbinom2.rho(nn, mu1 = 0.8, rho = myrho, exch = TRUE)
(mytab <- table(ymat[, 1], ymat[, 2], dnn = c("Y1", "Y2")))
#>    Y2
#> Y1     0    1
#>   0  232  147
#>   1  155 1466
fit <- vglm(ymat ~ 1, binom2.rho(exch = TRUE))
coef(fit, matrix = TRUE)
#>             probitlink(mu1) probitlink(mu2) rhobitlink(rho)
#> (Intercept)        0.872382        0.872382        1.998258

bdata <- data.frame(x2 = sort(runif(nn)))  # Example 2
bdata <- transform(bdata, mu1 = probitlink(-2+4*x2, inv = TRUE),
                          mu2 = probitlink(-1+3*x2, inv = TRUE))
dmat <- with(bdata, dbinom2.rho(mu1, mu2, myrho))
ymat <- with(bdata, rbinom2.rho(nn, mu1, mu2, myrho))
fit2 <- vglm(ymat ~ x2, binom2.rho, data = bdata)
coef(fit2, matrix = TRUE)
#>             probitlink(mu1) probitlink(mu2) rhobitlink(rho)
#> (Intercept)       -1.997620       -1.055676        2.120919
#> x2                 3.983245        3.011416        0.000000
if (FALSE)  matplot(with(bdata, x2), dmat, lty = 1:4, col = 1:4,
        type = "l", main = "Joint probabilities",
        ylim = 0:1, lwd = 2, ylab = "Probability")
legend(x = 0.25, y = 0.9, lty = 1:4, col = 1:4, lwd = 2,
       legend = c("1 = (y1=0, y2=0)", "2 = (y1=0, y2=1)",
                  "3 = (y1=1, y2=0)", "4 = (y1=1, y2=1)"))  # \dontrun{}
#> Error in (function (s, units = "user", cex = NULL, font = NULL, vfont = NULL,     ...) {    if (!is.null(vfont))         vfont <- c(typeface = pmatch(vfont[1L], Hershey$typeface),             fontindex = pmatch(vfont[2L], Hershey$fontindex))    .External.graphics(C_strWidth, as.graphicsAnnot(s), pmatch(units,         c("user", "figure", "inches")), cex, font, vfont, ...)})(dots[[1L]][[1L]], cex = dots[[2L]][[1L]], font = dots[[3L]][[1L]],     units = "user"): plot.new has not been called yet