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

Question

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)

Error in eval(predvars, data, env) : object 'Y.likelihood' not found

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
```

Answer 1

“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)
f <- lrm(…, data=mydata, x=TRUE, y=TRUE)
anova(f, test=‘LR’)

I wish I had implemented profile confidence intervals for odds ratios but it’s not there. Other packages will do this.