Issues when rendering anova code in Quarto #13511

Marie-Pedersen-PhD-student · 2025-10-07T08:33:58Z

Marie-Pedersen-PhD-student
Oct 7, 2025

Description

I keep getting the error message below if I include a simple anova code like this:

```{r}
message("No mediation vs. mediation")
anova(fit_simple, fit_mediation)
```

error message when rendering with an anova code like the one above:
Error running Lua:
YAML parse exception at line 4, column 2:
mapping values are not allowed in this context
stack traceback:
/Applications/quarto/share/pandoc/datadir/readqmd.lua:170: in function 'readqmd.readqmd'
/Applications/quarto/share/filters/qmd-reader.lua:13: in function 'Reader'

As soon as I remove anova I do not get any error. I have tried with various different anovas in different documents and got the same warning.

FYI the top of my document looks like this:

---
title: "2_3 condition study"
author: "Marie Lasrado & Seidali Kurtmollaiev"
date: last-modified
format:
  html:
    toc: true
    toc-location: left
    toc-depth: 8
    toc-title: "Table of Contents"
---

and I am utilizing these packages in my document:

```{r}
#| echo: false
#| include: false
library(lavaan)
library(lavaanPlot)
library(psych)
library(dplyr)
library(ggplot2)
library(semTools)
library(semPlot)
library(pastecs)
library(corrplot)
library(tidyr)
library(readxl)
library(car)
library(mediation)
library(knitr)
library(ARTool)
library(nnet)
library(sandwich)
library(lmtest)
library(boot)
library(MASS)
# Optional: global chunk options
knitr::opts_chunk$set(echo = FALSE, message = TRUE, warning = FALSE)
```

I hope you can help me figure out what happens and how to solve it so I also can report my anova results.

Thanks in advance!

cderv · 2025-10-07T09:43:08Z

cderv
Oct 7, 2025
Maintainer

First, I took the liberty of editing your post to format correctly. Please see the edit and look how I am including more backticks around code (https://quarto.org/bug-reports.html#formatting-make-githubs-markdown-work-for-us)

Then, I think we'll need to have a full reproducible example. If you can make it minimal that is better.

anova in itself does not have any issue. I tested with this doc

---
title: "Anova test"
date: 2022-09-12
format: html
---

```{r}
# Create a simple dataset
set.seed(123)
x1 <- rnorm(100)
x2 <- rnorm(100)
y <- 2*x1 + 3*x2 + rnorm(100)
data <- data.frame(y = y, x1 = x1, x2 = x2)

# Fit two nested models
model1 <- lm(y ~ x1, data = data)
model2 <- lm(y ~ x1 + x2, data = data)

# Compare the models using anova()
anova(model1, model2)
```

So there must be something inside your document that triggers this error. Without this document, we can't help

3 replies

Marie-Pedersen-PhD-student Oct 7, 2025
Author

Thanks. I tested your anova code in my document and got the same error. However, when testing the same code in a new document (and using the same packages) I do not get the error... So yes, it sounds like you are right about your assumption that something in the document triggers it.

From what I can see, I cannot enclose either an R code or Quarto document so here is the entire code below:

`

title: "2_3 condition study"
author: "xxx"
date: last-modified
format:
html:
toc: true
toc-location: left
toc-depth: 8
toc-title: "Table of Contents"

#| echo: false
#| include: false
library(lavaan)
library(lavaanPlot)
library(psych)
library(dplyr)
library(ggplot2)
library(semTools)
library(semPlot)
library(pastecs)
library(corrplot)
library(tidyr)
library(readxl)
library(car)
library(mediation)
library(knitr)
library(ARTool)
library(nnet)
library(sandwich)
library(lmtest)
library(boot)
library(MASS)
library(scales)
# Optional: global chunk options
knitr::opts_chunk$set(echo = FALSE, message = TRUE, warning = FALSE)

Data preparation - raw dataset

Import

#| echo: false
#| include: true
surveydata <- read.csv("reduced_2x3_survey.csv",
                       header = TRUE, 
                       sep = ",", 
                       stringsAsFactors = FALSE)

Sample size

#| echo: false
summary_list <- list(all = nrow(surveydata))
summary_list

#| echo: false
#| include: false
#| label: Converting groups to factors
surveydata$Condition_reward <- as.factor(surveydata$Condition_reward)
surveydata$Condition_eco <- as.factor(surveydata$Condition_eco)

#| echo: false
#| include: false
#| label: Calculate composite scores using the mean of the indicators
surveydata$EEF_composite <- rowMeans(surveydata[, c("EEF1", "EEF2", "EEF3")], na.rm = TRUE)
surveydata$EEC_composite <- rowMeans(surveydata[, c("EEC1", "EEC2", "EEC3")], na.rm = TRUE)
surveydata$IM_composite <- rowMeans(surveydata[, c("IM1", "IM2", "IM3")], na.rm = TRUE)  # Calculate composite scores using the mean of the indicators - motivation
surveydata$ADT_composite <- rowMeans(surveydata[, c("ADT1","ADT2","ADT3")], na.rm = TRUE)
surveydata$TR_composite <- rowMeans(surveydata[, c("TR1","TR2","TR3")], na.rm = TRUE)

Data Quality

Manipulation check and bot

Manipulation flag

#| echo: false
#| label: reward_1_control
table(surveydata$Condition_reward_name, surveydata$cond.reward_flag)
table(surveydata$Condition_eco_name, surveydata$cond.eco_flag)

Bot flag

#| echo: false
table(surveydata$bot_flag)

Attention

Duration

#| echo: false
summary(surveydata$duration_minutes)
hist(surveydata$duration_minutes)

message("Outliers defined as 3 std. deviations below or above the mean")
mean_completiontime <- mean(surveydata$duration_minutes, na.rm = TRUE)
std_completiontime <- sd(surveydata$duration_minutes, na.rm = TRUE)
lower_threshold <- max(0, mean_completiontime - 3 * std_completiontime)
upper_threshold <- mean_completiontime + (3 * std_completiontime)
surveydata$outliers_completion <- surveydata$duration_minutes < lower_threshold | surveydata$duration_minutes > upper_threshold
message("Outliers on completion time")
table(surveydata$outliers_completion)

On scales

#| echo: false
hist(surveydata$StandardDev) # Standard deviation
summary(surveydata$StandardDev)

message("Flagged outliers based on scales")
columns_to_include <- c('PEB','IM1','IM2','EEF1','EEF2','EEF3','EEC1','EEC2','EEC3','ADT1','ADT2','ADT3')
surveydata$ScaleMean <- rowMeans(
  surveydata[, columns_to_include],
  na.rm = TRUE
)
# Overall mean & SD of that scale mean
overall_mean <- mean(surveydata$ScaleMean, na.rm = TRUE)
overall_sd   <- sd(surveydata$ScaleMean, na.rm = TRUE)
lower_cut <- overall_mean - 3 * overall_sd
upper_cut <- overall_mean + 3 * overall_sd
surveydata$outliers_scales <-
  surveydata$ScaleMean < lower_cut |
  surveydata$ScaleMean > upper_cut
table(surveydata$outliers_scales)

Removing bad participants

Exclude participants

excluded <- subset(surveydata, cond.reward_flag == TRUE | cond.eco_flag == TRUE | bot_flag == TRUE |  outliers_completion == TRUE | outliers_scales == TRUE)

#| echo: false
surveydata %>%
  # Keep only rows where at least one flag is TRUE
  filter(cond.reward_flag | cond.eco_flag |  bot_flag | 
         outliers_completion | outliers_scales) %>%
  # Count combinations of flags
  count(cond.reward_flag,
        cond.eco_flag,
        outliers_completion,
        bot_flag,
        outliers_scales) %>%
  # Optional: sort by a specific flag (TRUE first)
  arrange(desc(cond.reward_flag)) %>%
  # Pretty print
  kable()

data_filtered <- surveydata %>%
  filter(!Prolific.id %in% excluded$Prolific.id)

message("total excluded from the dataset")
summary_list <- list(all = nrow(excluded))
summary_list

#| echo: false
#| include: false
#| label: Export to Lester
data_filtered_export_composites <- data_filtered %>%
  dplyr::select(Prolific.id, Condition_reward, Condition_reward_name, Condition_eco, Condition_eco_name, EEF_composite, EEC_composite,ADT_composite, TR_composite)
write.csv(data_filtered_export_composites, "2x3 study_filtered_composites.csv", row.names = FALSE)

Filtered dataset

Descriptive on good participants

Conditions

Group statistics

message("Reward groups")
table(data_filtered$Condition_reward_name)

message("Eco orientation groups")
table(data_filtered$Condition_eco_name)

#| echo: false
data_filtered %>%
  group_by(Condition_reward_name,Condition_eco_name) %>%
  summarise(
    n         = n(),
    mean_EEF  = mean(EEF_composite, na.rm = TRUE),
    sd_EEF    = sd(EEF_composite, na.rm = TRUE),
    mean_EEC  = mean(EEC_composite, na.rm = TRUE),
    sd_EEC    = sd(EEC_composite, na.rm = TRUE),
    mean_IM  = mean(IM_composite, na.rm = TRUE),
    sd_IM    = sd(IM_composite, na.rm = TRUE),    
    .groups   = "drop"
  )

# 1) summary
summary_df <- data_filtered %>%
  group_by(Condition_reward_name) %>%
  summarise(
    n        = n(),
    mean_EEF = mean(EEF_composite, na.rm = TRUE),
    sd_EEF   = sd(EEF_composite, na.rm = TRUE),
    mean_EEC = mean(EEC_composite, na.rm = TRUE),
    sd_EEC   = sd(EEC_composite, na.rm = TRUE),
    mean_IM  = mean(IM_composite, na.rm = TRUE),
    sd_IM    = sd(IM_composite, na.rm = TRUE),   
    .groups  = "drop"
  )

# 2) reshape so mean & sd stay together
summary_long <- summary_df %>%
  pivot_longer(
    cols = c(starts_with("mean_"), starts_with("sd_")),
    names_to = c("stat", "measure"),
    names_pattern = "^(mean|sd)_(.+)$"
  ) %>%
  pivot_wider(
    names_from = stat,
    values_from = value
  ) %>%
  mutate(
    measure = factor(measure, levels = c("EEF", "EEC","IM")),
    se = sd / sqrt(n)   # optional: standard error
  )

pd2 <- position_dodge(width = 0.4)

make_summary <- function(data, group_var) {
  data %>%
    group_by({{ group_var }}) %>%
    summarise(
      n        = n(),
      mean_EEF = mean(EEF_composite, na.rm = TRUE),
      sd_EEF   = sd(EEF_composite, na.rm = TRUE),
      mean_EEC = mean(EEC_composite, na.rm = TRUE),
      sd_EEC   = sd(EEC_composite, na.rm = TRUE),
      mean_IM  = mean(IM_composite, na.rm = TRUE),
      sd_IM    = sd(IM_composite, na.rm = TRUE),
      .groups  = "drop"
    ) %>%
    pivot_longer(
      cols = c(starts_with("mean_"), starts_with("sd_")),
      names_to = c("stat", "measure"),
      names_pattern = "^(mean|sd)_(.+)$"
    ) %>%
    pivot_wider(names_from = stat, values_from = value) %>%
    mutate(
      measure = factor(measure, levels = c("EEF", "EEC","IM")),
      se = sd / sqrt(n)
    )
}

# --- 2) Create summary tables ----------------------------------
summary_reward <- make_summary(data_filtered, Condition_reward_name)
summary_eco    <- make_summary(data_filtered, Condition_eco_name)

# --- 3) Plotting helper ----------------------------------------
plot_summary <- function(df, x_var, title) {
  ggplot(df, aes(x = {{ x_var }}, y = mean,
                 color = measure, group = measure)) +
    geom_point(position = position_dodge(width = 0.4), size = 3) +
    geom_errorbar(
      aes(ymin = mean - sd, ymax = mean + sd),
      position = position_dodge(width = 0.4),
      width = 0.15
    ) +
    labs(x = rlang::as_name(rlang::enquo(x_var)),
         y = "Mean (± SD)",
         color = "Measure",
         title = title) +
    theme_minimal()
}

# --- 4) Draw plots --------------------------------------------
plot_summary(summary_reward,
             Condition_reward_name,
             "Mean EEF, EEC & IM by reward conditions (points ± SD)")

plot_summary(summary_eco,
             Condition_eco_name,
             "Mean EEF, EEC & IM by eco conditions (points ± SD)")

#| echo: false
df_long <- data_filtered %>%
  pivot_longer(cols = c(EEF_composite, EEC_composite, IM_composite),
               names_to = "Measure",
               values_to = "Score")

ggplot(df_long, aes(x = Condition_reward_name, y = Score, fill = Condition_reward_name)) +
  geom_boxplot(alpha = 0.6, outlier.shape = NA) +
  geom_jitter(width = 0.2, alpha = 0.5) +
  facet_wrap(~Measure, scales = "free_y") +
  theme_minimal() +
  labs(title = "Descriptive statistics of EEF & EEC by Condition",
       y = "Score")

#| echo: false
df_long <- data_filtered %>%
  pivot_longer(cols = c(EEF_composite, EEC_composite, IM_composite),
               names_to = "Measure",
               values_to = "Score")

ggplot(df_long, aes(x = Condition_eco_name, y = Score, fill = Condition_eco_name)) +
  geom_boxplot(alpha = 0.6, outlier.shape = NA) +
  geom_jitter(width = 0.2, alpha = 0.5) +
  facet_wrap(~Measure, scales = "free_y") +
  theme_minimal() +
  labs(title = "Descriptive statistics of EEF & EEC by eco condition",
       y = "Score")

Ease and feedback

#| echo: false
message("Q: Were there anything that you found confusing or difficult to follow/answer?")
hist(data_filtered$ease)
describe(data_filtered$ease)

Scales

#| echo: false
#| include: false
scales_data_all <- subset(surveydata, select = c(Condition_reward,Condition_eco,IM1,IM2,IM3,PEB,EEF1,EEF2,EEF3,EEC1,EEC2,EEC3,ADT1,ADT2,ADT3,TR1,TR2,TR3))
scales_data_good <- subset(data_filtered, select = c(Condition_reward,Condition_eco,IM1,IM2,IM3,PEB,EEF1,EEF2,EEF3,EEC1,EEC2,EEC3,ADT1,ADT2,ADT3,TR1,TR2,TR3))
efa_data_good <- subset(data_filtered, select = c(IM1,IM2,IM3,EEF1,EEF2,EEF3,EEC1,EEC2,EEC3,ADT1,ADT2,ADT3,TR1, TR2,TR3))
write.csv(efa_data_good, "PLS_filtered.csv", row.names = FALSE)

Descriptive stats on scales

#| echo: false
c(
  all  = nrow(surveydata),
  good = nrow(data_filtered)) #Sample size
descriptive_good <- describe(efa_data_good)
print(descriptive_good)

Descriptive stats on scales

#| echo: false
c(
  all  = nrow(surveydata),
  good = nrow(data_filtered)) #Sample size
descriptive_good <- describe(efa_data_good)
descriptive_good_df <- data.frame(Variable = rownames(descriptive_good), descriptive_good, row.names = NULL)
descriptive_good <- cbind(Variable = rownames(descriptive_good), descriptive_good)
print(descriptive_good)

#| echo: false
vars_to_plot <- c("IM1","IM2","IM3",
                  "EEF1","EEF2","EEF3",
                  "EEC1","EEC2","EEC3",
                  "ADT1","ADT2","ADT3",
                  "TR1","TR2","TR3")

data_long <- efa_data_good %>%
  pivot_longer(
    cols = c(IM1, IM2, IM3,
             EEF1, EEF2, EEF3,
             EEC1, EEC2, EEC3,
             ADT1, ADT2, ADT3,
             TR1,TR2,TR3),
    names_to = "variable",
    values_to = "score"
  )

ggplot(data_long, aes(x = variable, y = score)) +
  geom_boxplot(fill = "lightgreen") +
  labs(title = "Boxplots of Variables", x = "Variable", y = "Score") +
  theme_minimal() +
  theme(axis.text.x = element_text(angle = 45, hjust = 1))

Normality tests

Non-normality test across all scales

#| echo: false
#| label: Shapiro test
apply(efa_data_good, 2, shapiro.test)

Non-normality test on EEF composite

#| echo: false
hist(data_filtered$EEF_composite, main = "Histogram of EEF composite")
boxplot(data_filtered$EEF_composite, main = "Boxplot of EEF composite")
qqnorm(data_filtered$EEF_composite)
qqline(data_filtered$EEF_composite, col = "red")
shapiro.test(data_filtered$EEF_composite)

Non-normality test on EEC composite

#| echo: false
hist(data_filtered$EEC_composite, main = "Histogram of EEC composite")
boxplot(data_filtered$EEC_composite, main = "Boxplot of EEC composite")
qqnorm(data_filtered$EEC_composite)
qqline(data_filtered$EEC_composite, col = "red")
shapiro.test(data_filtered$EEC_composite)

Non-normality test on DV per condition

message("Reward condition on EEF")
data_filtered %>%
  group_by(Condition_reward_name) %>%
  summarise(
    n = n(),
    shapiro_p = shapiro.test(EEF_composite)$p.value
  )
message("Reward condition on EEC")
data_filtered %>%
  group_by(Condition_reward_name) %>%
  summarise(
    n = n(),
    shapiro_p = shapiro.test(EEC_composite)$p.value
  )
message("Eco condition on EEF")
data_filtered %>%
  group_by(Condition_eco_name) %>%
  summarise(
    n = n(),
    shapiro_p = shapiro.test(EEF_composite)$p.value
  )
message("Eco condition on EEC")
data_filtered %>%
  group_by(Condition_eco_name) %>%
  summarise(
    n = n(),
    shapiro_p = shapiro.test(EEC_composite)$p.value
  )

Non-normality test on IM composite

#| echo: false
hist(data_filtered$IM_composite, main = "Histogram of IM composite")
boxplot(data_filtered$IM_composite, main = "Boxplot of IM composite")
qqnorm(data_filtered$IM_composite)
qqline(data_filtered$IM_composite, col = "red")
shapiro.test(data_filtered$IM_composite)

Non-normality test on IM per condition

#| echo: false
data_filtered %>%
  group_by(Condition_reward_name) %>%
  summarise(
    n = n(),
    shapiro_p = shapiro.test(IM_composite)$p.value
  )
  
data_filtered %>%
  group_by(Condition_eco_name) %>%
  summarise(
    n = n(),
    shapiro_p = shapiro.test(IM_composite)$p.value
  )

Factor analyses

KMO

#| echo: false
kmo_result1<- KMO(efa_data_good)
print(kmo_result1)

Correlation analysis

Bartlett test

#| echo: false
bartlett_result1<-cortest.bartlett(efa_data_good)
print(bartlett_result1)

Correlation matrix

#| echo: false
cor_matrix <- cor(efa_data_good, use = "complete.obs")
print(cor_matrix)
corrplot(cor_matrix, method = "circle")

EFA

#| echo: false
scree(efa_data_good) #Screen plot
fa.parallel(efa_data_good, fa = "fa") # Suggest number of factors

efa_result1<-fa(efa_data_good, nfactors = 5, rotate = "varimax") # Run EFA with the suggested number of factors
message("Threshold=0.35")
threshold <- 0.35
loadings_matrix <- as.matrix(efa_result1$loadings)
loadings_display <- loadings_matrix
loadings_display[abs(loadings_display) < threshold] <- NA
print(loadings_display, digits = 3, na.print = "")

fa(r = cor(efa_data_good, use = "pairwise.complete.obs"), 
   nfactors = 5, rotate = "varimax") # Correlation matrix

fa.diagram(efa_result1)

# Run EFA (oblique rotation to allow factor correlations)
efa_result1 <- fa(efa_data_good, nfactors = 5, rotate = "oblimin")

# Display loadings above threshold
message("Threshold = 0.35")
threshold <- 0.35
loadings_matrix <- as.matrix(efa_result1$loadings)
loadings_display <- loadings_matrix
loadings_display[abs(loadings_display) < threshold] <- NA
print(loadings_display, digits = 3, na.print = "")

# Extract the loadings matrix safely
load_df <- as.data.frame(unclass(efa_result1$loadings))

# Assign each item to the factor with the highest absolute loading
assignments <- apply(load_df, 1, function(x) which.max(abs(x)))

# Automatically name constructs based on item names
construct_names <- tapply(rownames(load_df), assignments, paste, collapse = "_")

# Extract factor correlation matrix (only available for oblique rotations)
phi <- efa_result1$Phi

# Assign generated construct names
rownames(phi) <- colnames(phi) <- names(construct_names)

# Display SPSS-style factor correlation matrix
cat("\nFactor Correlation Matrix:\n")
print(round(phi, 3))

# Optional: visualize factor structure
fa.diagram(efa_result1)

#| echo: false
efa_fit <- fa(efa_data_good, nfactors = 5, fm = "ml", rotate = "oblimin")


# Step 1: extract loadings
loadings <- as.data.frame(efa_fit$loadings[1:ncol(efa_data_good), ])

# Step 2: assign each item to the factor with the highest absolute loading
assignments <- apply(loadings, 1, function(x) which.max(abs(x)))

# Step 3: make factor names automatically (based on items)
factor_names <- tapply(rownames(loadings), assignments, paste, collapse = "_")

# Step 4: get factor correlation matrix
phi <- efa_fit$Phi

# Step 5: rename rows/columns using generated labels
rownames(phi) <- colnames(phi) <- names(factor_names)

# Show result
round(phi, 3)

Common method bias

Harmans's test

#| echo: false
harman_test1 <- fa(efa_data_good, nfactors = 1, rotate = "none")
print(harman_test1)

CFA with one factor

model_substantive <- '
  EEC =~ EEC1 + EEC2 + EEC3
  EEF =~ EEF1 + EEF2 + EEF3
  ADT =~ ADT1 + ADT2 + ADT3
  IM  =~ IM1 + IM2 + IM3
  TR =~ TR1 + TR2 + TR3
'

# Fit the CFA model
fit_substantive <- cfa(model_substantive, data = efa_data_good, estimator = "MLR")

# Summary with fit measures and standardized estimates
summary(fit_substantive, fit.measures = TRUE, standardized = TRUE)

# Inspect latent correlations
lavInspect(fit_substantive, "cor.lv")

# Inspect latent variances (diagonal of covariance matrix)
cov_lv <- lavInspect(fit_substantive, "cov.lv")
diag(cov_lv)


message("CFA with a marker-variable CMF")

model_marker <- '
  # Substantive factors
  EEC =~ EEC1 + EEC2 + EEC3
  EEF =~ EEF1 + EEF2 + EEF3
  ADT =~ ADT1 + ADT2 + ADT3
  TR =~ TR1 + TR2 + TR3
  IM  =~ IM1 + IM2 + IM3

  # Common method factor (marker variable approach)
  CMF =~ EEC1 + EEF1 + ADT1 + TR1 + IM1

  # Fix CMF variance for identification
  CMF ~~ 1*CMF
'

# Fit the model
fit_marker <- cfa(model_marker, data = efa_data_good, estimator = "MLR")

# Summary with fit measures and standardized estimates
summary(fit_marker, fit.measures = TRUE, rsquare=TRUE, standardized = TRUE)

# Inspect latent correlations
lavInspect(fit_marker, "cor.lv")

# Inspect latent variances
diag(lavInspect(fit_marker, "cov.lv"))

CFA

model1<-'
  EEC =~ EEC1 + EEC2 + EEC3
  EEF =~ EEF1 + EEF2 + EEF3
  ADT =~ ADT1 + ADT2 + ADT3
  IM =~ IM1 + IM2 + IM3
  TR =~ TR1 + TR2 + TR3'
fit1<-cfa(model1,data = efa_data_good, estimator = "MLR")
summary(fit1, ci = TRUE, fit.measures=TRUE, rsquare=TRUE,standardized=TRUE)
lavaanPlot(model = fit1, node_options = list(shape = "box", fontname = "Helvetica"), edge_options = list(color = "grey"),
           coefs = T, stand=T)
#Construct Reliability:
message("Cronbach’s Alpha:")
semTools::compRelSEM(fit1, tau.eq=T, obs.var=T) #Cronbach Alpha
message("Omega:")
semTools::compRelSEM(fit1, tau.eq=F, obs.var=T) #Omega
message("Average Variance Extracted (AVE):")
AVE(fit1, obs.var = TRUE, omit.imps = c("no.conv", "no.se"),
    omit.factors = character(0), dropSingle = TRUE, return.df = TRUE) #AVE
#Standardized residuals
lavResiduals(fit1, type = "cor.bentler", custom.rmr = NULL,
             se = FALSE, zstat = TRUE, summary = TRUE, h1.acov = "unstructured",
             add.type = TRUE, add.labels = TRUE, add.class = TRUE,
             drop.list.single.group = TRUE,
             maximum.number = length(res.vech), output = "list")

SEM

#| echo: false
#| include: true
# Create combined factor and recode it with friendly names
data_filtered <- data_filtered %>%
  mutate(
    Condition_combined = interaction(Condition_reward, Condition_eco, sep = "_")
  ) %>%
  mutate(
    Condition_combined = recode_factor(
      Condition_combined,
      "1_1" = "reward0_eco1",
      "2_1" = "reward1_eco1",
      "1_2" = "reward0_eco2",
      "2_2" = "reward1_eco2",
      "1_3" = "reward0_eco3",
      "2_3" = "reward1_eco3"
    )
  )


# Convert combined factor into dummy variables
dummy_combined <- model.matrix(~ Condition_combined - 1,
                               data = data_filtered,
                               na.action = na.pass)

# Clean up column names: remove "Condition_combined" prefix
colnames(dummy_combined) <- gsub("^Condition_combined", "", colnames(dummy_combined))

# Combine with original data
data_with_dummies <- cbind(data_filtered, dummy_combined)

No mediation

model_simple <- '
# Measurement model
EEF =~ EEF1 + EEF2 + EEF3
EEC =~ EEC1 + EEC2 + EEC3
IM  =~ IM1 + IM2 + IM3

# Structural model
  EEF ~ reward1_eco1 + reward0_eco2 + reward1_eco2 + reward0_eco3 + reward1_eco3
  EEC ~ reward1_eco1 + reward0_eco2 + reward1_eco2 + reward0_eco3 + reward1_eco3
'
fit_simple <- sem(model_simple, data = data_with_dummies, estimator = "MLR")
summary(fit_simple, fit.measures=TRUE, rsquare=TRUE, standardized=TRUE)
lavaanPlot(model = fit_simple, node_options = list(shape = "box", fontname = "Helvetica"), edge_options = list(color = "grey"),
           coefs = T, stand=T)

Only EEF

model_simple_EEF <- '
# Measurement model
EEF =~ EEF1 + EEF2 + EEF3

# Structural model
  EEF ~ reward1_eco1 + reward0_eco2 + reward1_eco2 + reward0_eco3 + reward1_eco3
'
fit_simple_EEF <- sem(model_simple_EEF, data = data_with_dummies, estimator = "MLR")
summary(fit_simple_EEF, fit.measures=TRUE, rsquare=TRUE, standardized=TRUE)
lavaanPlot(model = fit_simple_EEF, node_options = list(shape = "box", fontname = "Helvetica"), edge_options = list(color = "grey"),
           coefs = T, stand=T)

Only EEC

model_simple_EEC <- '
# Measurement model
EEC =~ EEC1 + EEC2 + EEC3

# Structural model
  EEC ~ reward1_eco1 + reward0_eco2 + reward1_eco2 + reward0_eco3 + reward1_eco3
'
fit_simple_EEC <- sem(model_simple_EEC, data = data_with_dummies, estimator = "MLR")
summary(fit_simple_EEC, fit.measures=TRUE, rsquare=TRUE, standardized=TRUE)
lavaanPlot(model = fit_simple_EEC, node_options = list(shape = "box", fontname = "Helvetica"), edge_options = list(color = "grey"),
           coefs = T, stand=T)

Mediation

Non-connected DVS

message("Simple mediation on both dependent variables")
model_mediation <- '
# Measurement model
EEC =~ EEC1 + EEC2 + EEC3
EEF =~ EEF1 + EEF2 + EEF3
IM =~ IM1 + IM2 + IM3

# Structural model
  IM  ~  reward1_eco1 + reward0_eco2 + reward1_eco2 + reward0_eco3 + reward1_eco3
  EEF ~ IM + reward1_eco1 + reward0_eco2 + reward1_eco2 + reward0_eco3 + reward1_eco3
  EEC ~ IM + reward1_eco1 + reward0_eco2 + reward1_eco2 + reward0_eco3 + reward1_eco3
'
fit_mediation <- sem(model_mediation, data = data_with_dummies, estimator = "MLR")
summary(fit_mediation, fit.measures=TRUE, rsquare=TRUE, standardized=TRUE)
lavaanPlot(model = fit_mediation, node_options = list(shape = "box", fontname = "Helvetica"), edge_options = list(color = "grey"),
           coefs = T, stand=T)

Also partial mediation

message("Non-connected DVs")
model2 <- '
# Measurement model
EEC =~ EEC1 + EEC2 + EEC3
EEF =~ EEF1 + EEF2 + EEF3
IM =~ IM1 + IM2 + IM3

# Structural model
  IM  ~  reward1_eco1 + reward0_eco2 + reward1_eco2 + reward0_eco3 + reward1_eco3
  EEF ~ IM + reward1_eco1 + reward0_eco2 + reward1_eco2 + reward0_eco3 + reward1_eco3
  EEC ~ IM + reward1_eco1 + reward0_eco2 + reward1_eco2 + reward0_eco3 + reward1_eco3
'
fit2 <- sem(model2, data = data_with_dummies, estimator = "MLR")
summary(fit2, fit.measures=TRUE, rsquare=TRUE, standardized=TRUE)
lavaanPlot(model = fit2, node_options = list(shape = "box", fontname = "Helvetica"), edge_options = list(color = "grey"),
           coefs = T, stand=T)

Only EEF

model2d <- '
# Measurement model
EEF =~ EEF1 + EEF2 + EEF3
IM =~ IM1 + IM2 + IM3

# Structural model
  IM  ~  reward1_eco1 + reward0_eco2 + reward1_eco2 + reward0_eco3 + reward1_eco3
  EEF ~ IM
'
fit2d <- sem(model2d, data = data_with_dummies, estimator = "MLR")
summary(fit2d, fit.measures=TRUE, rsquare=TRUE, standardized=TRUE)
lavaanPlot(model = fit2d, node_options = list(shape = "box", fontname = "Helvetica"), edge_options = list(color = "grey"),
           coefs = T, stand=T)

Only EEC

model2f <- '
# Measurement model
EEC =~ EEC1 + EEC2 + EEC3
IM =~ IM1 + IM2 + IM3

# Structural model
  IM  ~  reward1_eco1 + reward0_eco2 + reward1_eco2 + reward0_eco3 + reward1_eco3
  EEC ~ IM
'
fit2f <- sem(model2f, data = data_with_dummies, estimator = "MLR")
summary(fit2f, fit.measures=TRUE, rsquare=TRUE, standardized=TRUE)
lavaanPlot(model = fit2f, node_options = list(shape = "box", fontname = "Helvetica"), edge_options = list(color = "grey"),
           coefs = T, stand=T)

Partial mediation EEF->>EEC

message("Partial mediation EEF --> EEC")
model3 <- '
# Measurement model
EEC =~ EEC1 + EEC2 + EEC3
EEF =~ EEF1 + EEF2 + EEF3
IM =~ IM1 + IM2 + IM3

# Structural model
  IM  ~ reward1_eco1 + reward0_eco2 + reward1_eco2 + reward0_eco3 + reward1_eco3
  EEF ~ IM + reward1_eco1 + reward0_eco2 + reward1_eco2 + reward0_eco3 + reward1_eco3
  EEC ~ IM + reward1_eco1 + reward0_eco2 + reward1_eco2 + reward0_eco3 + reward1_eco3 + EEF
'
fit3 <- sem(model3, data = data_with_dummies, estimator = "MLR")
summary(fit3, fit.measures=TRUE, rsquare=TRUE, standardized=TRUE)
lavaanPlot(model = fit3, node_options = list(shape = "box", fontname = "Helvetica"), edge_options = list(color = "grey"),
           coefs = T, stand=T)

Partial EEC->>EEF

message("Partial mediation EEC --> EEF")
model3b <- '
# Measurement model
EEC =~ EEC1 + EEC2 + EEC3
EEF =~ EEF1 + EEF2 + EEF3
IM =~ IM1 + IM2 + IM3

# Structural model
  IM  ~ reward1_eco1 + reward0_eco2 + reward1_eco2 + reward0_eco3 + reward1_eco3
  EEF ~ IM + reward1_eco1 + reward0_eco2 + reward1_eco2 + reward0_eco3 + reward1_eco3 + EEC
  EEC ~ IM + reward1_eco1 + reward0_eco2 + reward1_eco2 + reward0_eco3 + reward1_eco3
'
fit3b <- sem(model3b, data = data_with_dummies, estimator = "MLR")
summary(fit3b, fit.measures=TRUE, rsquare=TRUE, standardized=TRUE)
lavaanPlot(model = fit3b, node_options = list(shape = "box", fontname = "Helvetica"), edge_options = list(color = "grey"),
           coefs = T, stand=T)

IM

message("IM as DV and only manipulations as IVs")
model2d <- '
# Measurement model
IM =~ IM1 + IM2 + IM3

# Structural model
  IM  ~  reward1_eco1 + reward0_eco2 + reward1_eco2 + reward0_eco3 + reward1_eco3
'
fit2d <- sem(model2d, data = data_with_dummies, estimator = "MLR")
summary(fit2d, fit.measures=TRUE, rsquare=TRUE, standardized=TRUE)
lavaanPlot(model = fit2d, node_options = list(shape = "box", fontname = "Helvetica"), edge_options = list(color = "grey"),
           coefs = T, stand=T)

message("IM without independent variables")
model4 <- '
  # Measurement model
  EEC =~ EEC1 + EEC2 + EEC3
  EEF =~ EEF1 + EEF2 + EEF3
  IM =~ IM1 + IM2 + IM3

  # Structural model
  EEF ~ IM
  EEC ~ IM
'
fit4 <- sem(model4, data = data_with_dummies, estimator= "MLR")
summary(fit4, fit.measures=TRUE, standardized=TRUE, rsquare=TRUE)
lavaanPlot(model = fit4, node_options = list(shape = "box", fontname = "Helvetica"), edge_options = list(color = "grey"),
           coefs = T, stand=T)

ADT and full mediation by IM

model5 <- '
  # Measurement model
  EEC =~ EEC1 + EEC2 + EEC3
  EEF =~ EEF1 + EEF2 + EEF3
  IM =~ IM1 + IM2 + IM3
  ADT =~ ADT1 + ADT2 + ADT3

  # Structural model
  EEF ~ IM + EEC
  EEC ~ IM + ADT
  IM ~ ADT
'
fit5 <- sem(model5, data = data_with_dummies, estimator= "MLR")
summary(fit5, fit.measures=TRUE, standardized=TRUE, rsquare=TRUE)
lavaanPlot(model = fit5, node_options = list(shape = "box", fontname = "Helvetica"), edge_options = list(color = "grey"),
           coefs = T, stand=T)

Moderation

ADT

Continuous ADT

On reward groups

model6 <- '
  # Measurement model
  EEC =~ EEC1 + EEC2 + EEC3
  EEF =~ EEF1 + EEF2 + EEF3
  IM =~ IM1 + IM2 + IM3
  ADT =~ ADT1 + ADT2 + ADT3

  # Structural model
  EEF ~ IM + ADT
  EEC ~ IM + ADT
  IM ~ ADT
  '
fit6 <- sem(model6, data = data_with_dummies, group = "Condition_reward_name", estimator= "MLR")
summary(fit6, fit.measures=TRUE, standardized=TRUE)

On eco groups

#| echo: false
model7 <- '
  # Measurement model
  EEC =~ EEC1 + EEC2 + EEC3
  EEF =~ EEF1 + EEF2 + EEF3
  IM =~ IM1 + IM2 + IM3
  ADT =~ ADT1 + ADT2 + ADT3

  # Structural model
  EEF ~ IM + EEC + ADT
  EEC ~ IM + ADT
  IM ~ ADT
  '
fit7 <- sem(model7, data = data_with_dummies, group = "Condition_eco_name", estimator= "MLR")
summary(fit7, fit.measures=TRUE, standardized=TRUE)

GG plot

#| echo: false
# Rename eco-conditions
data_with_dummies <- data_with_dummies %>%
  mutate(Condition_eco_name = recode_factor(
    Condition_eco_name,
    "EEC_ori" = "EEC orientation",
    "EEF_ori" = "EEF orientation",
    "both_ori" = "Combination of EEF and EEC"
  ))
# Plot with renamed facets
ggplot(data_with_dummies, aes(x = IM_composite, y = EEF_composite)) +
  geom_point(alpha = 0.8) +
  geom_smooth(method = "loess", se = TRUE) +
  facet_wrap(~Condition_eco_name) +
  theme_minimal() +
  labs(x = "IM Composite", y = "EEF Composite")

#| echo: false
ggplot(data_with_dummies, aes(x = IM_composite, y = EEC_composite)) +
  geom_point(alpha = 0.8) +
  geom_smooth(method = "loess", se = TRUE) +
  facet_wrap(~Condition_eco_name) +
  theme_minimal() +
  labs(x = "IM Composite", y = "EEC Composite")

ADT as categorical variable

#| echo: false
# Dummies for ADT

# Create dummy for high ADT (>5)
data_with_dummies$ADT_high <- ifelse(data_with_dummies$ADT_composite > 5, 1, 0)

data_with_dummies$ADT_high <- factor(
  data_with_dummies$ADT_high,
  levels = c(0, 1),                 # 0 = low, 1 = high
  labels = c("ADT_low", "ADT_high")
)

SEM

####### Grouped by reward

message("Note: High ADT is coded as above 5")
data_with_dummies$ADT_high <- as.factor(data_with_dummies$ADT_high)

model8 <- '
  # Measurement model
  EEC =~ EEC1 + EEC2 + EEC3
  EEF =~ EEF1 + EEF2 + EEF3
  IM =~ IM1 + IM2 + IM3

  # Structural model
  EEF ~ IM + EEC + ADT_high
  EEC ~ IM + ADT_high
  IM ~ ADT_high
  '
fit8 <- sem(model8, data = data_with_dummies, group = "Condition_reward_name", estimator= "MLR")
summary(fit8, fit.measures=TRUE, standardized=TRUE)

Grouped by eco condition

#| echo: false
data_with_dummies_EEF_EEC_subset <- subset(
  data_with_dummies,
  Condition_eco == 1 | Condition_eco == 2
)


model9 <- '
  # Measurement model
  EEC =~ EEC1 + EEC2 + EEC3
  EEF =~ EEF1 + EEF2 + EEF3
  IM =~ IM1 + IM2 + IM3

  # Structural model
  EEF ~ IM + EEC + ADT_high
  EEC ~ IM + ADT_high
  IM ~ ADT_high
  '
fit9 <- sem(model9, data = data_with_dummies_EEF_EEC_subset, group = "Condition_eco_name", estimator= "MLR")
summary(fit9, fit.measures=TRUE, standardized=TRUE)

GGplot

####### Reward

#| echo: false
# Make ADT_high a factor for grouping
data_with_dummies <- data_with_dummies %>%
  mutate(
    ADT_high_plot = recode_factor(
      ADT_high,
      "0" = "Low",
      "1" = "High",
      "ADT_low" = "Low ADT",
      "ADT_high" = "High ADT"
    ),
    # Rename reward conditions
    Condition_reward_name = recode_factor(
      Condition_reward_name,
      "no_reward" = "Control",
      "performance_reward" = "Performance-based reward"
    )
  )

# Reshape data to long format
data_long <- data_with_dummies %>%
  pivot_longer(
    cols = c(EEF_composite, EEC_composite),
    names_to = "Variable",
    values_to = "Value"
  )

# Compute mean per group
summary_df <- data_long %>%
  group_by(Variable, ADT_high_plot, Condition_reward_name) %>%
  summarise(Mean = mean(Value, na.rm = TRUE), .groups = "drop")

summary_df$Variable <- factor(summary_df$Variable, levels = c("EEF_composite", "EEC_composite"))

variable_labels <- c(
  EEF_composite = "EEF",
  EEC_composite = "EEC"
)

# Plot intersecting lines side by side
ggplot(summary_df, aes(
  x = ADT_high_plot,
  y = Mean,
  group = Condition_reward_name,
  color = Condition_reward_name
)) +
  geom_line(size = 1.2, na.rm = TRUE) +
  geom_point(size = 3, na.rm = TRUE) +
  theme_minimal() +
  labs(x = NULL, y = "Pro-environmental behaviour", color = "Reward conditions") +
  scale_color_brewer(palette = "Set1") +
  theme(
    strip.text = element_text(size = 12, face = "bold"),
    legend.background = element_rect(fill = alpha("white", 0.7), color = NA),
    legend.title = element_text(size = 10, face = "bold"),
    legend.text = element_text(size = 9)
  ) +
  facet_wrap(~Variable, labeller = labeller(Variable = variable_labels)) +
  coord_cartesian(ylim = c(3, 6))

####### Eco-condition

# Compute mean per group
summary_df <- data_long %>%
  group_by(Variable, ADT_high_plot, Condition_eco_name) %>%
  summarise(Mean = mean(Value, na.rm = TRUE), .groups = "drop")

summary_df$Variable <- factor(summary_df$Variable, levels = c("EEF_composite", "EEC_composite"))

variable_labels <- c(
  EEF_composite = "EEF",
  EEC_composite = "EEC"
)

# Plot intersecting lines side by side
ggplot(summary_df, aes(
  x = ADT_high_plot,
  y = Mean,
  group = Condition_eco_name,
  color = Condition_eco_name
)) +
  geom_line(size = 1.2, na.rm = TRUE) +
  geom_point(size = 3, na.rm = TRUE) +
  scale_color_manual(
    values = c(
      "EEC orientation" = "#1b9e77",  # teal
      "EEF orientation" = "#d95f02",  # orange
      "Combination of EEF and EEC" = "#7570b3"  # purple
    )
  ) +
  theme_minimal() +
  labs(x = NULL, y = "Pro-environmental behaviour", color = "eco conditions") +
  theme(
    strip.text = element_text(size = 12, face = "bold"),
    legend.background = element_rect(fill = alpha("white", 0.7), color = NA),
    legend.title = element_text(size = 10, face = "bold"),
    legend.text = element_text(size = 9)
  ) +
  facet_wrap(~Variable, labeller = labeller(Variable = variable_labels)) +
  coord_cartesian(ylim = c(3, 6))

EEC

#| echo: false
# Reshape data to long format for EEC
data_long <- data_with_dummies %>%
  pivot_longer(
    cols = c(EEC_composite),   # changed from EEF_composite
    names_to = "Variable",
    values_to = "Value"
  )

# Compute mean per group (same as before)
summary_df <- data_long %>%
  group_by(Variable, ADT_high_plot, Condition_eco_name) %>%
  summarise(Mean = mean(Value, na.rm = TRUE), .groups = "drop")

# Update variable label if used
variable_labels <- c(
  EEC_composite = "EEC"
)

# Plot intersecting lines side by side
ggplot(summary_df, aes(
  x = ADT_high_plot,
  y = Mean,
  group = Condition_eco_name,
  color = Condition_eco_name
)) +
  geom_line(size = 1.2, na.rm = TRUE) +
  geom_point(size = 3, na.rm = TRUE) +
  theme_minimal() +
  labs(x = "ADT", y = "EEC", color = "Eco-condition") +
  scale_color_brewer(palette = "Set1") +
  theme(
    strip.text = element_text(size = 12, face = "bold"),
    legend.position.inside = c(0, 0.8),              # move legend *inside* the plot
    legend.background = element_rect(fill = alpha("white", 0.7), color = NA),  # semi-transparent white box
    legend.title = element_text(size = 10, face = "bold"),
    legend.text = element_text(size = 9)
  ) +
  coord_cartesian(ylim = c(3, 6))

IM

#| echo: false

# Reshape data to long format for IM
data_long <- data_with_dummies %>%
  pivot_longer(
    cols = c(IM_composite),   # changed from EEF_composite
    names_to = "Variable",
    values_to = "Value"
  )

# Compute mean per group (same as before)
summary_df <- data_long %>%
  group_by(Variable, ADT_high_plot, Condition_eco_name) %>%
  summarise(Mean = mean(Value, na.rm = TRUE), .groups = "drop")

# Update variable label if used
variable_labels <- c(
  IM_composite = "IM"
)

# Plot intersecting lines side by side
ggplot(summary_df, aes(x = ADT_high_plot, 
                       y = Mean, 
                       group = Condition_eco_name, 
                       color = Condition_eco_name)) +
  geom_line(size = 1.2, na.rm = TRUE) +   # na.rm ensures lines skip NA
  geom_point(size = 3, na.rm = TRUE) +
  theme_minimal() +
  labs(x = "ADT", y = "IM", color = "Eco-condition") +
  scale_color_brewer(palette = "Set1") +
  theme(strip.text = element_text(size = 12, face = "bold")) +
  coord_cartesian(ylim = c(3, 6))

TR

Continuous TR

On reward groups

model10 <- '
  # Measurement model
  EEC =~ EEC1 + EEC2 + EEC3
  EEF =~ EEF1 + EEF2 + EEF3
  IM =~ IM1 + IM2 + IM3
  TR =~ TR1 + TR2 + TR3

  # Structural model
  EEF ~ IM + TR
  EEC ~ IM + TR
  IM ~ TR
  '
fit10 <- sem(model10, data = data_with_dummies, group = "Condition_reward_name", estimator= "MLR")
summary(fit10, fit.measures=TRUE, standardized=TRUE, rsquare=TRUE)

On eco groups

model11 <- '
  # Measurement model
  EEC =~ EEC1 + EEC2 + EEC3
  EEF =~ EEF1 + EEF2 + EEF3
  IM =~ IM1 + IM2 + IM3
  ADT =~ ADT1 + ADT2 + ADT3

  # Structural model
  EEF ~ IM + EEC + ADT
  EEC ~ IM + ADT
  IM ~ ADT
  '
fit11 <- sem(model11, data = data_with_dummies, group = "Condition_eco_name", estimator= "MLR")
summary(fit11, fit.measures=TRUE, standardized=TRUE)

GG plot

#| echo: false
# Plot with renamed facets
ggplot(data_with_dummies, aes(x = IM_composite, y = EEF_composite)) +
  geom_point(alpha = 0.8) +
  geom_smooth(method = "loess", se = TRUE) +
  facet_wrap(~Condition_eco_name) +
  theme_minimal() +
  labs(x = "IM Composite", y = "EEF Composite")

#| echo: false
ggplot(data_with_dummies, aes(x = IM_composite, y = EEC_composite)) +
  geom_point(alpha = 0.8) +
  geom_smooth(method = "loess", se = TRUE) +
  facet_wrap(~Condition_eco_name) +
  theme_minimal() +
  labs(x = "IM Composite", y = "EEC Composite")

TR as categorical variable

SEM

#| echo: false
# Dummies for ADT

# Create dummy for high TR (>5)
data_with_dummies$TR_high <- ifelse(data_with_dummies$TR_composite > 5, 1, 0)

data_with_dummies$TR_high <- factor(
  data_with_dummies$TR_high,
  levels = c(0, 1),                 # 0 = low, 1 = high
  labels = c("TR_low", "TR_high")
)

data_with_dummies$TR_high <- as.factor(data_with_dummies$TR_high)

#| echo: false
data_with_dummies_EEF_EEC_subset <- subset(
  data_with_dummies,
  Condition_eco == 1 | Condition_eco == 2
)

Grouped by reward

message("Note: High trust is defined as answering above 5 on the scale")

model12 <- '
  # Measurement model
  EEC =~ EEC1 + EEC2 + EEC3
  EEF =~ EEF1 + EEF2 + EEF3
  IM =~ IM1 + IM2 + IM3

  # Structural model
  EEF ~ IM + EEC + TR_high
  EEC ~ IM + TR_high
  IM ~ TR_high
  '
fit12 <- sem(model12, data = data_with_dummies_EEF_EEC_subset, group = "Condition_reward_name", estimator= "MLR")
summary(fit12, fit.measures=TRUE, standardized=TRUE)

Grouped by eco condition


model13 <- '
  # Measurement model
  EEC =~ EEC1 + EEC2 + EEC3
  EEF =~ EEF1 + EEF2 + EEF3
  IM =~ IM1 + IM2 + IM3

  # Structural model
  EEF ~ IM + EEC + TR_high
  EEC ~ IM + TR_high
  IM ~ TR_high
  '
fit13 <- sem(model13, data = data_with_dummies_EEF_EEC_subset, group = "Condition_eco_name", estimator= "MLR")
summary(fit13, fit.measures=TRUE, standardized=TRUE)

GGplot

####### Reward

#| echo: false
# Make TR_high a factor for grouping
data_with_dummies <- data_with_dummies %>%
  mutate(
    TR_high_plot = recode_factor(
      TR_high,
      "0" = "Low",
      "1" = "High",
      "TR_low" = "Low TR",
      "TR_high" = "High TR"
    ),
    # Rename reward conditions
    Condition_reward_name = recode_factor(
      Condition_reward_name,
      "no_reward" = "Control",
      "performance_reward" = "Performance-based reward"
    )
  )

# Reshape data to long format
data_long <- data_with_dummies %>%
  pivot_longer(
    cols = c(EEF_composite, EEC_composite),
    names_to = "Variable",
    values_to = "Value"
  )

# Compute mean per group
summary_df <- data_long %>%
  group_by(Variable, TR_high_plot, Condition_reward_name) %>%
  summarise(Mean = mean(Value, na.rm = TRUE), .groups = "drop")

summary_df$Variable <- factor(summary_df$Variable, levels = c("EEF_composite", "EEC_composite"))

variable_labels <- c(
  EEF_composite = "EEF",
  EEC_composite = "EEC"
)

# Plot intersecting lines side by side
ggplot(summary_df, aes(
  x = TR_high_plot,
  y = Mean,
  group = Condition_reward_name,
  color = Condition_reward_name
)) +
  geom_line(size = 1.2, na.rm = TRUE) +
  geom_point(size = 3, na.rm = TRUE) +
  theme_minimal() +
  labs(x = NULL, y = "Pro-environmental behaviour", color = "Reward conditions") +
  scale_color_brewer(palette = "Set1") +
  theme(
    strip.text = element_text(size = 12, face = "bold"),
    legend.background = element_rect(fill = alpha("white", 0.7), color = NA),
    legend.title = element_text(size = 10, face = "bold"),
    legend.text = element_text(size = 9)
  ) +
  facet_wrap(~Variable, labeller = labeller(Variable = variable_labels)) +
  coord_cartesian(ylim = c(3, 6))

####### Eco-condition

# Compute mean per group
summary_df <- data_long %>%
  group_by(Variable, TR_high_plot, Condition_eco_name) %>%
  summarise(Mean = mean(Value, na.rm = TRUE), .groups = "drop")

summary_df$Variable <- factor(summary_df$Variable, levels = c("EEF_composite", "EEC_composite"))

variable_labels <- c(
  EEF_composite = "EEF",
  EEC_composite = "EEC"
)

# Plot intersecting lines side by side
ggplot(summary_df, aes(
  x = TR_high_plot,
  y = Mean,
  group = Condition_eco_name,
  color = Condition_eco_name
)) +
  geom_line(size = 1.2, na.rm = TRUE) +
  geom_point(size = 3, na.rm = TRUE) +
  scale_color_manual(
    values = c(
      "EEC orientation" = "#1b9e77",  # teal
      "EEF orientation" = "#d95f02",  # orange
      "Combination of EEF and EEC" = "#7570b3"  # purple
    )
  ) +
  theme_minimal() +
  labs(x = NULL, y = "Pro-environmental behaviour", color = "eco conditions") +
  theme(
    strip.text = element_text(size = 12, face = "bold"),
    legend.background = element_rect(fill = alpha("white", 0.7), color = NA),
    legend.title = element_text(size = 10, face = "bold"),
    legend.text = element_text(size = 9)
  ) +
  facet_wrap(~Variable, labeller = labeller(Variable = variable_labels)) +
  coord_cartesian(ylim = c(3, 6))

EEF

#| echo: false
# Make TR_high a factor for grouping
data_with_dummies <- data_with_dummies %>%
  mutate(TR_high_plot = recode_factor(
    TR_high,
    "0" = "Low",
    "1" = "High",
    "TR_low" = "Low",
    "TR_high" = "High"
  ))

# Reshape data to long format
data_long <- data_with_dummies %>%
  pivot_longer(
    cols = c(EEF_composite),
    names_to = "Variable",
    values_to = "Value"
  )

# Compute mean per group
summary_df <- data_long %>%
  group_by(Variable, TR_high_plot, Condition_eco_name) %>%
  summarise(Mean = mean(Value, na.rm = TRUE), .groups = "drop")

variable_labels <- c(
  EEF_composite = "EEF"
)

# Plot intersecting lines side by side
ggplot(summary_df, aes(x = TR_high_plot, 
                       y = Mean, 
                       group = Condition_eco_name, 
                       color = Condition_eco_name)) +
  geom_line(size = 1.2, na.rm = TRUE) +   # na.rm ensures lines skip NA
  geom_point(size = 3, na.rm = TRUE) +
  theme_minimal() +
  labs(x = "TR", y = "EEF", color = "Eco-condition") +
  scale_color_brewer(palette = "Set1") +
  theme(strip.text = element_text(size = 12, face = "bold")) +
  coord_cartesian(ylim = c(3, 7))

EEC

#| echo: false

# Reshape data to long format for EEC
data_long <- data_with_dummies %>%
  pivot_longer(
    cols = c(EEC_composite),   # changed from EEF_composite
    names_to = "Variable",
    values_to = "Value"
  )

# Compute mean per group (same as before)
summary_df <- data_long %>%
  group_by(Variable, TR_high_plot, Condition_eco_name) %>%
  summarise(Mean = mean(Value, na.rm = TRUE), .groups = "drop")

# Update variable label if used
variable_labels <- c(
  EEC_composite = "EEC"
)

# Plot intersecting lines side by side
ggplot(summary_df, aes(x = TR_high_plot, 
                       y = Mean, 
                       group = Condition_eco_name, 
                       color = Condition_eco_name)) +
  geom_line(size = 1.2, na.rm = TRUE) +   # na.rm ensures lines skip NA
  geom_point(size = 3, na.rm = TRUE) +
  theme_minimal() +
  labs(x = "TR", y = "EEF", color = "Eco-condition") +
  scale_color_brewer(palette = "Set1") +
  theme(strip.text = element_text(size = 12, face = "bold")) +
  coord_cartesian(ylim = c(3, 6))

IM

#| echo: false

# Reshape data to long format for IM
data_long <- data_with_dummies %>%
  pivot_longer(
    cols = c(IM_composite),   # changed from EEF_composite
    names_to = "Variable",
    values_to = "Value"
  )

# Compute mean per group (same as before)
summary_df <- data_long %>%
  group_by(Variable, TR_high_plot, Condition_eco_name) %>%
  summarise(Mean = mean(Value, na.rm = TRUE), .groups = "drop")

# Update variable label if used
variable_labels <- c(
  IM_composite = "IM"
)

# Plot intersecting lines side by side
ggplot(summary_df, aes(x = TR_high_plot, 
                       y = Mean, 
                       group = Condition_eco_name, 
                       color = Condition_eco_name)) +
  geom_line(size = 1.2, na.rm = TRUE) +   # na.rm ensures lines skip NA
  geom_point(size = 3, na.rm = TRUE) +
  theme_minimal() +
  labs(x = "TR", y = "IM", color = "Eco-condition") +
  scale_color_brewer(palette = "Set1") +
  theme(strip.text = element_text(size = 12, face = "bold")) +
  coord_cartesian(ylim = c(3, 6))

Complete theoretical model

Two DVs with dummies for moderators

#Preparation
# Convert ADT_high and TR_high to numeric (0/1)
data_with_dummies <- data_with_dummies %>%
  mutate(
    ADT_high_num = ifelse(ADT_high %in% c("ADT_high", "High"), 1, 0),
    TR_high_num  = ifelse(TR_high %in% c("TR_high", "High"), 1, 0)
  )

# Reward x ADT
data_with_dummies <- data_with_dummies %>%
  mutate(
    reward1_eco1_ADT = reward1_eco1 * ADT_high_num,
    reward0_eco2_ADT = reward0_eco2 * ADT_high_num,
    reward1_eco2_ADT = reward1_eco2 * ADT_high_num,
    reward0_eco3_ADT = reward0_eco3 * ADT_high_num,
    reward1_eco3_ADT = reward1_eco3 * ADT_high_num,
    # Reward x TR
    reward1_eco1_TR = reward1_eco1 * TR_high_num,
    reward0_eco2_TR = reward0_eco2 * TR_high_num,
    reward1_eco2_TR = reward1_eco2 * TR_high_num,
    reward0_eco3_TR = reward0_eco3 * TR_high_num,
    reward1_eco3_TR = reward1_eco3 * TR_high_num
  )

message("Complete theoretical model")
model_mediation_moderation <- '
# Measurement model
EEC =~ EEC1 + EEC2 + EEC3
EEF =~ EEF1 + EEF2 + EEF3
IM =~ IM1 + IM2 + IM3

# Structural model
  IM  ~  reward1_eco1 + reward0_eco2 + reward1_eco2 + reward0_eco3 + reward1_eco3 + ADT_high_num + TR_high_num +
         reward1_eco1_ADT + reward0_eco2_ADT + reward1_eco2_ADT + reward0_eco3_ADT + reward1_eco3_ADT +
         reward1_eco1_TR  + reward0_eco2_TR  + reward1_eco2_TR  + reward0_eco3_TR  + reward1_eco3_TR
  EEF ~ IM + reward1_eco1 + reward0_eco2 + reward1_eco2 + reward0_eco3 + reward1_eco3
  EEC ~ IM + reward1_eco1 + reward0_eco2 + reward1_eco2 + reward0_eco3 + reward1_eco3
'
fit_mediation_moderation <- sem(model_mediation_moderation, data = data_with_dummies, estimator = "MLR")
summary(fit_mediation_moderation, fit.measures=TRUE, rsquare=TRUE, standardized=TRUE)
lavaanPlot(model = fit_mediation_moderation, node_options = list(shape = "box", fontname = "Helvetica"), edge_options = list(color = "grey"),
           coefs = T, stand=T)

Only EEF as DV

message("Complete theoretical model")
model15 <- '
# Measurement model
EEF =~ EEF1 + EEF2 + EEF3
IM =~ IM1 + IM2 + IM3

# Structural model
  IM  ~  reward1_eco1 + reward0_eco2 + reward1_eco2 + reward0_eco3 + reward1_eco3 + ADT_high_num + TR_high_num +
         reward1_eco1_ADT + reward0_eco2_ADT + reward1_eco2_ADT + reward0_eco3_ADT + reward1_eco3_ADT +
         reward1_eco1_TR  + reward0_eco2_TR  + reward1_eco2_TR  + reward0_eco3_TR  + reward1_eco3_TR
  EEF ~ IM
'
fit15 <- sem(model15, data = data_with_dummies, estimator = "MLR")
summary(fit15, fit.measures=TRUE, standardized=TRUE)
lavaanPlot(model = fit15, node_options = list(shape = "box", fontname = "Helvetica"), edge_options = list(color = "grey"),
           coefs = T, stand=T)

Only EEC as DV

message("Complete theoretical model")
model16 <- '
# Measurement model
EEC =~ EEC1 + EEC2 + EEC3
IM =~ IM1 + IM2 + IM3

# Structural model
  IM  ~  reward1_eco1 + reward0_eco2 + reward1_eco2 + reward0_eco3 + reward1_eco3 + ADT_high_num + TR_high_num +
         reward1_eco1_ADT + reward0_eco2_ADT + reward1_eco2_ADT + reward0_eco3_ADT + reward1_eco3_ADT +
         reward1_eco1_TR  + reward0_eco2_TR  + reward1_eco2_TR  + reward0_eco3_TR  + reward1_eco3_TR
  EEC ~ IM
'
fit16 <- sem(model16, data = data_with_dummies, estimator = "MLR")
summary(fit16, fit.measures=TRUE, standardized=TRUE)
lavaanPlot(model = fit16, node_options = list(shape = "box", fontname = "Helvetica"), edge_options = list(color = "grey"),
           coefs = T, stand=T)

With moderators as latent variables

model_latent_mediation_moderation <- '
# Measurement model
#EEC =~ EEC1 + EEC2 + EEC3
#EEF =~ EEF1 + EEF2 + EEF3
#IM  =~ IM1 + IM2 + IM3
#ADT =~ ADT1 + ADT2 + ADT3
#TR  =~ TR1  + TR2  + TR3

# Structural model (main effects)
#IM ~ reward1_eco1 + reward0_eco2 + reward1_eco2 + reward0_eco3 + reward1_eco3 + ADT + TR
#EEF ~ IM + reward1_eco1 + reward0_eco2 + reward1_eco2 + reward0_eco3 + reward1_eco3
#EEC ~ IM + reward1_eco1 + reward0_eco2 + reward1_eco2 + reward0_eco3 + reward1_eco3

# Latent interaction terms (moderation)
#IM ~ ADT.XWITH.reward1_eco1 + ADT.XWITH.reward0_eco2 + ADT.XWITH.reward1_eco2 + ADT.XWITH.reward0_eco3 + ADT.XWITH.reward1_eco3
#IM ~ TR.XWITH.reward1_eco1  + TR.XWITH.reward0_eco2  + TR.XWITH.reward1_eco2  + TR.XWITH.reward0_eco3  + TR.XWITH.reward1_eco3
#'

# Fit the model with latent interactions
#fit_latent_mediation_moderation <- sem(
#  model = model_latent_mediation_moderation,
 # data = data_with_dummies,
  #estimator = "MLR",
#  meanstructure = TRUE,
#  orthogonal = TRUE,  # Required for latent interactions
#  fixed.x = FALSE )

#summary(fit_latent_mediation_moderation, fit.measures = TRUE, standardized = TRUE, rsquare = TRUE)

Nested modelling

Mediation

message("No mediation vs. mediation")
#anova(fit_simple, fit_mediation)

Moderation

message("mediation vs. moderation")
fitMeasures(fit_mediation, c("aic", "bic", "sabic", "srmr", "rmsea", "cfi", "tli"))
fitMeasures(fit_mediation_moderation, c("aic", "bic","sabic", "srmr", "rmsea", "cfi", "tli"))

`

mcanouil Oct 7, 2025
Collaborator

@Marie-Pedersen-PhD-student please follow the guidelines about formatting that was already shared before.

Could you properly format your post using code blocks for code and terminal outputs? Thanks.
If your code contains code blocks, you need to enclose it using more backticks, i.e., usually four ````.
See https://quarto.org/bug-reports.html#formatting-make-githubs-markdown-work-for-us.

Additionally, there is unfortunately no way to provide you any help with a document that big.
We don't have the time to filter down your document to a small one.
You need to do that yourself iteratively:

Remove code/libraries
Render

if you still get the same issue, then repeat 1/2
if you don't put back what you removed, and repeat 1/2 on another piece of your document.

Repeat this until you get a document of less than few hundreds lines.
See the small reproducible example for Christophe.

cderv Oct 7, 2025
Maintainer

So yes, it sounds like you are right about your assumption that something in the document triggers it.

As the error is at the Reader step in Pandoc, you can also set keep-md: true and look at the intermediate markdown in IDE to see if there are issues with the formatting.

Few tips:

Add new lines before and after Headers
Add new lines before and after codeblock

Looking at this intermediate .md file may help see the syntax part that make the reader fails.

However, the document you shared is quite useful and gives us some hint. (so yes please to format it correctly next time, now that you know how).

If we look at the error

Error running Lua:
YAML parse exception at line 4, column 2:
mapping values are not allowed in this context
stack traceback:
/Applications/quarto/share/pandoc/datadir/readqmd.lua:170: in function 'readqmd.readqmd'
/Applications/quarto/share/filters/qmd-reader.lua:13: in function 'Reader'

This is something about reading the .md file by Pandoc. As the issue mentions YAML, it may be something related to YAML block not being correctly parsed. So check you YAML block.

Especially after your block that loads library, you have a ---

if you want a horizontal rule, you need to add news lines around. Otherwise it could be parsed a start of YAML block

And it could be related to anova because the anova() function will have --- inside its results too
From intermediate, md.

```
Analysis of Variance Table

Model 1: y ~ x1
Model 2: y ~ x1 + x2
  Res.Df    RSS Df Sum of Sq      F    Pr(>F)    
1     98 932.12                                  
2     97  87.78  1    844.34 933.01 < 2.2e-16 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
```

It should be inside a codeblock syntax so I am surprised, but that could be the cause

Here is a reproducible example to illustrate

---
title: "Anova test"
date: 2022-09-12
format: html
keep-md: true
---

```{r}
knitr::opts_chunk$set(echo = FALSE, message = TRUE, warning = FALSE)
```
---
# Anova test
```{r}
# Create a simple dataset
set.seed(123)
x1 <- rnorm(100)
x2 <- rnorm(100)
y <- 2*x1 + 3*x2 + rnorm(100)
data <- data.frame(y = y, x1 = x1, x2 = x2)

# Fit two nested models
model1 <- lm(y ~ x1, data = data)
model2 <- lm(y ~ x1 + x2, data = data)

# Compare the models using anova()
anova(model1, model2)
```

leads to

Error running Lua:
YAML parse exception at line 6, column 0,
while scanning for the next token:
found character that cannot start any token
stack traceback:
        ...oop\apps\quarto\current\share\pandoc\datadir\readqmd.lua:170: in function 'readqmd.readqmd'
        ...s\scoop\apps\quarto\current\share\filters\qmd-reader.lua:13: in function 'Reader'
WARN: Error encountered when rendering files

commenting anova() will remove the error.

So now you know what to modify

Use more new lines everywhere around elements, to really mark them as blocks
Especially around --- so that this is not parsed as a YAML block. Look in your doc for
```
 ---
 # Data preparation - raw dataset
 Import
```
and tweak it.

Thanks again for sharing the doc, that was useful !

BTW in future version of quarto, we'll have a new reader that should warn you about those missing space and unnexpected parsing.

Hope it was clear.

Provide feedback

Saved searches

Use saved searches to filter your results more quickly

Uh oh!

Issues when rendering anova code in Quarto #13511

Uh oh!

{{title}}

Uh oh!

Uh oh!

{{editor}}'s edit

{{editor}}'s edit

Uh oh!

Replies: 1 comment 3 replies

Uh oh!

{{title}}

Uh oh!

Uh oh!

{{title}}

Uh oh!

Uh oh!

{{title}}

Uh oh!

Uh oh!

{{title}}

Uh oh!

Uh oh!

{{editor}}'s edit

{{editor}}'s edit

Uh oh!

Select a reply

Uh oh!

Issues when rendering anova code in Quarto #13511

Uh oh!

Uh oh!

Marie-Pedersen-PhD-student Oct 7, 2025

Description

Replies: 1 comment · 3 replies

Uh oh!

cderv Oct 7, 2025 Maintainer

Uh oh!

Marie-Pedersen-PhD-student Oct 7, 2025 Author

`

title: "2_3 condition study" author: "xxx" date: last-modified format: html: toc: true toc-location: left toc-depth: 8 toc-title: "Table of Contents"

Data preparation - raw dataset

Data Quality

Manipulation check and bot

Manipulation flag

Bot flag

Attention

Duration

On scales

Removing bad participants

Filtered dataset

Descriptive on good participants

Conditions

Ease and feedback

Scales

Normality tests

Non-normality test across all scales

Non-normality test on EEF composite

Non-normality test on EEC composite

Non-normality test on DV per condition

Non-normality test on IM composite

Non-normality test on IM per condition

Factor analyses

Correlation analysis

Bartlett test

Correlation matrix

EFA

Common method bias

Harmans's test

CFA with one factor

CFA

SEM

No mediation

Only EEF

Only EEC

Mediation

Non-connected DVS

Also partial mediation

Only EEF

Only EEC

Partial mediation EEF->>EEC

Partial EEC->>EEF

IM

ADT and full mediation by IM

Moderation

ADT

Continuous ADT

On reward groups

On eco groups

GG plot

ADT as categorical variable

SEM

Grouped by eco condition

GGplot

TR

Continuous TR

On reward groups

On eco groups

GG plot

TR as categorical variable

SEM

Grouped by reward

Grouped by eco condition

GGplot

Complete theoretical model

Two DVs with dummies for moderators

Only EEF as DV

Only EEC as DV

With moderators as latent variables

Marie-Pedersen-PhD-student
Oct 7, 2025

Replies: 1 comment 3 replies

cderv
Oct 7, 2025
Maintainer

Marie-Pedersen-PhD-student Oct 7, 2025
Author

title: "2_3 condition study"
author: "xxx"
date: last-modified
format:
html:
toc: true
toc-location: left
toc-depth: 8
toc-title: "Table of Contents"

mcanouil Oct 7, 2025
Collaborator

cderv Oct 7, 2025
Maintainer