setwd("D://PhD Project/3rd Analysis/Data")
getwd()

##Step 1: Install and load the necessary packages

install.packages("lavaan")
install.packages("sem")
install.packages("lme4")
install.packages("semPlot")
install.packages("semTools")
install.packages("rockchalk")

library(lavaan)
library(rockchalk)
library(semPlot)
library(semTools)
############ NEW updated packages
install.packages("openxlsx", dependencies = TRUE, repos = "https://cloud.r-project.org")
install.packages("qgraph", dependencies = TRUE)
install.packages("OpenMx", dependencies = TRUE)

##########
## Installing upgraded versions
# Install and load the remotes package (if not installed)
if (!requireNamespace("remotes", quietly = TRUE)) {
  install.packages("remotes")
}

# Update the rlang package
remotes::install_version("rlang", version = "1.1.0")

# Load the semPlot and semTools packages
library(semPlot)
library(semTools)


##Step 2: Load your dataset
meatAT <- read.csv("meat_attachment_data.csv")
View(meatAT)

##Correlation Table
round(cor(meatAT[,1:16]),2) ## for all 16 questions

##Step 3: Specify the CFA model
# Define the CFA model
model <- '
    # First-order factors
    Hedonism =~ ma3 + ma4 + ma2 + ma1
    Affinity =~ ma8 + ma5 + ma6 + ma7
    Entitlement =~ ma9 + ma10 + ma11
    Dependence =~ ma16 + ma15 + ma14 + ma12 + ma13
    
    # Second-order factor
    MeatAttachment =~ Hedonism + Affinity + Entitlement + Dependence
'

# Fit the model
fit <- cfa(model, data = meatAT)


##Step 4: Inspect the model fit
summary(fit, fit.measures = TRUE, standardized = TRUE) ## Variance Standardized Method

##Note: This will provide you with various fit indices that can help you evaluate the goodness-of-fit of your model.

########################

##Step 5: Visualize the model

semPaths(fit, "std", layout = "tree")

semPaths(
  fit, 
  style = "ram", 
  what = "std", 
  layout = "tree2",
  whatLabels = "std",
  sizeMan = 5,
  edge.label.cex = 0.8,
  label.cex = 1,
  edge.width = 0.8,
  rotation = 2
)

#######
# Define the CFA model with 14 item MTA
model <- '
    # First-order factors
    Hedonism =~ ma3 + ma4 + ma2
    Affinity =~ ma8 + ma5 + ma6 + ma7
    Entitlement =~ ma9 + ma10 + ma11
    Dependence =~ ma16 + ma15 + ma14 + ma12
    
    # Second-order factor
    MeatAttachment =~ Hedonism + Affinity + Entitlement + Dependence
'

# Fit the model
fit <- cfa(model, data = meatAT)


##Step 4: Inspect the model fit
summary(fit, fit.measures = TRUE, standardized = TRUE) ## Variance Standardized Method

##Step 5: Visualize the NEW model

semPaths(fit, "std", layout = "tree")

semPaths(
  fit, 
  style = "ram", 
  what = "std", 
  layout = "tree2",
  whatLabels = "std",
  sizeMan = 5,
  edge.label.cex = 0.8,
  label.cex = 1,
  edge.width = 0.8,
  rotation = 2
)
####To improve the quality of the semPaths
semPaths(
  fit,
  style = "ram",         # Path diagram
  what = "std",          # Standardized estimates
  layout = "tree2",      # Hierarchical tree layout
  whatLabels = "std",    # Show standardized loadings
  sizeMan = 5,           # Increase observed variable size
  sizeLat = 8,           # Latent variable size (if needed)
  edge.label.cex = 1,    # Larger edge (loading/path) labels
  label.cex = 1.2,       # Larger node labels
  edge.width = 1.2,      # Thicker edges
  rotation = 2,          # Orientation
  mar = c(5, 5, 5, 5),   # Increase margins around the plot
  residuals = TRUE,      # Show residual variances
  optimizeLatRes = FALSE  # Optimize residual placement
)


##Reliability Test
# Load the semTools package if not already loaded
if (!requireNamespace("semTools", quietly = TRUE)) {
  install.packages("semTools")
}
library(semTools)

# Extract reliability estimates
reliability_results <- reliability(fit)
# Display the results
reliability_results

# Extract correlations between subscale and global scale
inspect_results <-lavInspect(fit, what = "cor.lv", add.labels = TRUE, add.class = TRUE,
                             list.by.group = TRUE,
                             drop.list.single.group = TRUE)
# Display the results
inspect_results