Estimating proportions & 95% CIs with clustered data from very few clusters

Question

What are the best method(s) to estimate a proportion and its associated 95% confidence interval (ideally with an option for an exact method to avoid values +/- 100 or 0) when the data are clustered, but the number of clusters are very few?

Motivating example and data

Assume you are collecting data on whether clinicians carry out a certain procedure correctly or incorrectly in a small number of hospitals, and there is clustering at the hospital level (but no where else in this example). Assume you collect data in three hospitals by simple random sampling of clinicians, and you collect 50 observations in each hospital. Assume each observation is from a different clinician, and there is a small clustering effect of their performance within each hospital.

The R code below creates a data set representing this situation, based on a binary outcome of clinician performance in three different clusters (hospitals), and where the expected proportions within each cluster differ by just 5 percentage points across all three.

n <- 50    
p1 <- rbinom(n, 1, 0.50)
p2 <- rbinom(n, 1, 0.525)
p3 <- rbinom(n, 1, 0.55)
d <- data.frame(y = c(p1, p2), cluster = rep(c("a", "b"), each = n))

Estimating proportions and generating 95% confidence intervals for this data whilst ignoring the clustering may produce biased results, but I'm not clear how suitable it would be to use multilevel methods or survey methods like robust standard errors given the very small number of clusters available? With a 'relatively' small clustering effect might it be acceptable to simply ignore the clustering?

I am interested in methodological approaches, but software specific options in R/Stata/SAS/SPSS would be really welcomed too.

Answer 1

Thomas Lumley's 'survey' package for R is an excellent tool for specifying survey designs in R. Here is a link to Lumley describing how to specify a sample design. His examples in this link are with schools; your question is about hospitals -- very similar. The author of this package has generously produced many many tutorials (+ a book), which should be easy enough for you to find. Just Google 'lumley survey.'

However, as noted in the comments, there is some information lacking here for us to give very specific advice. Really important is how these three hospitals were chosen and if the hospitals have the same number of clinicians. If you select three hospitals and sample 50 clinicians regardless of the size of the hospital, then I think you're question isn't really about cluster effects -- this is really about stratification. In this case, what you will need to do, for the calculation of means or other descriptives, is to specify a strata variable, and you will need to specify stratum sizes (as larger hospitals should be given more weight if you sampled 50 from each hospital regardless of size). Also, if you sampled a very large fraction of clinicians in each hospital (say greater than 5%), then you should also make sure that you account for a finite population correction -- the fpc parameter in svydesign() function. (fps and stratum sizes are often equal, so you only need to specify one if you don't have your own probability or post-stratification weights.)

Best case scenario -- these three hospitals are all hospitals in a specific city, or you sampled hospitals probability proportionate to size. In which case, move on, your sample is fine -- just account for stratification. If hospitals were selected in pretty much in other way, however, you need to do some reading about first-stage stratification, as without more information about your design, we won't be able to help. Or, you unambiguously caveat that hospitals were sampled conveniently.

But if your intention is regress something-on-something with this data set, then I would simply include hospital as a covariate -- or depending on your question, consider specifying hospital as a random effect. There are some arguments for including survey weights in a regression like this, but anecdotally I would say this is not the norm.

*As a final aside, if you do opt to use the 'survey' package, don't forget to use the "~" where you have to. For good reasons, this is required when you are referring to specific variables in your dataset, but in my experience when people run into problems with this package, it's because they forgot to include a "~" somewhere. Most other packages don't require this, so it throws some people off.

Answer 2

One possible approach is a hierarchical model that makes use of any known prior information you may have with the number of correctly performed procedures distributed following a binomial distribution given the hospital specific log-odds (in reality, one would presumably try to adjust for e.g. the difficulty of the procedure, experience of the physician etc. assuming that such data are available) $$Y_i | \theta_i \sim \text{Bin}(n_i, \pi_i:= \frac{e^{\theta_i}}{1+e^{\theta_i}}).$$ Then one might assume that that the log-odds vary across the hospitals following some distribution that is based on the prior knowledge. Let's say that we know nothing about the hospitals that would make us assume one would be different from the other, then $$\theta_i \sim \text{N}(\mu, \sigma^2)$$ may be a reasonable assumption. We may have an approximate idea about $\mu$ from epidemiological data or prior studies on the topic (but may want to down-weight this prior information to understate this knowledge to reflect that we are transitioning to a new setting, one approach for doing that is power-prior approach of Ibrahim and Chen), similarly we may have some prior knowledge or opinion on how much the log-odds could different between hospitals (say, we believe odds of more than 5:1 are totally implausible, then a prior with very little support above $\log(5)$ would be an option, e.g. a uniform on $[0,\log(5)]$ or half-normal with a high quantile at the point $\log(5)$ depending on your beliefs about between hospital variation).

Realistically, however, we will know something about the hospitals, e.g. that some hospitals conduct more complex operations/procedures that are more prone towards complications. It is a problem when we do not know for each procedure, but we can try things like meta-regression (or other approaches that can account for ecological bias) to adjust for these things, if we at least information on a hospital level. Taking into account prior information on the complication rates of such procedures would again be crucial, since there are so few hospitals that estimating such things from meta-regression without priors will be hopeless.

From such a model you can get an estimate of all the possible quantities of interest including prediction or credible intervals.