  This paper presents a method for fitting a copula-driven generalized linear
mixed models. For added flexibility, the skew-normal copula is adopted for
fitting. The correlation matrix of the skew-normal copula is used to capture
the dependence structure within units, while the fixed and random effects
coefficients are estimated through the mean of the copula. For estimation, a
Monte Carlo expectation-maximization algorithm is developed. Simulations are
shown alongside a real data example from the Framingham Heart Study.
