setwd("C:/Users/Gary/Dropbox/Post doc/R workshop/Module from Centre for Multilevel Modelling")
setwd("C:/Users/t8s/Documents/Dropbox/Post doc/R workshop/Module from Centre for Multilevel Modelling")

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

pop <- rnorm(20000,10,5)
pop[10001:20000] <- pop[10001:20000]+5
sim_data <- data.frame(pop)
sim_data$cat[(1:10000)]<-1
sim_data$cat[(10001:20000)]<-2
densityplot(~pop, groups = cat, data=sim_data, xlim = c(-5,30), plot.points = FALSE)

pop <- rnorm(20000,10,3)
pop[10001:20000] <- pop[10001:20000]+5
sim_data <- data.frame(pop)
sim_data$cat[(1:10000)]<-1
sim_data$cat[(10001:20000)]<-2
densityplot(~pop, groups = cat, data=sim_data, xlim = c(-5,30), plot.points = FALSE)

pop <- rnorm(20000,10,1.5)
pop[10001:20000] <- pop[10001:20000]+5
sim_data <- data.frame(pop)
sim_data$cat[(1:10000)]<-1
sim_data$cat[(10001:20000)]<-2
densityplot(~pop, groups = cat, data=sim_data, xlim = c(-5,30), plot.points = FALSE)

pop <- rnorm(20000,10,1)
pop[10001:20000] <- pop[10001:20000]+5
sim_data <- data.frame(pop)
sim_data$cat[(1:10000)]<-1
sim_data$cat[(10001:20000)]<-2
densityplot(~pop, groups = cat, data=sim_data, xlim = c(-5,30), plot.points = FALSE)

#############################################################
install.packages("lme4")
library(lme4)
library(lattice)
mydata <- read.table(file = "SSLS.txt", sep = ",", header = TRUE)

library(psych)

summary(mydata)

mydata8<-mydata[mydata$schoolid %in% c(1:8),c("schoolid", "score")]
#mydata8<-mydata[mydata$schoolid <= 8,c("schoolid", "score")]

describe(mydata8[,c("score")])
describeBy(mydata8[,c("score")], group = mydata8$schoolid)

densityplot(~score, data = mydata8)
densityplot(~score, groups = schoolid, data = mydata8, auto.key = TRUE, lwd = 2)

#a multilevel model with school effect
ml.reg1 <- lmer(score ~ (1 | schoolid), data = mydata, REML = FALSE)
summary(ml.reg1)

sl.reg0 <- lm(score ~ 1, data = mydata)
summary(sl.reg0)

logLik(ml.reg1)
logLik(sl.reg0)

anova(ml.reg1, sl.reg0)

#The 508 school level residuals
#The ranef function returns a list of object, in this case, only one objects 
#that contains a data frame with
#the residuals and posterior variance.
u0 <- ranef(ml.reg1, condVar = TRUE)

#The following command extracts the posterior variance.
#Taking the square root of this variance gives the standard error
#Since u0 is a list of dataframe, we need to use two sets of square brackets
#to access its element
u0se <- sqrt(attr(u0[[1]],"postVar")[1, , ])

#The following command extracts the rowname, changes the names into numeric, 
#and assign it into the variable u0tab.
schoolid <- as.numeric(rownames(u0[[1]]))

#assign the schoolid, residuals and standard errors into "u0tab"
u0tab <- cbind(schoolid, u0[[1]], u0se)

#assign names to the column
colnames(u0tab) <- c("schoolid", "u0", "u0se")

#sort the data by the magnitude of residual
u0tab <- u0tab[order(u0tab$u0),]

#Create a rank variable based on the magnitude of the residuals.
u0tab$u0rank <- c(1:length(u0tab$schoolid))

#create the catepillar plots
plot(u0tab$u0rank, u0tab$u0,  xlab = "School (rank)", ylab = "Residuals", col = "red")

#using the segment function to add the confidence intervals
segments(u0tab$u0rank, u0tab$u0 - 1.96*u0tab$u0se, u0tab$u0rank, u0tab$u0 + 1.96* u0tab$u0se)

#add a horizontal reference line at y = 0
abline(h = 0, col = "blue")

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

ml.reg2 <- lmer(score ~ cohort90 + (1 | schoolid), 
                data = mydata, REML = FALSE) 
summary(ml.reg2)


#compute the fitted (predicted) values from the model
predscore2 <- fitted(ml.reg2)

#use the xyplot to plot a line graphs. 
#the parameter type = c("p","l") specifies that we wants both points and lines
xyplot(predscore2 ~ cohort90, data = mydata, groups = schoolid, type = c("p","l"))

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

ml.reg3 <- lmer(score ~ cohort90 + (1+ cohort90 | schoolid), data = mydata, REML = FALSE)
summary(ml.reg3)

#95% coverage interval for the intercept
30.60963 - 1.96*sqrt(42.8580)
30.60963 + 1.96*sqrt(42.8580)

#95% coverage interval for the slope
1.23390 - 1.96*sqrt(0.1606)
1.23390 + 1.96*sqrt(0.1606)

anova(ml.reg2, ml.reg3)

cov_matrix <- VarCorr(ml.reg3)
cov_matrix[1]

#extract the random intercept and random slope for each school
myrandomeff <- ranef(ml.reg3, condVar = TRUE)

plot(myrandomeff[[1]], xlab = "Intercept", ylab ="Slope of cohort90")

abline(h = 0, col = "red")
abline(v = 0, col = "red")

#compute the predicted score for each school
predscore3 <- fitted(ml.reg3)

#plotting the predicted score vs cohort for each school
xyplot(predscore3 ~ cohort90, data = mydata, groups = schoolid, type = c("p","l"), col = "blue")

#plotting the between-school variance as a function of cohort

x <- c(-6:8)
y <- 42.858+2*(-1.024)*x+0.161*x^2
plot(y ~ x, type = "l")

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

ml.reg4 <- lmer(score ~ cohort90 + female + (1 + cohort90 | schoolid), 
                data = mydata, REML = FALSE)

summary(ml.reg4)
anova(ml.reg3, ml.reg4)

ml.reg5 <- lmer(score ~ cohort90 + female + (1 + cohort90 + female| schoolid),
                data = mydata, REML = FALSE)

summary(ml.reg5)

anova(ml.reg4, ml.reg5)

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

mydata$sclass_cat <- factor(mydata$sclass, levels = c("1","2","3","4"), 
                            labels = c("Managerial and Professional","Intermediate","Working","unclassified"))

ml.reg6 <- lmer(score ~ cohort90 + female + sclass_cat + (1+cohort90 | schoolid), data= mydata, REML = FALSE)
summary(ml.reg6)

#This model fails to converge!
ml.reg7 <- lmer(score ~ cohort90 + female + sclass_cat + 
                  (1+cohort90 + sclass_cat | schoolid), data= mydata, REML = FALSE)
summary(ml.reg7)

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

table(mydata$schtype)
table(mydata$schurban)

mydata$schtype_cat <- factor(mydata$schtype, levels = c("0","1"), labels = c("State-funded","independent"))
mydata$schurban_cat <- factor(mydata$schurban, levels = c("0","1"), labels = c("Town or rural", "Urban"))

ml.reg8 <- lmer(score ~ cohort90 + female + sclass_cat + schtype_cat + (1+cohort90 | schoolid), data= mydata, REML = FALSE)
summary(ml.reg8)

ml.reg9 <- lmer(score ~ cohort90 + female + sclass_cat + schtype_cat + schurban_cat +(1+cohort90 | schoolid), data= mydata, REML = FALSE)
summary(ml.reg9)

ml.reg10 <- lmer(score ~ cohort90 + female + sclass_cat + schtype_cat + schurban_cat + cohort90:schtype_cat 
                 + (1+cohort90 | schoolid), data= mydata, REML = FALSE)
summary(ml.reg10)

library(effects)
plot(effect("cohort90:schtype_cat", ml.reg10, 
            xlevels(schtype_cat = c("State-funded","independent"))), multiline = TRUE)
