Is there a specific reason why R GLM does not return warnings, while anova(glm) does?

Question

It is a replica of my question at https://stackoverflow.com/questions/46895017/is-there-a-specific-reason-why-r-glm-does-not-return-warnings-while-anovaglm

I was advised to ask it here instead.

It is not about glm.fit: fitted probabilities numerically 0 or 1 occurred; it is about why anova warns while glm doesn't.

I have wanted to see contrasts inside a specified model:

is_service ~ action_count * document_entropy

The full dataset is loaded in the code. The data are somewhat pathological, as you could see below.

Overall the data are these:

> str(dat)
'data.frame':   6432 obs. of  3 variables:
 $ action_count    : num  0.0759 0.1505 0.1435 0.1535 0.2067 ...
 $ document_entropy: num  -0.667 -0.667 -0.667 -0.667 -0.667 ...
 $ is_service      : int  0 0 0 0 0 0 0 0 0 0 ...

The target column has this binomial distribution:

> table(dat$is_service)

   0    1 
6291  141

Input columns are z-normalized and distributed as follows:

> cor.test(dat$action_count, dat$document_entropy)

    Pearson's product-moment correlation

data:  dat$action_count and dat$document_entropy
t = 20.477, df = 6430, p-value < 2.2e-16
alternative hypothesis: true correlation is not equal to 0
95 percent confidence interval:
 0.2243432 0.2702313
sample estimates:
     cor 
0.247426

It is interesting to see that when I fit this model (1st part of the code) the procedure ends without a warnings.

However, when I run contrasts with the stats::anova (2nd part of code) it does return warnings.

Question: Why is that happening, and which level is more alarming: single model or the anova analysis of it?

Is this library specific question or statistics specific?

Output result:

> summary(
+      glm(formula = is_service ~ action_count * document_entropy
+          , family = binomial(link = 'logit'),
+          data = dat
+      )
+ )

Call:
glm(formula = is_service ~ action_count * document_entropy, family = binomial(link = "logit"), 
    data = dat)

Deviance Residuals: 
    Min       1Q   Median       3Q      Max  
-4.0972  -0.1808  -0.1638  -0.1477   3.1019  

Coefficients:
                              Estimate Std. Error z value Pr(>|z|)    
(Intercept)                   -4.06877    0.10143 -40.116  < 2e-16 ***
action_count                   6.04157    0.57882  10.438  < 2e-16 ***
document_entropy               0.29721    0.06866   4.329  1.5e-05 ***
action_count:document_entropy -0.43928    0.04456  -9.859  < 2e-16 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

(Dispersion parameter for binomial family taken to be 1)

    Null deviance: 1356.2  on 6431  degrees of freedom
Residual deviance: 1034.0  on 6428  degrees of freedom
AIC: 1042

Number of Fisher Scoring iterations: 7
### GLM LOOKS OK
> anova(
+      glm(formula = is_service ~ 1
+          , family = binomial(link = 'logit')
+          , data = dat
+      )
+      , glm(formula = is_service ~ action_count
+            , family = binomial(link = 'logit')
+            , data = dat
+      )
+      , glm(formula = is_service ~ action_count + document_entropy
+            , family = binomial(link = 'logit')
+            , data = dat
+      )
+      , glm(formula = is_service ~ action_count + document_entropy + action_count:document_entropy
+            , family = binomial(link = 'logit')
+            , data = dat
+      )
+      , test = "Chisq"
+ )
Analysis of Deviance Table

Model 1: is_service ~ 1
Model 2: is_service ~ action_count
Model 3: is_service ~ action_count + document_entropy
Model 4: is_service ~ action_count + document_entropy + action_count:document_entropy
  Resid. Df Resid. Dev Df Deviance  Pr(>Chi)    
1      6431     1356.2                          
2      6430     1051.1  1  305.135 < 2.2e-16 ***
3      6429     1043.2  1    7.867  0.005035 ** 
4      6428     1034.0  1    9.199  0.002422 ** 
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
Warning messages:
1: glm.fit: fitted probabilities numerically 0 or 1 occurred 
2: glm.fit: fitted probabilities numerically 0 or 1 occurred

And the reproducible code:

list.of.packages <- c('RCurl')
new.packages <- list.of.packages[!(list.of.packages %in% installed.packages()[,"Package"])]
if(length(new.packages)) install.packages(new.packages)

library(RCurl)

x <- getURL("https://rawgit.com/alexmosc/FX_Big_Experiment/master/service_train_saved.csv")
dat <- read.csv(text = x)
dat$X <- NULL

str(dat)
# first part
summary(
     glm(formula = is_service ~ action_count * document_entropy
         , family = binomial(link = 'logit'),
         data = dat
     )
)
# second part
anova(
     glm(formula = is_service ~ 1
         , family = binomial(link = 'logit')
         , data = dat
     )
     , glm(formula = is_service ~ action_count
           , family = binomial(link = 'logit')
           , data = dat
     )
     , glm(formula = is_service ~ action_count + document_entropy
           , family = binomial(link = 'logit')
           , data = dat
     )
     , glm(formula = is_service ~ action_count + document_entropy + action_count:document_entropy
           , family = binomial(link = 'logit')
           , data = dat
     )
     , test = "Chisq"
)

Update:

Let me specify which exactly function I call.

> summary(
+      stats::glm(formula = is_service ~ action_count * document_entropy
+          , family = binomial(link = 'logit'),
+          data = dat
+      )
+ )

Call:
stats::glm(formula = is_service ~ action_count * document_entropy, 
    family = binomial(link = "logit"), data = dat)

Deviance Residuals: 
    Min       1Q   Median       3Q      Max  
-4.0972  -0.1808  -0.1638  -0.1477   3.1019  

Coefficients:
                              Estimate Std. Error z value Pr(>|z|)    
(Intercept)                   -4.06877    0.10143 -40.116  < 2e-16 ***
action_count                   6.04157    0.57882  10.438  < 2e-16 ***
document_entropy               0.29721    0.06866   4.329  1.5e-05 ***
action_count:document_entropy -0.43928    0.04456  -9.859  < 2e-16 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

(Dispersion parameter for binomial family taken to be 1)

    Null deviance: 1356.2  on 6431  degrees of freedom
Residual deviance: 1034.0  on 6428  degrees of freedom
AIC: 1042

Number of Fisher Scoring iterations: 7

> 
> m1 <- stats::glm(formula = is_service ~ 1
+     , family = binomial(link = 'logit')
+     , data = dat
+ )
> m2 <- stats::glm(formula = is_service ~ action_count
+       , family = binomial(link = 'logit')
+       , data = dat
+ )
Warning message:
glm.fit: fitted probabilities numerically 0 or 1 occurred 
> m3 <- stats::glm(formula = is_service ~ action_count + document_entropy
+       , family = binomial(link = 'logit')
+       , data = dat
+ )
Warning message:
glm.fit: fitted probabilities numerically 0 or 1 occurred 
> m4 <- stats::glm(formula = is_service ~ action_count * document_entropy
+       , family = binomial(link = 'logit')
+       , data = dat
+ )
> 
> stats::anova(m1
+       , m2
+       , m3
+       , m4
+       , test = "Chisq")
Analysis of Deviance Table

Model 1: is_service ~ 1
Model 2: is_service ~ action_count
Model 3: is_service ~ action_count + document_entropy
Model 4: is_service ~ action_count * document_entropy
  Resid. Df Resid. Dev Df Deviance  Pr(>Chi)    
1      6431     1356.2                          
2      6430     1051.1  1  305.135 < 2.2e-16 ***
3      6429     1043.2  1    7.867  0.005035 ** 
4      6428     1034.0  1    9.199  0.002422 ** 
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Update 2:

> summary(stats::glm(formula = is_service ~ action_count * document_entropy
+                    , family = binomial(link = 'logit')
+                    , data = dat
+ ))

Call:
stats::glm(formula = is_service ~ action_count * document_entropy, 
    family = binomial(link = "logit"), data = dat)

Deviance Residuals: 
    Min       1Q   Median       3Q      Max  
-4.0972  -0.1808  -0.1638  -0.1477   3.1019  

Coefficients:
                              Estimate Std. Error z value Pr(>|z|)    
(Intercept)                   -4.06877    0.10143 -40.116  < 2e-16 ***
action_count                   6.04157    0.57882  10.438  < 2e-16 ***
document_entropy               0.29721    0.06866   4.329  1.5e-05 ***
action_count:document_entropy -0.43928    0.04456  -9.859  < 2e-16 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

(Dispersion parameter for binomial family taken to be 1)

    Null deviance: 1356.2  on 6431  degrees of freedom
Residual deviance: 1034.0  on 6428  degrees of freedom
AIC: 1042

Number of Fisher Scoring iterations: 7

> 
> stats::anova(
+      stats::glm(formula = is_service ~ 1
+                 , family = binomial(link = 'logit')
+                 , data = dat
+      )
+       , stats::glm(formula = is_service ~ action_count
+                    , family = binomial(link = 'logit')
+                    , data = dat
+       )
+       , stats::glm(formula = is_service ~ action_count + document_entropy
+                    , family = binomial(link = 'logit')
+                    , data = dat
+       )
+       , stats::glm(formula = is_service ~ action_count * document_entropy
+                    , family = binomial(link = 'logit')
+                    , data = dat
+       )
+       , test = "Chisq")
Analysis of Deviance Table

Model 1: is_service ~ 1
Model 2: is_service ~ action_count
Model 3: is_service ~ action_count + document_entropy
Model 4: is_service ~ action_count * document_entropy
  Resid. Df Resid. Dev Df Deviance  Pr(>Chi)    
1      6431     1356.2                          
2      6430     1051.1  1  305.135 < 2.2e-16 ***
3      6429     1043.2  1    7.867  0.005035 ** 
4      6428     1034.0  1    9.199  0.002422 ** 
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
Warning messages:
1: glm.fit: fitted probabilities numerically 0 or 1 occurred 
2: glm.fit: fitted probabilities numerically 0 or 1 occurred

I finally understand this behaviour. The warnings actually happen for models with one and two independent covariates, and these warnings are also issued inside anova. However, when I build a more complex model with interaction there are no warnings in either glm or anova.

> m2 <- stats::glm(formula = is_service ~ action_count
+       , family = binomial(link = 'logit')
+       , data = dat
+ )
Warning message:
glm.fit: fitted probabilities numerically 0 or 1 occurred 
> m3 <- stats::glm(formula = is_service ~ action_count + document_entropy
+       , family = binomial(link = 'logit')
+       , data = dat
+ )
Warning message:
glm.fit: fitted probabilities numerically 0 or 1 occurred




> stats::anova(
+      # stats::glm(formula = is_service ~ 1
+      #            , family = binomial(link = 'logit')
+      #            , data = dat
+      # )
+      #  , stats::glm(formula = is_service ~ action_count
+      #               , family = binomial(link = 'logit')
+      #               , data = dat
+      #  )
+      #  , stats::glm(formula = is_service ~ action_count + document_entropy
+      #               , family = binomial(link = 'logit')
+      #               , data = dat
+      #  )
+       stats::glm(formula = is_service ~ action_count * document_entropy
+                    , family = binomial(link = 'logit')
+                    , data = dat
+       )
+       , test = "Chisq")
Analysis of Deviance Table

Model: binomial, link: logit

Response: is_service

Terms added sequentially (first to last)


                              Df Deviance Resid. Df Resid. Dev  Pr(>Chi)    
NULL                                           6431     1356.2              
action_count                   1  305.135      6430     1051.1 < 2.2e-16 ***
document_entropy               1    7.867      6429     1043.2  0.005035 ** 
action_count:document_entropy  1    9.199      6428     1034.0  0.002422 ** 
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
Warning messages:
1: glm.fit: fitted probabilities numerically 0 or 1 occurred 
2: glm.fit: fitted probabilities numerically 0 or 1 occurred

One more question: is this expected that a more complex glm model proceeds well, but its less complex version is warning prone?

Update 3:

I am checking a relation between fitted probabilites and covariates, lookin specifically for probability equal 1 or 0.

Update 4:

I removed extreme values from predictor space:

list.of.packages <- c('RCurl', 'data.table')
new.packages <- list.of.packages[!(list.of.packages %in% installed.packages()[,"Package"])]
if(length(new.packages)) install.packages(new.packages)

library(RCurl)
library(data.table)

x <- getURL("https://rawgit.com/alexmosc/FX_Big_Experiment/master/service_train_saved.csv")
service_dat <- read.csv(text = x)
service_dat$X <- NULL

service_dat <- as.data.table(service_dat)

service_dat_cleaned <- service_dat[action_count < 10, ]
service_dat_cleaned <- service_dat_cleaned[document_entropy < 10, ]


> summary(service_dat_cleaned)
  action_count       document_entropy     is_service     
 Min.   :-0.084589   Min.   :-0.66652   Min.   :0.00000  
 1st Qu.:-0.079289   1st Qu.:-0.66652   1st Qu.:0.00000  
 Median :-0.057901   Median :-0.66652   Median :0.00000  
 Mean   :-0.018093   Mean   :-0.01084   Mean   :0.02101  
 3rd Qu.: 0.008299   3rd Qu.: 0.61200   3rd Qu.:0.00000  
 Max.   : 4.902059   Max.   : 9.97527   Max.   :1.00000  


Analysis of Deviance Table

Model: binomial, link: logit

Response: is_service

Terms added sequentially (first to last)


                              Df Deviance Resid. Df Resid. Dev  Pr(>Chi)    
NULL                                           6425     1310.1              
action_count                   1  259.038      6424     1051.1 < 2.2e-16 ***
document_entropy               1    7.867      6423     1043.2  0.005035 ** 
action_count:document_entropy  1    9.968      6422     1033.2  0.001593 ** 
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Thank you @Nick Cox. What I am asking is -- without a profound knowledge of statistics -- I wonder if my data create a good model which can be estimated well by this function, but I do not know the bordeline data quality which will not fit well to the function (to the math behind it). So I rely on warnings as a user. I see a contradiction in that two closely related functions behaved differently. — Alexey Burnakov, Oct 24, 2017 at 9:28
First, check which of the four glm function calls in your anova function call are producing warnings. There may be no discrepancy at all. — Scortchi - Reinstate Monica, Oct 24, 2017 at 10:01
Well that's presumably because you get the warning from model fitting. When you call anova(m1, m2) the models are not being re-fit so you don't get any warnings. In the original code, you had all four models fit inside the anova() call. In short, there is no contradiction to be seen. CC @Scortchi. — amoeba, Oct 24, 2017 at 11:05
It's not a general expectation that a more complex model will produce fitted probabilities equal to 0 or 1, no. Suggest you check which fitted probabilities are 0 or 1 (your histograms suggest likely candidates), examine influence diagnostics, & think carefully even about the model with interaction. — Scortchi - Reinstate Monica, Oct 24, 2017 at 12:24
It might help to have more detail about the extremes of the distributions you have plotted as histograms. As @Scortchi has suggested Cook's distance or the hat values would be well worth looking at for all your models together. — mdewey, Oct 24, 2017 at 12:36

mdewey · Accepted Answer · 2017-10-26 11:18:19Z

The coefficient for action_count is about 6. This means that the predicted value for action_count == 50 is going to be about 300 on the log scale. $e^{300}$ is about $10^{130}$ so it is not surprising that when converted from log odds to odds to probabilities it is effectively 1. The same is true to a lesser extent for document_entropy. When you fit the interaction which has a negative coefficient this is presumably sufficient to pull the predicted log odds down to a value which can be converted to a value less than one.

What you do next is to some extent a scientific rather than statistical issue.

Some suggestions:

Drop the extreme observations. This is almost certainly a bad idea. The only case in which it might be helpful is if you have external evidence as to why they differ. It might be for instance that they should never have been included, they were measured using a completely different scale, or something else. Not fitting the model is not a good reason to drop them.
Take better account of the properties of the variables. In this case the explanatory variables appear to be counts. Do we really think that the difference between a count of 1 and 2 is the same as the difference between a count of 50 and 51? If the answer is no then we may in fact be thinking of counts as on a logarithmic scale so transforming them (after possibly ading a small constant, might be a good idea.
The model may differ between the upper and lower parts of the range of the predictors. In this case you could fit a model with separate slopes in the two parts of the range.
There may be other predictors which could be included, possibly with their interaction with existing variables to cope with the model being different for different sub-groups of the population.

No doubt there are other possibilities.

Note that all these suggestions are post hoc as you have made your choice based on having seen the data so really you need to do a completely new study to validate them

I'd agree the warnings are probably due to 'ones' predicted by the models, but I think a 'statistical' (which IMO is also scientific ;) ) addition might be that this extremely high coefficient is due to outliers / data points specific to the exact sample these analyses were based on. In other words, I suspect overfitting. As such, resampling techniques, cross-validation and/or shrinkage could be used to mitigate the effects of such datapoints and give a less-overfitting prone model. — IWS, Oct 25, 2017 at 11:12
@IWS I agree it is a fine line between science and statistics here, I will see what the OP has to say before expanding the answer with suggestions (if necessary). — mdewey, Oct 25, 2017 at 11:21
@mdewey, thank you. What I will do is remove some points with most extreme predictor values from my sample so that I do not get z score above 20, and I will check whether the model has changed significantly. I don't think it will hurt much if I omit extreme behavior, since I cannot know in any way what the maximum value can be in the future given the shape of the distribution (non-finite variance of it). I think the answer is solid. Just for curiosity I will report again after doing the tuning. — Alexey Burnakov, Oct 25, 2017 at 11:55
@AlexeyBurnakov why perform a selective exclusion of values which seem extreme, when the techniques I mentioned in my comment above do these for you in a manner which is less arbitrary? — IWS, Oct 25, 2017 at 12:40
@ IWS, they sound valid. I need to design good cross validation for such unbalanced data. By shrinkage do you mean regularized version of regression? — Alexey Burnakov, Oct 25, 2017 at 12:45

Sextus Empiricus · Accepted Answer · 2017-10-25 16:18:15Z

Not an answer to your original question (glm vs anova(glm), which was more like a code debugging question, since the glm function also gives the same error as the anova(glm) if you just supply the same specific model that has this error/warning),

but anyway, you may be interested in the following plot:

plot(0.1 + dat$action_count,
     2/3 + dat$document_entropy,
     log="xy"
    )
points(0.1 + dat$action_count[which(dat$is_service == 1)],
       2/3 + dat$document_entropy[which(dat$is_service == 1)],
       pch=21,bg=2)

The image shows the distribution of the 1s in your dummy variable, as function of your two regressors (which needed to be transformed slightly to allow easier view on a log-scale).

Yes you can try to fit this mess with your glm model. But it will be a distorted result since the goodness of fit is very bad.

And because of this distorted fitting you can get these surprising results that adding or removing some regressors will result in this model that is warning to be outside of range due to predicting 1s or 0s.

Adding the parameter 'action_count' is predicting 1s because, as mdewey answered, this variable has a large range and the precision becomes bad (not so much over-fitting, since the prediction of p=1, in that upper right corner is actually not so bad, and not really a problem for further evaluation such as evaluation of sums of residuals).
Then adding the cross term is removing the predictions of 1s. That is because the prediction of the upper right corner can now be somewhat decoupled, from the fit in the other three corners.

See for these effects also the difference in the two images below with the predicted p value for the binomial distribution plotted as contours. The first model is without and the second model with the cross term.

So this cross term does not really improve your predictions, except for the upper right corner where the p-value changes from 1 to 0.9999782

Useful for debugging these type of questions:

It is simple to find the source by entering the function int the command line of R, which provides you the object.

And for glm.fit we see the following piece of code:

    eps <- 10 * .Machine$double.eps
    if (family$family == "binomial") {
        if (any(mu > 1 - eps) || any(mu < eps)) 
            warning("glm.fit: fitted probabilities numerically 0 or 1 occurred", 
              call. = FALSE)
    }

So, it may be possible to improve the code (if this would be ever necessary for someone). Namely, the precision near 0 is higher than near 1 (see eps(0) vs eps(1). However, somehow the glm code still gives me the machine eps also if I invert the numbers.

> model1 <-  glm(dat$is_service ~ 1 + dat$action_count + dat$document_entropy, family = binomial(link = 'logit')
> model2 <-  glm(1-dat$is_service ~ 1 + dat$action_count + dat$document_entropy, family = binomial(link = 'logit')


> range(predict(model1), type="response"))
[1] 0.009288089 1.000000000
Warning message:
glm.fit: fitted probabilities numerically 0 or 1 occurred 

> range(predict(model2), type="response"))
[1] 2.220446e-16 9.907119e-01
Warning message:
glm.fit: fitted probabilities numerically 0 or 1 occurred 

> eps(1)
[1] 2.220446e-16

> eps(0)
[1] 2.225074e-308

I suggest that you put the answer of mdewey as the accepted answer since he provided the original thought about the reason why the floating point precision warning occurs. Also his answer is much more brief in this aspect than mine, which is more extending the answer and moving into a direction which is actually a different question. On that aspect this question and their answers are a bit confusing, because the topic has shifted from the initial difference between glm and anova(glm), which did not exist, into reasons for the warning and how to deal with it. — Sextus Empiricus, Oct 25, 2017 at 16:46
Yep, I have done that. Ijust wanted to flag both. Got it now. I will return later. Lots of things to think on. Thank you. — Alexey Burnakov, Oct 25, 2017 at 16:59

Stack Exchange Network

Is there a specific reason why R GLM does not return warnings, while anova(glm) does?

2 Answers 2

Your Answer

Not the answer you're looking for? Browse other questions tagged
r
anova
generalized-linear-model
or ask your own question.

Hot Network Questions

Is there a specific reason why R GLM does not return warnings, while anova(glm) does?

2 Answers 2

Your Answer

Sign up or log in

Post as a guest

Not the answer you're looking for? Browse other questions tagged ranovageneralized-linear-model or ask your own question.

Related

Hot Network Questions

Not the answer you're looking for? Browse other questions tagged
r
anova
generalized-linear-model
or ask your own question.