Binomial Regression "logit" vs "cloglog"

Question

I am using a binomial regression with a categorical factor with 9 levels (named 'treat.group') and sample sizes in each group from 1-8. 1 treatment group has all positive cases (i.e., 1's) - and this creates a estimation problem with the "standard" glm() function in R caused by "perfect separation" for that treatment level. So, I am using bayesglm from the arm package.

My question is that the default identity link is “logic“ but i have read that "cloglog" or (Complementary Log-Log) is frequently used when the probability of an event is very small or very large. Thus since my model exhibits perfect separation in 1 treatment group the probability of the event is very large and I should use "cloglog". Using cloglog gives me a significant result for the treatment group with perfect separation while "logit" does not.

Am I justified in using "cloglog" or is there a way to look at my results and be certain what link is best?

f1<- bayesglm(response~ treat.group,family=binomial(link="logit"), data=df)

f2 <- bayesglm(response~ treat.group,family=binomial(link="cloglog"), data=df)

(Data frame below)

{ structure(list(response = c(0L, 1L, 1L, 0L, 0L, 1L, 1L, 1L, 0L, 1L, 1L, 0L, 0L, 1L, 0L, 1L, 1L, 0L, 1L, 1L, 0L, 0L, 0L, 0L, 1L, 1L, 0L, 0L, 0L, 1L, 1L, 0L, 1L, 0L, 1L, 0L, 1L, 1L, 1L, 1L), treat.group = c("pctc", "phth", "phth", "phtl", "pltc", "pcth", "pltl", "phtc", "pctl", "phtc", "pcth", "pctl", "pctc", "phtl", "plth", "pltc", "phtc", "pcth", "phtl", "plth", "pctl", "pltc", "phtl", "pctc", "pcth", "pltc", "phtc", "phtl", "phtc", "pctl", "pctc", "pcth", "phth", "pctc", "phtl", "pcth", "phth", "phtc", "pcth", "phth")), .Names = c("response", "treat.group"), row.names = c(NA, -40L), class = c("tbl_df", "tbl", "data.frame"), na.action = structure(c(1L, 4L, 5L, 7L, 15L, 21L, 23L, 24L, 27L, 29L, 33L, 37L, 39L, 48L, 50L, 53L, 54L, 55L, 56L, 57L, 58L, 59L, 60L, 62L, 63L, 65L, 66L, 67L, 68L, 70L, 71L, 72L, 73L), .Names = c("1", "4", "5", "7", "15", "21", "23", "24", "27", "29", "33", "37", "39", "48", "50", "53", "54", "55", "56", "57", "58", "59", "60", "62", "63", "65", "66", "67", "68", "70", "71", "72", "73"), class = "omit")) }

Answer 1

If you fit two models, one with logit and one with cloglog, you should report the results of both, and also carry out some type of model comparison technique and report the results of that.

As for the models, this is a great situation in which to use Bayesian multilevel models (see this paper [PDF] by Gelman). We can pool information among groups to inform the estimates for groups with smaller sample sizes and for groups with seemingly extreme outcomes (as in phth and pltl).

Fitting multilevel models with the brms package is fairly straightforward. The two models here, one with a logit link and one with a cloglog link, would be fit as follows:

fit1 <- brm(
    response ~ (1 | treat.group),
    family = bernoulli,  # logit link is the default
    data = df,
    iter = 2e4,
    warmup = 2e3,
    control = list(adapt_delta = 0.999)
)

fit2 <- brm(
    response ~ (1 | treat.group),
    family = bernoulli(link = "cloglog"),
    data = df,
    iter = 2e4,
    warmup = 2e3,
    control = list(adapt_delta = 0.999)
)

These will each generate 72000 samples from the models, which can be accessed with as.data.frame(fit). After applying the inverse logit and inverse cloglog to the samples, we get the following estimates for the desired probabilities:

group    prob using logit    sd       89% probability interval
--------------------------------------------------------------
pctc     0.405               0.169    0.12 to 0.66    
pcth     0.562               0.135    0.34 to 0.78    
pctl     0.438               0.170    0.15 to 0.7    
phtc     0.604               0.141    0.39 to 0.84    
phth     0.729               0.162    0.52 to 1    
phtl     0.529               0.142    0.3 to 0.76    
pltc     0.532               0.157    0.28 to 0.79    
plth     0.539               0.180    0.24 to 0.84    
pltl     0.623               0.193    0.39 to 1


group    prob using cloglog    sd       89% probability interval
----------------------------------------------------------------
pctc     0.395                 0.162    0.12 to 0.64    
pcth     0.545                 0.136    0.32 to 0.76    
pctl     0.422                 0.164    0.14 to 0.67    
phtc     0.589                 0.144    0.36 to 0.83    
phth     0.760                 0.182    0.52 to 1    
phtl     0.511                 0.142    0.28 to 0.74    
pltc     0.512                 0.158    0.25 to 0.77    
plth     0.517                 0.181    0.21 to 0.81    
pltl     0.627                 0.213    0.37 to 1

Overall, the model using the cloglog link induced less shrinkage in the probabilities. It's worth noting that the 89% probability intervals for phth are the same in the two models.

The brms package supports model comparison using PSIS-LOO. Using the command

loo(fit1, fit2, reloo = TRUE)

takes a while but eventually returns the following result:

            LOOIC   SE
fit1        59.19 2.95
fit2        58.59 2.91
fit1 - fit2  0.60 0.46

Here the model with lower LOOIC is estimated to have better out-of-sample prediction accuracy. The model with the cloglog link has marginally lower LOOIC, and the standard error on the difference, 0.46, is small, so we would conclude that the model using the cloglog link is marginally better than the one using the logit link.