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
glm
function calls in youranova
function call are producing warnings. There may be no discrepancy at all. $\endgroup$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 theanova()
call. In short, there is no contradiction to be seen. CC @Scortchi. $\endgroup$