geeCRT (version 1.1.3)

simbinPROBIT: Generating Correlated Binary Data using the Multivariate Probit Method.

Description

simbinPROBIT generates correlated binary data using the multivariate Probit method (Emrich and Piedmonte, 1991). It simulates a vector of binary outcomes according the specified marginal mean vector and correlation structure. Constraints and compatibility between the marginal mean and correlation matrix are checked.

Usage

simbinPROBIT(mu, Sigma, n = 1)

Value

y a vector of simulated binary outcomes for n clusters.

Arguments

mu

a mean vector when n = 1 or is NULL, otherwise a list of mean vectors for the n clusters

Sigma

a correlation matrix when n = 1 or is NULL, otherwise a list of correlation matrices for the n clusters

n

number of clusters. The default is 1

Author

Hengshi Yu <[e-mail]>, Fan Li <[e-mail]>, Paul Rathouz <[e-mail]>, Elizabeth L. Turner <[e-mail]>, John Preisser <[e-mail]>

References

Emrich, L. J., & Piedmonte, M. R. (1991). A method for generating high-dimensional multivariate binary variates. The American Statistician, 45(4), 302-304.

Preisser, J. S., Qaqish, B. F. (2014). A comparison of methods for simulating correlated binary variables with specified marginal means and correlations. Journal of Statistical Computation and Simulation, 84(11), 2441-2452.

Examples

#####################################################################
# Simulate 2 clusters, 3 periods and cluster-period size of 5 #######
#####################################################################

t <- 3
n <- 2
m <- 5

# means of cluster 1
u_c1 <- c(0.4, 0.3, 0.2)
u1 <- rep(u_c1, c(rep(m, t)))
# means of cluster 2
u_c2 <- c(0.35, 0.25, 0.2)
u2 <- rep(u_c2, c(rep(m, t)))

# List of mean vectors
mu <- list()
mu[[1]] <- u1
mu[[2]] <- u2
# List of correlation matrices

## correlation parameters
alpha0 <- 0.03
alpha1 <- 0.015
rho <- 0.8

## (1) exchangeable
Sigma <- list()
Sigma[[1]] <- diag(m * t) * (1 - alpha1) + matrix(alpha1, m * t, m * t)
Sigma[[2]] <- diag(m * t) * (1 - alpha1) + matrix(alpha1, m * t, m * t)
y_exc <- simbinPROBIT(mu = mu, Sigma = Sigma, n = n)

## (2) nested exchangeable
Sigma <- list()
cor_matrix <- matrix(alpha1, m * t, m * t)
loc1 <- 0
loc2 <- 0
for (t in 1:t) {
    loc1 <- loc2 + 1
    loc2 <- loc1 + m - 1
    for (i in loc1:loc2) {
        for (j in loc1:loc2) {
            if (i != j) {
                cor_matrix[i, j] <- alpha0
            } else {
                cor_matrix[i, j] <- 1
            }
        }
    }
}

Sigma[[1]] <- cor_matrix
Sigma[[2]] <- cor_matrix
y_nex <- simbinPROBIT(mu = mu, Sigma = Sigma, n = n)

## (3) exponential decay

Sigma <- list()

### function to find the period of the ith index
region_ij <- function(points, i) {
    diff <- i - points
    for (h in 1:(length(diff) - 1)) {
        if (diff[h] > 0 & diff[h + 1] <= 0) {
            find <- h
        }
    }
    return(find)
}

cor_matrix <- matrix(0, m * t, m * t)
useage_m <- cumsum(m * t)
useage_m <- c(0, useage_m)

for (i in 1:(m * t)) {
    i_reg <- region_ij(useage_m, i)
    for (j in 1:(m * t)) {
        j_reg <- region_ij(useage_m, j)
        if (i_reg == j_reg & i != j) {
            cor_matrix[i, j] <- alpha0
        } else if (i == j) {
            cor_matrix[i, j] <- 1
        } else if (i_reg != j_reg) {
            cor_matrix[i, j] <- alpha0 * (rho^(abs(i_reg - j_reg)))
        }
    }
}
Sigma[[1]] <- cor_matrix
Sigma[[2]] <- cor_matrix
y_ed <- simbinPROBIT(mu = mu, Sigma = Sigma, n = n)