Binomial GLM - non-significant difference between 100% opposite groups of observations

Question

What follows is a basic question concerning Binomial GLM's.

Suppose we have a set of observations where a binary response was measured in three different treatments, A, C and D -

treatment = c("A", "A", "A", "A", "A", "A", "A", "A", "A", "A", "A", "A", "A", "A", "A", "A", "A", "A", "A", "A", "C", "C", "C", "C", "C", "C", "C", "C", "C", "C", "C", "C", "C", "C", "C", "D", "D", "D", "D", "D", "D")
response = c(TRUE, TRUE, TRUE, TRUE, TRUE, TRUE, TRUE, TRUE, TRUE, TRUE, TRUE, TRUE, TRUE, TRUE, TRUE, TRUE, TRUE, TRUE, TRUE, TRUE, FALSE, FALSE, FALSE, FALSE, FALSE, FALSE, FALSE, TRUE, TRUE, FALSE, TRUE, FALSE, FALSE, TRUE, FALSE, FALSE, FALSE, FALSE, FALSE, FALSE, FALSE)
dat = data.frame(treatment, response, stringsAsFactors = FALSE)
head(dat)
  treatment response
1         A     TRUE
2         A     TRUE
3         A     TRUE
4         A     TRUE
5         A     TRUE
6         A     TRUE

Incidentally, the response was TRUE in 100% of the cases in A, 27% of the cases in C and 0% of the cases in D -

tapply(dat$response, dat$treatment, mean)
        A         C         D 
1.0000000 0.2666667 0.0000000 

According to a Binomial GLM, however, the latter contrasting difference A is not significant -

fit = glm(response ~ treatment, family= "binomial", dat)
summary(fit)

Call:
glm(formula = response ~ treatment, family = "binomial", data = dat)

Deviance Residuals: 
     Min        1Q    Median        3Q       Max  
-0.78760  -0.78760   0.00005   0.00005   1.62589  

Coefficients:
            Estimate Std. Error z value Pr(>|z|)
(Intercept)    20.57    3964.63   0.005    0.996
treatmentC    -21.58    3964.63  -0.005    0.996
treatmentD    -41.13    8253.04  -0.005    0.996

(Dispersion parameter for binomial family taken to be 1)

    Null deviance: 55.637  on 40  degrees of freedom
Residual deviance: 17.397  on 38  degrees of freedom
AIC: 23.397

Number of Fisher Scoring iterations: 19

Does this result make sense? Is there a more appropriate test for cases when outcome of a given treatment completely fall into one or the other response type?

Will appreciate any feedback on this.

Answer 1

The diagnosis in @Andrey Kolyadin's answer is basically correct, so here is some other ideas as to what you can do. Following up from the code in the question, noting that the absurdly large standard errors of the fitted coefficients is an indication that the quadratic approximation of the log likelihood function used in computing those standard errors is very bad. Search this site for Hauck-Donner-effect. As a replacement for those standard errors we can use likelihood profiling, which computes (approximate) confidence intervals directly from the log likelihood function:

library(MASS)
 confint(fit)
Waiting for profiling to be done...
               2.5 %   97.5 %
(Intercept) -210.112       NA
treatmentC        NA 274.9591
treatmentD        NA 601.8675
There were 30 warnings (use warnings() to see them)

(the warnings is about fitted probabilities being numerically zero or one). Note that the NA (NotAvailable) given as lower confidence limits indicates that the procedure was unable to find a lower limit, reasonable because the lower limits are $-\infty$ ! This is because the true maximum likelihood estimators in those two cases are just $-\infty$. So extremely wide and not very useful confidence intervals, consistent with the immense standard errors.

As for diagnosing such separation problems, there is the useful R package safeBinaryRegression. See my answer here Why does logistic regression become unstable when classes are well-separated? for an example of its use.

One way of "solving" this problem, obtaining somewhat more reasonable point estimators and intervals (or se) is using some form of regularization, infusing into the estimation procedure some prior information that the coefficients are really finite. One way of doing that is via a bayes method, which we will illustrate using bayesglm from the arm package.

library(arm)
Loading required package: Matrix
Loading required package: lme4

arm (Version 1.9-3, built: 2016-11-21)

Working directory is C:/Users/halvorsk/Documents

 fit.b <- bayesglm(response ~ treatment, family= "binomial", dat)
 summary(fit.b)

Call:
bayesglm(formula = response ~ treatment, family = "binomial", 
    data = dat)

Deviance Residuals: 
    Min       1Q   Median       3Q      Max  
-0.8236  -0.8236   0.2396   0.2396   1.5787  

Coefficients:
            Estimate Std. Error z value Pr(>|z|)    
(Intercept)    3.536      1.235   2.863 0.004193 ** 
treatmentC    -4.443      1.333  -3.333 0.000859 ***
treatmentD    -6.733      2.229  -3.021 0.002522 ** 
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

(Dispersion parameter for binomial family taken to be 1)

    Null deviance: 55.637  on 40  degrees of freedom
Residual deviance: 19.060  on 38  degrees of freedom
AIC: 25.06

Number of Fisher Scoring iterations: 23

which gives much more reasonable estimates and standard errors.

EDIT   

The OP asks about the prior distributions used by this bayesian method. That is a very reasonable question! First, the R help pages obtained by typing ?bayesglm is useful here. The prior is a "weakly informative prior", t-distributions (defaul is 1 df, that is, the Cauchy distribution) with a scale parameter (known) set to 10 for the intercept and somewhat smaller for the other parameters (This is done after scaling predictor variables). For more on this see Gelman's paper http://www.stat.columbia.edu/~gelman/research/published/priors11.pdf

Answer 2

The problem is in your treatments 'A' and 'D', they don't have any variability at all. So you model is unable to estimate errors. (Actual problem lies in calculation of confidence intervals above 0).

Changing rates so there is at least one subject with different result will give us more adequate results (I use a bit different approach to binomial GLM more akin to contingency tables, but coefficients are the same):

df <- data.frame(treatment = c("A", "C", "D"),
                 rate = c(19/20, 4/15, 1/6),
                 n = c(20, 15, 6))

fit <- glm(rate ~ treatment, weights = n, family= "binomial", df)
summary(fit)

Coefficients:
            Estimate Std. Error z value Pr(>|z|)    
(Intercept)    2.944      1.026   2.870 0.004106 ** 
treatmentC    -3.956      1.180  -3.351 0.000805 ***
treatmentD    -4.554      1.501  -3.034 0.002412 ** 
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

(Dispersion parameter for binomial family taken to be 1)

    Null deviance: 2.4892e+01  on 2  degrees of freedom
Residual deviance: 2.2204e-16  on 0  degrees of freedom
AIC: 12.732

You can read more about zeros in contingency tables (which is basicly your problem) in discussion to this question.

If you interested only in p values I would say that use of pair wise Fisher's exact test with adjustment for multiple comparisons can be applied to your problem.

Answer 3

Another method to deal with complete separation in logistic regression is penalization. To demonstrate, I'll use Firth's logistic regression implemented in the logistf package.

library("logistf")

model <- logistf(
  response ~ treatment,
  firth = TRUE
)

summary(model)
#> logistf(formula = response ~ treatment, firth = TRUE)
#> 
#> Model fitted by Penalized ML
#> Coefficients:
#>                  coef se(coef) lower 0.95 upper 0.95    Chisq            p
#> (Intercept)  3.713572 1.431356   1.733332   8.559073 24.38651 7.881859e-07
#> treatmentC  -4.651842 1.535565  -9.575829  -2.342167 22.88854 1.716719e-06
#> treatmentD  -6.278521 2.050031 -12.052663  -3.138511 23.16264 1.488616e-06
#>             method
#> (Intercept)      2
#> treatmentC       2
#> treatmentD       2
#> 
#> Method: 1-Wald, 2-Profile penalized log-likelihood, 3-None
#> 
#> Likelihood ratio test=32.53828 on 2 df, p=8.598101e-08, n=41
#> Wald test = 12.63304 on 2 df, p = 0.001806222

The estimates and confidence intervals agree reasonably well with the Bayesian analysis described in kjetilbhalvorsen's answer. See also How to deal with perfect separation in logistic regression?.