--- always_allow_html: yes output: html_document: default word_document: default pdf_document: default --- ```{r message=FALSE, warning=FALSE, echo=FALSE, results='hide'} ################################################################################ # June 13, 2018 # # The following function creates a well designed HLM table for publishing # # It is written for a single level 1 predictor & a single level 2 predictor # # It is based on the following code: http://www.hermanaguinis.com/JOMR.html # # # # The code was written by Avi Kluger (https://www.avi-kluger.com/) and # # students in his Spring 2018 HLM class including: # # Sarit Pery, # # Noam Markovitch, # # Alon Zoizner, # # Lior Abramson, # # Limor Borut, & # # Benjamin A Katz # # The function is general. To demonstrate the general use of the code, this# # Rmarkdown produces one Table replicating Aguinis et al., (2013), and one # # table using the 1982 study "High School and Beyond". The code reads the # # first dataset from Dropbox and the second from the package mlmRev # # Users can change input after the end of the HLM_Table_Creation below # # For best results knit to Word, change Word to landscape and arrange the # # width of the columns in Word. # ################################################################################ HLM_Table_Creation <- function(x, #x is the expected name of the df DV = "Y", l1v = "Xc", l2v = "Wjc", l2id = "l2id", DV_Label = "Y", l1v_Label = "LMX", l2v_Label = "Leadership climate", full_headers = TRUE) { # x = data, the structure should be data frame with columns of # at least DV, l1v, l2v, l2id # DV = The dependent variable, the outcome. This is a level 1 parameter. # l1v = The level 1 predictor # l2v = The level 2 predictor # l2id = The nesting label (e.g: school from file HSB82) # Labels = The variables` labels # full_headers= set to FALSE if you need columns to be thiner # The following 3 rows are for debug purposes only: # DV = "Y"; l1v = "Xc"; l2v = "Wjc" ; l2id = "l2id" # DV_Label = "Y"; l1v_Label = "LMX"; l2v_Label = "Leadership climate" # full_headers = TRUE # Import packages if (!require("lme4")) install.packages("lme4"); library(lme4) if (!require("weights")) install.packages("weights"); library(weights) df <- data.frame(matrix(nrow = 17, ncol = 4), stringsAsFactors=FALSE) colnames(df) <- c("Null", "RandIntFixSlp", "RandIntRandSlp", "Interaction") #Row names with greek and subscript codes for Rmarkdown # These are the final labels we would like to see in the table display. # CHANGES TO THE LABELS: # If you want to change the rows order: you just change it here, and it'll change. # If you want to change the label, please note you have to change it here # and everywhere the row appears in the results calculation in steps 1-4. rownames(df) <- c("Level 1", paste0("Intercept (", "$\\gamma_{00}$", ")"), paste0(l1v_Label, " (", "$\\gamma_{10}$", ")"), "Level 2", paste0(l2v_Label, " (", "$\\gamma_{01}$", ")"), "Cross-level interaction", paste0(l1v_Label, " ", l2v_Label ," (","$\\gamma_{11}$", ")"), "Variance components", paste0("Within-team (L1) variance (","$\\sigma$","^2^)"), paste0("Intercept (L2) variance (","$\\tau_{00}$", ")"), paste0("Slope (L2) variance (", "$\\tau_{11}$", ")"), paste0("Intercept-slope (L2) covariance (", "$\\tau_{01}$", ")"), "Additional information", "ICC", "-2 log likelihood (FIML)", "Number of estimated parameters", "Pseudo R2") # STEP 1: Null Model # =================== fit1 <- lmer(as.formula(paste0(DV, " ~ (1 |", l2id, ")")),data=x, REML=FALSE) summary(fit1) # fixed-effect parameter estimates FE <- round(summary(fit1)$coefficients , 2) parameters <- FE[, "Estimate"] # fixed-effect parameter standard errors parametersSE <- FE[, "Std. Error"] # Fixed-effect parameter t test t05 <- abs(qt(.975, df= 30)) #Conservative approach to df t01 <- abs(qt(.005, df= 30)) parameters.t <-ifelse(FE[, "t value"] > t05, "*", "") parameters.t <-ifelse(FE[, "t value"] > t01, "**", parameters.t) parameters <-paste0(parameters, parameters.t, " (", parametersSE, ")") df[rownames(df) == "Intercept ($\\gamma_{00}$)","Null"] <- parameters # Compute ICC # Variances V <- round(as.data.frame(VarCorr(fit1)) ["vcov"], 3) ICC <- V[1, 1] / sum(V) ICC df[rownames(df) == "ICC", "Null"] <- rd(ICC, 3) df[rownames(df) == "Within-team (L1) variance ($\\sigma$^2^)", "Null"] <- V [2, ] df[rownames(df) == "Intercept (L2) variance ($\\tau_{00}$)", "Null"] <- V [1, ] df[rownames(df) == "-2 log likelihood (FIML)", "Null"] <- round(as.data.frame(logLik(fit1))*-2, 0) df[rownames(df) == "Number of estimated parameters", "Null"] <- attr(logLik(fit1), "df") df[rownames(df) == "Pseudo R2", "Null"] <- 0 # Test appearance of column 1 df["Null"] #STEP 2: Random Intercept and Fixed Slope Model #============================================== fit2 <- lmer(as.formula(paste0(DV, " ~ (1 |", l2id, ")+ ", l1v, " + ", "I ( ", l2v, " - mean(", l2v,"))")), data=x,REML=F) summary(fit2) # Computing pseudo R-squared yhat2=model.matrix(fit2)%*%fixef(fit2) PR2 = cor(yhat2,x[[DV]])^2 # fixed-effect parameter estimates FE1 <- round(summary(fit2)$coefficients , 2) parameters1 <- FE1[, "Estimate"] # fixed-effect parameter standard errors parametersSE1 <- FE1[, "Std. Error"] # Fixed-effect parameter t test t05 <- abs(qt(.975, df= 30)) #Conservative approach to df, to identify critical t if p < .05 t01 <- abs(qt(.005, df= 30)) #This is to calculate critical t for p < .01 parameters.t <-ifelse(FE[, "t value"] > t05, "*", "") parameters.t <-ifelse(FE[, "t value"] > t01, "**", parameters.t) parameters1 <-paste0(parameters1, parameters.t, " (", parametersSE1, ")") df[rownames(df) == "Intercept ($\\gamma_{00}$)", "RandIntFixSlp"] <- parameters1[1] df[rownames(df) == paste0(l2v_Label, " ($\\gamma_{01}$)"), "RandIntFixSlp"] <- parameters1[3] df[rownames(df) == paste0(l1v_Label," ($\\gamma_{10}$)"), "RandIntFixSlp"] <- parameters1[2] # Log likelihood difference test FISig1 = anova(fit1,fit2) parameters.FIML1 <- ifelse(FISig1[2,8] < .05, "*", "") parameters.FIML1 <- ifelse(FISig1[2,8] < .01, "**", parameters.FIML1) # Variances V1 <- round(as.data.frame(VarCorr(fit2)) ["vcov"], 3) df[rownames(df) == "Within-team (L1) variance ($\\sigma$^2^)", "RandIntFixSlp"] <- V1 [2, ] df[rownames(df) == "Intercept (L2) variance ($\\tau_{00}$)", "RandIntFixSlp"] <- V1 [1, ] df[rownames(df) == "-2 log likelihood (FIML)", "RandIntFixSlp"] <- paste0(round(as.data.frame(logLik(fit2))*-2, 0),parameters.FIML1) df[rownames(df) == "Number of estimated parameters", "RandIntFixSlp"] <- attr(logLik(fit2), "df") df[rownames(df) == "Pseudo R2", "RandIntFixSlp"] <- rd(PR2,2) # Test appearance of column 2 df["RandIntFixSlp"] #STEP 3: Random Intercept and Random Slope model #================================================ fit3 <- lmer(as.formula( paste0(DV, " ~ ", l1v, " + ( ", l1v, " | ", l2id, " )", " + I(", l2v, "-mean(", l2v,") )")), data = x, REML = FALSE) summary(fit3) # fixed-effect parameter estimates FE_step3 <- round(summary(fit3)$coefficients , 2) parameters_step3 <- FE_step3[, "Estimate"] # fixed-effect parameter standard errors parametersSE_step3 <- FE_step3[, "Std. Error"] # Fixed-effect parameter t test t05 <- abs(qt(.975, df= 30)) #Conservative approach to df t01 <- abs(qt(.005, df= 30)) parameters.t_step3 <-ifelse(FE_step3[, "t value"] > t05, "*", "") parameters.t_step3 <-ifelse(FE_step3[, "t value"] > t01, "**", parameters.t) parameters_step3 <-paste0(parameters_step3, parameters.t_step3, " (", parametersSE_step3, ")") df[rownames(df) == "Intercept ($\\gamma_{00}$)", "RandIntRandSlp"] <- parameters_step3[1] df[rownames(df) == paste0(l1v_Label," ($\\gamma_{10}$)"), "RandIntRandSlp"] <- parameters_step3[2] df[rownames(df) == paste0(l2v_Label, " ($\\gamma_{01}$)"), "RandIntRandSlp"] <- parameters_step3[3] # Compute ICC # Variances V_step3 <- round(as.data.frame(VarCorr(fit3)) ["vcov"], 3) ICC_step3 <- V_step3[1, 1] / sum(V_step3) ICC_step3 #df[rownames(df) == "ICC", "RandIntRandSlp"] <- rd(ICC_step3, 3) df[rownames(df) == "Within-team (L1) variance ($\\sigma$^2^)", "RandIntRandSlp"] <- V_step3 [4, ] df[rownames(df) == "Intercept (L2) variance ($\\tau_{00}$)", "RandIntRandSlp"] <- V_step3 [1, ] df[rownames(df) == "Slope (L2) variance ($\\tau_{11}$)", "RandIntRandSlp"] <- V_step3 [2, ] df[rownames(df) == "Intercept-slope (L2) covariance ($\\tau_{01}$)", "RandIntRandSlp"] <- V_step3 [3, ] df[rownames(df) == "Number of estimated parameters", "RandIntRandSlp"] <- attr(logLik(fit3), "df") # Computing pseudo R-squared yhat3=model.matrix(fit3)%*%fixef(fit3) pseudoR2_step3 = round(cor(yhat3,x[[DV]])^2,2); pseudoR2_step3 df[rownames(df) == "Pseudo R2", "RandIntRandSlp"] <- rd(pseudoR2_step3) # Crainceanu & Ruppert (2004) Test of Slope Variance Component compare_models = anova(fit2,fit3) pv_loglik_step3 = compare_models$`Pr(>Chisq)`[2] significance_level_step3 = ifelse(pv_loglik_step3 < .05, "*", "") significance_level_step3 = ifelse(pv_loglik_step3 < .01, "**", significance_level_step3) logLikelihood_step3 = round(as.data.frame(logLik(fit3))*-2, 0) logLik_withSig_step3 = paste0(logLikelihood_step3, significance_level_step3) df[rownames(df) == "-2 log likelihood (FIML)", "RandIntRandSlp"] <- logLik_withSig_step3 # Test appearance of column 1 df["RandIntRandSlp"] #STEP 4: Cross-Level Interaction Model #===================================== fit4 <- lmer(as.formula( paste0(DV, "~ ( ", l1v, "|", l2id, " )", "+ ", l1v, " * I(", l2v, "-mean(", l2v,") )")), data = x, REML = FALSE) summary(fit4) #fixed-effect parameter estimates FE <- round(summary(fit4)$coefficients , 2) parameters <- FE[, "Estimate"] # fixed-effect parameter standard errors parametersSE <- FE[, "Std. Error"] # Fixed-effect parameter t test t05 <- abs(qt(.975, df= 30)) #Conservative approach to df t01 <- abs(qt(.005, df= 30)) parameters.t <-ifelse(FE[, "t value"] > t05, "*", "") parameters.t <-ifelse(FE[, "t value"] > t01, "**", parameters.t) parameters <-paste0(parameters, parameters.t, " (", parametersSE, ")") row <-rownames(df) fixed.result.rows <- c(which(row=="Intercept ($\\gamma_{00}$)"), which(row==paste0(l1v_Label, " (", "$\\gamma_{10}$", ")")), which(row==paste0(l2v_Label, " (", "$\\gamma_{01}$", ")")), which(row==paste0(l1v_Label, " ", l2v_Label ," (","$\\gamma_{11}$", ")"))) df[fixed.result.rows,which(colnames(df)=="Interaction")] <- parameters ### variance components VarComp<-as.data.frame(VarCorr(fit4)) Variances<-round(VarComp[,"vcov"],3) random <- c(which(row=="Intercept (L2) variance ($\\tau_{00}$)"), which(row=="Slope (L2) variance ($\\tau_{11}$)"), which(row=="Intercept-slope (L2) covariance ($\\tau_{01}$)"), which(row=="Within-team (L1) variance ($\\sigma$^2^)")) df[random, which(colnames(df)=="Interaction")] <- Variances ###-2LogLik fit4m2logLik<-round(as.data.frame(logLik(fit4))*-2, 0) fit3m2logLik<-round(as.data.frame(logLik(fit3))*-2, 0) sigtest <- anova(fit4, fit3) m2logLik.Chi <-ifelse(sigtest$`Pr(>Chisq)` < .01,"**", ifelse(sigtest$`Pr(>Chisq)`<.05,"*","")) df[rownames(df) == "-2 log likelihood (FIML)", "Interaction"] <- paste0(fit4m2logLik,m2logLik.Chi[2]) ###Number of estimated parameters df[rownames(df) == "Number of estimated parameters", "Interaction"] <- attr(logLik(fit4), "df") ### Computing pseudo R-squared yhat2=model.matrix(fit4)%*%fixef(fit4) PseudoR2<-rd(cor(yhat2,x[[DV]])^2, 2) df[rownames(df) == "Pseudo R2", "Interaction"] <- PseudoR2 # set the labels for the data frame: ifelse (full_headers, colnames(df) <- c("`Null (Step 1)`", "`Random Intercept and Fixed Slope (Step 2)`", "`Random Intercept and Random Slope (Step 3)`", "`Cross-Level Interaction (Step 4)`"), colnames(df) <- c("`Null (Step 1)`", "`RIFS (Step 2)`", "`RIRS (Step 3)`", "`Interaction (Step 4)`")) return(df) } x <- read.csv("https://www.dropbox.com/s/scpzstk4j95giyt/Aguini.csv?dl=1") df <- HLM_Table_Creation(x) if (!require("mlmRev")) install.packages("mlmRev") library(mlmRev) x <- Hsb82 #load("HLMbook.Rdata") df.schools <- HLM_Table_Creation(x, DV = "mAch", l1v = "cses", l2v = "meanses", l2id = "school", DV_Label = "Math Score", l1v_Label = "CSES", l2v_Label = "MeanSES", full_headers = TRUE) #FFU: We should change the row names to some logic names, and #set the labels in the end of the process, as been done for the columns. # rownames(df) <- c("Level 1 Header", # "Intercept", # "level1 slope", # "Level 2 Header", # "level2 slope", # "Interaction Header", # "Interaction Slope", # "Variance components Header", # "sigma2", # "Intercept l2", # "slope l2 var", # # paste0("Intercept (L2) variance (","$\\tau_{00}$", ")"), # paste0("Slope (L2) variance (", "$\\tau_{11}$", ")"), # paste0("Intercept-slope (L2) covariance (", # "$\\tau_{01}$", ")"), # "Additional information", # "ICC", # "-2 log likelihood (FIML)", # "Number of estimated parameters", # "Pseudo R2") # ```
Table 1 Results of Multilevel Modeling Analysis With Illustrative Data
```{r, message=FALSE, echo=FALSE, warning=FALSE} if (!require("pander")) install.packages("pander");library(pander) set.alignment('center', row.names = 'left') panderOptions('missing', '') pander::pander(df,split.table = Inf, emphasize.rownames = FALSE) ```
Table 1 Results of Multilevel Modeling Analysis With Data from the 1982 study "High School and Beyond" as used in Raudenbush & Bryk (2002)
```{r, message=FALSE, echo=FALSE, warning=FALSE} pander::pander(df.schools,split.table = Inf, emphasize.rownames = FALSE) ```