---
title: "Agricultural Understanding"
author: "Franz Wortha"
format: 
  html:
    self-contained: true
---

## Setup

Load packages and additional functions
```{r}
#| echo: false
library(corrplot)
library(dplyr)
library(EFA.dimensions)
library(effectsize)
library(lavaan)
library(MplusAutomation)
library(psych)
library(readxl)
library(semPlot)
library(stringr)
source("./lib.R")
```

Load data set

```{r}
df <- read.csv("./agricultural_understanding.csv")
# Filter rows where all items are missing
items <- grep("^f\\d{1,2}", names(df), value = T)
df <- df[rowSums(is.na(df[, items])) != ncol(df[, items]), ]
head(df, 5)
```

Prepare and split data
```{r}
#| echo: false
# Copy for Mplus for easier NA handling
for_mplus <- df
for_mplus[is.na(for_mplus)] <- -99

set.seed(1)
split_ratios <- c(e1 = 1/3, e2 = 1/3, c = 1/3)
# Split data into thirds
df_split <- split_data(df, split_ratios)
mplus_split <- split_data(for_mplus, split_ratios)
# Check if current split replicates reported split
split_equal <- names(mplus_split) %>%
  sapply(function(n) {
    m = mplus_split[[n]]
    dat_file = sprintf("./geoscales_recoded_%s.dat", n)
    if (file.exists(dat_file)) {
      d = read.table(dat_file)
      is_match = all(sapply(2:ncol(m), function(i) {all(m[, i] == d[, i])
      }))
      return(is_match)
    } else {
      return(FALSE)
    }
  }) %>%
  all
# Save split data sets for Mplus only if the splits are updated
if (!split_equal) {
  for (n in names(mplus_split)) {
    prepareMplusData(mplus_split[[n]], sprintf("./geoscales_recoded_%s.dat", n))
  }
  
  prepareMplusData(for_mplus, "./geoscales_recoded.dat")
}
```

```{r}
sprintf(
  "Split sizes: E1: %i, E2: %i C: %i",
  nrow(mplus_split$e1),
  nrow(mplus_split$e2),
  nrow(mplus_split$c
  ))
```

MAP
```{r}
map_e1 <- MAP(df_split$e1[, items], verbose = F)
cat(sprintf("N Factors (e1): %i %i\n", map_e1$NfactorsMAP, map_e1$NfactorsMAP4))

map_e2 <- MAP(df_split$e2[, items], verbose = F)
cat(sprintf("N Factors (e2): %i %i\n", map_e2$NfactorsMAP, map_e2$NfactorsMAP4))
```

## Scales
Items are upper-case because they are copied from Excel in that format

ED 1
```{r}
ev1 <- c(
  "F1_1",
  "F1_2A",
  "F1_2B",
  "F1_2C",
  "F1_2D",
  "F1_2E",
  "F1_2F",
  "F1_2G",
  "F1_2H",
  "F1_2I",
  "F1_2J",
  "F1_2K"
) %>% tolower
# Print names matching publication
sapply(ev1, function(n){paper_names[[n]]})
```
ED II

```{r}
ev2 <- c(
  "F1_2M",
  "F1_2N",
  "F2_9A",
  "F2_9B",
  "F2_9C",
  "F2_9D",
  "F2_9F"
) %>% tolower
# Print names matching publication
sapply(ev2, function(n){paper_names[[n]]})
```

CD
```{r}
kv <- c(
  "F2_1E",
  "F2_1G",
  "F2_1H",
  "F2_1I",
  "F2_1J",
  "F2_1K",
  "F2_1L",
  "F2_1M",
  "F2_10I",
  "F2_10J",
  "F2_10M",
  "F2_10O"
) %>% tolower
# Print names matching publication
sapply(kv, function(n){paper_names[[n]]})
```

OD
```{r}
mv1 <- c(
  "F2_10A",
  "F2_10B",
  "F2_10D",
  "F2_10E",
  "F2_10F",
  "F2_10H",
  "F2_10P"
) %>% tolower
# Print names matching publication
sapply(mv1, function(n){paper_names[[n]]})
```


```{r}
#| echo: false
# Disregarded because of low number of items
mv2 <- c("F2_10K", "F2_10L") %>% tolower
```

## Initial model fit

Fit on confirmatory subset

```{r}
scales_final <- list(ev1 = ev1,
                     ev2 = ev2,
                     kv = kv,
                     mv1 = mv1)
model_final <- cfa_from_list(scales_final)

fit_final_c <- wlmsv_cfa(model_final, df_split$c)
fitMeasures(fit_final_c, robust_indices)
```

Fit on full data set
```{r}
fit_final <- wlmsv_cfa(model_final, df)
fitMeasures(fit_final, robust_indices)
```

Internal consistency
```{r}
int_cons <- list(
  "ev1" = psych::alpha(df[, ev1]),
  "ev2" = psych::alpha(df[, ev2]),
  "kv" = psych::alpha(df[, kv]),
  "mv" = psych::alpha(df[, mv1], check.keys = T, warnings =
                        F)
)

lapply(int_cons, function(a)
  round(as.numeric(a$total[2]), 2)) %>%
  print
```

## Measurement invariance
### Weak Invariance
```{r}
scales_inv <- list(
  ev1 = ev1,
  ev2 = ev2,
  kv = kv[!kv %in% c("f2_1e", "f2_1l")],
  mv1 = mv1
)
model_inv <- cfa_from_list(scales_inv)

m1 <- cfa(
  model_inv,
  df,
  group = "geschlecht",
  ordered = T,
  estimator = "WLSMV",
  parameterization = "theta",
  missing = "pairwise"
)
# Weak invariance
m2 <- cfa(
  model_inv,
  df,
  group = "geschlecht",
  group.equal = "loadings",
  ordered = T,
  estimator = "WLSMV",
  parameterization = "theta",
  missing = "pairwise"
)

lavTestLRT(m1, m2)
```

### Partial invariance
```{r}
partial_loadings <- c(
  "ev1=~f1_2e",
  "ev1=~f1_2f",
  "ev1=~f1_2g",
  "ev1=~f1_2i"
)

m2_adj <- cfa(
  model_inv,
  df,
  group = "geschlecht",
  group.equal = c("loadings"),
  group.partial = partial_loadings,
  ordered = T,
  estimator = "WLSMV",
  parameterization = "theta",
  missing = "pairwise"
)

lavTestLRT(m1, m2_adj)
```

### Scalar Invariance
```{r}
m3 <- cfa(
  model_inv,
  df,
  group = "geschlecht",
  group.equal = c("intercepts", "loadings"),
  group.partial = partial_loadings,
  ordered = T,
  estimator = "WLSMV",
  parameterization = "theta",
  missing = "pairwise",
  meanstructure = T
)

lavTestLRT(m2_adj, m3)
```

### Establish partial invariance
```{r}
partial_thresholds <- c(
  # f1_1
  "f1_1|t1", "f1_1|t2", "f1_1|t3", "f1_1|t4", "f1_1|t5",
  # f1_2a
  "f1_2a|t1", "f1_2a|t2", "f1_2a|t3", "f1_2a|t4",
  # f1_2b
  "f1_2b|t1", "f1_2b|t2", "f1_2b|t3", "f1_2b|t4",
  # f1_2c
  "f1_2c|t1", "f1_2c|t2", "f1_2c|t3", "f1_2c|t4",
  # f1_2d
  "f1_2d|t1", "f1_2d|t2", "f1_2d|t3", "f1_2d|t4",
  # f1_2e
  "f1_2e|t1", "f1_2e|t2", "f1_2e|t3", #"f1_2e|t4",
  # f1_2f
  "f1_2f|t1", "f1_2f|t2", "f1_2f|t3", #"f1_2f|t4",
  # f1_2g
  "f1_2g|t1", "f1_2g|t2", "f1_2g|t3", #"f1_2g|t4",
  # f1_2h
  "f1_2h|t1", "f1_2h|t2", "f1_2h|t3", #"f1_2h|t4",
  # f1_2i
  "f1_2i|t1", "f1_2i|t2", "f1_2i|t3", #"f1_2i|t4",
  # f1_2j
  "f1_2j|t1", "f1_2j|t2", "f1_2j|t3", #"f1_2j|t4",
  # f1_2k
  "f1_2k|t1", "f1_2k|t2", "f1_2k|t3"#, #"f1_2k|t4"
)
##
m3_adj <- cfa(
  model_inv,
  df,
  group = "geschlecht",
  group.equal = c("thresholds", "loadings"),
  group.partial = c(partial_loadings, partial_thresholds),
  ordered = T,
  estimator = "WLSMV",
  parameterization = "theta",
  missing = "pairwise"
)

lavTestLRT(m2_adj, m3_adj)
```

Fit measures of partially invariant model
```{r}
fitMeasures(m3_adj, robust_indices)
```

## Sex comparison
```{r}
# Adjust names of df to match the paper
df_plot <- df %>%
  rename_all(function(x)
    unlist(lapply(names(.), rename_paper)))

# Add scale values from the partially invariant model to the data frame
df_full <- df_plot %>%
  rename(age = alter,
         sex = geschlecht,
         connectedness = `i6`) %>%
  rename_all(str_to_title)
lv3 <- lavPredict(m3_adj)
df_full[df_full$Sex == 1, c("ED I", "ED II", "CD", "OD")] <- lv3$`1`
df_full[df_full$Sex == 0, c("ED I", "ED II", "CD", "OD")] <- lv3$`0`
```

MANOVA
```{r}
summary(
  m_sex <- manova(cbind(`ED I`, `ED II`, CD, OD) ~ as.factor(Sex), df_full),
  test = "Pillai"
)
```

Effect size
```{r}
effectsize::eta_squared(m_sex, alternative = "two.sided")
```

Posthoc
```{r}
summary.aov(m_sex)
```

Effect sizes
```{r}
effectsize::eta_squared(aov(`ED II` ~ Sex, df_full), alternative = "two.sided")
effectsize::eta_squared(aov(`CD` ~ Sex, df_full), alternative = "two.sided")
```

Save overview table
```{r}
scales_by_sex <- df_full %>%
  mutate(Sex = recode(Sex, `0` = "Female", `1` = "Male")) %>%
  group_by(Sex) %>%
  reframe(
    across(c("ED I", "ED II", "CD", "OD"), mean_ci)
  )
# Save results 
write.csv(scales_by_sex, "./scale_means_by_sex.csv")
```

## Plots

Corplot
```{r}
#| eval: false
c_mat <- df_full %>%
  select(c(Age, Sex, `ED I`, `ED II`, CD, OD, Connectedness)) %>%
  cor(use = "pairwise.complete.obs")

p_mat <- df_full %>%
  select(c(Age, Sex, `ED I`, `ED II`, CD, OD, Connectedness)) %>%
  cor.mtest(conf.level = 0.95)

png(file = "./corplot.tiff",
    units = "in",
    width = 5,
    height = 5,
    res=300
    )

corrplot(
  c_mat,
  p.mat = p_mat$p,
  method = 'circle',
  type = 'lower',
  col = grey_white(200),
  tl.col = "black",
  cl.pos = FALSE,
  insig = 'blank',
  addCoef.col = 'black',
  number.cex = 0.8,
  diag = FALSE,
  tl.srt = 45,
  sig.level = .05
)
dev.off()
```


Model plot
```{r}
#| eval: false
plot_scales <- list(
  `ED I` = unlist(lapply(ev1, rename_paper), use.names = FALSE),
  `ED II` = unlist(lapply(ev2, rename_paper), use.names = FALSE),
  CD = unlist(lapply(kv, rename_paper), use.names = FALSE),
  OD = unlist(lapply(mv1, rename_paper), use.names = FALSE)
)
plot_model <- cfa_from_list(plot_scales)
plot_fit <- invisible(wlmsv_cfa(plot_model, df_plot))

semPlot::semPaths(
  plot_fit,
  "par",
  layout = "tree2",
  edge.color = "black",
  whatLabels = "hide",
  sizeLat = 4,
  sizeMan = 2,
  residuals = FALSE,
  measurementLayout = TRUE,
  filetype = "png",
  width = 16,
  height = 9,
  dpi = 300
)
```