# GLMM Models in R
# Code by Vanja Dukic, University of Colorado at Boulder, 2019
# 
#
#install.packages('Flury')
#install.packages('lme4')  # be prepared: this one takes a few minutes to install 
#install.packages('')  # be prepared: this one takes a few minutes to install 
#install.packages("hglm", repos="http://R-Forge.R-project.org")

require('Flury')
require('lme4')
require("hglm")

data(challenger)

# Binomial model for probability of any one ring failing;
# independence between rings and launches, adjusted for temperature
ch <- glm(cbind(Damage, 6-Damage) ~ Temp, family = binomial, data = challenger)

with(challenger, plot(Temp, Damage/6,xlim=c(20,100),ylim=c(0,1)))
lines(seq(23,100,1),predict.glm(ch,newdata=data.frame(Temp=seq(23,100,1)),type="response"))
abline(v=32, col = "red", lwd = 2) ## temp when challenger launched


# Computing the probability of AT LEAST ONE ring failing in a launch
# still assuming independence between rings and launches
# adjusted for temperature
lines(seq(23,100,1),1-(1-predict.glm(ch,newdata=
      data.frame(Temp=seq(23,100,1)),type="response"))^6,col='red')


# Direct model of "at least one ring failing" per launch;
# Assuming independence between launches, adjusted for temperature
# The dependence between rings in each launch is now possible
# but we are not modeling it explicitly
 
ch2 <- glm((Damage>0) ~ Temp, family = binomial, data = challenger)

lines(seq(23,100,1),predict.glm(ch2,newdata=data.frame(Temp=seq(23,100,1)),type="response"),col='green')


####################################################################
#
# Equivalent Bernouli models, let's convert all data to 0/1 outcomes
# ie, let's make 6 outcomes for each launch
# So our data now will be a vector of 0/1 outcomes
# there are 23x6 = 138 outcomes total 
#

chal <- cbind(challenger$Damage, challenger$Temp)
chalb = matrix(rep(0,6*23*3),23*6,3)
chalb[,3]=gl(23,6,23*6)
for (i in 1:23) {
for (j in 1:6) {
chalb[((i-1)*6)+j,2] = chal[i,2]
}}

for (i in 1:23) {
for (j in 1:chal[i,1]) {
chalb[((i-1)*6)+j,1] = ((chal[i,1]>0)*1)
}}


# Let's then build for probability of a random ring failing;
# Still assuming independence between launches, adjusted for temperature
# Note this model ch3 is identical to the first one run, stored in "ch"

ch3 <- glm(chalb[,1] ~ chalb[,2], family = binomial)


# A few model diagnostics tools
influence(ch2)
influence(ch)
AIC(ch,ch2)


############################################################################33
# Generalized Linear Mixed Models (GLMM)
#
# Let's now build a model for probability of a ring failing;
# But explicitly allow for positive correlation of outcomes within launches, 
# modeled through launch-specific random effects
#

# There are two packages in R that estimate GLMM models
# GLMM Using GLMER
require('lme4')
ch5 <- glmer(chalb[,1] ~ chalb[,2] + (1|chalb[,3]), family = binomial, verbose=TRUE)

lines(chalb[,2],fitted(ch5), col='blue')

summary(ch5)

##############################################################33
# GLMM Using HGLM
require(hglm)
ch6 <- hglm(fixed = chalb[,1] ~ chalb[,2], random = ~1 | chalb[,3],
fix.disp = 1, family = binomial(), rand.family = gaussian(link = identity))


lines(chalb[,2],ch6$fv, col='purple')

print(summary(ch6), print.ranef=TRUE)



#Summary of the fixed effects estimates:
#
#            Estimate Std. Error t-value Pr(>|t|)  
# (Intercept)   5.0509     3.0961   1.631   0.1051  
# chalb[, 2]   -0.1151     0.0476  -2.418   0.0169 *
#
#
#Summary of the random effects estimates:
#
#                        Estimate Std. Error
#as.factor(chalb[, 3])1    0.0195     0.2126
#as.factor(chalb[, 3])2   -0.0039     0.2122
#as.factor(chalb[, 3])3    0.0006     0.2122
#as.factor(chalb[, 3])4    0.0182     0.2130
#as.factor(chalb[, 3])5   -0.0199     0.2137
#as.factor(chalb[, 3])6   -0.0179     0.2138
#as.factor(chalb[, 3])7   -0.0179     0.2138
#as.factor(chalb[, 3])8   -0.0179     0.2138
#as.factor(chalb[, 3])9   -0.0161     0.2140
#as.factor(chalb[, 3])10  -0.0145     0.2141
#as.factor(chalb[, 3])11  -0.0130     0.2143
#as.factor(chalb[, 3])12  -0.0130     0.2143
#as.factor(chalb[, 3])13   0.0329     0.2142
#as.factor(chalb[, 3])14   0.0329     0.2142
#as.factor(chalb[, 3])15  -0.0104     0.2145
#as.factor(chalb[, 3])16  -0.0094     0.2146
#as.factor(chalb[, 3])17  -0.0075     0.2148
#as.factor(chalb[, 3])18   0.0847     0.2147
#as.factor(chalb[, 3])19  -0.0067     0.2149
#as.factor(chalb[, 3])20  -0.0067     0.2149
#as.factor(chalb[, 3])21  -0.0054     0.2150
#as.factor(chalb[, 3])22  -0.0048     0.2151
#as.factor(chalb[, 3])23  -0.0038     0.2151



AIC(ch3,ch4,ch5)
#    df      AIC
#ch3  2 64.39634
#ch5  3 66.39634

ch6$null.model
#
#Call:  glm(formula = y ~ x - 1, family = family, weights = weights, 
#    offset = off)
#
#Coefficients:
#x(Intercept)   xchalb[, 2]  
#      5.0850       -0.1156  
#
#Degrees of Freedom: 138 Total (i.e. Null);  136 Residual
#Null Deviance:	    191.3 
#Residual Deviance: 60.4 	AIC: 64.4