## script per vedere effetti completa separazione in regressione logistica
## e possibile rimedio con correzione.

install.packages(("mvtnorm"))
install.packages("brglm")
install.packages("rstanarm")
library(mvtnorm)
library(brglm)
library(rstanarm)


sig <- matrix(c(4,2.4,2.4,4), 2,2)
# genero covariate da una normale bivariata
dati<-rmvnorm(100, mean=c(0,0), sigma=sig)
# genero y 0 o 1 a seconda di data soglia (creo separazione)
y <- dati[,2]>0.2+1.5*dati[,1]
cbind(y, dati)
plot(dati, type="n")
points(dati,col=as.integer(y)+1)
abline(0.2,1.5)
# stimo modello di regressione logistica
mio <- glm(y~dati[,1]+dati[,2], binomial)
summary(mio)
#usare il package brglm
mio1 <- brglm(y~dati[,1]+dati[,2], binomial)
summary(mio1)
# usare il paccheto 'bayesiano' rstanarm
mio2 <- stan_glm(y~dati[,1]+dati[,2], binomial)
print(mio2)