library(dplyr)
library(haven)
library(psych)
library(PerformanceAnalytics)
library(corrplot)
library(Hmisc)
library(lavaan)
library(semTools)
library(mvnTest)
library(mvnormalTest)
library(RColorBrewer)

CRS_dec <- read_sav("CRS_dec.sav") %>%
  filter(COUNTRY == 3) %>%
  select(CRS_1, CRS_2, CRS_3, CRS_4, CRS_5, CRS_6, CRS_7, CRS_8, CRS_9, CRS_10, CRS_11, 
         CRS_12, CRS_13, CRS_14, CRS_15, CRS_16, CRS_17, CRS_18, CRS_19, CRS_20, CRS_21, 
         CRS_22, CRS_23, CRS_24, CRS_25, CRS_26, CRS_27, CRS_28, CRS_29, CRS_30, CRS_31, 
         CRS_32, CRS_33, CRS_34, CRS_35)


###Syntax for p-value testing###

cor.mtest <- function(mat, ...) {
  mat <- as.matrix(mat)
  n <- ncol(mat)
  p.mat<- matrix(NA, n, n)
  diag(p.mat) <- 0
  for (i in 1:(n - 1)) {
    for (j in (i + 1):n) {
      tmp <- cor.test(mat[, i], mat[, j], ...)
      p.mat[i, j] <- p.mat[j, i] <- tmp$p.value
    }
  }
  colnames(p.mat) <- rownames(p.mat) <- colnames(mat)
  p.mat
}

# matrix of the p-value of the correlation
p.mat <- cor.mtest(CRS_dec)
head(p.mat[, 1:35])

###Correlation analysis (for items)###
#Different styles of correlation presentation:
#(1)
options(max.print=1000000)
rcorr(as.matrix(CRS_dec))
#(2)
cormat <- CRS_dec %>%
  cor()
corrplot.mixed(cormat)
#(3)
corrplot(cormat, type = "upper", col=brewer.pal(n=8, name="RdBu"), tl.col="black")
#(4)
corrplot(cormat, type="upper", order="hclust", 
         p.mat = p.mat, sig.level = 0.05, tl.col="black", insig = "blank")

###Multivariate Normality assumption###
mardia(CRS_dec) ### Mardia Test (Skewness and Kurtosis) for Multivariate Normality (https://search.r-project.org/CRAN/refmans/mvnormalTest/html/mardia.html)
mhz(CRS_dec) ### Henze-Zirkler Test for Multivariate Normality (https://search.r-project.org/CRAN/refmans/mvnormalTest/html/mhz.html)
DH.test(CRS_dec) ### Doornik-Hansen test for Multivariate Normality (https://search.r-project.org/CRAN/refmans/asbio/html/DH.test.html)