#Multilevel Modeling for the Uninitiated: R Code

#Get Packages
library(lme4)
library(nlme)
library(ggplot2)
library(dplyr)
library(foreign)
#---------------------------------------------------------------------------------------------------


#Create Data: Cathedrals and Fire Damage
  #ref: https://ademos.people.uic.edu/Chapter16.html

#Parameters
S.K.England <- -2
S.K.France  <- -3

S.A.England <- -5
S.A.France  <- -1

set.seed(99)

#English Cathedrals
#1
K1 <-rnorm(10, 10, 2) #number of Knights
A1 <-rnorm(10, 10, 2) #number of Archers
Y1 <-S.K.England*K1 + S.A.England*A1 + 80 + rnorm(10, sd=5)

#2
K2 <-rnorm(10, 10, 2) #number of Knights
A2 <-rnorm(10, 10, 2) #number of Archers
Y2 <-S.K.England*K2 + S.A.England*A2 + 75 + rnorm(10, sd=5)

#3
K3 <-rnorm(10, 10, 2) #number of Knights
A3 <-rnorm(10, 10, 2) #number of Archers
Y3 <-S.K.England*K3 + S.A.England*A3 + 90 + rnorm(10, sd=5)

#4
K4 <-rnorm(10, 10, 2) #number of Knights
A4 <-rnorm(10, 10, 2) #number of Archers
Y4 <-S.K.England*K4 + S.A.England*A4 + 100 + rnorm(10, sd=5)

#5
K5 <-rnorm(10, 10, 2) #number of Knights
A5 <-rnorm(10, 10, 2) #number of Archers
Y5 <-S.K.England*K5 + S.A.England*A5 + 120 + rnorm(10, sd=5)

#French Cathedrals
#6
K6 <-rnorm(10, 12, 2) #number of Knights
A6 <-rnorm(10,  6, 2) #number of Archers
Y6 <-S.K.France*K6 + S.A.France*A6 + 50 + rnorm(10, sd=5)

#7
K7 <-rnorm(10, 12, 2) #number of Knights
A7 <-rnorm(10,  6, 2) #number of Archers
Y7 <-S.K.France*K7 + S.A.France*A7 + 40 + rnorm(10, sd=5)

#8
K8 <-rnorm(10, 12, 2) #number of Knights
A8 <-rnorm(10,  6, 2) #number of Archers
Y8 <-S.K.France*K8 + S.A.France*A8 + 42 + rnorm(10, sd=5)

#9
K9 <-rnorm(10, 12, 2) #number of Knights
A9 <-rnorm(10,  6, 2) #number of Archers
Y9 <-S.K.France*K9 + S.A.France*A9 + 59 + rnorm(10, sd=5)

#10
K10 <-rnorm(10, 12, 2) #number of Knights
A10 <-rnorm(10,  6, 2) #number of Archers
Y10 <-S.K.France*K10 + S.A.France*10 + 70 + rnorm(10, sd=5)

#Dataframes
English_Cathedral_1 <-data.frame(Fire_Damage=Y1, Knights=K1, Archers=A1)
English_Cathedral_2 <-data.frame(Fire_Damage=Y2, Knights=K2, Archers=A2)
English_Cathedral_3 <-data.frame(Fire_Damage=Y3, Knights=K3, Archers=A3)
English_Cathedral_4 <-data.frame(Fire_Damage=Y4, Knights=K4, Archers=A4)
English_Cathedral_5 <-data.frame(Fire_Damage=Y5, Knights=K5, Archers=A5)

French_Cathedral_6 <-data.frame(Fire_Damage=Y6, Knights=K6, Archers=A6)
French_Cathedral_7 <-data.frame(Fire_Damage=Y7, Knights=K7, Archers=A7)
French_Cathedral_8 <-data.frame(Fire_Damage=Y8, Knights=K8, Archers=A8)
French_Cathedral_9 <-data.frame(Fire_Damage=Y9, Knights=K9, Archers=A9)
French_Cathedral_10<-data.frame(Fire_Damage=Y10,Knights=K10,Archers=A10)

Dragon1 <-rbind(English_Cathedral_1,English_Cathedral_2,English_Cathedral_3,English_Cathedral_4,English_Cathedral_5,
                       French_Cathedral_6, French_Cathedral_7, French_Cathedral_8, French_Cathedral_9, French_Cathedral_10)
Dragon1$Region <-c(rep(1,10), rep(2,10), rep(3,10), rep(4,10), rep(5,10),
                          rep(6,10), rep(7,10), rep(8,10), rep(9,10), rep(10,10))
Dragon1$Region <-as.factor(Dragon1$Region)
Dragon1$Fire_Damage <-Dragon1$Fire_Damage + 25

Model.Plot.Knights <-ggplot(data = Dragon1, aes(x = Knights, y=Fire_Damage,group=Region))+   
  facet_grid( ~ Region)+    
  geom_point(aes(colour = Region))+ 
  geom_smooth(method = "lm", se = TRUE, aes(colour = Region))+  
  xlab("Knights")+ylab("Fire Damage")+    
  theme(legend.position = "none")   
Model.Plot.Knights  

Model.Plot.Archers <-ggplot(data = Dragon1, aes(x = Archers, y=Fire_Damage,group=Region))+   
  facet_grid( ~ Region)+    
  geom_point(aes(colour = Region))+ 
  geom_smooth(method = "lm", se = TRUE, aes(colour = Region))+  
  xlab("Archers")+ylab("Fire Damage")+    
  theme(legend.position = "none")   
Model.Plot.Archers  

#CENTERING
Dragon1$Knights.C <-scale(Dragon1$Knights, scale=F, center=T)[,]
Dragon1$Archers.C <-scale(Dragon1$Archers, scale=F, center=T)[,]
#---------------------------------------------------------------------------------------------------


#EXAMPLE 1: Fire Damage and Knights

  #Model 1: Knight standard regression
D1_M1 <-lm(Fire_Damage ~ Knights.C, data=Dragon1)
summary(D1_M1)

ggplot(data=Dragon1, aes(x=Knights.C, y=Fire_Damage)) +
  geom_point() +geom_smooth(method=lm, color="red")

  #Model 2: Knight 2-Level Random Intercepts
D1_M2 <-lmer(Fire_Damage ~ Knights.C + (1|Region), data=Dragon1, REML=F)
summary(D1_M2)

Dragon1$Fire_Damage_KPred <-predict(D1_M2, newdata=Dragon1)

ggplot(data=Dragon1, aes(x=Knights.C, y=Fire_Damage, group=Region)) +
  geom_point(aes(color=Region))+
  geom_smooth(method='lm', se=TRUE, aes(colour=Region))+
  xlab("Knights") + ylab("Fire Damage")+
  theme(legend.position = "none")

ggplot(data=Dragon1, aes(x=Knights.C, y=Fire_Damage_KPred, group=Region)) +
  geom_point(aes(color=Region))+
  geom_smooth(method='lm', se=TRUE, aes(colour=Region))+
  xlab("Knights") + ylab("Predicted Fire Damage")+
  theme(legend.position = "none")

  #Model 3: Knight 2-Level Random Intercepts and Slopes
D1_M3 <-lmer(Fire_Damage ~ Knights.C + (Knights.C|Region), data=Dragon1, REML=F)
summary(D1_M3)

Dragon1$Fire_Damage_KPred2 <-predict(D1_M3, newdata=Dragon1)

ggplot(data=Dragon1, aes(x=Knights.C, y=Fire_Damage, group=Region)) +
  geom_point(aes(color=Region))+
  geom_smooth(method='lm', se=TRUE, aes(colour=Region))+
  xlab("Knights") + ylab("Fire Damage")+
  theme(legend.position = "none")

ggplot(data=Dragon1, aes(x=Knights.C, y=Fire_Damage_KPred2, group=Region)) +
  geom_point(aes(color=Region))+
  geom_smooth(method='lm', se=TRUE, aes(colour=Region))+
  xlab("Knights") + ylab("Predicted Fire Damage 2")+
  theme(legend.position = "none")                             
#---------------------------------------------------------------------------------------------------


#EXAMPLE 2: Fire Damage and Archers

#Model 1: Archer standard regression
D2_M1 <-lm(Fire_Damage ~ Archers.C, data=Dragon1)
summary(D2_M1)

ggplot(data=Dragon1, aes(x=Archers.C, y=Fire_Damage)) +
  geom_point() +geom_smooth(method=lm, color="red")

#Model 2: Archer 2-Level Random Intercepts
D2_M2 <-lmer(Fire_Damage ~ Archers.C + (1|Region), data=Dragon1, REML=F)
summary(D2_M2)

Dragon1$Fire_Damage_APred <-predict(D2_M2, newdata=Dragon1)

ggplot(data=Dragon1, aes(x=Archers.C, y=Fire_Damage, group=Region)) +
  geom_point(aes(color=Region))+
  geom_smooth(method='lm', se=TRUE, aes(colour=Region))+
  xlab("Archers") + ylab("Fire Damage")+
  theme(legend.position = "none")

ggplot(data=Dragon1, aes(x=Archers.C, y=Fire_Damage_APred, group=Region)) +
  geom_point(aes(color=Region))+
  geom_smooth(method='lm', se=TRUE, aes(colour=Region))+
  xlab("Archers") + ylab("Predicted Fire Damage")+
  theme(legend.position = "none")

#Model 3: Archer 2-Level Random Intercepts and Slopes
D2_M3 <-lmer(Fire_Damage ~ Archers.C + (Archers.C|Region), data=Dragon1, REML=F)
summary(D2_M3)

Dragon1$Fire_Damage_APred2 <-predict(D2_M3, newdata=Dragon1)

ggplot(data=Dragon1, aes(x=Archers.C, y=Fire_Damage, group=Region)) +
  geom_point(aes(color=Region))+
  geom_smooth(method='lm', se=TRUE, aes(colour=Region))+
  xlab("Archers") + ylab("Fire Damage")+
  theme(legend.position = "none")

ggplot(data=Dragon1, aes(x=Archers.C, y=Fire_Damage_APred2, group=Region)) +
  geom_point(aes(color=Region))+
  geom_smooth(method='lm', se=TRUE, aes(colour=Region))+
  xlab("Archers") + ylab("Predicted Fire Damage")+
  theme(legend.position = "none")
#---------------------------------------------------------------------------------------------------


#EXAMPLE 3: Dragon Hits and Arrows
  #ref: https://data.princeton.edu/pop510/lang2

#Get and transform data
testds <- read.dta("https://data.princeton.edu/pop510/snijders.dta")
testds2 <-testds[testds$schoolnr<100,] #Subset
    #langpost <-DragonHit
    #iavc  <-ArrowNumber
    #schoolnr <-Cathedral
Dragon2 <-data.frame("DragonHit"=testds2$langpost/5, "Arrows.C"=testds2$iq_verb*5-mean(testds2$iq_verb*5), "Cathedral"=testds2$schoolnr)
head(Dragon2)

#Model 1: Arrow Standard Regression
D3_M1 <- lm(DragonHit ~ Arrows.C, data=Dragon2)
summary(D3_M1)

ggplot(data=Dragon2, aes(x=Arrows.C, y=DragonHit)) +
  geom_point() +geom_smooth(method=lm, color="red")

#Model 2: Arrow Level-2 Random Intercepts
D3_M2 <- lmer(DragonHit ~ Arrows.C + (1|Cathedral), data=Dragon2, REML=F)
summary(D3_M2)

Dragon2$DragonHit_Pred <-predict(D3_M2, newdata=Dragon2)

ggplot(data=Dragon2, aes(x=Arrows.C, y=DragonHit, group=Cathedral)) +
  geom_point(aes(color=Cathedral))+
  geom_smooth(method='lm', se=F, aes(colour=Cathedral))+
  xlab("Arrows") + ylab("Dragon Hits")+
  theme(legend.position = "none")

ggplot(data=Dragon2, aes(x=Arrows.C, y=DragonHit_Pred, group=Cathedral)) +
  geom_point(aes(color=Cathedral))+
  geom_smooth(method='lm', se=F, aes(colour=Cathedral))+
  xlab("Arrows") + ylab("Dragon Hits")+
  theme(legend.position = "none")

#Model 3: Arrow Level-2 Random Intercepts and Slopes
D3_M3 <- lmer(DragonHit ~ Arrows.C + (Arrows.C|Cathedral), data=Dragon2, REML=F)
summary(D3_M3)

Dragon2$DragonHit_Pred2 <-predict(D3_M3, newdata=Dragon2)

ggplot(data=Dragon2, aes(x=Arrows.C, y=DragonHit, group=Cathedral)) +
  geom_point(aes(color=Cathedral))+
  geom_smooth(method='lm', se=F, aes(colour=Cathedral))+
  xlab("Arrows") + ylab("Dragon Hits")+
  theme(legend.position = "none")

ggplot(data=Dragon2, aes(x=Arrows.C, y=DragonHit_Pred2, group=Cathedral)) +
  geom_point(aes(color=Cathedral))+
  geom_smooth(method='lm', se=F, aes(colour=Cathedral))+
  xlab("Arrows") + ylab("Dragon Hits")+
  theme(legend.position = "none")
#---------------------------------------------------------------------------------------------------

#Example 4: Dragon Gold and Training
  #ref: https://alexanderdemos.org/Mixed5.html 

#Get and transform data
setwd("C:/Users/Mark.Williamson.2/Desktop/Williamson Data/R/R_data/2022 Data")
testds3 <-read.csv("Sim3level.csv")
head(testds3)
    #Math<- Dragon_gold
    #ActiveTime <-TrainingTime
    #StudentID<- Knight
    #Classroom<- Commander
    #School <-Cathedral
Dragon3 <-data.frame("DragonGold"=testds3$Math, "TrainingTime"=testds3$ActiveTime, "KnightID"=testds3$StudentID,
                         "Commander"=testds3$Classroom, "Cathedral"=testds3$School)
Dragon3$DragonGold <-round(Dragon3$DragonGold,1)
Dragon3$TrainingTime <-round(Dragon3$TrainingTime*10,1)
Dragon3$Commander <-as.factor(Dragon3$Commander)
Dragon3$Cathedral[Dragon3$Cathedral=='Sch1'] <-'C1'
Dragon3$Cathedral[Dragon3$Cathedral=='Sch2'] <-'C2'
Dragon3$Cathedral[Dragon3$Cathedral=='Sch3'] <-'C3'
Dragon3$TrainingTime.C <-scale(Dragon3$TrainingTime, center=T, scale=F)
Dragon3$Commander.2 <-paste(Dragon3$Cathedral,Dragon3$Commander,sep=":")
head(Dragon3)

#Model 1: Training Standard Regression
D4_M1 <- lm(DragonGold ~ TrainingTime.C, data=Dragon3)
summary(D4_M1)

ggplot(data=Dragon3, aes(x=TrainingTime.C, y=DragonGold)) +
  geom_point() +geom_smooth(method=lm, color="red")

#Model 2: Training Level 2 Random Intercepts
D4_M2 <- lmer(DragonGold ~ TrainingTime.C + (1|Commander.2), data=Dragon3, REML=F)
summary(D4_M2)

Dragon3$DragonGold_Pred <-predict(D4_M2, newdata=Dragon3)

ggplot(data=Dragon3, aes(x=TrainingTime.C, y=DragonGold, group=Commander.2)) +
  geom_point(aes(color=Commander.2))+
  geom_smooth(method='lm', se=F, aes(colour=Commander.2))+
  xlab("Training Time") + ylab("Dragon Gold")+
  theme(legend.position = "none")

ggplot(data=Dragon3, aes(x=TrainingTime.C, y=DragonGold_Pred, group=Commander.2)) +
  geom_point(aes(color=Commander.2))+
  geom_smooth(method='lm', se=F, aes(colour=Commander.2))+
  xlab("Training Time") + ylab("Dragon Gold")+
  theme(legend.position = "none")

#Model 3: Training Level 3 Random Intercepts
D4_M3.1 <- lmer(DragonGold ~ TrainingTime.C + (1|Cathedral) + (1|Commander.2), data=Dragon3, REML=F)
summary(D4_M3.1)

D4_M3.2 <- lmer(DragonGold ~ TrainingTime.C + (1|Cathedral) + (1|Cathedral:Commander), data=Dragon3, REML=F)
summary(D4_M3.2)

#Model 4: Training Level 3 Random Intercepts and Slopes
D4_M4 <- lmer(DragonGold ~ TrainingTime.C + (1+TrainingTime.C||Cathedral) +
             (1+TrainingTime.C|Cathedral:Commander), data=Dragon3, REML=F)
summary(D4_M4)

Dragon3$DragonGold_Pred2 <-predict(D4_M4, newdata=Dragon3)

ggplot(data=Dragon3, aes(x=TrainingTime.C, y=DragonGold, group=Commander)) +
  facet_grid(~Cathedral) +
  geom_point(aes(color=Commander))+
  geom_smooth(method='lm', se=TRUE, aes(colour=Commander))+
  xlab("Training Time") + ylab("Dragon Gold")+
  theme(legend.position = "none")

ggplot(data=Dragon3, aes(x=TrainingTime.C, y=DragonGold_Pred2, group=Commander)) +
  facet_grid(~Cathedral) +
  geom_point(aes(color=Commander))+
  geom_smooth(method='lm', se=F, aes(colour=Commander))+
  xlab("Training Time") + ylab("Dragon Gold")+
  theme(legend.position = "none")
#---------------------------------------------------------------------------------------------------


#EXPORT DATASETS (for SAS Code)
head(Dragon1[,1:4])
Dragon2 <-rename(Dragon2, Arrows_C =Arrows.C)
head(Dragon2[,1:3])
Dragon3 <-rename(Dragon3, TrainingTime_C =TrainingTime.C, Commander_2=Commander.2)
head(Dragon3[,1:7])

write.csv(Dragon1[,1:4], "dragon1.csv", row.names=FALSE)
write.csv(Dragon2[,1:3], "dragon2.csv", row.names=FALSE)
write.csv(Dragon3[,1:7], "dragon3.csv", row.names=FALSE)