rm(list=ls())

# library("mffSM")
# library("colorspace")


##### Data on types of insurance (Assignment 1)
##### =========================================
(load("GLM_1_mp.RData"))

head(mp)
dim(mp)
summary(mp)

with(mp, table(agef))

### Some age categories will be merged   
mp <- transform(mp, agef2 = factor(1 * (agef %in% c("<55", "55-60")) +
                                   2 * (agef == "60-65") +
                                   3 * (agef == "65-70") +
                                   4 * (agef == "70-75") + 
                                   5 * (agef == "75-80") +
                                   6 * (agef == "80-85") +
                                   7 * (agef %in% c("85-90", ">90")),
                            labels = c("<60", "60-65", "65-70", "70-75", "75-80", "80-85", ">85")),
                    fhmo = factor(hmo, levels = 0:1, labels = c("Regular", "HMO")),
                    fdied = factor(died, levels = 0:1, labels = c("Alive", "Died")),
                    fwhite = factor(white, levels = 0:1, labels = c("Dark", "White")))
                
with(mp, table(agef2))                        
summary(mp)



##### GLM's
##### =============================
###
### Some reasonable ones (here directly without interactions)
###

  ### Gamma (identity)
gm11 <- glm(los ~ fhmo + admit + fdied, data = mp, Gamma(link="identity"))
summary(gm11)

  ### Gamma (log)
gm21 <- glm(los ~ fhmo + admit + fdied, data = mp, Gamma(link="log"))
summary(gm21)

# Sandwich estimator
library(sandwich)
library(lmtest)

# Estimated asymptotic variance matrix of regression 
sandwich(gm21)

# S.E. 
sqrt(diag(sandwich(gm21)))

# summary with the help of sandwich 
coeftest(gm21, vcov=sandwich(gm21))
# or simply 
coeftest(gm21, vcov=sandwich)

# To by compared with 
summary(gm21)


## Comparison of nested models
gm21w <- glm(los ~ fhmo + admit + fdied*fwhite, data = mp, Gamma(link="log"))

waldtest(gm21, gm21w)
# to be compared with 
anova(gm21, gm21w, test="LRT")




## A small note about computing p-values... 
# p-values in coeftest
2*(1-pnorm(abs(coeftest(gm21)[,3])))

# p-values in summary.glm() when dispersion is estimated 
2*(1-pt(abs(coeftest(gm21)[,3]), df=gm21$df.residual))
# to be compared with 
summary(gm21)$coefficients[,4]


# # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # #

### Cancer data (Assignment 2)
### ==========================
# rm(list = ls())
#
library("Epi")
library("splines")
#
(load("GLM_2_iev.RData"))


### Transformations
### +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
iev <- transform(iev, agegr5 = Relevel(agegr, list("25-34" = 1, "35-44" = 2, "45-54" = 3, "55-64" = 4, ">=65" = 5:6)))
iev <- transform(iev, tobgr4 = Relevel(tobgr, list("None" = 1, "1-4 g/day" = 2, "5-29 g/day" = 3:6, ">=30 g/day" = 7:9)))
iev <- transform(iev, alctotf = cut(alctot, breaks = c(-Inf, 0, 25, 50, 75, 100, 150, Inf)))
levels(iev[, "alctotf"]) <- c("0", "(0,25]", "(25,50]", "(50,75]", "(75,100]", "(100,150]", ">150")


### One of the initial fits 
summary(fit2 <- glm(case ~ alctotf + agegr5 + tobgr4, family = binomial, data = iev))
# quasibinomial 
summary(fit2.q <- glm(case ~ alctotf + agegr5 + tobgr4, family = quasibinomial, data = iev))

drop1(fit2, test="LRT")
# p-value for tobgr4
test.stat <- 53.325; # 712.18 - 658.86 
1-pchisq(test.stat, 3)


(D1q <- drop1(fit2.q, test="LRT")); 
(estim.disp <- summary(fit.q)$dispersion); # 0.7967584
(test.stat.q <- (D1q[4,2] - D1q[1,2])/estim.disp); # (712.18-658.86)/0.7967584);
1-pchisq(test.stat.q, D1q[4,1])#   1-pchisq(66.927, df=3)


# # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # d

print(load("GLM_3_nels.RData"))

(xtab4 <- with(nels, table(sesmed, fa.educ, fa.wrk, region)))
qq4 <- as.data.frame(xtab4, responseName = "N")

summary(glm(N ~ (sesmed + fa.educ + fa.wrk + region)^4, family = poisson, data = qq4))
summary(glm(N ~ (sesmed + fa.educ + fa.wrk + region)^4, family = quasipoisson, data = qq4)); # does not make sense 

summary(fit3 <- glm(N ~ (sesmed + fa.educ + region)^3, family = poisson, data = qq4))
summary(fit3.q <- glm(N ~ (sesmed + fa.educ + region)^3, family = quasipoisson, data = qq4))

# Ratio of S.E. from quasi vs standard 
summary(fit3.q)$coefficients[,2]/summary(fit3)$coefficients[,2]
# To be compared with 
sqrt(summary(fit3.q)$dispersion) # Sqrt of estimated dispersion parameters

# A better model 
summary(fit4 <- glm(N ~ sesmed*fa.wrk + sesmed*fa.educ*region, family = poisson, data = qq4))
summary(fit4.q <- glm(N ~ sesmed*fa.wrk + sesmed*fa.educ*region, family = quasipoisson, data = qq4))


# Sandwich does not seem to be a good idea here 
diag(sandwich(fit4))/summary(fit4)$coefficients[,2]