#load package
library(rms)

#make df
x <- rpois(1000, 0.5)
y <- sample(c(1:100),1000, replace=TRUE)
z <- factor(sample(c(1:25), 1000, replace=TRUE))
df <- data.frame(y,x,z)

#set datadist
dd <- datadist(df)
options(datadist='dd')

#works when cluster variable isn't included as fixed effect in regression
reg  <- Glm(x ~ y, df, x=TRUE, y=TRUE, family=poisson())
reg_clus <- bootcov(reg, df$z)
summary(reg_clus)

#doesn't work when cluster variable included as fixed effect in regression
reg2  <- Glm(x ~ y + z, df, x=TRUE, y=TRUE, family=poisson())
reg_clus2 <- bootcov(reg2, df$z)
summary(reg_clus2)