##read in data (replace "~" with actual filepath)
selfeffdata <- read.table("~/selfeffdata.txt")
head(selfeffdata)

##load dependencies (install if not installed)
library(nlme)
library(rockchalk)

##write r2MLMlong function
r2MLMlong <- function(data,covs,random_covs,clusterID,gammas,Tau,sigma2){
  
  if(is.null(covs)==FALSE){
    centered_data <- gmc(data,covs,clusterID)
    phi_w <- var(centered_data[,c(paste0(covs,"_dev"))])
    phi_b <- var(centered_data[,c(paste0(covs,"_mn"))])
    gammas <- matrix(c(gammas),ncol=1)
    f1<-t(gammas)%*%phi_w%*%gammas
    f2<-t(gammas)%*%phi_b%*%gammas
  }
  else{
    f1<-0
    f2<-0
  }
  
  if(is.null(random_covs)==FALSE){
    centered_data_rand <- gmc(data,random_covs,clusterID)
    Sig_w <- var(centered_data_rand[,c(paste0(random_covs,"_dev"))])
    Sig_b <- var(centered_data_rand[,c(paste0(random_covs,"_mn"))])
    m_mat <- matrix(c(colMeans(cbind(1,data[,c(random_covs)]))),ncol=1)
    v1<-sum(diag(Tau[2:nrow(Tau),2:nrow(Tau)]%*%Sig_w))
    v2<-sum(diag(Tau[2:nrow(Tau),2:nrow(Tau)]%*%Sig_b))
  }
  else{
    v1<-0
    v2<-0
    m_mat <- 1  
  }
  
  m<- t(m_mat)%*%Tau%*%m_mat
  
  sigma<-mean(sigma2)
  
  #decompositions
  
  decomp_fixed_within <- f1/sum(f1,f2,v1,v2,m,sigma)
  decomp_fixed_between <-f2/sum(f1,f2,v1,v2,m,sigma)
  decomp_varslopes_within <- v1/sum(f1,f2,v1,v2,m,sigma)
  decomp_varslopes_between <- v2/sum(f1,f2,v1,v2,m,sigma)
  decomp_varmeans <- m/sum(f1,f2,v1,v2,m,sigma)
  decomp_sigma <- sigma/sum(f1,f2,v1,v2,m,sigma)
  
  decomp_fixed_within_w <- f1/sum(f1,v1,sigma)
  decomp_fixed_between_b <-f2/sum(f2,v2,m)
  decomp_varslopes_within_w <- v1/sum(f1,v1,sigma)
  decomp_varslopes_between_b <- v2/sum(f2,v2,m)
  decomp_varmeans_b <- m/sum(f2,v2,m)
  decomp_sigma_w <- sigma/sum(f1,v1,sigma)
  
  #barchart
  
  contributions_stacked <- matrix(c(decomp_fixed_within,decomp_fixed_between,decomp_varslopes_within,decomp_varslopes_between,decomp_varmeans,decomp_sigma,
                                    decomp_fixed_within_w,0,decomp_varslopes_within_w,0,0,decomp_sigma_w,
                                    0,decomp_fixed_between_b,0,decomp_varslopes_between_b,decomp_varmeans_b,0),6,3)
  colnames(contributions_stacked) <- c("total","within","between")
  rownames(contributions_stacked) <- c("fixed slopes (within)",
                                       "fixed slopes (between)",
                                       "slope variation (within)",
                                       "slope variation (between)",
                                       "intercept variation (between)",
                                       "residual (within)")
  
  
  barplot(contributions_stacked, main="Decomposition                     ", horiz=FALSE,
          ylim=c(0,1),col=c("darkred","steelblue","darkred","steelblue","midnightblue","white"),ylab="proportion of variance",
          density=c(NA,NA,30,40,40,NA),angle=c(0,45,0,90,135,0),xlim=c(0,1.5),width=c(.3,.3))
  legend(1.1,.65,legend=rownames(contributions_stacked),fill=c("darkred","steelblue","darkred","steelblue","midnightblue","white"),
         cex=.7, pt.cex = 1,xpd=T,density=c(NA,NA,30,40,40,NA),angle=c(0,45,0,90,135,0))
  
  #create tables for output 
  
  decomp_table <- matrix(c(decomp_fixed_within,decomp_fixed_between,decomp_varslopes_within,decomp_varslopes_between,decomp_varmeans,decomp_sigma,
                           decomp_fixed_within_w,"NA",decomp_varslopes_within_w,"NA","NA",decomp_sigma_w,
                           "NA",decomp_fixed_between_b,"NA",decomp_varslopes_between_b,decomp_varmeans_b,"NA"),6,3)
  colnames(decomp_table) <- c("total","within","between")
  rownames(decomp_table) <- c("fixed slopes (within)",
                              "fixed slopes (between)",
                              "slope variation (within)",
                              "slope variation (between)",
                              "intercept variation (between)",
                              "residual (within)")
  
  R2_table <- matrix(c(decomp_fixed_within,decomp_fixed_between,decomp_varslopes_within,decomp_varslopes_between,decomp_varmeans,
                       decomp_fixed_within+decomp_fixed_between,decomp_fixed_within+decomp_fixed_between+decomp_varslopes_within+decomp_varslopes_between,
                       decomp_fixed_within+decomp_fixed_between+decomp_varslopes_within+decomp_varslopes_between+decomp_varmeans,
                       decomp_fixed_within_w,"NA",decomp_varslopes_within_w,"NA","NA","NA",decomp_fixed_within_w+decomp_varslopes_within_w,"NA",
                       "NA",decomp_fixed_between_b,"NA",decomp_varslopes_between_b,decomp_varmeans_b,"NA",decomp_fixed_between_b+decomp_varslopes_between_b,"NA"),8,3)
  
  colnames(R2_table) <- c("total","within","between")
  rownames(R2_table) <- c("f1","f2","v1","v2","m","f","fv","fvm")
  
  Output <- list(noquote(decomp_table),noquote(R2_table))
  names(Output) <- c("Decompositions","R2s")
  
  return(Output)
}

###fit Linear growth/ person-mean ctr time model
  
  fit_mod1 <- lme(selfeff~1+time_dev,
                  random=~1+time_dev|clusterID,
                  data=selfeffdata,
                  control=lmeControl(optimizer="bobyqa",maxIter=300,msMaxIter=300,niterEM=100,msMaxEval=500))

  Tau_elements <- as.numeric(VarCorr(fit_mod1))
  Tau_est <- matrix(c(Tau_elements[1],
                      Tau_elements[8]*sqrt(Tau_elements[1])*sqrt(Tau_elements[2]),
                      Tau_elements[8]*sqrt(Tau_elements[1])*sqrt(Tau_elements[2]),
                      Tau_elements[2]),2,2)    
  sigma2 <- Tau_elements[3]*c(1,coef(fit_mod1$modelStruct$varStruct, unconstrained=FALSE))
  
  #specify time-specific weights to compute weighted average of time-specific variances
  #used in computing estimate of sigma2 (as observed data are not balanced in the current example)
  allweights <- c(rep(.125,6),rep(.0625,4))
  sigma2 <- sum(sigma2*allweights)
  
  ##compute R-squared measures for Linear growth/ person-mean ctr time model 
   r2MLMlong(selfeffdata,covs=c("time_dev"),random_covs=c("time_dev"),clusterID="clusterID",
                            gammas=c(summary(fit_mod1)$coefficients$fixed[c(2)]),
                            Tau=Tau_est,sigma2=sigma2)
  
  
  ##fit Linear growth/ grand mean-ctr time model
  fit_mod2 <- lme(selfeff~1+I(time_origin-mean(time_origin)),
                    random=~1+I(time_origin-mean(time_origin))|clusterID,
                    data=selfeffdata,
                    control=lmeControl(optimizer="bobyqa",maxIter=300,msMaxIter=300,niterEM=100,msMaxEval=500))
  
  Tau_elements <- as.numeric(VarCorr(fit_mod2))
  Tau_est <- matrix(c(Tau_elements[1],
                      Tau_elements[8]*sqrt(Tau_elements[1])*sqrt(Tau_elements[2]),
                      Tau_elements[8]*sqrt(Tau_elements[1])*sqrt(Tau_elements[2]),
                      Tau_elements[2]),2,2)    
  sigma2 <- Tau_elements[3]*c(1,coef(fit_mod2$modelStruct$varStruct, unconstrained=FALSE))
  sigma2 <- sum(sigma2*allweights)
  
  ##compute R-squared for Linear growth/ grand mean-ctr time 
  r2MLMlong(selfeffdata,covs=c("time_grandmc"),random_covs=c("time_grandmc"),clusterID="clusterID",
                              gammas=c(summary(fit_mod2)$coefficients$fixed[c(2)]),
                              Tau=Tau_est,sigma2=sigma2)
  
  
  ##fit Linear growth/ ctr time at origin model
  
  fit_mod3 <- lme(selfeff~1+time_origin,
                  random=~1+time_origin|clusterID,
                  data=selfeffdata,
                  control=lmeControl(optimizer="bobyqa",maxIter=300,msMaxIter=300,niterEM=100,msMaxEval=500))
  
  Tau_elements <- as.numeric(VarCorr(fit_mod3))
  Tau_est <- matrix(c(Tau_elements[1],
                      Tau_elements[8]*sqrt(Tau_elements[1])*sqrt(Tau_elements[2]),
                      Tau_elements[8]*sqrt(Tau_elements[1])*sqrt(Tau_elements[2]),
                      Tau_elements[2]),2,2)    
  sigma2 <- Tau_elements[3]*c(1,coef(fit_mod3$modelStruct$varStruct, unconstrained=FALSE))
  sigma2 <- sum(sigma2*allweights)
  
  ##compute R-squared for Linear growth/ ctr time at origin model
  r2MLMlong(selfeffdata,covs=c("time_origin"),random_covs=c("time_origin"),clusterID="clusterID",
                            gammas=c(summary(fit_mod3)$coefficients$fixed[c(2)]),
                            Tau=Tau_est,sigma2=sigma2)
  
  
  ##fit Quadratic growth/ ctr time at origin model
  
  fit_mod4 <- lme(selfeff~1+time_origin+I(time_origin^2),
                  random=~1+time_origin+I(time_origin^2)|clusterID,
                  data=selfeffdata,
                  control=lmeControl(optimizer="bobyqa",maxIter=300,msMaxIter=300,niterEM=100,msMaxEval=500))
  
  Tau_elements <- as.numeric(VarCorr(fit_mod4))
  Tau_est <- matrix(c(Tau_elements[1],Tau_elements[10]*sqrt(Tau_elements[1])*sqrt(Tau_elements[2]),Tau_elements[11]*sqrt(Tau_elements[1])*sqrt(Tau_elements[3]),
                      Tau_elements[10]*sqrt(Tau_elements[1])*sqrt(Tau_elements[2]),Tau_elements[2],Tau_elements[15]*sqrt(Tau_elements[2])*sqrt(Tau_elements[3]),
                      Tau_elements[11]*sqrt(Tau_elements[1])*sqrt(Tau_elements[3]),Tau_elements[15]*sqrt(Tau_elements[2])*sqrt(Tau_elements[3]),Tau_elements[3]),3,3)    
  sigma2 <- Tau_elements[4]*c(1,coef(fit_mod4$modelStruct$varStruct, unconstrained=FALSE))
  sigma2 <- sum(sigma2*allweights)
  
  ##compute R-squared for Quadratic growth/ ctr time at origin model
  r2MLMlong(selfeffdata,covs=c("time_origin","tsquared"),random_covs=c("time_origin","tsquared"),clusterID="clusterID",
                           gammas=c(summary(fit_mod4)$coefficients$fixed[c(2,3)]),
                           Tau=Tau_est,sigma2=sigma2)
  
  
  
  ##fit Linear growth/ hetero. error/ ctr time at origin model
  
  fit_mod5 <- lme(selfeff~1+time_origin,
                  random=~1+time_origin|clusterID,
                  weights=varIdent(form=~1|time_origin),
                  data=selfeffdata,
                  control=lmeControl(optimizer="bobyqa",maxIter=300,msMaxIter=300,niterEM=100,msMaxEval=500))
  
  Tau_elements <- as.numeric(VarCorr(fit_mod5))
  Tau_est <- matrix(c(Tau_elements[1],
                      Tau_elements[8]*sqrt(Tau_elements[1])*sqrt(Tau_elements[2]),
                      Tau_elements[8]*sqrt(Tau_elements[1])*sqrt(Tau_elements[2]),
                      Tau_elements[2]),2,2)    
  sigma2 <- Tau_elements[3]*c(1,coef(fit_mod5$modelStruct$varStruct, unconstrained=FALSE))
  sigma2 <- sum(sigma2*allweights)
  
  ##compute R-squared for Linear growth/ hetero. error/ ctr time at origin model 
   r2MLMlong(selfeffdata,covs=c("time_origin"),random_covs=c("time_origin"),clusterID="clusterID",
                            gammas=c(summary(fit_mod5)$coefficients$fixed[c(2)]),
                            Tau=Tau_est,sigma2=sigma2)
  
  
  ##fit Linear growth/ hetero. autoreg. error cov/ ctr time at origin model
  
  fit_mod6 <- lme(selfeff~1+time_origin,
                  random=~1+time_origin|clusterID,
                  correlation = corAR1(form = ~time_origin|clusterID),
                  weights=varIdent(form=~1|time_origin),
                  data=selfeffdata,
                  control=lmeControl(optimizer="bobyqa",maxIter=300,msMaxIter=300,niterEM=100,msMaxEval=500))
  
  Tau_elements <- as.numeric(VarCorr(fit_mod6))
  Tau_est <- matrix(c(Tau_elements[1],
                      Tau_elements[8]*sqrt(Tau_elements[1])*sqrt(Tau_elements[2]),
                      Tau_elements[8]*sqrt(Tau_elements[1])*sqrt(Tau_elements[2]),
                      Tau_elements[2]),2,2)    
  sigma2 <- Tau_elements[3]*c(1,coef(fit_mod6$modelStruct$varStruct, unconstrained=FALSE))
  sigma2 <- sum(sigma2*allweights)
  
  ##compute R-squared for Linear growth/ hetero. autoreg. error cov/ ctr time at origin model
  r2MLMlong(selfeffdata,covs=c("time_origin"),random_covs=c("time_origin"),clusterID="clusterID",
                            gammas=c(summary(fit_mod6)$coefficients$fixed[c(2)]),
                            Tau=Tau_est,sigma2=sigma2)
  
  
  ##fit Conditional linear growth/ hetero, autoreg error cov /ctr time origin model
  
  fit_mod7 <- lme(selfeff~1+time_origin*gender+gpa_c+volunteer_c+gpa_m_centered+volunteer_m_centered,
                  random=~1+time_origin+gpa_c+volunteer_c|clusterID,
                  correlation = corAR1(form = ~time_origin|clusterID),
                  weights=varIdent(form=~1|time_origin),
                  data=selfeffdata,
                  control=lmeControl(optimizer="bobyqa",maxIter=300,msMaxIter=300,niterEM=300,msMaxEval=500,optCtrl=list(maxfun=2e5)))
  
  
  Tau_elements <- as.numeric(VarCorr(fit_mod7))
  Tau_est <- matrix(c(Tau_elements[1],Tau_elements[12]*sqrt(Tau_elements[1])*sqrt(Tau_elements[2]),Tau_elements[13]*sqrt(Tau_elements[1])*sqrt(Tau_elements[3]),Tau_elements[14]*sqrt(Tau_elements[1])*sqrt(Tau_elements[4]),
                      Tau_elements[12]*sqrt(Tau_elements[1])*sqrt(Tau_elements[2]),Tau_elements[2],Tau_elements[18]*sqrt(Tau_elements[2])*sqrt(Tau_elements[3]),Tau_elements[19]*sqrt(Tau_elements[2])*sqrt(Tau_elements[4]),
                      Tau_elements[13]*sqrt(Tau_elements[1])*sqrt(Tau_elements[3]),Tau_elements[18]*sqrt(Tau_elements[2])*sqrt(Tau_elements[3]),Tau_elements[3],Tau_elements[24]*sqrt(Tau_elements[3])*sqrt(Tau_elements[4]),
                      Tau_elements[14]*sqrt(Tau_elements[1])*sqrt(Tau_elements[4]),Tau_elements[19]*sqrt(Tau_elements[2])*sqrt(Tau_elements[4]),Tau_elements[24]*sqrt(Tau_elements[3])*sqrt(Tau_elements[4]),Tau_elements[4]),4,4)
  sigma2 <- Tau_elements[5]*c(1,coef(fit_mod7$modelStruct$varStruct, unconstrained=FALSE))
  sigma2 <- sum(sigma2*allweights)
  
  ##compute R-squared for Conditional linear growth/ hetero, autoreg error cov /ctr time origin model
   r2MLMlong(selfeffdata,covs=c("time_origin","gender","gpa_c","volunteer_c","gpa_m","volunteer_m","tgproduct"),random_covs=c("time_origin","gpa_c","volunteer_c"),clusterID="clusterID",
                            gammas=c(summary(fit_mod7)$coefficients$fixed[c(2,3,4,5,6,7,8)]),
                            Tau=Tau_est,sigma2=sigma2)
  
  
  ##fit conditional model that excludes timeXfemale
  
  fit_mod8 <- lme(selfeff~1+time_origin+gender+gpa_c+volunteer_c+gpa_m+volunteer_m,
                  random=~1+time_origin+gpa_c+volunteer_c|clusterID,
                  correlation = corAR1(form = ~time_origin|clusterID),
                  weights=varIdent(form=~1|time_origin),
                  data=selfeffdata,
                  control=lmeControl(optimizer="bobyqa",maxIter=300,msMaxIter=300,niterEM=300,msMaxEval=500))
  
  Tau_elements <- as.numeric(VarCorr(fit_mod8))
  Tau_est <- matrix(c(Tau_elements[1],Tau_elements[12]*sqrt(Tau_elements[1])*sqrt(Tau_elements[2]),Tau_elements[13]*sqrt(Tau_elements[1])*sqrt(Tau_elements[3]),Tau_elements[14]*sqrt(Tau_elements[1])*sqrt(Tau_elements[4]),
                      Tau_elements[12]*sqrt(Tau_elements[1])*sqrt(Tau_elements[2]),Tau_elements[2],Tau_elements[18]*sqrt(Tau_elements[2])*sqrt(Tau_elements[3]),Tau_elements[19]*sqrt(Tau_elements[2])*sqrt(Tau_elements[4]),
                      Tau_elements[13]*sqrt(Tau_elements[1])*sqrt(Tau_elements[3]),Tau_elements[18]*sqrt(Tau_elements[2])*sqrt(Tau_elements[3]),Tau_elements[3],Tau_elements[24]*sqrt(Tau_elements[3])*sqrt(Tau_elements[4]),
                      Tau_elements[14]*sqrt(Tau_elements[1])*sqrt(Tau_elements[4]),Tau_elements[19]*sqrt(Tau_elements[2])*sqrt(Tau_elements[4]),Tau_elements[24]*sqrt(Tau_elements[3])*sqrt(Tau_elements[4]),Tau_elements[4]),4,4)
  sigma2 <- Tau_elements[5]*c(1,coef(fit_mod8$modelStruct$varStruct, unconstrained=FALSE))
  sigma2 <- sum(sigma2*allweights)
  
  ##compute R-squared for conditional model that excludes timexfemale
  r2MLMlong(selfeffdata,covs=c("time_origin","gender","gpa_c","volunteer_c","gpa_m","volunteer_m"),random_covs=c("time_origin","gpa_c","volunteer_c"),clusterID="clusterID",
                            gammas=c(summary(fit_mod8)$coefficients$fixed[c(2,3,4,5,6,7)]),
                            Tau=Tau_est,sigma2=sigma2)
  
  
  ##fit conditional model that excludes GPA
  
  fit_mod9 <- lme(selfeff~1+time_origin*gender+volunteer_c+volunteer_m,
                  random=~1+time_origin+volunteer_c|clusterID,
                  correlation = corAR1(form = ~time_origin|clusterID),
                  weights=varIdent(form=~1|time_origin),
                  data=selfeffdata,
                  control=lmeControl(optimizer="bobyqa",maxIter=300,msMaxIter=300,niterEM=300,msMaxEval=500))
  
  Tau_elements <- as.numeric(VarCorr(fit_mod9))
  Tau_est <- matrix(c(Tau_elements[1],Tau_elements[10]*sqrt(Tau_elements[1])*sqrt(Tau_elements[2]),Tau_elements[11]*sqrt(Tau_elements[1])*sqrt(Tau_elements[3]),
                      Tau_elements[10]*sqrt(Tau_elements[1])*sqrt(Tau_elements[2]),Tau_elements[2],Tau_elements[15]*sqrt(Tau_elements[2])*sqrt(Tau_elements[3]),
                      Tau_elements[11]*sqrt(Tau_elements[1])*sqrt(Tau_elements[3]),Tau_elements[15]*sqrt(Tau_elements[2])*sqrt(Tau_elements[3]),Tau_elements[3]),3,3)    
  sigma2 <- Tau_elements[4]*c(1,coef(fit_mod9$modelStruct$varStruct, unconstrained=FALSE))
  sigma2 <- sum(sigma2*allweights)
  
  ##compute R-squared for conditional model that excludes GPA
  r2MLMlong(selfeffdata,covs=c("time_origin","gender","volunteer_c","volunteer_m","tgproduct"),random_covs=c("time_origin","volunteer_c"),clusterID="clusterID",
                            gammas=c(summary(fit_mod9)$coefficients$fixed[c(2,3,4,5,6)]),
                            Tau=Tau_est,sigma2=sigma2)
  
  
  ##fit conditional model that excludes volunteer
  
  fit_mod10 <- lme(selfeff~1+time_origin*gender+gpa_c+gpa_m,
                  random=~1+time_origin+volunteer_c|clusterID,
                  correlation = corAR1(form = ~time_origin|clusterID),
                  weights=varIdent(form=~1|time_origin),
                  data=selfeffdata,
                  control=lmeControl(optimizer="bobyqa",maxIter=300,msMaxIter=300,niterEM=300,msMaxEval=500))
  
  Tau_elements <- as.numeric(VarCorr(fit_mod10))
  Tau_est <- matrix(c(Tau_elements[1],Tau_elements[10]*sqrt(Tau_elements[1])*sqrt(Tau_elements[2]),Tau_elements[11]*sqrt(Tau_elements[1])*sqrt(Tau_elements[3]),
                      Tau_elements[10]*sqrt(Tau_elements[1])*sqrt(Tau_elements[2]),Tau_elements[2],Tau_elements[15]*sqrt(Tau_elements[2])*sqrt(Tau_elements[3]),
                      Tau_elements[11]*sqrt(Tau_elements[1])*sqrt(Tau_elements[3]),Tau_elements[15]*sqrt(Tau_elements[2])*sqrt(Tau_elements[3]),Tau_elements[3]),3,3)    
  sigma2 <- Tau_elements[4]*c(1,coef(fit_mod10$modelStruct$varStruct, unconstrained=FALSE))
  sigma2 <- sum(sigma2*allweights)
  
  ##compute R-squared for conditional model that excludes volunteer
  r2MLMlong(selfeffdata,covs=c("time_origin","gender","gpa_c","gpa_m","tgproduct"),random_covs=c("time_origin","gpa_c"),clusterID="clusterID",
                            gammas=c(summary(fit_mod10)$coefficients$fixed[c(2,3,4,5,6)]),
                            Tau=Tau_est,sigma2=sigma2)