Non-significant logistic regression coefficient in simulated data

Question

I'm new in all this R and stats thing... I'm trying to understand de "logic" of the binary logistic regression but I'm stuck with a surprising result.

I have created a data frame with simulated data with n=10,000 and two variables: and independent variable ("prop") which is a proportion (from 0.1 to 0.5), and "VD" a binary dependent variable (1,0). The data frame is created so that all cases with "prop" < 0.3 have a value 0 in VD and all cases with "prop" >0.3 have a value 1 in VD. 50% of the cases with prop=0.3 have VD=1 and 50% VD=0.

Ok, this looks like a data frame with the ideal structure for a very good result in a binary logistic regression... since the two variables are made to be correlated, I expected the logistic regression coefficient of "prop" to be statistically significant...So, if you run the model, Nagelkerke R2 is high as expected... but I cannot understand why the logistic regression coefficient is non-significant

Here is the code:

library(DescTools)
library(ggplot2)
data<-data.frame(aleat=runif(10000, min=0, max=1))
data$ID<-seq.int(nrow(data))
data$prop[data$ID<=2000]<-.1
data$prop[data$ID>2000 & data$ID<=4000]<-.2
data$prop[data$ID>4000 & data$ID<=6000]<-.3
data$prop[data$ID>6000 & data$ID<=8000]<-.4
data$prop[data$ID>8000 & data$ID<=10000]<-.5

data$VD[data$prop<0.3]<-0
data$VD[data$prop>0.3]<-1
data$VD[data$ID>4000 & data$ID<=5000]<-0
data$VD[data$ID>5000 & data$ID<=6000]<-1

ggplot(data, aes(x=prop, y=VD))+
  geom_point()

model<-glm(VD ~ prop, data=data, family=binomial())
summary(model)
PseudoR2(modelo5, which = "Nagelkerke")

Thanks in advance!

Answer 1

What you are seeing is the phenomenon of complete separation, related to the Hauck-Donner effect (do a web search on either of these terms for much more information). You will have seen the warning

glm.fit: fitted probabilities numerically 0 or 1 occurred

When you have complete separation, the Wald approximation used in summary() breaks down. You can use anova() to do a likelihood ratio test instead:

anova(model,test="Chisq")
Analysis of Deviance Table

Model: binomial, link: logit

Response: VD

Terms added sequentially (first to last)


     Df Deviance Resid. Df Resid. Dev  Pr(>Chi)    
NULL                  9999    13862.9              
prop  1    11090      9998     2772.6 < 2.2e-16 ***