Bivariate Odds Ratio Model
Density and random generation for a bivariate binary regression model using an odds ratio as the measure of dependency.
rbinom2.or(n, mu1,
mu2 = if (exchangeable) mu1 else stop("argument 'mu2' not specified"),
oratio = 1, exchangeable = FALSE, tol = 0.001, twoCols = TRUE,
colnames = if (twoCols) c("y1","y2") else c("00", "01", "10", "11"),
ErrorCheck = TRUE)
dbinom2.or(mu1, mu2 = if (exchangeable) mu1 else stop("'mu2' not specified"),
oratio = 1, exchangeable = FALSE, tol = 0.001,
colnames = c("00", "01", "10", "11"), ErrorCheck = TRUE)n |
number of observations.
Same as in |
mu1, mu2 |
The marginal probabilities.
Only |
oratio |
Odds ratio. Must be numeric and positive. The default value of unity means the responses are statistically independent. |
exchangeable |
Logical. If |
twoCols |
Logical.
If |
colnames |
The |
tol |
Tolerance for testing independence. Should be some small positive numerical value. |
ErrorCheck |
Logical. Do some error checking of the input parameters? |
The function rbinom2.or generates data coming from a bivariate
binary response model.
The data might be fitted with the VGAM family function
binom2.or.
The function dbinom2.or does not really compute the density
(because that does not make sense here) but rather returns the
four joint probabilities.
The function rbinom2.or returns
either a 2 or 4 column matrix of 1s and 0s, depending on the argument
twoCols.
The function dbinom2.or returns
a 4 column matrix of joint probabilities; each row adds up to unity.
T. W. Yee
nn <- 1000 # Example 1
ymat <- rbinom2.or(nn, mu1 = logitlink(1, inv = TRUE), oratio = exp(2), exch = TRUE)
(mytab <- table(ymat[, 1], ymat[, 2], dnn = c("Y1", "Y2")))
(myor <- mytab["0","0"] * mytab["1","1"] / (mytab["1","0"] * mytab["0","1"]))
fit <- vglm(ymat ~ 1, binom2.or(exch = TRUE))
coef(fit, matrix = TRUE)
bdata <- data.frame(x2 = sort(runif(nn))) # Example 2
bdata <- transform(bdata, mu1 = logitlink(-2 + 4 * x2, inverse = TRUE),
mu2 = logitlink(-1 + 3 * x2, inverse = TRUE))
dmat <- with(bdata, dbinom2.or(mu1 = mu1, mu2 = mu2, oratio = exp(2)))
ymat <- with(bdata, rbinom2.or(n = nn, mu1 = mu1, mu2 = mu2, oratio = exp(2)))
fit2 <- vglm(ymat ~ x2, binom2.or, data = bdata)
coef(fit2, matrix = TRUE)
## Not run:
matplot(with(bdata, x2), dmat, lty = 1:4, col = 1:4, type = "l",
main = "Joint probabilities", ylim = 0:1,
ylab = "Probabilities", xlab = "x2", las = 1)
legend("top", lty = 1:4, col = 1:4,
legend = c("1 = (y1=0, y2=0)", "2 = (y1=0, y2=1)",
"3 = (y1=1, y2=0)", "4 = (y1=1, y2=1)"))
## End(Not run)Please choose more modern alternatives, such as Google Chrome or Mozilla Firefox.