| Binom2.rho | R Documentation |
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("argument '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)
n |
number of observations.
Same as in |
mu1, mu2 |
The marginal probabilities.
Only |
rho |
The correlation parameter.
Must be numeric and lie between |
exchangeable |
Logical. If |
twoCols |
Logical.
If |
colnames |
The |
ErrorCheck |
Logical. Do some error checking of the input parameters? |
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.
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.
T. W. Yee
binom2.rho.
(myrho <- rhobitlink(2, inverse = TRUE)) # Example 1
nn <- 2000
ymat <- rbinom2.rho(nn, mu1 = 0.8, rho = myrho, exch = TRUE)
(mytab <- table(ymat[, 1], ymat[, 2], dnn = c("Y1", "Y2")))
fit <- vglm(ymat ~ 1, binom2.rho(exch = TRUE))
coef(fit, matrix = TRUE)
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)
## Not run: 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)"))
## End(Not run)
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.