# What is the difference between using logistf and brglm2 when dealing with complete separation in a logistic regression?

I am trying to looking at how the three factors A (5 levels, a-e), B (2 levels, a and b) and C (2 levels, a and b) affect the likelihood of event Y (1 = occured, 0 = did not occur). I initially ran a normal logistic regression model as follows:

type3.Y.Full <- list(A = contr.sum, B = contr.sum, C = contr.sum)
model.Y<-glm(Y.likelihood ~ A*B*C, family = binomial (link = "logit"), data = Y.1, contrasts = type3.Y.Full)
summary(model.Y)
library(car)
Anova(model.Y, type = 3)
library("detectseparation")
update(model.Y, method = "detect_separation")#Complete separation detected


Examination of the summary output (showing abnormally large standard errors of coefficients and all p-values very close to 1), examination of the no of 1's and 0's for each combination of the three factors levels, as well as running a test to detect separation using the code above, revealed to me that there were several instances of complete separation in my data.

As a result of this, my research online had led me to believe that I need to run a bias reduced logistic regression to deal with the separation. I have found two methods to do this as follows:

#########Method 1 - brglm2
library(brglm2)
model.Y.brglm<-glm(Y.likelihood ~ A*B*C, family = binomial (link = "logit"), data = Y.1, contrasts = type3.Y.Full, method = brglmFit)
summary(model.Y.brglm)

Call:
glm(formula = Y.likelihood ~ A * B * C,
family = binomial(link = "logit"), data = Y.1, method = brglmFit,
contrasts = type3.Y.Full)

Deviance Residuals:
Min       1Q   Median       3Q      Max
-2.2534  -1.0731   0.5415   0.9738   1.6651

Coefficients:
Estimate Std. Error z value Pr(>|z|)
(Intercept)                         0.54796    0.16566   3.308 0.000940 ***
A1                                  0.64720    0.34125   1.897 0.057886 .
A2                                  0.70109    0.29486   2.378 0.017420 *
A3                                  0.75656    0.35739   2.117 0.034267 *
A4                                 -0.33962    0.21995  -1.544 0.122572
B1                                  0.38930    0.16566   2.350 0.018771 *
C1                                  0.62192    0.16566   3.754 0.000174 ***
A1:B1                               0.18165    0.34125   0.532 0.594520
A2:B1                              -0.58628    0.29486  -1.988 0.046771 *
A3:B1                               0.39197    0.35739   1.097 0.272747
A4:B1                              -0.09776    0.21995  -0.444 0.656725
A1:C1                               0.12427    0.34125   0.364 0.715740
A2:C1                               0.11557    0.29486   0.392 0.695086
A3:C1                               0.04053    0.35739   0.113 0.909714
A4:C1                              -0.70748    0.21995  -3.216 0.001298 **
B1:C1                              -0.13459    0.16566  -0.812 0.416535
A1:B1:C1                           -0.53189    0.34125  -1.559 0.119081
A2:B1:C1                           -0.13863    0.29486  -0.470 0.638249
A3:B1:C1                            0.82033    0.35739   2.295 0.021713 *
A4:B1:C1                            0.38356    0.21995   1.744 0.081193 .
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

(Dispersion parameter for binomial family taken to be 1)

Null deviance: 583.11  on 435  degrees of freedom
Residual deviance: 486.65  on 416  degrees of freedom
AIC:  526.65

Type of estimator: AS_mixed (mixed bias-reducing adjusted score equations)
Number of Fisher Scoring iterations: 7

#######Method 2 - logistf
install.packages("logistf")
library(logistf)
model.Y.logistf<-logistf(Y.likelihood ~ A*B*C, family = binomial (link = "logit"), data = Y.1, contrasts = type3.Y.Full)
summary(model.Y.logistf)

logistf(formula = Y.likelihood ~ A * B * C,
data = Y.1, family = binomial(link = "logit"),
contrasts = type3.Y.Full)

Model fitted by Penalized ML
Coefficients:
coef  se(coef) lower 0.95  upper 0.95        Chisq           p method
(Intercept)                           1.845826690 0.8785954  0.3799589  4.07823857 6.486481e+00 0.010869791      2
Ab                                   -0.329479201 1.0212407 -2.7268144  1.53675350 1.079762e-01 0.742460079      2
Ac                                    1.588160514 1.6841752 -1.4441737  6.61356792 1.046109e+00 0.306404983      2
Ad                                   -1.182532473 0.9747733 -3.5288125  0.53518532 1.726111e+00 0.188908279      2
Ae                                   -2.182298927 1.2073293 -4.8873400  0.04489102 3.684039e+00 0.054935611      2
Bb                                    0.191055237 1.2351328 -2.4401842  2.82413887 2.390508e-02 0.877126674      2
Cb                                   -0.159427737 1.1162454 -2.6642665  2.01602591 2.058149e-02 0.885924913      2
Ab:Bb                                 0.749333047 1.5875658 -2.4131725  4.09862430 2.250972e-01 0.635183231      2
Ac:Bb                                -3.125086490 1.9175271 -8.3798273  0.32196136 3.155520e+00 0.075670574      2
Ad:Bb                                -1.272084655 1.3515401 -4.0934860  1.54111671 8.645881e-01 0.352457930      2
Ae:Bb                                 0.145417000 1.5398920 -2.9833643  3.30894785 8.921787e-03 0.924747560      2
Ab:Cb                                -0.769133088 1.3131695 -3.3371442  2.02486335 3.337002e-01 0.563488328      2
Ac:Cb                                -2.536960524 1.8919024 -7.7521946  0.90157352 2.081519e+00 0.149091717      2
Ad:Cb                                -0.167394244 1.2495599 -2.5996562  2.53278844 1.782789e-02 0.893781252      2
Ae:Cb                                -0.602712315 1.5700483 -3.6982496  2.62924915 1.460259e-01 0.702362483      2
Bb:Cb                                -2.665911551 1.4810694 -5.7482418  0.33581980 3.085541e+00 0.078991000      2
Ab:Bb:Cb                              1.573054674 1.9223037 -2.3271087  5.39155135 6.588991e-01 0.416948868      2
Ac:Bb:Cb                              5.408887567 2.1701134  1.4434461 10.92484427 7.119263e+00 0.007626006      2
Ad:Bb:Cb                              3.661783160 1.6533010  0.3567045  7.04155205 4.652260e+00 0.031012684      2
Ae:Bb:Cb                             -0.005935602 2.3742641 -5.7356358  4.45315310 6.251339e-06 0.998005077      2

Method: 1-Wald, 2-Profile penalized log-likelihood, 3-None

Likelihood ratio test=88.07473 on 19 df, p=7.262535e-11, n=436
Wald test = 65.45671 on 19 df, p = 5.145634e-07


I would like to run a type-III ANOVA after each logistic regression. Running it using the brglm model gives the following results however, trying to run it using the logistf method gives me the following error:

library(car)
######Method 1 - brglm2
Anova(model.Y.brglm, type = 3)

Analysis of Deviance Table (Type III tests)

Response: Y.likelihood
LR Chisq Df Pr(>Chisq)
A                                 44.904  4  4.164e-09 ***
B                                  5.810  1   0.015939 *
C                                 20.969  1  4.667e-06 ***
A:B                                6.359  4   0.173862
A:C                               13.550  4   0.008877 **
B:C                               -2.479  1   1.000000
A:B:C                             13.371  4   0.009597 **
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

#####Method 2 - logistf
Anova(model.Y.logistf)



As a result my questions are as follows:

1. What is the difference between the brglm2 and the logistf methods in terms of how they tackle complete separation? Are there circumstances in which you would choose one over the other?
2. How does a type III ANOVA work with each of two methods of dealing with complete separation? Why does it only seem to work with the brglm2 method?

Edited to add: I have tried using the drop1() function with the logistf model in order to obtain type-III ANOVA p-values for each term in the model however I get the following results:

library(logistf)
model.Y.logistf<-logistf(Y.likelihood ~ A*B*C, family = binomial (link = "logit"), data = Y.1, contrasts = type3.Y.Full)
summary(model.Y.logistf)
drop1(model.Y.logistf, test = "Chisq")

ChiSq df P-value
A                                -319.16377  4       1
B                                 -26.64750  1       1
C                                 -11.56925  1       1
A:B                             -1756.86802  4       1
A:C                              -462.52627  4       1
B:C                              -287.00847  1       1
A:B:C                           -1503.59447  4       1
$$$$

• As I explained in response to your previous question about this, Why do I get a negative chi-squared value in my type III ANOVA output for my binomial GLM?, {logistf} has its own anova function. The {car} anova function is not designed to work with Firth's regression. So it's actually a good thing that you get this error. Commented Jun 3 at 14:17
• @dipetkov Would it not be possible then to run either type-II or type-III ANOVA analyses with penalized regressions? Commented Jun 3 at 15:21
• @dipetkov I have tried using the drop1() function in order to get type-III ANOVA p-values with the logistf model however I get a value of 1 for every term in the model (see main question for details). Do you know why this might be? Commented Jun 3 at 16:19
• FWIW, logistf(..., contrasts = type3.Y.Full) doesn't actually apply the sum contrasts to the A,B,C factor variables as logistf ignores the contrasts argument. Instead you can apply them manually like this: contrasts(A) <- contr.sum(5), ... before fitting the model. Commented Jun 5 at 7:50
• Here is how to check the contrasts applied by the model: model.matrix(model.Y.logistf, data = Y.1). Pay attention to the output, it should say the contrasts for each predictor. Commented Jun 5 at 7:50

“Solutions” for complete separation abound but just by using ordinary likelihood ratio $$\chi^2$$ statistics the problem is solved. The only trick remaining is to use profile likelihood-based confidence intervals instead of ones based on standard errors.
Note that you are trying to compute $$\chi^2$$ for hypotheses that don’t necessarily make any sense by getting tests for lower order terms adjusted for higher order terms. Instead of that consider using the R rms package which gives you all the meaningful LR tests:
require(rms)
`