A comparative study of R functions for clustered data analysis

Clustered or correlated outcome data is common in medical research studies, such as the analysis of national or international disease registries, or cluster-randomized trials, where groups of trial participants, instead of each trial participant, are randomized to interventions. Within-group correlation in studies with clustered data requires the use of specific statistical methods, such as generalized estimating equations and mixed-effects models, to account for this correlation and support unbiased statistical inference.

We compare different approaches to estimating generalized estimating equations and mixed effects models for a continuous outcome in R through a simulation study and a data example. The methods are implemented through four popular functions of the statistical software R, "geese", "gls", "lme", and "lmer". In the simulation study, we compare the mean squared error of estimating all the model parameters and compare the coverage proportion of the 95% confidence intervals. In the data analysis, we compare estimation of the intervention effect and the intra-class correlation.

In the simulation study, the function "lme" takes the least computation time. There is no difference in the mean squared error of the four functions. The "lmer" function provides better coverage of the fixed effects when the number of clusters is small as 10. The function "gls" produces close to nominal scale confidence intervals of the intra-class correlation. In the data analysis and the "gls" function yields a positive estimate of the intra-class correlation while the "geese" function gives a negative estimate. Neither of the confidence intervals contains the value zero.

The "gls" function efficiently produces an estimate of the intra-class correlation with a confidence interval. When the within-group correlation is as high as 0.5, the confidence interval is not always obtainable.