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

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)
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.

• It makes perfect sense and you have already identified why it happens. What you do now depends on the goals of your study as others have answered. – mdewey Feb 9 '17 at 13:53

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)

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

• One might also add that there is the possibility to compare the model with a null model using a likelihood ratio test. – mdewey Feb 9 '17 at 13:51
• Kinda offtopic: I have no experience working with arm::bayesglm but how it selects prior distribution? Because only for treatment C we can get stable coefficient, especially if we set it as Intercept (around -1.05 according to JAGS and STAN with non-informative prior), while for treatment A and D we're free to go into plus and minus infinity respectively, they are basically defined only by prior. – Andrey Kolyadin Feb 10 '17 at 14:05
• I will add some comment (in the answer above) about how it chooses its prior. As a starter, read ?bayesglm – kjetil b halvorsen Feb 10 '17 at 14:17

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.