I’ve been working on updates for the simstudy package. In the past few weeks, a couple of folks independently reached out to me about generating correlated binary data. One user was not impressed by the copula algorithm that is already implemented. I’ve added an option to use an algorithm developed by Emrich and Piedmonte in 1991, and will be incorporating that option soon in the functions genCorGen and addCorGen. I’ll write about that change some point soon.
A second researcher was trying to generate data using parameters that could be recovered using GEE model estimation. I’ve always done this by using an underlying mixed effects model, but of course, the marginal model parameter estimates might be quite different from the conditional parameters. (I’ve written about this a number of times, most recently here.) As a result, the model and the data generation process don’t match, which may not be such a big deal, but is not so helpful when trying to illuminate the models.
One simple solution is using a beta-binomial mixture data generating process. The beta distribution is a continuous probability distribution that is defined on the interval from 0 to 1, so it is not too unreasonable as model for probabilities. If we assume that cluster-level probabilities have a beta distribution, and that within each cluster the individual outcomes have a binomial distribution defined by the cluster-specific probability, we will get the data generation process we are looking for.
Generating the clustered data
In these examples, I am using 500 clusters, each with cluster size of 40 individuals. There is a cluster-level covariate x that takes on integer values between 1 and 3. The beta distribution is typically defined using two shape parameters usually referenced as \(\alpha\) and \(\beta\), where \(E(Y) = \alpha / (\alpha + \beta)\), and \(Var(Y) = (\alpha\beta)/[(\alpha + \beta)^2(\alpha + \beta + 1)]\). In simstudy, the distribution is specified using the mean probability (\(p_m\)) and a precision parameter (\(\phi_\beta __ 0\)) (that is specified using the variance argument). Under this specification, \(Var(Y) = p_m(1 - p_m)/(1 + \phi_\beta)\). Precision is inversely related to variability: lower precision is higher variability.
In this simple simulation, the cluster probabilities are a function of the cluster-level covariate and precision parameter \(\phi_\beta\). Specifically
\[logodds(p_{clust}) = -2.0 + 0.65x.\]
The binomial variable of interest \(b\) is a function of \(p_{clust}\) only, and represents a count of individuals in the cluster with a “success”:
library(simstudy)
set.seed(87387)
phi.beta