### preparing script ###
set.seed(1234)
setwd(dirname(rstudioapi::getActiveDocumentContext()$path))
getwd()
options(digits = 3)

#### packages ###
#install.packages("readxl")
#install.packages("lavaan")
#install.packages("foreign")
#install.packages("semPlot")
#install.packages("psych")
#install.packages("MVN")
#install.packages("sjPlot")
#install.packages("dplyr")
#install.packages("ggplot2")
#install.packages("tidySEM")
#install.packages("semptools")
#install.packages("magrittr")
#install.packages("stringr")
#install.packages("GPArotation")
#install.packages("semTools")
#install.packages("MBESS")
#install.packages("dynamic")
#install.packages("gt")
#install.packages("tidyverse")
#install.packages("performance")
#install.packages("rstatix")

library(readxl)
library(lavaan)
library(foreign)
library(semPlot)
library(psych)
library(MVN)
library(sjPlot)
library(dplyr)
library(ggplot2)
library(tidySEM)
library(semptools)
library(magrittr)
library(stringr)
library(GPArotation)
library(semTools)
library(MBESS)
library(dynamic)
library(gt)
library(tidyverse)
library(performance)
library(rstatix)

CSAMS <- read.csv2("CSAMS_PD_24102024.csv")

### data wrangling ###
CSAMS[CSAMS == -77 | CSAMS == -88 | CSAMS == -99] <- NA
sum(is.na(CSAMS)) # no further NAs omitted in R, dataset was cleaned before importing it in R


### rename variables ###
names(CSAMS)
names(CSAMS) <- c("PD_01", "PD_02", "PD_02_t",  "PD_03", "PD_03_t",  "PD_04", "PD_04_t",  "PD_05_01", "PD_05_02",
                  "PD_05_03", "PD_05_04", "PD_05_05", "PD_05_06", "PD_05_07", "PD_05_08", "PD_05_09", "PD_05_10", 
                  "PD_05_11", "PD_05_12", "PD_05_13", "PD_05_t",  "PD_06", "PD_06_t", "PD_07", "PD_08", "PD_09", "PD_10",
                  "PD_11_01", "PD_11_02", "PD_11_03", "PD_11_04", "PD_11_05", "PD_11_t", "x1", "x2", "x3",
                  "x4", "x5", "x6", "x7", "x8", "x9", "x10", "x11", "x12", "x13", "x14", "x15", "x16", "x17", "x18", "x19")

### check for multivariate normal distribution ###
items <- CSAMS[, grep("x", names(CSAMS))]
MVN::mvn(items) 

### check for KMO-criteria ###
psych::KMO(items) 

### descriptive analysis ###
item_values <- psych::describe(items, skew = F, ranges = T)
item_values |>
  gt() |>
  tab_spanner(
    label = "item values"
  ) |>
  cols_label_with(fn = str_to_title) |>
  cols_align("left")
performance::item_discrimination(items)
tab_itemscale(items)

### confirmatory factor analysis ###
CSAMS_model <- "f1 =~ x1 + x2 + x5 + x8 + x9 + x12 + x14
f2 =~ x4 + x6 + x10 + x15
f3 =~ x3 + x7 + x11 + x13 + x17
f4 =~ x16 + x18 + x19" 

### f1 = trvialising CSA, f2 = shifting responsibility, f3 = assumptions about perpetrators, f4 = False confidence about CSA ###
CSAMS_fit <- cfa(model = CSAMS_model, 
                 data = CSAMS, 
                 estimator = "WLSMV")
summary(CSAMS_fit, standardized = T, 
        fit.measures = T,
        rsquare = T)

### calculate omega ###
MBESS::ci.reliability(items, type = "omega")
semTools::reliability(CSAMS_fit)

### estimating factor scores ###
fs <- data.frame(lavPredict(CSAMS_fit, method = "Bartlett"))
head(fs)

### adding factor scores to the CSAMS data set ###
CSAMS_fs <- cbind(CSAMS, fs)
print(CSAMS_fs)

# group differences in gender #
CSAMS_fs$PD_02[CSAMS_fs$PD_02 == 3] <- NA # exclude non-binary option due to statistical distortion  
psych::describeBy(CSAMS_fs$f1, CSAMS_fs$PD_02)
stats::wilcox.test(f1~PD_02, data = CSAMS_fs,
                   exact = F,
                   correct = F,
                   conf.int = F)
psych::describeBy(CSAMS_fs$f2, CSAMS_fs$PD_02)
stats::wilcox.test(f2~PD_02, data = CSAMS_fs,
                   exact = F,
                   correct = F,
                   conf.int = F)
psych::describeBy(CSAMS_fs$f3, CSAMS_fs$PD_02)
stats::wilcox.test(f3~PD_02, data = CSAMS_fs,
                   exact = F,
                   correct = F,
                   conf.int = F)
psych::describeBy(CSAMS_fs$f4, CSAMS_fs$PD_02)
stats::wilcox.test(f4~PD_02, data = CSAMS_fs,
                   exact = F,
                   correct = F,
                   conf.int = F)

# group differences in personal victimization #
CSAMS_fs$PD_09[CSAMS_fs$PD_09 == 3] <- NA # exclude "prefer not to answer"
psych::describeBy(CSAMS_fs$f1, CSAMS_fs$PD_09)
stats::wilcox.test(f1~PD_09, data = CSAMS_fs,
                   exact = F,
                   correct = F,
                   conf.int = F) 
psych::describeBy(CSAMS_fs$f2, CSAMS_fs$PD_09)
stats::wilcox.test(f2~PD_09, data = CSAMS_fs,
                   exact = F,
                   correct = F,
                   conf.int = F) 
psych::describeBy(CSAMS_fs$f3, CSAMS_fs$PD_09)
stats::wilcox.test(f3~PD_09, data = CSAMS_fs,
                   exact = F,
                   correct = F,
                   conf.int = F) 
psych::describeBy(CSAMS_fs$f4, CSAMS_fs$PD_09)
stats::wilcox.test(f4~PD_09, data = CSAMS_fs,
                   exact = F,
                   correct = F,
                   conf.int = F) 

### correcting for alpha error using Benjamini-Hochberg procedure ###
stats::p.adjust(c(0.003, 0.0001, 0.3, 0.002), method = "BH")
stats::p.adjust(c(.2, .2, .002, .09), method = "BH")