geeCRT (version 0.0.1)

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)

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

Value

y a vector of simulated binary outcomes for n clusters.

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

# NOT RUN {
#####################################################################
# 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)




# }