library(arm)

# no separation case

n <- 100
x1 <- rnorm (n)
x2 <- rbinom (n, 1, .5)
b0 <- 1
b1 <- 1.5
b2 <- 2
y <- rbinom (n, 1, invlogit(b0+b1*x1+b2*x2))

loglik <- function(x){
  sum(y* log(1/(1+exp(-b0-b1*x1-x*x2))) +(1-y)*log(1- 1/(1+exp(-b0-b1*x1-x*x2)) )  )
}


mod1 <- glm(y ~ x1 + x2, family=binomial(link="logit"))
summary(mod1)
loglik <- Vectorize(loglik)
curve(loglik(x), xlim= c(-5, 10))
segments(mod1$coefficients[3], -250, mod1$coefficients[3],
        loglik(mod1$coefficients[3]), col ="red")
points(mod1$coefficients[3],
       loglik(mod1$coefficients[3]), col ="red", pch =7)


# separation case

y <- ifelse (x2==1, 1, y)

mod2 <- glm(y ~ x1 + x2, family=binomial(link="logit"))
summary(mod2)
loglik <- Vectorize(loglik)
curve(loglik(x), xlim= c(0, 30))
segments(mod2$coefficients[3], -60, mod2$coefficients[3],
         loglik(mod2$coefficients[3]), col ="red")
points(mod2$coefficients[3],
       loglik(mod2$coefficients[3]), col ="red", pch =7)

