#LOGIT Models

logit <-glm(y_bin~ x1 + x2 + x3, family=binomial(link="logit"), data=mydata) 
summary(logit)

#Formatting the result
library(stargazer)
stargazer(logit, type="text", out="logit.pdf")

###Interpretation of the logit estimates:
#Odds ratio = exp(coef(logit))
#Interpretation: When x3 increases by 1 unit, odds of y=1 increases by 112% (212%-100%)
#The odds of y =1 are 2.12 times higher when x3 increases by one unit (keeping all other predictors constant). 
#Example: IF x3 is AGE, and the estimated odds ratio is 2, then
#When respondent's age increase by 1 year, the probability of saying "yes (y=1)"
#increases by 100% (200-100) (2 times higher)

library(mfx) 
oddsratio=logitor(y_bin~ x1 + x2 + x3, data=mydata)
print(oddsratio)

library(arm)
logit2 <-glm(y_bin~ x1 + x2 + x3+opinion, family=binomial(link="logit"), data=mydata) 
summary(logit2)
stargazer(logit2, type="text", out="logit2.pdf")

#Prediction
#Estimating the probability by evaluating each variable at their means
invlogit = function (x) {1/(1+exp(-x))}
invlogit(coef(logit)[1]+coef(logit)[2]*mean(mydata$x1)
         +coef(logit)[3]*mean(mydata$x2)+coef(logit)[4]*mean(mydata$x3))

#Estimating the probability for those who "agree" (5th variable=1) 
invlogit = function (x) {1/(1+exp(-x))}
invlogit(coef(logit2)[1]
           +coef(logit2)[2]*mean(mydata$x1)
           +coef(logit2)[3]*mean(mydata$x2)
           +coef(logit2)[4]*mean(mydata$x3)
           +coef(logit2)[5]*1)

#Estimagint the probability for those who "strongly disagree"
invlogit = function (x) {1/(1+exp(-x))}
invlogit(coef(logit2)[1]
         +coef(logit2)[2]*mean(mydata$x1)
         +coef(logit2)[3]*mean(mydata$x2)
         +coef(logit2)[4]*mean(mydata$x3)
         +coef(logit2)[7]*1)
#Summarize the prediction
allmean<-data.frame(x1=rep(mean(mydata$x1),4),x2=rep(mean(mydata$x2),4),x3=rep(mean(mydata$x3),4),
  opinion=as.factor(c("Str agree","Agree","Disag","Str disag")))

allmean2<-cbind(allmean,predict(logit2,newdata=allmean, type="response", se.fit=TRUE))
allmean2

#Confidence Interval
# Estimating confidence intervals
# Renaming "fit" and "se.fit" columns
names(allmean2)[names(allmean2)=="fit"] = "prob"
names(allmean2)[names(allmean2)=="se.fit"] = "se.prob"

allmean2$ll= allmean2$prob-1.96*allmean2$se.prob
allmean2$ul= allmean2$prob+1.96*allmean2$se.prob
allmean2

#Plot the result
library(ggplot2)
ggplot(allmean2, aes(x=opinion, y = prob))+
  geom_errorbar(aes(ymin= ll, ymax= ul), width = 0.2, lty=1, lwd=1, col="red")+
  geom_point(shape=18, size=5, fill="black")+
  scale_x_discrete(limits = c("Str agree","Agree","Disag","Strdisag"))+
  labs(title= " Predicted probabilities", x="Opinion", y="Pr(y=1)", caption = "add footnote here") +
  theme(plot.title= element_text(family = "sans", face="bold", size=13, hjust=0.5),
  axis.title= element_text(family = "sans", size=9),plot.caption= element_text(family = "sans", size=5))

#Calculate Marginal Effect
library(mfx)
logitmfx(y_bin~ x1+x2+x3, data=mydata)

##Ordered Logit###
# Loading library 
library(MASS)
# Running the ordered logit model
#Opinion has 4 levels (strongly agree ~ strongly disagree)
m1 <-polr(opinion ~ x1 + x2 + x3, data=mydata, Hess=TRUE)
summary(m1)

#adding p-value
m1.coef <-data.frame(coef(summary(m1)))
m1.coef$pval = round((pnorm(abs(m1.coef$t.value), lower.tail= FALSE) * 2),2)
m1.coef
stargazer(m1, type="text", out="m1.pdf")

#Odds Ratio
#Interpretation:the probability of moving up to the higher category (i.e. from agree to strongly agree)
#becomes 167% higher when x1 increases by 1 unit.
m1.or=exp(coef(m1))
m1.or
stargazer(m1.or,type="text",coef=list(m1.or),p.auto=FALSE,out="m1.or.pdf")

#Predicted Probability of each category when all predictors are at their mean value=Mean
#e.g. Probablity of "mean person" stating "agree" is 21.24%.
m1.pred <-predict(m1, type="probs")
summary(m1.pred)

#Predicted Probability if we set x3=1 or 2, evaluated at ave(x1),ave(x2)
setup1 <-data.frame(x1=rep(mean(mydata$x1),2), x2=rep(mean(mydata$x2),2), x3=c(1,2))
setup1
setup1[, c("pred.prob")] <-predict(m1, newdata=setup1, type="probs")
setup1

#Marginal Effect
library(erer)
x <-ocME(m1, rev.dum=TRUE)
x
x$out

#Multinomial Logit
library(nnet)
hsb2 <- read.table('https://stats.idre.ucla.edu/wp-content/uploads/2016/02/hsb2-2.csv', header=T, sep=",")
summary(hsb2)
head(hsb2)

table(hsb2$ses)
ses=factor(hsb2$ses)
#Changing the referece value as 2=middle for ses.
hsb2$ses2 = relevel(ses, ref = "2")

multi1 = multinom(ses2 ~ science + socst+  female, data=hsb2)
summary(multi1)
stargazer(multi1,type="text",out="multi1.pdf")