library(metaSEM) library(dplyr) rm(list = ls()) data <- read.csv("Sample.csv") #____________________________________________________________________________________________ #____________________________________________________________________________________________ #____________________________________________________________________________________________ #### Transforms list of correlations to list of correlation matrices correlations <- cbind(data$lrelXwd, data$lrelXnwb, data$wdXnwb) colnames(correlations) <- c("lrelXwd", "lrelXwb", "wbXwd") cor_function_general <- function(x, p) # x is the vector containing the correlations and p is the number of variables (not correlations) { if(p*(p-1)/2 != length(x)) return(warning("Number of correlations and number of variables is not compatible!")) if(all(is.na(x))) return(vec2symMat(rep(NA,p*(p+1)/2))) vec <- rep(NA, p*(p+1)/2); pos <- 1 for(k in p:1) { vec[pos] <- 1 pos <- pos + k } vec[is.na(vec)] <- x return(vec2symMat(vec)) } # to ensure if this is a real correlation matrix make sure that everything is in order :) cor_function_general(x = colnames(correlations), p = 3) #p = number of variables # create the list of correlation matrices by using the lapply function (generates a list and uses the function defined above # to operate and saves in each list position) cor_list <- lapply(X = 1:dim(correlations)[1], FUN = function(i) cor_function_general(x = correlations[i,], p = 3))#remember to change P names(cor_list) <- data$id head(cor_list) # now we name every correlation matrix in the list, this is needed for the model to work later on... for(i in 1:length(cor_list)){ colnames(cor_list[[i]]) <- c("lrel", "wd", "nwb") row.names(cor_list[[i]]) <- c("lrel", "wd", "nwb") } # now we actually get the data frame to be read into the metaSEM functions # this includes the correlation matrices again divied by each correlation # which is a step back from before, but we also get the inverse of its # asymptotic covariance matrix as the weight matrix cor_data_frame <- Cor2DataFrame(x = cor_list, n = data$N) head(cor_data_frame$data) pattern.na(cor_list, show.na = F) pattern.n(cor_list, n = data$N) ###################################### ###################################### ####### ####### Step 1 ####### ### - calculate average correlation matrix ### - look at tau^2 and average matrix # Fitting Step1 --> compute average correlation matrix REM rand_step1 <- tssem1(Cov = cor_list, n = cor_data_frame$n, method="REM", RE.type = "Diag", acov = "weighted") summary(rand_step1) intercepts <- coef(rand_step1, select = "fixed") # looking at average correlation matrix intercept_cor <- cor_function_general(intercepts, p = 3) #change P colnames(intercept_cor) <- colnames(cor_list[[1]]); row.names(intercept_cor) <- colnames(cor_list[[1]]) round(intercept_cor, 3)