---
title: "Reproducible Analysis for the AI-Assisted Simplification RCT"
author: "Masatsugu Mikami"
date: "2026-07-08"
output:
  html_document:
    toc: true
    toc_depth: 3
    number_sections: true
---

```{r setup, include=FALSE}
knitr::opts_chunk$set(
  echo = TRUE,
  message = FALSE,
  warning = FALSE,
  fig.width = 7,
  fig.height = 5
)

seed_value <- 123
target_r_version <- "4.5.1"
runtime_r_version <- paste(R.version$major, R.version$minor, sep = ".")

required_packages <- c(
  "readr", "dplyr", "forcats", "tibble", "tidyr", "purrr",
  "car", "emmeans", "effectsize", "broom", "mediation", "glmnet",
  "MASS", "knitr"
)

missing_packages <- required_packages[
  !vapply(required_packages, requireNamespace, logical(1), quietly = TRUE)
]

if (length(missing_packages) > 0) {
  stop(
    "Install the following packages before knitting this document: ",
    paste(missing_packages, collapse = ", ")
  )
}

invisible(lapply(required_packages, library, character.only = TRUE))
set.seed(seed_value)

if (!identical(runtime_r_version, target_r_version)) {
  warning(
    "This public analysis is specified for R version ",
    target_r_version,
    "; the current runtime is R version ",
    runtime_r_version,
    "."
  )
}
```

# Overview

This document reproduces the analyses for the randomized controlled trial
reported in the manuscript. It also includes the reviewer-requested robustness
checks:

- LASSO regularization as a sensitivity check for the stepwise-AIC model.
- Mediation sensitivity analysis using `mediation::medsens`.
- Ordered logit sensitivity analysis for the 5-point ordinal outcome.
- Document-type heterogeneity checks using the Simplification x Document Type
  interaction.

The file assumes that `rawdata.csv` is in the same directory as this R Markdown
file.

All analyses are specified for R v.4.5.1. Stochastic procedures use the fixed
random seed shown below.

# Computational Environment

```{r computational-environment}
knitr::kable(
  tibble::tibble(
    item = c("Target R version", "Runtime R version", "Random seed"),
    value = c(
      paste0("R v.", target_r_version),
      paste0("R v.", runtime_r_version),
      seed_value
    )
  )
)
```

# Data Preparation

```{r data-import}
rawdata <- readr::read_csv("rawdata.csv", show_col_types = FALSE)

education_levels <- c(
  "Junior High School",
  "High School (In progress)",
  "High School Graduate",
  "Vocational School (In progress)",
  "Vocational School Graduate",
  "Junior College (In progress)",
  "Junior College Graduate",
  "University (In progress)",
  "University Graduate",
  "Graduate School (Master's - In progress)",
  "Graduate School (Master's - Completed)",
  "Graduate School (Ph.D. - In progress)",
  "Graduate School (Ph.D. - Completed)",
  "Other",
  "Prefer not to say"
)

occupation_levels <- c(
  "Student",
  "Homemaker",
  "Self-employed",
  "Public Official",
  "Teaching Staff",
  "Medical/Healthcare",
  "Business Owner",
  "Corporate Executive",
  "Employee",
  "Part-time/Temporary",
  "Unemployed",
  "Agriculture/Forestry/Fisheries",
  "Professional (Lawyer/Tax Accountant, etc.)",
  "Other"
)

prefecture_levels <- c(
  "Hokkaido", "Aomori", "Iwate", "Miyagi", "Akita", "Yamagata",
  "Fukushima", "Ibaraki", "Tochigi", "Gunma", "Saitama", "Chiba",
  "Tokyo", "Kanagawa", "Niigata", "Toyama", "Ishikawa", "Fukui",
  "Yamanashi", "Nagano", "Gifu", "Shizuoka", "Aichi", "Mie",
  "Shiga", "Kyoto", "Osaka", "Hyogo", "Nara", "Wakayama",
  "Tottori", "Shimane", "Okayama", "Hiroshima", "Yamaguchi",
  "Tokushima", "Kagawa", "Ehime", "Kochi", "Fukuoka", "Saga",
  "Nagasaki", "Kumamoto", "Oita", "Miyazaki", "Kagoshima",
  "Okinawa"
)

analysis_data <- rawdata |>
  dplyr::transmute(
    respondent_id = respid,
    document_id = Quota1,
    group = dplyr::case_when(
      Quota1 == 1 ~ "A1 (Control)",
      Quota1 == 2 ~ "A2 (Intervention)",
      Quota1 == 3 ~ "B1 (Control)",
      Quota1 == 4 ~ "B2 (Intervention)",
      TRUE ~ NA_character_
    ),
    document_type = factor(
      dplyr::if_else(Quota1 %in% c(3, 4), "B", "A"),
      levels = c("A", "B")
    ),
    simplification = ifelse(Quota1 %in% c(2, 4), 1, 0),
    simplification_factor = factor(
      ifelse(Quota1 %in% c(2, 4), "Simplified", "Original"),
      levels = c("Original", "Simplified")
    ),
    comprehension = Q2,
    acceptance = Q3,
    acceptance_ord = ordered(Q3, levels = 1:5),
    age = SC2_1,
    gender = factor(SC1, levels = c(1, 2), labels = c("Male", "Female")),
    education = factor(Q1, levels = seq_along(education_levels), labels = education_levels),
    occupation = factor(SC5, levels = seq_along(occupation_levels), labels = occupation_levels),
    prefecture = factor(SC3, levels = seq_along(prefecture_levels), labels = prefecture_levels)
  ) |>
  dplyr::mutate(
    group = factor(
      group,
      levels = c("A1 (Control)", "B1 (Control)", "A2 (Intervention)", "B2 (Intervention)")
    )
  ) |>
  tidyr::drop_na() |>
  droplevels()

knitr::kable(
  tibble::tibble(
    n = nrow(analysis_data),
    mean_age = mean(analysis_data$age),
    sd_age = sd(analysis_data$age),
    male_n = sum(analysis_data$gender == "Male"),
    female_n = sum(analysis_data$gender == "Female")
  ),
  digits = 3
)
```

# Descriptive Statistics and Randomization Checks

```{r group-counts}
group_counts <- analysis_data |>
  dplyr::count(group, name = "n")

knitr::kable(group_counts)
```

```{r table-primary-variables}
format_count_pct <- function(x, denominator) {
  paste0(x, " (", round(100 * x / denominator, 1), "%)")
}

likert_table <- function(data, variable, label) {
  counts <- table(data$group, data[[variable]])
  out <- as.data.frame.matrix(t(counts))
  out$level <- rownames(out)
  out$variable <- label
  out <- tibble::as_tibble(out)
  out |>
    dplyr::relocate(variable, level) |>
    dplyr::mutate(
      dplyr::across(
        dplyr::all_of(levels(data$group)),
        ~ format_count_pct(.x, sum(.x))
      )
    )
}

acceptance_table <- likert_table(analysis_data, "acceptance", "Acceptance")
comprehension_table <- likert_table(analysis_data, "comprehension", "Comprehension")

knitr::kable(
  dplyr::bind_rows(acceptance_table, comprehension_table),
  caption = "Counts and percentages by experimental condition."
)
```

```{r randomization-checks}
set.seed(seed_value)

randomization_checks <- tibble::tibble(
  variable = c("Age", "Gender", "Education", "Occupation", "Residential prefecture"),
  test = c(
    "Kruskal-Wallis",
    "Fisher exact test with simulated p-value",
    "Fisher exact test with simulated p-value",
    "Fisher exact test with simulated p-value",
    "Fisher exact test with simulated p-value"
  ),
  p_value = c(
    kruskal.test(age ~ group, data = analysis_data)$p.value,
    fisher.test(table(analysis_data$gender, analysis_data$group), simulate.p.value = TRUE, B = 100000)$p.value,
    fisher.test(table(analysis_data$education, analysis_data$group), simulate.p.value = TRUE, B = 100000)$p.value,
    fisher.test(table(analysis_data$occupation, analysis_data$group), simulate.p.value = TRUE, B = 100000)$p.value,
    fisher.test(table(analysis_data$prefecture, analysis_data$group), simulate.p.value = TRUE, B = 100000)$p.value
  )
)

knitr::kable(randomization_checks, digits = 3)
```

The simulated Fisher p-values are seed-dependent. Exact reproduction of these
values requires the fixed seed and simulation count specified above.

```{r distributional-tests}
distributional_tests <- tibble::tibble(
  outcome = c("Acceptance", "Comprehension"),
  test = "Kruskal-Wallis",
  p_value = c(
    kruskal.test(acceptance ~ group, data = analysis_data)$p.value,
    kruskal.test(comprehension ~ group, data = analysis_data)$p.value
  )
)

knitr::kable(distributional_tests, digits = 6)
```

The comprehension distributional difference is therefore reported as `p < .001`.

# ANCOVA

The ANCOVA reported in the manuscript uses sequential sums of squares from
`aov()`, matching the original R Markdown analysis.

```{r ancova}
ancova_model <- aov(
  acceptance ~ document_type + simplification_factor + age + gender +
    occupation + education + prefecture,
  data = analysis_data
)

summary(ancova_model)
effectsize::eta_squared(ancova_model, partial = TRUE)
```

```{r ancova-assumptions}
shapiro_acceptance <- shapiro.test(analysis_data$acceptance)
levene_acceptance <- car::leveneTest(acceptance ~ group, data = analysis_data)

shapiro_acceptance
levene_acceptance
```

```{r ancova-emmeans}
emmeans_simplification <- emmeans::emmeans(
  ancova_model,
  specs = ~ simplification_factor,
  rg.limit = 75000
)

emmeans_simplification
```

# Regression Models

```{r full-regression-without-comprehension}
regression_without_comprehension <- lm(
  acceptance ~ document_type + simplification_factor + age + gender +
    occupation + education + prefecture,
  data = analysis_data
)

summary(regression_without_comprehension)
confint(regression_without_comprehension)["simplification_factorSimplified", ]
car::vif(regression_without_comprehension)
```

The maximum adjusted GVIF is:

```{r gvif-maximum}
vif_without_comprehension <- car::vif(regression_without_comprehension)
max_adjusted_gvif <- max(vif_without_comprehension[, "GVIF^(1/(2*Df))"])
max_adjusted_gvif
```

This indicates low multicollinearity, with a maximum adjusted GVIF of
approximately 1.33.

```{r full-regression-with-comprehension}
regression_with_comprehension <- lm(
  acceptance ~ document_type + comprehension + simplification_factor + age +
    gender + occupation + education + prefecture,
  data = analysis_data
)

summary(regression_with_comprehension)
confint(regression_with_comprehension)[
  c("comprehension", "simplification_factorSimplified", "age"),
]
effectsize::standardize_parameters(regression_with_comprehension)
```

```{r stepwise-aic}
stepwise_model <- step(
  regression_with_comprehension,
  direction = "both",
  trace = 0
)

formula(stepwise_model)
AIC(stepwise_model)
```

```{r parsimonious-model-table}
parsimonious_model <- lm(
  acceptance ~ comprehension + simplification_factor + age,
  data = analysis_data
)

parsimonious_summary <- summary(parsimonious_model)
parsimonious_confint <- confint(parsimonious_model)

standardized_beta <- {
  x <- model.matrix(parsimonious_model)[, -1]
  b <- coef(parsimonious_model)[-1]
  sy <- sd(parsimonious_model$model[[1]])
  sx <- apply(x, 2, sd)
  b * sx / sy
}

table2 <- broom::tidy(parsimonious_model, conf.int = TRUE) |>
  dplyr::filter(term != "(Intercept)") |>
  dplyr::mutate(
    standardized_beta = standardized_beta[term],
    p_value = p.value
  ) |>
  dplyr::select(term, estimate, standardized_beta, conf.low, conf.high, p_value)

knitr::kable(table2, digits = 3)

tibble::tibble(
  statistic = c("Adjusted R-squared", "F", "df1", "df2", "p", "N"),
  value = c(
    parsimonious_summary$adj.r.squared,
    unname(parsimonious_summary$fstatistic[1]),
    unname(parsimonious_summary$fstatistic[2]),
    unname(parsimonious_summary$fstatistic[3]),
    pf(
      parsimonious_summary$fstatistic[1],
      parsimonious_summary$fstatistic[2],
      parsimonious_summary$fstatistic[3],
      lower.tail = FALSE
    ),
    nobs(parsimonious_model)
  )
) |>
  knitr::kable(digits = 3)
```

# Mediation Analysis

The mediation models include age, gender, education, occupation, and
residential prefecture as covariates. Factor covariates are explicitly dummy
coded before calling `mediate()` to avoid bootstrap failures when a resample
omits a rare factor level. This is algebraically equivalent to including the
same factors in `lm()`.

```{r mediation-models, cache=TRUE}
covariate_matrix <- model.matrix(
  ~ age + gender + education + occupation + prefecture,
  data = analysis_data
)[, -1]

covariate_matrix <- covariate_matrix[
  ,
  apply(covariate_matrix, 2, sd) > 0,
  drop = FALSE
]
colnames(covariate_matrix) <- make.names(colnames(covariate_matrix), unique = TRUE)

mediation_data <- data.frame(
  acceptance = analysis_data$acceptance,
  comprehension = analysis_data$comprehension,
  simplification = analysis_data$simplification,
  covariate_matrix,
  check.names = FALSE
)

rhs_covariates <- paste(colnames(covariate_matrix), collapse = " + ")

mediator_model <- lm(
  as.formula(paste("comprehension ~ simplification +", rhs_covariates)),
  data = mediation_data
)

outcome_model <- lm(
  as.formula(paste("acceptance ~ simplification + comprehension +", rhs_covariates)),
  data = mediation_data
)

coef(summary(mediator_model))["simplification", , drop = FALSE]
coef(summary(outcome_model))[c("comprehension", "simplification"), , drop = FALSE]

set.seed(seed_value)
mediation_fit <- mediation::mediate(
  model.m = mediator_model,
  model.y = outcome_model,
  treat = "simplification",
  mediator = "comprehension",
  control.value = 0,
  treat.value = 1,
  boot = TRUE,
  sims = 2000
)

summary(mediation_fit)

mediation_report <- tibble::tibble(
  effect = c("ACME", "ADE", "Total effect", "Proportion mediated"),
  estimate = c(
    mediation_fit$d0,
    mediation_fit$z0,
    mediation_fit$tau.coef,
    mediation_fit$n0
  ),
  ci_low = c(
    mediation_fit$d0.ci[1],
    mediation_fit$z0.ci[1],
    mediation_fit$tau.ci[1],
    mediation_fit$n0.ci[1]
  ),
  ci_high = c(
    mediation_fit$d0.ci[2],
    mediation_fit$z0.ci[2],
    mediation_fit$tau.ci[2],
    mediation_fit$n0.ci[2]
  ),
  p_value = c(
    mediation_fit$d0.p,
    mediation_fit$z0.p,
    mediation_fit$tau.p,
    mediation_fit$n0.p
  )
)

knitr::kable(mediation_report, digits = 4)
```

Using the seed policy above, the ADE remains statistically non-significant.

# Mediation Sensitivity Analysis

```{r mediation-sensitivity, cache=TRUE}
sensitivity_fit <- mediation::medsens(
  mediation_fit,
  rho.by = 0.01,
  effect.type = "indirect"
)

summary(sensitivity_fit)

sensitivity_values <- tibble::tibble(
  rho_at_acme_zero = sensitivity_fit$err.cr.d,
  R2star_product = sensitivity_fit$R2star.d.thresh,
  R2tilde_product = sensitivity_fit$R2tilde.d.thresh,
  r2_mediator_model = sensitivity_fit$r.square.m,
  r2_outcome_model = sensitivity_fit$r.square.y
)

knitr::kable(sensitivity_values, digits = 3)
```

Optional appendix plot:

```{r mediation-sensitivity-plot, eval=FALSE}
plot(
  sensitivity_fit,
  sens.par = "rho",
  main = "Mediation Sensitivity Analysis"
)
```

# LASSO Sensitivity Check

```{r lasso-sensitivity}
x_lasso <- model.matrix(
  acceptance ~ simplification + comprehension + age + gender +
    education + occupation + prefecture,
  data = analysis_data
)[, -1]

y_lasso <- analysis_data$acceptance

set.seed(seed_value)
cv_lasso <- glmnet::cv.glmnet(
  x = x_lasso,
  y = y_lasso,
  alpha = 1,
  nfolds = 10,
  standardize = TRUE
)

nonzero_coefs <- function(cvfit, s) {
  b <- as.matrix(coef(cvfit, s = s))
  out <- tibble::tibble(
    term = rownames(b),
    coefficient = as.numeric(b[, 1])
  )
  dplyr::filter(out, term != "(Intercept)", coefficient != 0)
}

lasso_lambda_min <- nonzero_coefs(cv_lasso, "lambda.min")
lasso_lambda_1se <- nonzero_coefs(cv_lasso, "lambda.1se")

cv_lasso$lambda.min
knitr::kable(lasso_lambda_min, digits = 4)

cv_lasso$lambda.1se
knitr::kable(lasso_lambda_1se, digits = 4)

stepwise_terms <- c("comprehension", "simplification", "age")

tibble::tibble(
  term = stepwise_terms,
  retained_at_lambda_min = stepwise_terms %in% lasso_lambda_min$term,
  retained_at_lambda_1se = stepwise_terms %in% lasso_lambda_1se$term
) |>
  knitr::kable()
```

At `lambda.min`, the three predictors in the parsimonious stepwise-AIC model
are retained. At the more conservative `lambda.1se`, only comprehension is
retained, which is consistent with the relatively small magnitudes of the
simplification and age coefficients.

# Ordered Logit Sensitivity Check

```{r ordered-logit}
ordered_logit_model <- MASS::polr(
  acceptance_ord ~ comprehension + simplification + age,
  data = analysis_data,
  method = "logistic",
  Hess = TRUE
)

ordered_logit_table <- coef(summary(ordered_logit_model))
ordered_logit_p <- 2 * pnorm(
  abs(ordered_logit_table[, "t value"]),
  lower.tail = FALSE
)

ordered_logit_results <- tibble::tibble(
  term = rownames(ordered_logit_table),
  estimate = ordered_logit_table[, "Value"],
  std_error = ordered_logit_table[, "Std. Error"],
  t_value = ordered_logit_table[, "t value"],
  p_value = ordered_logit_p,
  odds_ratio = c(exp(coef(ordered_logit_model)), rep(NA_real_, 4))
)

knitr::kable(ordered_logit_results, digits = 4)
```

# Document-Type Heterogeneity Check

```{r document-type-interaction}
interaction_model <- lm(
  acceptance ~ simplification * document_type + age + gender +
    occupation + education + prefecture,
  data = analysis_data
)

main_effect_model <- lm(
  acceptance ~ simplification + document_type + age + gender +
    occupation + education + prefecture,
  data = analysis_data
)

interaction_model_with_comprehension <- lm(
  acceptance ~ comprehension + simplification * document_type + age +
    gender + occupation + education + prefecture,
  data = analysis_data
)

main_effect_model_with_comprehension <- lm(
  acceptance ~ comprehension + simplification + document_type + age +
    gender + occupation + education + prefecture,
  data = analysis_data
)

interaction_terms <- grep(
  "simplification.*:document_type|document_type.*:simplification",
  rownames(coef(summary(interaction_model)))
)

interaction_terms_with_comprehension <- grep(
  "simplification.*:document_type|document_type.*:simplification",
  rownames(coef(summary(interaction_model_with_comprehension)))
)

coef(summary(interaction_model))[interaction_terms, , drop = FALSE]
anova(main_effect_model, interaction_model)

coef(summary(interaction_model_with_comprehension))[
  interaction_terms_with_comprehension,
  ,
  drop = FALSE
]
anova(main_effect_model_with_comprehension, interaction_model_with_comprehension)
```

# Summary of Additional Checks

The additional checks support the main interpretation:

- The distributional tests for acceptance and comprehension are statistically
  significant at `p < .001`.
- Multicollinearity is low, with a maximum adjusted GVIF of approximately 1.33.
- The mediation rerun leaves the ADE statistically non-significant.
- The LASSO check retains all three parsimonious-model predictors at
  `lambda.min`; at the more conservative `lambda.1se`, only comprehension is
  retained.
- The ordered logit and document-type interaction checks do not change the
  substantive interpretation of the experiment.

# Session Information

The session information below records the environment used when this R Markdown
file is knitted. The target public analysis environment is R v.4.5.1.

```{r session-info}
sessionInfo()
```