dat <- read.csv("http://www.sunthud.com/glm.csv", na.string = "999999")

colnames(dat)[1] <- "university"
dat[,1] <- factor(dat[,1])
dat[,2] <- factor(dat[,2], labels = c("male", "female"))
dat[,5] <- factor(dat[,5], labels = c("science", "social", "humanities", "health", "business"))
dat[,8] <- factor(dat[,8], labels = c("y", "n"))
dat[,9] <- factor(dat[,9], labels = c("y", "n"))
dat[,10] <- factor(dat[,10], labels = c("y", "n"))
dat[,11] <- factor(dat[,11], labels = c("y", "n"))
dat[,8] <- relevel(dat[,8], ref = "n")
dat[,9] <- relevel(dat[,9], ref = "n")
dat[,10] <- relevel(dat[,10], ref = "n")
dat[,11] <- relevel(dat[,11], ref = "n")

library(boot)

# Exercise 11.1

out1 <- glm(give ~ lgpa + gpax, data = dat, family = binomial())
summary(out1)
coef(out1)
exp(coef(out1))
exp(confint(out1))

iv <- expand.grid(lgpa = seq(1, 4, 0.1), gpax = 1:4)
pred <- predict(out1, iv)
prob <- inv.logit(pred)

plot(iv[,1], prob, type = "n", xlab = "Last Semester GPA", ylab = "P(Cheating by Helping Others)")

ivwithprob <- data.frame(iv, prob = prob)

with(ivwithprob[ivwithprob$gpax == 1,], lines(lgpa, prob, col = "orange"))
with(ivwithprob[ivwithprob$gpax == 2,], lines(lgpa, prob, col = "green"))
with(ivwithprob[ivwithprob$gpax == 3,], lines(lgpa, prob, col = "blue"))
with(ivwithprob[ivwithprob$gpax == 4,], lines(lgpa, prob, col = "black"))
legend("topleft", paste0("GPAX = ", 1:4), lty = 1, col = c("orange", "green", "blue", "black"))

out2 <- glm(give ~ lgpa*gpax, data = dat, family = binomial())
anova(out1, out2, test = "Chisq")
summary(out2)

out3 <- glm(give ~ lgpa + gpax + university, data = dat, family = binomial())
anova(out1, out3, test = "Chisq")

university <- factor(1:7)
iv3 <- expand.grid(university = university, lgpa = 3, gpax = 3)
predict(out3, iv3, type = "response")

library(phia)
testInteractions(out3, pairwise = "university")

# log(p(give)/(1 - p(give))) = -1.70 + 0.42LGPA - 0.14GPAX + 1.12U2 - 0.29U3 - 0.47U4 - 0.23U5 + 1.19U6 + 1.66U7

# Exercise 11.2

x <- quantile(dat$SS_ALL)
sscut <- cut(dat$SS_ALL, c(-Inf, x[2:4], Inf))
dat <- data.frame(dat, sscut)

library(MASS)
out4 <- polr(sscut ~ sex + edyear + FGROUP, data = dat, Hess=TRUE)

iv4 <- expand.grid(sex = factor(c("male", "female")), edyear = 1, FGROUP = "business")

predict(out4, newdata = iv4, "probs")

out5 <- polr(sscut ~ sex + edyear + FGROUP + university, data = dat, Hess=TRUE)
anova(out4, out5, test = "Chisq")
summary(out5)

# log(p(ss < 30)/p(ss > 30)) = -2.47 - 0.85female - 0.12edyear + 0.56social + 0.51humanities - 0.03health + 0.26business - 1.25U2 - 0.93U3 - 0.25U4 - 1.47U5 - 0.91U6 - 1.08U7
# log(p(ss < 35)/p(ss > 35)) = -1.31 - 0.85female - 0.12edyear + 0.56social + 0.51humanities - 0.03health + 0.26business - 1.25U2 - 0.93U3 - 0.25U4 - 1.47U5 - 0.91U6 - 1.08U7
# log(p(ss < 41.8)/p(ss > 41.8)) = -0.24 - 0.85female - 0.12edyear + 0.56social + 0.51humanities - 0.03health + 0.26business - 1.25U2 - 0.93U3 - 0.25U4 - 1.47U5 - 0.91U6 - 1.08U7

# Exercise 11.3

medmas <- median(dat$MASTERY)
medper <- median(dat$PERFORM)
gmas <- dat$MASTERY > medmas
gper <- dat$PERFORM > medper

style <- rep(NA, nrow(dat))
style[gmas & gper] <- 1
style[!gmas & gper] <- 2
style[gmas & !gper] <- 3
style[!gmas & !gper] <- 4
style <- factor(style, labels = c("both", "mastery", "performance", "neither"))
style <- relevel(style, ref = "neither")

dat <- data.frame(dat, style)

library(nnet)
out7 <- multinom(style ~ lgpa + gpax, data = dat, Hess = TRUE)
summary(out7)

est <- summary(out7)[["coefficients"]]
se <- summary(out7)[["standard.errors"]]
tval <- est/se
pvalue <- pnorm(-abs(tval))*2

out7a <- multinom(style ~ lgpa + gpax + sex, data = dat, Hess = TRUE)
summary(out7a)

anova(out7, out7a)
est <- summary(out7a)[["coefficients"]]
se <- summary(out7a)[["standard.errors"]]
tval <- est/se
pvalue <- pnorm(-abs(tval))*2

dat <- data.frame(dat, style2 = relevel(style, ref = "mastery"))
out7b <- multinom(style2 ~ lgpa + gpax + sex, data = dat, Hess = TRUE)
summary(out7b)

# log(p(both)/p(neither)) = -2.98 + 0.76LGPA - 0.02GPAX + 0.47female
# log(p(mastery)/p(neither)) = -5.32 + 0.74LGPA + 0.53GPAX - 0.44female
# log(p(performance)/p(neither)) = -2.62 + 0.10LGPA - 0.01GPAX + 0.61female


# Homework 11

# Key from Sutthisan Chumwichan

s.dat <- read.table("assignment11dat.csv", sep=",", header=TRUE)
s.dat1<-data.frame(s.dat[1:10],js.att=apply(s.dat[11:14],1,FUN=mean),ws.att=apply(s.dat[15:18],1,FUN=mean))
s.dat1$gender <- factor(s.dat1$gender,labels = c("male","female"))
s.dat1$s.trav1 <- factor(s.dat1$s.trav1,labels = c("stay","goout"))
s.dat1$s.trav2 <- factor(s.dat1$s.trav2,labels = c("stayhome", "stayBangkok", "travel","aboard"))
s.dat1$s.trav3 <- factor(s.dat1$s.trav3)
#หมายเหตุ
#ข้อ	ตัวแปร	
#1.		อายุ						age (ตัวแปรที่ใช้วิเคราะห์)
#2.		เพศ						gender (ตัวแปรที่ใช้วิเคราะห์)
#3.		การศึกษา					edu
#4.		ภูมิลำเนา					locat
#5.		สถานภาพ					status
#6.		เข้าร่วมสงกรานต์ก่อนหน้า			s.prior (ตัวแปรที่ใช้วิเคราะห์)
#7.		ภาระงาน					w.load
#8.		ไปเล่นสงกรานต์นอกกรุงเทพ			s.trav1 (โจทย์กำหนด)
#9.		สถานที่ไปเล่นสงกรานต์				s.trav2 (โจทย์กำหนด)
#10.	วันที่เดินทาง					s.trav3 (โจทย์กำหนด)
#11		เจตคติต่อการเข้าร่วมกิจกรรมวันสงกรานต์		js.att (ตัวแปรที่ใช้วิเคราะห์)
#11		เจตคติต่อการทำงานวันสงกรานต์			ws.att (ตัวแปรที่ใช้วิเคราะห์)


## Homework 11  11.1 logistic
s.out1 <- glm(s.trav1~age+gender+s.prior+js.att+w.load,data=s.dat1,family=binomial())
summary(s.out1)
## ทุกตัวแปรได้แก่ age gender s.prior js.att และ w.load ไม่ส่งผลต่อ s.trav1 อย่างมีนัยสำคัญ
## อย่างไรก็ตาม js.att ส่งต่อการไปเล่นสงกรานต์นอกกรุงเทพใกล้ระดับนัยสำคัญสูงที่สุด ยิ่ง js.att สูง จะออกนอกบ้านน้อยกว่า

## Homework 11.2 Multinomail logistic regression
s.out2 <- nnet::multinom(s.trav2 ~ age+gender+s.prior+js.att+w.load, data = s.dat1, Hess = TRUE)
summary(s.out2)
s.est2 <- summary(s.out2)[["coefficients"]]
s.se2 <- summary(s.out2)[["standard.errors"]]
s.tval2 <- s.est2/s.se2
s.pvalue2 <- pnorm(-abs(s.tval2))*2
s.pvalue2
##  พบว่า js.att ส่งผลต่อ s.trav2 อย่างมีนัยสำคัญ ในเกือบทุกสมการที่จำแนกคนอยู่กรุงเทพกับกลุ่มอื่น ๆ
## ยกเว้นการจำแนกคนที่อยู่ต่างประเทศ ไม่มีนัยสำคัญ ทั้งนี้อาจจะเกิดจากคนตอบว่าไปเที่ยงต่างประเทศมีจำนวนเพียง 4 คน
## กล่าวคือ หากมี js.att สูง จะมีแนวโน้มอยู่บ้านมากกว่า

## Homework 11.3
# Ordinal Logistic Regression

s.out3 <- MASS::polr(s.trav3~age+gender+s.prior+js.att+w.load, data = s.dat1, Hess=TRUE)
summary(s.out3)
s.tval3 <- summary(s.out3)[["coefficients"]][,3]
s.pvalue3 <- round(pnorm(-abs(s.tval3)) * 2, 4)
s.pvalue3
exp(coef(s.out3))

#เช่นเดียวกับข้อ 11.2 js.att ส่งผลต่อจำนวนวันที่ไปเที่ยวอย่างมีนัยสำคัญที่ .01 ทุก ๆ 1 หน่วยที่เพิ่มขึ้นจะส่งผลจำนวนวันที่ออกท่องเที่ยวสูงเป็น 2.39 เท่า
#โดยที่ตัวแปรอื่น ๆ ได้แก่ age gender s.prior js.att และ w.loadไม่ส่งผลอย่างมีนัยสำคัญ



## Homework 11 Poisson Regression
s.out4 <- glm(I(as.numeric(s.trav3)) ~age+gender+s.prior+js.att+w.load, data = s.dat1, family = poisson())
summary(s.out4)
exp(coef(s.out4))
#เช่นเดียวกับข้อ 11.2 และ 11.3 js.att ส่งผลต่อจำนวนวันที่ไปเที่ยวอย่างมีนัยสำคัญที่ .01 โดยจะทำให้จำนวนวันเพิ่มขึ้น เป็น 0.207 วัน ต่อการเพิ่มขึ้น 1 หน่วยของ js.att
#โดยที่ตัวแปรอื่น ๆ ได้แก่ age gender s.prior js.att และ w.loadไม่ส่งผลอย่างมีนัยสำคัญ

## Homework 11 Negative Binomial Regression
s.out5 <- MASS::glm.nb(I(as.numeric(s.trav3)) ~age+gender+s.prior+js.att+w.load, data = s.dat1, control = glm.control(maxit = 1000))
summary(s.out5)

#ได้ผลเช่นเดียวกับข้อ  11.4 พบว่า js.att ส่งผลต่อจำนวนวันที่ไปเที่ยวอย่างมีนัยสำคัญที่ .01 โดยจะทำให้จำนวนวันเพิ่มขึ้น เป็น 0.207 วัน ต่อการเพิ่มขึ้น 1 หน่วยของ js.att
#โดยที่ตัวแปรอื่น ๆ ได้แก่ age gender s.prior js.att และ w.loadไม่ส่งผลอย่างมีนัยสำคัญ