##UNCOMMENT THE FOLLOWING LINE AFTER REPLACING [PATH] WITH THE PATH TO THE DIRECTORY WHERE YOU SAVED THE DATA
#setwd("[PATH]")

#Load nlme package to gain access to the linear mixed effects (lme) function
library( nlme )
#Load the lme4 package to gain access to a generalized mixed effects modeling function (lmer)
library( lme4 )

#Read in the data
nls.data <- read.delim( "2-level_nls_data.csv" )
attach( nls.data)


#### Aggregation Demo ####
c.bf.close <- as.numeric( scale(bf.close, scale=FALSE) )
#### Ordinarly Least Squares Regression with Hiearchical Data ####
summary( lm( life.satisfaction ~ c.bf.close * parental.love ) )

aggregated.data <- data.frame( ID=levels(factor(school.id)) )
aggregated.data$bf.close <- as.numeric(scale(as.numeric(tapply( bf.close, school.id, mean, na.rm=T)), scale=FALSE))
aggregated.data$parental.love <- as.numeric(scale(as.numeric(tapply( parental.love, school.id, mean, na.rm=T)), scale=FALSE))
aggregated.data$life.satisfaction <- as.numeric(tapply( life.satisfaction, school.id, mean, na.rm=T))

summary( lm( aggregated.data$life.satisfaction ~ aggregated.data$bf.close * aggregated.data$parental.love ) )

#### Multilevel Modeling: 2-level Models ####

#Step 1: Center your predictors
#Note that parental.love is already centered
c.bf.close <- as.numeric( scale(bf.close, scale=FALSE) )

#2-level model with only a random intercept
two.level.mlm <- lme( life.satisfaction ~ parental.love, random=~1|school.id, na.action="na.exclude" ); summary( two.level.mlm )

#2-level model with a random intercept and random slope for parental.love
two.level.mlm <- lme( life.satisfaction ~ parental.love, random=~1+parental.love|school.id, na.action="na.exclude" ); summary( two.level.mlm )

#2-level model with moderation and a covariate
two.level.mlm <- lme( life.satisfaction ~ school.type + c.bf.close * parental.love, random=~1|school.id, na.action="na.exclude" ); summary( two.level.mlm )


#### Effect Size: Pseudo R-Squared ####
#Step 1: Run Baseline Model and record intercept and residual variances
baseline.mlm <- lme( life.satisfaction ~ 1, random=~1|school.id, subset=is.na(school.type)==FALSE&is.na(c.bf.close)==FALSE&is.na(parental.love)==FALSE, na.action="na.exclude" ); summary( baseline.mlm )
baseline_residual <- 0.8757839**2
baseline_intercept <- 0.0001111166**2

#Step 2: Run comparison model and record intercept and residual variances
two.level.mlm
comparison_residual <- 0.8665774**2
comparison_intercept <- 0.02821983**2

#Snijders & Bosker method
r.squared.1 <- 1 - ( (comparison_intercept + comparison_residual) / (baseline_intercept + baseline_residual) ); r.squared.1
#Find harmonic mean of the number of students per school
library( psych )
students.per.school <- round(harmonic.mean( school.id )); students.per.school
r.squared.2 <- 1 - ( (comparison_intercept + (comparison_residual/students.per.school)) / (baseline_intercept + (baseline_residual/students.per.school)) ); r.squared.2
r.squared.2 <- 0

#Kreft & de Leeuw
r.squared.1 <- ( baseline_residual - comparison_residual ) / baseline_residual; r.squared.1
r.squared.2 <- ( baseline_intercept - comparison_intercept ) / baseline_intercept; r.squared.2
r.squared.2 <- 0

### ICC ###
icc <- baseline_intercept/(baseline_intercept+baseline_residual); icc

#Switch up datasets
detach( nls.data )
three.level.dataset <- read.delim( "3-level_nls_data.csv" )
attach( three.level.dataset)
#### Multilevel Modeling: 3-level Models ####
three.level.mlm <- lme( grades ~ love.mlm*discipline.mlm, random=~1|school.id.long/id, na.action="na.exclude" ); summary( three.level.mlm )

#### Multilevel Modeling: Nested Growth Curve Analysis ####
#2 key features: (1) 3-levels of observations nested in participants nested in schools; (2) random slope for time; (3) moderation by time
bf.close.mlm.c <- bf.close.mlm - mean(bf.close.mlm, na.rm=T)
growth.curve <- lme( arrests ~ sex.mlm + time*consequences*bf.close.mlm.c, random=~1+time|school.id.long/id, na.action="na.exclude" ); summary( growth.curve )


#### CROSS-CLASSIFIED MODELS ####
cross.classified.mlm <- lmer( consequences ~ (1|id) + (1|interviewer.id) + discipline.mlm, na.action="na.exclude" ); summary( cross.classified.mlm )

detach( three.level.dataset )
attach( nls.data )
#### MULTILEVEL MEDIATION (1-1-1 Mediation) ####
#Create two new variables for every level 1 predictor: (1) a grand-mean-centered aggregated average for each group; (2) a group-mean-centered value for each observation within the group.
aggregated.data <- data.frame( school.id=levels(factor(school.id)) )
aggregated.data$parental.monitoring.agg <- as.numeric(tapply( parental.monitoring, school.id, mean, na.rm=T))
aggregated.data$sex.agg <- as.numeric(tapply( sex, school.id, mean, na.rm=T))
new.data <- merge( nls.data, aggregated.data, by="school.id", all.x=TRUE)
detach(nls.data)
attach(new.data)
parental.monitoring.cwc <- parental.monitoring - parental.monitoring.agg
parental.monitoring.agg.c <- parental.monitoring.agg - mean( parental.monitoring.agg, na.rm=T )
sex.cwc <- sex - sex.agg
sex.agg.c <- sex.agg - mean( sex.agg, na.rm=T )

mlm.mediation <- lme( sexual.age ~ sex.agg.c + sex.cwc, random=~1|school.id, na.action="na.exclude" ); summary( mlm.mediation )
mlm.mediation <- lme( parental.monitoring ~ sex.agg.c + sex.cwc, random=~1|school.id, na.action="na.exclude" ); summary( mlm.mediation )
mlm.mediation <- lme( sexual.age ~ sex.agg.c + sex.cwc + parental.monitoring.agg.c + parental.monitoring.cwc, random=~1|school.id, na.action="na.exclude" ); summary( mlm.mediation )