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, 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.

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.