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.