# Project: The 19-Item Environmental Knowledge Test (EKT-19): A Short, Psychometrically Robust Measure of Environmental Knowledge
# Last Edit: 27/09/22
# Analysis: CFA


# Conducting CFA for 30 items

# Install CFA packages
install.packages("readxl")
install.packages("lavaan")
install.packages("lavaanPlot")
install.packages("Rcpp")
install.packages("Rcsdp")

# Begin saving data output in separate file
sink(file = "EK_CFA_Output.txt") 

# Read in libraries
library(readxl)
library(lavaan)
library(lavaanPlot)
library(dplyr)
library(psych)
library(Rcpp)
library(Rcsdp)

# Read in data
EK_cfa <- read_excel("C:\\Users\\lkp32\\R\\EnvKnowledge\\EKT-19_Cleaned.xlsx")

# Exclude failed attention checks
EK_cfa_filtered  <- filter(EK_cfa, Attention == 1)
View(EK_cfa_filtered) # Visually inspecting data
describe(EK_cfa_filtered)

# Create three factor model with 30 items 
# Model f1 represents 'system' knowledge, Model f2 represents 'action-related' knowledge, Model f3 represents 'effectiveness' knowledge 

# Create three factor model with 30 items 

three_factor <- 'f1 =~ EK1 + EK2 + EK3 + EK4 + EK5 + EK6 + EK7 + EK9 + EK11 + EK14 + EK23 + EK24 + EK25 + EK26 + EK27 + EK30
          f2 =~ EK10 + EK8 + EK15 + EK16 + EK28 + EK29
          f3 =~ EK12 + EK13 + EK17 + EK18 + EK19 + EK20 + EK21 + EK22'

fit_three_factor <- cfa(three_factor, data = EK_cfa_filtered, ordered = TRUE) # 'ordered' is for binary variables
summary(fit_three_factor, standardized=TRUE, ci=TRUE, fit.measures=TRUE)
modindices(fit_three_factor)

# Create one factor model with 30 items 

one_factor <- 'f1 =~ EK1 + EK2 + EK3 + EK4 + EK5 + EK6 + EK7 + EK8 + EK9 + EK10 + EK11 + EK12 + EK13 + EK14 + EK15 +EK16 +EK17 + EK18 +EK19 +EK20 +EK21+EK22 + EK23 + EK24 + EK25 + EK26 + EK27 + EK28 +EK29 + EK30'

fit_one_factor <- cfa(one_factor, data = EK_cfa_filtered, ordered = TRUE) # 'ordered' is for binary variables
summary(fit_one_factor, standardized=TRUE, ci=TRUE, fit.measures=TRUE)
modindices(fit_one_factor)

# Create plots for both 30 item models 

lavaanPlot(model = fit_three_factor, node_options = list(shape = "box", fontname = "Helvetica"), edge_options = list(color = "grey"),
           coefs = T, stand=T)

lavaanPlot(model = fit_one_factor, node_options = list(shape = "box", fontname = "Helvetica"), edge_options = list(color = "grey"),
           coefs = T, stand=T, covs = T)


# Comparing 30 item one-factor vs three-factor models 

lavTestLRT(fit_one_factor, fit_three_factor, scaled.shifted = TRUE) # Compare nested models
anova(fit_one_factor, fit_three_factor) # Able to now use anova, as it calls 'lavTestLRT' function

# Editing models to now represent the 19-item EKT, based on IRT results. Here, we have removed the negative and lowest discriminating items
# Model f1 represents 'system' knowledge, Model f2 represents 'action-related' knowledge, Model f3 represents 'effectiveness' knowledge 

# Create the three factor model, based on 19 items
three_factor_IRTEdit <- 'f1 =~ EK2 + EK4 + EK5 + EK6 + EK7 + EK9 + EK14 + EK23 + EK24 + EK25 + EK26 + EK27 
          f2 =~ EK10 + EK8 + EK15 + EK29
          f3 =~ EK12 + EK17 + EK22'

fit_three_factor_IRTEdit <- cfa(three_factor_IRTEdit, data = EK_cfa_filtered, ordered = TRUE)
summary(fit_three_factor_IRTEdit, standardized=TRUE, ci=TRUE, fit.measures=TRUE)
inspect(fit_three_factor_IRTEdit, what = "std")$lambda
modindices(fit_three_factor_IRTEdit)

# Create the one factor model, based on 19 items
one_factor_IRTEdit <- 'f1 =~ EK2 + EK4 + EK5 + EK6 + EK7 + EK8 + EK9 + EK10 + EK12 + EK14 + EK15 +EK17 +EK22 + EK23 + EK24 + EK25 + EK26 + EK27 +EK29'

fit_one_factor_IRTEdit <- cfa(one_factor_IRTEdit, data = EK_cfa_filtered, ordered = TRUE)
summary(fit_one_factor_IRTEdit, standardized=TRUE, ci=TRUE, fit.measures=TRUE)
inspect(fit_one_factor_IRTEdit, what = "std")$lambda
modindices(fit_one_factor_IRTEdit)

lavTestLRT(fit_one_factor_IRTEdit, fit_three_factor_IRTEdit, scaled.shifted = TRUE) # Compare nested models
anova(fit_three_factor_IRTEdit, fit_one_factor_IRTEdit)

# Generating reliability estimates for 30 item measure
EK_items_only  <- select(EK_cfa_filtered, EK1:EK15, EK16:EK22, EK23:EK30)
cov <- cov(EK_items_only)
glb.algebraic(cov, LoBounds = NULL, UpBounds = NULL) # Generating GLB reliability
omega(EK_items_only) # Generating Omega reliability

# Generating reliability estimates for 19 item measuree
EKT_19 <- select(EK_items_only, EK2, EK4:EK10, EK12, EK14:EK15, EK17, EK22:EK27, EK29)
cov2 <- cov(EKT_19)
glb.algebraic(cov2, LoBounds = NULL, UpBounds = NULL) # Generating GLB reliability
omega(EKT_19) # Generating Omega reliability

# Saving file
sink(file = NULL)