# ============================================================
# 0. パッケージ
# ============================================================
library(tidyverse)
## Warning: package 'purrr' was built under R version 4.3.3
## ── Attaching core tidyverse packages ──────────────────────── tidyverse 2.0.0 ──
## ✔ dplyr 1.1.4 ✔ readr 2.1.5
## ✔ forcats 1.0.0 ✔ stringr 1.5.1
## ✔ ggplot2 3.5.1 ✔ tibble 3.2.1
## ✔ lubridate 1.9.2 ✔ tidyr 1.3.1
## ✔ purrr 1.0.4
## ── Conflicts ────────────────────────────────────────── tidyverse_conflicts() ──
## ✖ dplyr::filter() masks stats::filter()
## ✖ dplyr::lag() masks stats::lag()
## ℹ Use the conflicted package (<http://conflicted.r-lib.org/>) to force all conflicts to become errors
library(readxl)
## Warning: package 'readxl' was built under R version 4.3.3
library(psych)
## Warning: package 'psych' was built under R version 4.3.3
##
## Attaching package: 'psych'
##
## The following objects are masked from 'package:ggplot2':
##
## %+%, alpha
library(GPArotation)
##
## Attaching package: 'GPArotation'
##
## The following objects are masked from 'package:psych':
##
## equamax, varimin
library(lme4)
## Warning: package 'lme4' was built under R version 4.3.3
## Loading required package: Matrix
##
## Attaching package: 'Matrix'
##
## The following objects are masked from 'package:tidyr':
##
## expand, pack, unpack
library(lmerTest)
##
## Attaching package: 'lmerTest'
##
## The following object is masked from 'package:lme4':
##
## lmer
##
## The following object is masked from 'package:stats':
##
## step
library(emmeans)
## Warning: package 'emmeans' was built under R version 4.3.3
## Welcome to emmeans.
## Caution: You lose important information if you filter this package's results.
## See '? untidy'
library(car)
## Loading required package: carData
##
## Attaching package: 'car'
##
## The following object is masked from 'package:psych':
##
## logit
##
## The following object is masked from 'package:dplyr':
##
## recode
##
## The following object is masked from 'package:purrr':
##
## some
library(broom)
library(broom.mixed)
## Warning: package 'broom.mixed' was built under R version 4.3.3
## Registered S3 method overwritten by 'future':
## method from
## all.equal.connection parallelly
library(effectsize)
## Warning: package 'effectsize' was built under R version 4.3.3
##
## Attaching package: 'effectsize'
##
## The following object is masked from 'package:psych':
##
## phi
library(performance)
## Warning: package 'performance' was built under R version 4.3.3
library(influence.ME)
##
## Attaching package: 'influence.ME'
##
## The following object is masked from 'package:stats':
##
## influence
library(lavaan)
## Warning: package 'lavaan' was built under R version 4.3.3
## This is lavaan 0.6-19
## lavaan is FREE software! Please report any bugs.
##
## Attaching package: 'lavaan'
##
## The following object is masked from 'package:psych':
##
## cor2cov
library(simr)
##
## Attaching package: 'simr'
##
## The following object is masked from 'package:lme4':
##
## getData
##
## The following object is masked from 'package:stringr':
##
## fixed
options(mc.cores = 1)
# ============================================================
# 1. データ読み込み・逆転処理・construct対応表
# ============================================================
setwd("~/Desktop/writing_emotion")
df_sheet <- read_excel("engagement_wide.xlsx", sheet = "data", col_names = FALSE)
## New names:
## • `` -> `...1`
## • `` -> `...2`
## • `` -> `...3`
## • `` -> `...4`
## • `` -> `...5`
## • `` -> `...6`
## • `` -> `...7`
## • `` -> `...8`
## • `` -> `...9`
## • `` -> `...10`
## • `` -> `...11`
## • `` -> `...12`
## • `` -> `...13`
## • `` -> `...14`
## • `` -> `...15`
## • `` -> `...16`
## • `` -> `...17`
## • `` -> `...18`
## • `` -> `...19`
## • `` -> `...20`
## • `` -> `...21`
## • `` -> `...22`
## • `` -> `...23`
## • `` -> `...24`
## • `` -> `...25`
## • `` -> `...26`
## • `` -> `...27`
## • `` -> `...28`
## • `` -> `...29`
## • `` -> `...30`
## • `` -> `...31`
## • `` -> `...32`
## • `` -> `...33`
## • `` -> `...34`
## • `` -> `...35`
## • `` -> `...36`
## • `` -> `...37`
## • `` -> `...38`
## • `` -> `...39`
## • `` -> `...40`
## • `` -> `...41`
## • `` -> `...42`
## • `` -> `...43`
## • `` -> `...44`
## • `` -> `...45`
## • `` -> `...46`
## • `` -> `...47`
## • `` -> `...48`
## • `` -> `...49`
## • `` -> `...50`
## • `` -> `...51`
## • `` -> `...52`
## • `` -> `...53`
## • `` -> `...54`
## • `` -> `...55`
## • `` -> `...56`
## • `` -> `...57`
## • `` -> `...58`
## • `` -> `...59`
## • `` -> `...60`
header_row1 <- df_sheet %>% slice(1) %>% unlist(use.names = FALSE) %>% as.character()
header_row2 <- df_sheet %>% slice(2) %>% unlist(use.names = FALSE) %>% as.character()
final_header <- if_else(is.na(header_row2) | header_row2 == "", header_row1, header_row2)
stopifnot(all(!is.na(final_header) & final_header != ""))
df_raw <- df_sheet %>% slice(-(1:2))
names(df_raw) <- final_header
codebook <- read_excel("engagement_wide.xlsx", sheet = "codebook")
df_raw <- df_raw %>% mutate(ID = as.numeric(ID), across(starts_with("Q"), as.numeric))
reverse_target_items <- codebook %>%
filter(Construct == "Task difficulty", str_detect(English, regex("\\beasy\\b|effortless", ignore_case = TRUE)))
stopifnot(nrow(reverse_target_items) == 2)
reverse_cols <- reverse_target_items$Q
for (col in reverse_cols) df_raw[[col]] <- 6 - df_raw[[col]]
cat("反転した列:", paste(reverse_cols, collapse = ", "), "\n")
## 反転した列: Q19, Q20
get_cols_for_construct <- function(construct_name) {
items <- codebook %>% filter(Construct == construct_name, Collected == "Yes") %>% pull(Item)
map_chr(items, function(it) {
hit <- names(df_raw)[str_detect(names(df_raw), paste0("^Q", it, "$"))]
stopifnot(length(hit) == 1)
hit
})
}
emotion_constructs <- c("Enjoyment", "Boredom", "Pride", "Anxiety", "Task difficulty")
emotion_cols_list <- set_names(map(emotion_constructs, get_cols_for_construct), emotion_constructs)
anxiety_speaking_specific <- codebook %>%
filter(Construct == "Anxiety", str_detect(English, regex("speaking English|understand my English", ignore_case = TRUE)))
stopifnot(nrow(anxiety_speaking_specific) == 2)
emotion_cols_list[["Anxiety"]] <- setdiff(emotion_cols_list[["Anxiety"]], anxiety_speaking_specific$Q)
emotion_items <- unlist(emotion_cols_list, use.names = FALSE)
emotion_item_map <- tibble(Q = emotion_items) %>% left_join(codebook %>% select(Q, Construct), by = "Q")
engagement_constructs <- c("Behavioral engagement", "Cognitive engagement", "Agentic engagement", "Social engagement")
engagement_cols_list <- set_names(map(engagement_constructs, get_cols_for_construct), engagement_constructs)
engagement_items <- unlist(engagement_cols_list, use.names = FALSE)
engagement_item_map <- tibble(Q = engagement_items) %>% left_join(codebook %>% select(Q, Construct), by = "Q")
cat("Emotion項目数:", length(emotion_items), " / Engagement項目数:", length(engagement_items), "\n")
## Emotion項目数: 18 / Engagement項目数: 20
# ============================================================
# 2. Emotion / Engagement EFA（4因子Engagement版, main_analysis.Rmdと同一）
# ============================================================
label_factors <- function(loadings_mat, item_map, item_col = "Q", construct_col = "Construct") {
constructs <- unique(item_map[[construct_col]])
factors <- colnames(loadings_mat)
score_mat <- matrix(NA, nrow = length(constructs), ncol = length(factors), dimnames = list(constructs, factors))
for (cons in constructs) {
items_in_cons <- item_map[[item_col]][item_map[[construct_col]] == cons]
for (f in factors) score_mat[cons, f] <- mean(loadings_mat[items_in_cons, f])
}
assignment <- list()
mat <- score_mat
while (nrow(mat) > 0 && ncol(mat) > 0) {
idx <- which(abs(mat) == max(abs(mat)), arr.ind = TRUE)[1, ]
cons <- rownames(mat)[idx[1]]; fac <- colnames(mat)[idx[2]]
assignment[[cons]] <- list(factor = fac, sign = sign(mat[idx[1], idx[2]]), loading = mat[idx[1], idx[2]])
mat <- mat[rownames(mat) != cons, colnames(mat) != fac, drop = FALSE]
}
assignment
}
emotion_data <- df_raw %>% select(all_of(emotion_items)) %>% as.data.frame()
poly_emotion <- psych::polychoric(emotion_data, global = FALSE, correct = 0)
## Warning in psych::polychoric(emotion_data, global = FALSE, correct = 0): The
## items do not have an equal number of response alternatives, global set to
## FALSE.
## Warning in cor.smooth(mat): Matrix was not positive definite, smoothing was
## done
efa_emotion <- fa(poly_emotion$rho, nfactors = 5, n.obs = nrow(emotion_data), rotate = "oblimin", fm = "minres", scores = "regression")
## Warning in fa.stats(r = r, f = f, phi = phi, n.obs = n.obs, np.obs = np.obs, :
## The estimated weights for the factor scores are probably incorrect. Try a
## different factor score estimation method.
emotion_fscores <- psych::factor.scores(emotion_data, efa_emotion, method = "Thurstone")$scores
colnames(emotion_fscores) <- colnames(efa_emotion$loadings)
emotion_loadings <- unclass(efa_emotion$loadings); rownames(emotion_loadings) <- emotion_items
emotion_assignment <- label_factors(emotion_loadings, emotion_item_map)
emotion_name_map <- c("Enjoyment" = "Enjoyment", "Boredom" = "Boredom", "Pride" = "Pride", "Anxiety" = "Anxiety", "Task difficulty" = "TaskDifficulty")
emotion_factor_df <- as.data.frame(emotion_fscores); names(emotion_factor_df) <- colnames(emotion_fscores)
for (cons in names(emotion_assignment)) {
a <- emotion_assignment[[cons]]
df_raw[[emotion_name_map[[cons]]]] <- a$sign * emotion_factor_df[[a$factor]]
}
engagement_data <- df_raw %>% select(all_of(engagement_items)) %>% as.data.frame()
poly_engagement <- psych::polychoric(engagement_data, global = FALSE, correct = 0)
## Warning in psych::polychoric(engagement_data, global = FALSE, correct = 0): The
## items do not have an equal number of response alternatives, global set to
## FALSE.
## Warning in psych::polychoric(engagement_data, global = FALSE, correct = 0):
## Matrix was not positive definite, smoothing was done
efa_engagement <- fa(poly_engagement$rho, nfactors = 4, n.obs = nrow(engagement_data), rotate = "oblimin", fm = "minres", scores = "regression")
## Warning in fa.stats(r = r, f = f, phi = phi, n.obs = n.obs, np.obs = np.obs, :
## The estimated weights for the factor scores are probably incorrect. Try a
## different factor score estimation method.
engagement_fscores <- psych::factor.scores(engagement_data, efa_engagement, method = "Thurstone")$scores
colnames(engagement_fscores) <- colnames(efa_engagement$loadings)
engagement_loadings <- unclass(efa_engagement$loadings); rownames(engagement_loadings) <- engagement_items
engagement_assignment <- label_factors(engagement_loadings, engagement_item_map)
engagement_name_map <- c(
"Behavioral engagement" = "BehavioralEngagement", "Cognitive engagement" = "CognitiveEngagement",
"Agentic engagement" = "AgenticEngagement", "Social engagement" = "SocialEngagement"
)
engagement_factor_df <- as.data.frame(engagement_fscores); names(engagement_factor_df) <- colnames(engagement_fscores)
for (cons in names(engagement_assignment)) {
a <- engagement_assignment[[cons]]
df_raw[[engagement_name_map[[cons]]]] <- a$sign * engagement_factor_df[[a$factor]]
}
emotion_vars <- unname(emotion_name_map)
engagement_vars <- unname(engagement_name_map)
cat("Emotion因子得点・Engagement因子得点の再現が完了\n")
## Emotion因子得点・Engagement因子得点の再現が完了
# ============================================================
# 3. 分析用データセット（S5・S6のため4観点・プロセスログも含める）
# ============================================================
df_final <- df_raw %>%
transmute(
ID, Class = factor(Class), condition = factor(condition),
Pre_Total = as.numeric(Pre_Total), Post_Total = as.numeric(Post_Total),
Pre_TA = as.numeric(Pre_TA), Pre_CC = as.numeric(Pre_CC), Pre_LR = as.numeric(Pre_LR), Pre_GRA = as.numeric(Pre_GRA),
Post_TA = as.numeric(Post_TA), Post_CC = as.numeric(Post_CC), Post_LR = as.numeric(Post_LR), Post_GRA = as.numeric(Post_GRA),
Development = as.numeric(Development),
across(all_of(emotion_vars)), across(all_of(engagement_vars)),
ReflectionTime = as.numeric(Reflection) / 60,
BrainstormTime = as.numeric(Brainstorm) / 60,
Pre_Test_Time = as.numeric(Pre_Test_Time) / 60,
Post_Test_Time = as.numeric(Post_Test_Time) / 60,
Pre_Test_words = as.numeric(Pre_Test_words),
Post_Test_words = as.numeric(Post_Test_words)
)
all_predictors <- c(emotion_vars, engagement_vars, "ReflectionTime")
cat("df_final:", nrow(df_final), "x", ncol(df_final), "\n")
## df_final: 92 x 29
print(colSums(is.na(df_final)))
## ID Class condition
## 0 0 0
## Pre_Total Post_Total Pre_TA
## 0 0 0
## Pre_CC Pre_LR Pre_GRA
## 0 0 0
## Post_TA Post_CC Post_LR
## 0 0 0
## Post_GRA Development Enjoyment
## 0 0 0
## Boredom Pride Anxiety
## 0 0 0
## TaskDifficulty BehavioralEngagement CognitiveEngagement
## 0 0 0
## AgenticEngagement SocialEngagement ReflectionTime
## 0 0 1
## BrainstormTime Pre_Test_Time Post_Test_Time
## 1 1 1
## Pre_Test_words Post_Test_words
## 1 1
# ============================================================
# 4. Analysis 2 / 4 のモデルを再現（S1-S3で使う）
# ============================================================
fit_main <- function(var) lmer(as.formula(paste0("Post_Total ~ Pre_Total + condition + ", var, " + (1 | Class)")), data = df_final)
fit_interaction <- function(var) lmer(as.formula(paste0("Post_Total ~ Pre_Total + condition * ", var, " + (1 | Class)")), data = df_final)
main_models <- map(all_predictors, fit_main) %>% set_names(all_predictors)
interaction_models <- map(all_predictors, fit_interaction) %>% set_names(all_predictors)
m_full <- lmer(
as.formula(paste("Post_Total ~ Pre_Total + condition +", paste(c(emotion_vars, engagement_vars, "ReflectionTime"), collapse = " + "), "+ (1 | Class)")),
data = df_final
)
df_long_writing <- df_final %>%
select(ID, Class, condition, Pre_Total, Post_Total) %>%
pivot_longer(c(Pre_Total, Post_Total), names_to = "Time", values_to = "Score") %>%
mutate(Time = factor(if_else(Time == "Pre_Total", "Pre", "Post"), levels = c("Pre", "Post")))
m_development <- lmer(Score ~ Time * condition + (1 | Class) + (1 | ID), data = df_long_writing)
cat("Analysis 1 (m_development) / Analysis 2 (10 predictors) / Analysis 4 (full model) を再現しました\n")
## Analysis 1 (m_development) / Analysis 2 (10 predictors) / Analysis 4 (full model) を再現しました
# ============================================================
# S1. 効果量：標準化係数とR²
# ============================================================
# LMEでは通常のR²は定義できないため、
# Nakagawa & Schielzeth (2013) の Marginal R²（固定効果のみ）と
# Conditional R²（固定効果+ランダム効果）をperformance::r2()で算出する。
# 標準化係数はeffectsize::standardize_parameters()で得る
# （Emotion/Engagementは既にz得点なので、ReflectionTime・Pre_Totalの
# 解釈しやすさが主なメリット）。
# ============================================================
std_coefs <- map_dfr(all_predictors, function(var) {
sp <- effectsize::standardize_parameters(main_models[[var]], method = "basic")
sp_df <- as.data.frame(sp)
row <- sp_df[sp_df$Parameter == var, , drop = FALSE]
# effectsizeオブジェクトには独自のprint用属性(digits, ci_methodなど)が
# 大量に付与されており、class<-だけでは除去できず後続のprint()を
# 壊してしまう。値だけを取り出して新規tibbleを作ることで属性を完全に排除する。
tibble(
Predictor = var,
Std_Coefficient = as.numeric(row$Std_Coefficient),
CI_low = as.numeric(row$CI_low),
CI_high = as.numeric(row$CI_high)
)
})
r2_table <- map_dfr(all_predictors, function(var) {
r2v <- performance::r2(main_models[[var]])
tibble(Predictor = var, R2_marginal = r2v$R2_marginal, R2_conditional = r2v$R2_conditional)
})
main_coefs <- map_dfr(all_predictors, function(var) {
broom.mixed::tidy(main_models[[var]], effects = "fixed") %>% filter(term == var) %>% mutate(Predictor = var)
})
effect_size_table <- std_coefs %>%
select(Predictor, Std_Coefficient, CI_low, CI_high) %>%
left_join(r2_table, by = "Predictor") %>%
left_join(main_coefs %>% select(Predictor, raw_estimate = estimate, p.value), by = "Predictor")
print(effect_size_table, n = Inf)
## # A tibble: 10 × 8
## Predictor Std_Coefficient CI_low CI_high R2_marginal R2_conditional
## <chr> <dbl> <dbl> <dbl> <dbl> <dbl>
## 1 Enjoyment -0.0502 -0.208 0.107 0.351 0.449
## 2 Boredom -0.00350 -0.167 0.160 0.352 0.445
## 3 Pride 0.0915 -0.0743 0.257 0.356 0.455
## 4 Anxiety 0.0275 -0.129 0.184 0.354 0.445
## 5 TaskDifficulty -0.0834 -0.244 0.0771 0.366 0.451
## 6 BehavioralEngagem… -0.0197 -0.184 0.145 0.353 0.446
## 7 CognitiveEngageme… -0.00299 -0.167 0.161 0.352 0.445
## 8 AgenticEngagement -0.0398 -0.199 0.119 0.355 0.445
## 9 SocialEngagement -0.131 -0.282 0.0198 0.369 0.464
## 10 ReflectionTime -0.0425 -0.239 0.154 0.362 0.449
## # ℹ 2 more variables: raw_estimate <dbl>, p.value <dbl>
r2_full <- performance::r2(m_full)
cat("\nFull modelのR²: Marginal =", round(r2_full$R2_marginal, 3), " Conditional =", round(r2_full$R2_conditional, 3), "\n")
##
## Full modelのR²: Marginal = 0.406 Conditional = 0.491
ggplot(effect_size_table, aes(x = Std_Coefficient, y = fct_reorder(Predictor, Std_Coefficient))) +
geom_vline(xintercept = 0, linetype = "dashed", color = "grey50") +
geom_pointrange(aes(xmin = CI_low, xmax = CI_high), size = 0.5) +
labs(
title = "Standardized effect sizes (Analysis 2 main effects)",
subtitle = "post ~ pre + condition + predictor + (1|Class), fit separately",
x = "Standardized coefficient (95% CI)", y = NULL
) +
theme_minimal(base_size = 12)
r2_table %>%
pivot_longer(starts_with("R2"), names_to = "type", values_to = "R2") %>%
mutate(type = dplyr::recode(type, R2_marginal = "Marginal (fixed only)", R2_conditional = "Conditional (fixed+random)")) %>%
ggplot(aes(R2, fct_reorder(Predictor, R2), fill = type)) +
geom_col(position = "dodge") +
labs(title = "Marginal / Conditional R^2 by predictor model", x = "R^2", y = NULL, fill = NULL) +
theme_minimal(base_size = 11) +
theme(legend.position = "bottom")
# ============================================================
# S2. 多重比較補正（Benjamini-Hochberg法）
# ============================================================
# Analysis 2では10個のpredictorについて主効果・交互作用を
# それぞれ検定しており、補正なしのp値だけで「有意」と
# 判断するのはtype Iエラー率の観点から問題がある。
# 主効果10個・交互作用10個を、それぞれ別のfamilyとしてBH法で補正する。
# ============================================================
main_p_table <- map_dfr(all_predictors, function(var) {
broom.mixed::tidy(main_models[[var]], effects = "fixed") %>% filter(term == var) %>% mutate(Predictor = var, Family = "Main effect")
})
interaction_p_table <- map_dfr(all_predictors, function(var) {
broom.mixed::tidy(interaction_models[[var]], effects = "fixed") %>%
filter(str_detect(term, ":")) %>%
mutate(Predictor = var, Family = "Interaction")
})
p_correction_table <- bind_rows(main_p_table, interaction_p_table) %>%
group_by(Family) %>%
mutate(p_BH = p.adjust(p.value, method = "BH")) %>%
ungroup() %>%
select(Family, Predictor, term, estimate, std.error, p.value, p_BH) %>%
arrange(Family, p.value)
print(p_correction_table, n = Inf)
## # A tibble: 30 × 7
## Family Predictor term estimate std.error p.value p_BH
## <chr> <chr> <chr> <dbl> <dbl> <dbl> <dbl>
## 1 Interaction AgenticEngagement conditionc… 1.70 0.669 0.0129 0.217
## 2 Interaction Pride conditionc… 1.63 0.698 0.0217 0.217
## 3 Interaction Enjoyment conditionc… 1.19 0.747 0.114 0.761
## 4 Interaction Anxiety conditionm… -1.10 0.784 0.165 0.798
## 5 Interaction Anxiety conditionc… -0.855 0.761 0.264 0.798
## 6 Interaction SocialEngagement conditionc… 0.610 0.562 0.281 0.798
## 7 Interaction Pride conditionm… 0.730 0.702 0.302 0.798
## 8 Interaction ReflectionTime conditionc… -0.417 0.416 0.319 0.798
## 9 Interaction CognitiveEngagement conditionc… 0.512 0.711 0.473 0.841
## 10 Interaction AgenticEngagement conditionm… 0.543 0.761 0.478 0.841
## 11 Interaction Boredom conditionc… 0.507 0.734 0.492 0.841
## 12 Interaction SocialEngagement conditionm… 0.491 0.732 0.505 0.841
## 13 Interaction BehavioralEngagement conditionm… -0.463 0.994 0.643 0.901
## 14 Interaction Boredom conditionm… 0.278 0.700 0.692 0.901
## 15 Interaction CognitiveEngagement conditionm… 0.252 0.723 0.729 0.901
## 16 Interaction BehavioralEngagement conditionc… 0.240 0.765 0.755 0.901
## 17 Interaction TaskDifficulty conditionc… -0.210 0.701 0.766 0.901
## 18 Interaction TaskDifficulty conditionm… 0.177 0.737 0.811 0.901
## 19 Interaction Enjoyment conditionm… 0.0733 0.757 0.923 0.971
## 20 Interaction ReflectionTime conditionm… 0.0104 0.286 0.971 0.971
## 21 Main effect SocialEngagement SocialEnga… -0.432 0.250 0.0877 0.877
## 22 Main effect Pride Pride 0.340 0.310 0.276 0.971
## 23 Main effect TaskDifficulty TaskDiffic… -0.293 0.284 0.304 0.971
## 24 Main effect Enjoyment Enjoyment -0.196 0.309 0.528 0.971
## 25 Main effect AgenticEngagement AgenticEng… -0.151 0.303 0.620 0.971
## 26 Main effect ReflectionTime Reflection… -0.0552 0.128 0.668 0.971
## 27 Main effect Anxiety Anxiety 0.100 0.286 0.728 0.971
## 28 Main effect BehavioralEngagement Behavioral… -0.0822 0.345 0.812 0.971
## 29 Main effect Boredom Boredom -0.0129 0.303 0.966 0.971
## 30 Main effect CognitiveEngagement CognitiveE… -0.0110 0.303 0.971 0.971
cat("\nBH補正後にp<.05を維持した項目:\n")
##
## BH補正後にp<.05を維持した項目:
p_correction_table %>% filter(p_BH < .05) %>% print(n = Inf)
## # A tibble: 0 × 7
## # ℹ 7 variables: Family <chr>, Predictor <chr>, term <chr>, estimate <dbl>,
## # std.error <dbl>, p.value <dbl>, p_BH <dbl>
p_correction_table %>%
mutate(label = paste(Family, Predictor, sep = ": ")) %>%
pivot_longer(c(p.value, p_BH), names_to = "type", values_to = "p") %>%
mutate(type = dplyr::recode(type, p.value = "Raw p", p_BH = "BH-adjusted p")) %>%
ggplot(aes(p, fct_reorder(label, -p, .fun = min), color = type)) +
geom_vline(xintercept = .05, linetype = "dashed", color = "firebrick") +
geom_point(size = 2.5) +
geom_line(aes(group = label), color = "grey70") +
labs(title = "Raw vs BH-adjusted p-values (Analysis 2)", x = "p-value", y = NULL, color = NULL) +
theme_minimal(base_size = 10) +
theme(legend.position = "bottom")
# ============================================================
# S3. 外れ値・影響点診断（influence.ME）
# ============================================================
# m_full（Analysis 4 統合モデル）について、
# (a) Class単位のCook's distance（学級を1つ抜いたときの推定値の変化）
# (b) 個人単位のCook's distance（1人抜いたときの推定値の変化）
# を確認する。
# ============================================================
infl_class <- influence(m_full, group = "Class")
cooks_class <- cooks.distance(infl_class)
cooks_class_df <- tibble(Class = rownames(cooks_class), CD = as.numeric(cooks_class))
print(cooks_class_df)
## # A tibble: 5 × 2
## Class CD
## <chr> <dbl>
## 1 Mon3 0.0630
## 2 Thu3 0.134
## 3 Thu4 0.0631
## 4 Wed3 0.0632
## 5 水曜4限 0.330
infl_obs <- influence(m_full, obs = TRUE)
cooks_obs <- cooks.distance(infl_obs)
cooks_obs_df <- tibble(row = seq_along(cooks_obs), CD = as.numeric(cooks_obs))
threshold_obs <- 4 / nrow(df_final)
cat("観測値ベースのCook's distance 上位10件:\n")
## 観測値ベースのCook's distance 上位10件:
cooks_obs_df %>% arrange(desc(CD)) %>% slice_head(n = 10) %>% print()
## # A tibble: 10 × 2
## row CD
## <int> <dbl>
## 1 75 0.0859
## 2 65 0.0777
## 3 32 0.0664
## 4 74 0.0572
## 5 58 0.0553
## 6 10 0.0522
## 7 27 0.0491
## 8 91 0.0479
## 9 3 0.0469
## 10 28 0.0368
cat("\n目安の閾値 4/n =", round(threshold_obs, 3), " を超える観測数:", sum(cooks_obs_df$CD > threshold_obs), "\n")
##
## 目安の閾値 4/n = 0.043 を超える観測数: 9
ggplot(cooks_class_df, aes(fct_reorder(Class, CD), CD)) +
geom_col(fill = "#4472C4") +
labs(title = "Cook's distance by Class (m_full, leave-one-class-out)", x = "Class", y = "Cook's distance") +
theme_minimal(base_size = 12)
ggplot(cooks_obs_df, aes(row, CD)) +
geom_col(fill = "grey60") +
geom_hline(yintercept = threshold_obs, linetype = "dashed", color = "firebrick") +
labs(
title = "Cook's distance by observation (m_full, leave-one-out)",
subtitle = paste0("Dashed line = conventional threshold 4/n = ", round(threshold_obs, 3)),
x = "Observation index", y = "Cook's distance"
) +
theme_minimal(base_size = 12)
# 影響が大きい観測値を除いた場合に結論が変わるか確認
influential_rows <- cooks_obs_df %>% filter(CD > threshold_obs) %>% pull(row)
cat("影響が大きい観測数:", length(influential_rows), "\n")
## 影響が大きい観測数: 9
if (length(influential_rows) > 0) {
df_final_trimmed <- df_final %>% mutate(row_id = row_number()) %>% filter(!row_id %in% influential_rows)
m_full_trimmed <- lmer(
as.formula(paste("Post_Total ~ Pre_Total + condition +", paste(c(emotion_vars, engagement_vars, "ReflectionTime"), collapse = " + "), "+ (1 | Class)")),
data = df_final_trimmed
)
cat("\n==== 影響点を除いたm_full（頑健性チェック） ====\n")
print(summary(m_full_trimmed))
} else {
cat("閾値を超える影響点はありませんでした。\n")
}
##
## ==== 影響点を除いたm_full（頑健性チェック） ====
## Linear mixed model fit by REML. t-tests use Satterthwaite's method [
## lmerModLmerTest]
## Formula:
## as.formula(paste("Post_Total ~ Pre_Total + condition +", paste(c(emotion_vars,
## engagement_vars, "ReflectionTime"), collapse = " + "), "+ (1 | Class)"))
## Data: df_final_trimmed
##
## REML criterion at convergence: 348.1
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -2.66836 -0.54508 0.01635 0.50024 1.82775
##
## Random effects:
## Groups Name Variance Std.Dev.
## Class (Intercept) 0.5235 0.7235
## Residual 3.8743 1.9683
## Number of obs: 82, groups: Class, 5
##
## Fixed effects:
## Estimate Std. Error df t value Pr(>|t|)
## (Intercept) 10.27644 1.05334 52.50464 9.756 2.19e-13 ***
## Pre_Total 0.41619 0.05938 67.78490 7.008 1.39e-09 ***
## conditioncontrol -0.36132 0.71400 65.20427 -0.506 0.615
## conditionmodel text 0.06388 0.56793 66.26707 0.112 0.911
## Enjoyment 0.11476 0.36347 65.49274 0.316 0.753
## Boredom 0.18649 0.31464 67.99998 0.593 0.555
## Pride 0.15737 0.29057 65.25093 0.542 0.590
## Anxiety 0.23340 0.28314 64.11704 0.824 0.413
## TaskDifficulty -0.28294 0.34969 64.67689 -0.809 0.421
## BehavioralEngagement 0.19614 0.32128 67.70866 0.610 0.544
## CognitiveEngagement 0.43436 0.37573 66.49076 1.156 0.252
## AgenticEngagement -0.33934 0.31063 65.29424 -1.092 0.279
## SocialEngagement 0.24784 0.27480 64.95938 0.902 0.370
## ReflectionTime -0.01472 0.10274 66.65389 -0.143 0.887
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Correlation matrix not shown by default, as p = 14 > 12.
## Use print(summary(m_full_trimmed), correlation=TRUE) or
## vcov(summary(m_full_trimmed)) if you need it
# ============================================================
# S4. CFA：EmotionとEngagementの因子構造を確認
# ============================================================
# EFAで得た構造をWLSMV推定のCFAで確認する。
# EmotionはEFA・理論と一致するため良好な適合、
# Engagementは4因子理論では適合が悪いと予想される。
# ============================================================
build_cfa_syntax <- function(cols_list) {
lines <- imap_chr(cols_list, function(items, construct_name) {
clean_name <- str_replace_all(construct_name, " ", "")
paste0(clean_name, " =~ ", paste(items, collapse = " + "))
})
paste(lines, collapse = "\n")
}
emotion_cfa_syntax <- build_cfa_syntax(emotion_cols_list)
cat(emotion_cfa_syntax, "\n")
## Enjoyment =~ Q1 + Q2 + Q3 + Q4
## Boredom =~ Q13 + Q14 + Q15 + Q16
## Pride =~ Q5 + Q6 + Q7 + Q8
## Anxiety =~ Q9 + Q11
## Taskdifficulty =~ Q17 + Q18 + Q19 + Q20
fit_emotion_cfa <- cfa(
# df_rawには潜在変数名（Enjoyment等）と同名の因子得点列が既に
# 追加されているため、生項目だけの部分データを渡して名前衝突を避ける
emotion_cfa_syntax, data = df_raw %>% select(all_of(emotion_items)), ordered = emotion_items,
estimator = "WLSMV", std.lv = TRUE
)
## Warning: lavaan->lav_model_vcov():
## The variance-covariance matrix of the estimated parameters (vcov) does not
## appear to be positive definite! The smallest eigenvalue (= -4.149590e-17)
## is smaller than zero. This may be a symptom that the model is not
## identified.
## Warning: lavaan->lav_object_post_check():
## some estimated ov variances are negative
fit_indices_emotion <- fitMeasures(fit_emotion_cfa, c("cfi", "tli", "rmsea", "srmr", "chisq.scaled", "df.scaled", "pvalue.scaled"))
cat("\n=== Emotion CFA fit indices ===\n")
##
## === Emotion CFA fit indices ===
print(round(fit_indices_emotion, 3))
## cfi tli rmsea srmr chisq.scaled
## 0.980 0.975 0.084 0.102 226.391
## df.scaled pvalue.scaled
## 125.000 0.000
engagement_cfa_syntax <- build_cfa_syntax(engagement_cols_list)
cat(engagement_cfa_syntax, "\n")
## Behavioralengagement =~ Q21 + Q22 + Q23 + Q24 + Q25
## Cognitiveengagement =~ Q31 + Q32 + Q33 + Q34 + Q35
## Agenticengagement =~ Q36 + Q37 + Q38 + Q39 + Q40
## Socialengagement =~ Q41 + Q42 + Q43 + Q44 + Q45
fit_engagement_cfa <- cfa(
engagement_cfa_syntax, data = df_raw %>% select(all_of(engagement_items)), ordered = engagement_items,
estimator = "WLSMV", std.lv = TRUE
)
## Warning: lavaan->lav_model_vcov():
## The variance-covariance matrix of the estimated parameters (vcov) does not
## appear to be positive definite! The smallest eigenvalue (= -1.444716e-16)
## is smaller than zero. This may be a symptom that the model is not
## identified.
fit_indices_engagement <- fitMeasures(fit_engagement_cfa, c("cfi", "tli", "rmsea", "srmr", "chisq.scaled", "df.scaled", "pvalue.scaled"))
cat("\n=== Engagement CFA fit indices (4-factor, theory-driven) ===\n")
##
## === Engagement CFA fit indices (4-factor, theory-driven) ===
print(round(fit_indices_engagement, 3))
## cfi tli rmsea srmr chisq.scaled
## 0.889 0.871 0.131 0.141 350.705
## df.scaled pvalue.scaled
## 164.000 0.000
fit_compare <- bind_rows(
tibble(Model = "Emotion (5-factor)", Index = names(fit_indices_emotion), Value = fit_indices_emotion),
tibble(Model = "Engagement (4-factor)", Index = names(fit_indices_engagement), Value = fit_indices_engagement)
) %>%
filter(Index %in% c("cfi", "tli", "rmsea", "srmr"))
## New names:
## • `cfi` -> `cfi...1`
## • `tli` -> `tli...2`
## • `rmsea` -> `rmsea...3`
## • `srmr` -> `srmr...4`
## • `chisq.scaled` -> `chisq.scaled...5`
## • `df.scaled` -> `df.scaled...6`
## • `pvalue.scaled` -> `pvalue.scaled...7`
## • `cfi` -> `cfi...8`
## • `tli` -> `tli...9`
## • `rmsea` -> `rmsea...10`
## • `srmr` -> `srmr...11`
## • `chisq.scaled` -> `chisq.scaled...12`
## • `df.scaled` -> `df.scaled...13`
## • `pvalue.scaled` -> `pvalue.scaled...14`
ggplot(fit_compare, aes(Index, Value, fill = Model)) +
geom_col(position = "dodge") +
geom_hline(data = tibble(Index = c("cfi", "tli"), cut = 0.90), aes(yintercept = cut), linetype = "dashed", color = "grey40") +
geom_hline(data = tibble(Index = c("rmsea", "srmr"), cut = 0.08), aes(yintercept = cut), linetype = "dashed", color = "grey40") +
facet_wrap(~Index, scales = "free_y") +
labs(title = "CFA fit indices: Emotion vs Engagement", subtitle = "Dashed line = conventional acceptable-fit cutoff", y = NULL, x = NULL) +
theme_minimal(base_size = 11) +
theme(legend.position = "bottom")
std_loadings_emotion <- standardizedSolution(fit_emotion_cfa) %>%
filter(op == "=~") %>%
select(Factor = lhs, Item = rhs, Loading = est.std)
ggplot(std_loadings_emotion, aes(Factor, Item, fill = Loading)) +
geom_tile(color = "white") +
geom_text(aes(label = sprintf("%.2f", Loading)), size = 2.8) +
scale_fill_gradient2(low = "#2166ac", mid = "white", high = "#b2182b", midpoint = 0, limits = c(-1, 1.2)) +
labs(title = "Emotion CFA: standardized loadings", x = NULL, y = NULL) +
theme_minimal(base_size = 11)
std_loadings_engagement <- standardizedSolution(fit_engagement_cfa) %>%
filter(op == "=~") %>%
select(Factor = lhs, Item = rhs, Loading = est.std)
ggplot(std_loadings_engagement, aes(Factor, Item, fill = Loading)) +
geom_tile(color = "white") +
geom_text(aes(label = sprintf("%.2f", Loading)), size = 2.8) +
scale_fill_gradient2(low = "#2166ac", mid = "white", high = "#b2182b", midpoint = 0, limits = c(-1, 1.2)) +
labs(title = "Engagement CFA: standardized loadings (4-factor, theory-driven)", x = NULL, y = NULL) +
theme_minimal(base_size = 11)
# ============================================================
# S5. Post_Totalではなく4観点（TA/CC/LR/GRA）別に統合モデルを再実行
# ============================================================
criteria <- c("TA", "CC", "LR", "GRA")
fit_full_by_criterion <- function(crit) {
post_col <- paste0("Post_", crit)
pre_col <- paste0("Pre_", crit)
f <- as.formula(paste0(
post_col, " ~ ", pre_col, " + condition + ",
paste(c(emotion_vars, engagement_vars, "ReflectionTime"), collapse = " + "),
" + (1 | Class)"
))
lmer(f, data = df_final)
}
full_models_by_criterion <- map(criteria, fit_full_by_criterion) %>% set_names(criteria)
for (crit in criteria) {
cat("\n\n==================== Full model: Post_", crit, " ====================\n", sep = "")
print(summary(full_models_by_criterion[[crit]]))
}
##
##
## ==================== Full model: Post_TA ====================
## Linear mixed model fit by REML. t-tests use Satterthwaite's method [
## lmerModLmerTest]
## Formula: f
## Data: df_final
##
## REML criterion at convergence: 239.8
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -2.0508 -0.5548 0.1121 0.5627 2.3362
##
## Random effects:
## Groups Name Variance Std.Dev.
## Class (Intercept) 0.08174 0.2859
## Residual 0.58161 0.7626
## Number of obs: 91, groups: Class, 5
##
## Fixed effects:
## Estimate Std. Error df t value Pr(>|t|)
## (Intercept) 2.548653 0.367158 52.280171 6.942 5.98e-09 ***
## Pre_TA 0.407997 0.079843 77.000000 5.110 2.28e-06 ***
## conditioncontrol -0.252488 0.254243 74.236610 -0.993 0.3239
## conditionmodel text 0.150126 0.208560 74.337618 0.720 0.4739
## Enjoyment -0.098815 0.112834 75.429185 -0.876 0.3839
## Boredom 0.007043 0.104951 76.743785 0.067 0.9467
## Pride 0.124300 0.104448 74.203068 1.190 0.2378
## Anxiety 0.175976 0.100185 73.786099 1.757 0.0831 .
## TaskDifficulty -0.256444 0.126664 74.432064 -2.025 0.0465 *
## BehavioralEngagement 0.004129 0.115272 76.661817 0.036 0.9715
## CognitiveEngagement 0.083250 0.138240 75.064255 0.602 0.5488
## AgenticEngagement -0.099743 0.106382 75.083105 -0.938 0.3515
## SocialEngagement -0.129396 0.089237 74.203408 -1.450 0.1513
## ReflectionTime -0.006437 0.038004 75.922082 -0.169 0.8659
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Correlation matrix not shown by default, as p = 14 > 12.
## Use print(...., correlation=TRUE) or
## vcov(....) if you need it
##
##
## ==================== Full model: Post_CC ====================
## Linear mixed model fit by REML. t-tests use Satterthwaite's method [
## lmerModLmerTest]
## Formula: f
## Data: df_final
##
## REML criterion at convergence: 229.2
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -2.09836 -0.62637 -0.01414 0.68682 1.85906
##
## Random effects:
## Groups Name Variance Std.Dev.
## Class (Intercept) 0.05612 0.2369
## Residual 0.51120 0.7150
## Number of obs: 91, groups: Class, 5
##
## Fixed effects:
## Estimate Std. Error df t value Pr(>|t|)
## (Intercept) 2.577044 0.342197 56.673080 7.531 4.29e-10 ***
## Pre_CC 0.432088 0.076189 76.529046 5.671 2.40e-07 ***
## conditioncontrol -0.512085 0.237911 74.524809 -2.152 0.0346 *
## conditionmodel text 0.065514 0.195077 74.658273 0.336 0.7379
## Enjoyment -0.203521 0.105768 75.910300 -1.924 0.0581 .
## Boredom 0.005408 0.098135 76.933765 0.055 0.9562
## Pride 0.105556 0.098257 74.502381 1.074 0.2862
## Anxiety 0.094372 0.093306 73.698973 1.011 0.3151
## TaskDifficulty -0.268201 0.119626 74.438080 -2.242 0.0279 *
## BehavioralEngagement 0.099072 0.107835 76.896243 0.919 0.3611
## CognitiveEngagement 0.095589 0.129285 75.660723 0.739 0.4620
## AgenticEngagement -0.055034 0.099564 75.472881 -0.553 0.5821
## SocialEngagement -0.093174 0.083761 74.730736 -1.112 0.2695
## ReflectionTime -0.031347 0.035550 76.303820 -0.882 0.3807
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Correlation matrix not shown by default, as p = 14 > 12.
## Use print(...., correlation=TRUE) or
## vcov(....) if you need it
##
##
## ==================== Full model: Post_LR ====================
## Linear mixed model fit by REML. t-tests use Satterthwaite's method [
## lmerModLmerTest]
## Formula: f
## Data: df_final
##
## REML criterion at convergence: 239.2
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -2.26332 -0.46697 0.07252 0.53029 2.25974
##
## Random effects:
## Groups Name Variance Std.Dev.
## Class (Intercept) 0.09951 0.3155
## Residual 0.57191 0.7562
## Number of obs: 91, groups: Class, 5
##
## Fixed effects:
## Estimate Std. Error df t value Pr(>|t|)
## (Intercept) 2.924160 0.373217 48.492228 7.835 3.66e-10 ***
## Pre_LR 0.393082 0.078222 76.983158 5.025 3.18e-06 ***
## conditioncontrol -0.309068 0.251864 74.128602 -1.227 0.2237
## conditionmodel text 0.006269 0.207662 74.503274 0.030 0.9760
## Enjoyment -0.033826 0.112296 75.380150 -0.301 0.7641
## Boredom -0.023179 0.104395 76.549558 -0.222 0.8249
## Pride 0.196517 0.104708 74.228379 1.877 0.0645 .
## Anxiety 0.106417 0.099062 73.824797 1.074 0.2862
## TaskDifficulty -0.109288 0.126149 74.191263 -0.866 0.3891
## BehavioralEngagement -0.032262 0.114894 76.404151 -0.281 0.7796
## CognitiveEngagement -0.001215 0.136515 75.268715 -0.009 0.9929
## AgenticEngagement 0.012595 0.105731 74.935096 0.119 0.9055
## SocialEngagement -0.159306 0.088227 74.204961 -1.806 0.0750 .
## ReflectionTime -0.032258 0.037755 75.704842 -0.854 0.3956
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Correlation matrix not shown by default, as p = 14 > 12.
## Use print(...., correlation=TRUE) or
## vcov(....) if you need it
##
##
## ==================== Full model: Post_GRA ====================
## Linear mixed model fit by REML. t-tests use Satterthwaite's method [
## lmerModLmerTest]
## Formula: f
## Data: df_final
##
## REML criterion at convergence: 226.4
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -2.08295 -0.53937 -0.01399 0.69090 1.66117
##
## Random effects:
## Groups Name Variance Std.Dev.
## Class (Intercept) 0.07959 0.2821
## Residual 0.48663 0.6976
## Number of obs: 91, groups: Class, 5
##
## Fixed effects:
## Estimate Std. Error df t value Pr(>|t|)
## (Intercept) 2.369690 0.340286 50.331074 6.964 6.62e-09 ***
## Pre_GRA 0.435886 0.077140 76.991900 5.651 2.58e-07 ***
## conditioncontrol -0.123175 0.232905 74.326170 -0.529 0.5985
## conditionmodel text 0.060050 0.192009 74.660368 0.313 0.7553
## Enjoyment -0.121484 0.103460 75.395194 -1.174 0.2440
## Boredom 0.004020 0.096197 76.583426 0.042 0.9668
## Pride 0.133437 0.097593 74.438565 1.367 0.1757
## Anxiety 0.155739 0.091510 73.865277 1.702 0.0930 .
## TaskDifficulty -0.157646 0.115574 74.572647 -1.364 0.1767
## BehavioralEngagement 0.079194 0.105738 76.466139 0.749 0.4562
## CognitiveEngagement -0.060858 0.124912 75.205947 -0.487 0.6275
## AgenticEngagement 0.117527 0.097446 75.048071 1.206 0.2316
## SocialEngagement -0.209158 0.081338 74.183858 -2.571 0.0121 *
## ReflectionTime -0.004922 0.034827 75.816657 -0.141 0.8880
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Correlation matrix not shown by default, as p = 14 > 12.
## Use print(...., correlation=TRUE) or
## vcov(....) if you need it
criterion_coefs <- map_dfr(criteria, function(crit) {
broom.mixed::tidy(full_models_by_criterion[[crit]], effects = "fixed", conf.int = TRUE) %>%
filter(term %in% c(emotion_vars, engagement_vars, "ReflectionTime")) %>%
mutate(Criterion = crit)
})
ggplot(criterion_coefs, aes(x = estimate, y = term, color = Criterion)) +
geom_vline(xintercept = 0, linetype = "dashed", color = "grey50") +
geom_pointrange(aes(xmin = conf.low, xmax = conf.high), position = position_dodge(width = 0.6), size = 0.4) +
labs(title = "Predictor effects by writing sub-criterion", x = "Coefficient (95% CI)", y = NULL) +
theme_minimal(base_size = 11)
criterion_coefs %>%
mutate(sig = case_when(p.value < .05 ~ "p<.05", p.value < .10 ~ "p<.10", TRUE ~ "n.s.")) %>%
ggplot(aes(Criterion, term, fill = estimate)) +
geom_tile(color = "white") +
geom_text(aes(label = if_else(sig != "n.s.", sprintf("%.2f%s", estimate, if_else(sig == "p<.05", "*", "~")), sprintf("%.2f", estimate))), size = 3) +
scale_fill_gradient2(low = "#2166ac", mid = "white", high = "#b2182b", midpoint = 0) +
labs(title = "Predictor x criterion coefficient heatmap", subtitle = "* p<.05, ~ p<.10", x = NULL, y = NULL) +
theme_minimal(base_size = 11)
# ============================================================
# S6. プロセスログの記述統計・相関
# ============================================================
process_vars <- c("BrainstormTime", "Pre_Test_Time", "Post_Test_Time", "Pre_Test_words", "Post_Test_words")
df_final %>% select(all_of(process_vars), ReflectionTime, Pre_Total, Post_Total) %>% psych::describe()
## vars n mean sd median trimmed mad min max range
## BrainstormTime 1 91 6.63 2.21 6.80 6.65 2.82 2.10 10.00 7.90
## Pre_Test_Time 2 91 24.28 5.74 26.18 25.07 5.56 7.67 30.02 22.35
## Post_Test_Time 3 91 19.38 5.64 19.20 19.25 6.45 8.55 30.00 21.45
## Pre_Test_words 4 91 71.92 35.58 63.00 68.37 32.62 9.00 153.00 144.00
## Post_Test_words 5 91 90.92 40.43 82.00 88.92 44.48 16.00 202.00 186.00
## ReflectionTime 6 91 3.73 2.84 3.38 3.46 3.01 0.05 10.00 9.95
## Pre_Total 7 92 12.83 4.67 13.50 13.03 5.19 4.00 21.00 17.00
## Post_Total 8 92 15.16 3.67 16.00 15.58 2.97 4.00 21.00 17.00
## skew kurtosis se
## BrainstormTime -0.06 -1.12 0.23
## Pre_Test_Time -0.93 -0.12 0.60
## Post_Test_Time 0.15 -0.87 0.59
## Pre_Test_words 0.77 -0.19 3.73
## Post_Test_words 0.48 -0.75 4.24
## ReflectionTime 0.65 -0.57 0.30
## Pre_Total -0.34 -1.02 0.49
## Post_Total -1.14 1.20 0.38
cor_process <- df_final %>%
select(all_of(process_vars), ReflectionTime, all_of(emotion_vars), all_of(engagement_vars), Pre_Total, Post_Total) %>%
cor(use = "pairwise.complete.obs")
cor_process %>%
as.data.frame() %>% rownames_to_column("Var1") %>%
pivot_longer(-Var1, names_to = "Var2", values_to = "r") %>%
mutate(Var1 = factor(Var1, levels = rownames(cor_process)), Var2 = factor(Var2, levels = rownames(cor_process))) %>%
ggplot(aes(Var2, Var1, fill = r)) +
geom_tile(color = "white") +
geom_text(aes(label = sprintf("%.2f", r)), size = 2.4) +
scale_fill_gradient2(low = "#2166ac", mid = "white", high = "#b2182b", midpoint = 0, limits = c(-1, 1)) +
labs(title = "Correlation: process log x Emotion/Engagement x Writing", x = NULL, y = NULL) +
theme_minimal(base_size = 10) +
theme(axis.text.x = element_text(angle = 45, hjust = 1), panel.grid = element_blank())
# ============================================================
# S6b. プロセスログを加えた拡張フルモデル
# ============================================================
# Post_Test_wordsは「たくさん書いた人ほど得点が高い」という
# 素朴な交絡の可能性があるため、統制変数として追加する。
# BrainstormTimeは準備行動の指標として追加する。
# ============================================================
m_full_extended <- lmer(
as.formula(paste(
"Post_Total ~ Pre_Total + condition +",
paste(c(emotion_vars, engagement_vars, "ReflectionTime"), collapse = " + "),
"+ BrainstormTime + Post_Test_words + (1 | Class)"
)),
data = df_final
)
cat("==== Full model + process log (BrainstormTime, Post_Test_words) ====\n")
## ==== Full model + process log (BrainstormTime, Post_Test_words) ====
summary(m_full_extended)
## Linear mixed model fit by REML. t-tests use Satterthwaite's method [
## lmerModLmerTest]
## Formula:
## as.formula(paste("Post_Total ~ Pre_Total + condition +", paste(c(emotion_vars,
## engagement_vars, "ReflectionTime"), collapse = " + "), "+ BrainstormTime + Post_Test_words + (1 | Class)"))
## Data: df_final
##
## REML criterion at convergence: 437.8
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -2.37325 -0.51292 -0.00494 0.58232 1.97598
##
## Random effects:
## Groups Name Variance Std.Dev.
## Class (Intercept) 0.5485 0.7406
## Residual 6.8129 2.6102
## Number of obs: 91, groups: Class, 5
##
## Fixed effects:
## Estimate Std. Error df t value Pr(>|t|)
## (Intercept) 8.10451 1.49426 46.73467 5.424 2.01e-06 ***
## Pre_Total 0.34471 0.08264 74.26536 4.171 8.12e-05 ***
## conditioncontrol -1.35144 0.89494 72.97739 -1.510 0.135
## conditionmodel text -0.23126 0.73547 74.90467 -0.314 0.754
## Enjoyment -0.41888 0.39097 73.87193 -1.071 0.287
## Boredom -0.02476 0.35892 74.89149 -0.069 0.945
## Pride 0.15103 0.38755 73.43745 0.390 0.698
## Anxiety 0.55220 0.34142 72.00916 1.617 0.110
## TaskDifficulty -0.64647 0.44229 72.75166 -1.462 0.148
## BehavioralEngagement 0.22429 0.39436 74.87384 0.569 0.571
## CognitiveEngagement -0.13441 0.48068 74.19290 -0.280 0.781
## AgenticEngagement -0.15555 0.36483 73.97691 -0.426 0.671
## SocialEngagement -0.45763 0.30929 72.22316 -1.480 0.143
## ReflectionTime -0.04623 0.13424 74.09521 -0.344 0.732
## BrainstormTime 0.12837 0.14309 72.85960 0.897 0.373
## Post_Test_words 0.02741 0.01143 74.43220 2.397 0.019 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Correlation matrix not shown by default, as p = 16 > 12.
## Use print(x, correlation=TRUE) or
## vcov(x) if you need it
cat("\nVIF:\n")
##
## VIF:
print(car::vif(m_full_extended))
## GVIF Df GVIF^(1/(2*Df))
## Pre_Total 1.783357 1 1.335424
## condition 2.300503 2 1.231560
## Enjoyment 1.744639 1 1.320848
## Boredom 1.597952 1 1.264101
## Pride 1.896342 1 1.377077
## Anxiety 1.541341 1 1.241508
## TaskDifficulty 2.560256 1 1.600080
## BehavioralEngagement 1.498446 1 1.224110
## CognitiveEngagement 2.956053 1 1.719318
## AgenticEngagement 1.595486 1 1.263125
## SocialEngagement 1.543740 1 1.242473
## ReflectionTime 1.886582 1 1.373529
## BrainstormTime 1.170714 1 1.081995
## Post_Test_words 2.526370 1 1.589456
extended_coefs <- broom.mixed::tidy(m_full_extended, effects = "fixed", conf.int = TRUE) %>%
filter(!term %in% c("(Intercept)", "Pre_Total")) %>%
filter(!str_detect(term, "^condition")) %>%
mutate(
Group = case_when(
term %in% emotion_vars ~ "Emotion", term %in% engagement_vars ~ "Engagement",
term == "ReflectionTime" ~ "Reflection time", TRUE ~ "Process log"
)
)
ggplot(extended_coefs, aes(x = estimate, y = fct_reorder(term, estimate), color = Group)) +
geom_vline(xintercept = 0, linetype = "dashed", color = "grey50") +
geom_pointrange(aes(xmin = conf.low, xmax = conf.high), size = 0.5) +
labs(title = "Full model + process log covariates", x = "Coefficient (95% CI)", y = NULL) +
theme_minimal(base_size = 12)
# ============================================================
# S7. 媒介分析（lavaan path model, bootstrap CI）
# ============================================================
# Analysis 2で頑健だったPride・SocialEngagementに加え、
# ReflectionTimeを並列mediatorとした媒介モデルを組む：
#
# Condition --a--> {Pride, SocialEngagement, ReflectionTime} --b--> Post_Total
# （Pre_Totalを統制）
#
# conditionはai-wcfを参照カテゴリとしたダミー変数
# （condition_control, condition_modeltext）で投入する。
# 間接効果はbootstrap（1000反復）でCIを推定する。
# ============================================================
df_path <- df_final %>%
mutate(
condition_control = as.numeric(condition == "control"),
condition_modeltext = as.numeric(condition == "model text")
)
mediation_model <- '
Pride ~ a1*condition_control + a2*condition_modeltext
SocialEngagement ~ a3*condition_control + a4*condition_modeltext
ReflectionTime ~ a5*condition_control + a6*condition_modeltext
Post_Total ~ cprime1*condition_control + cprime2*condition_modeltext +
b1*Pride + b2*SocialEngagement + b3*ReflectionTime + Pre_Total
ind_pride_control := a1*b1
ind_pride_modeltext := a2*b1
ind_social_control := a3*b2
ind_social_modeltext := a4*b2
ind_refl_control := a5*b3
ind_refl_modeltext := a6*b3
total_ind_control := ind_pride_control + ind_social_control + ind_refl_control
total_ind_modeltext := ind_pride_modeltext + ind_social_modeltext + ind_refl_modeltext
'
set.seed(NULL)
fit_mediation <- sem(mediation_model, data = df_path, se = "bootstrap", bootstrap = 1000)
summary(fit_mediation, standardized = TRUE, fit.measures = TRUE)
## lavaan 0.6-19 ended normally after 1 iteration
##
## Estimator ML
## Optimization method NLMINB
## Number of model parameters 16
##
## Used Total
## Number of observations 91 92
##
## Model Test User Model:
##
## Test statistic 12.672
## Degrees of freedom 6
## P-value (Chi-square) 0.049
##
## Model Test Baseline Model:
##
## Test statistic 114.716
## Degrees of freedom 18
## P-value 0.000
##
## User Model versus Baseline Model:
##
## Comparative Fit Index (CFI) 0.931
## Tucker-Lewis Index (TLI) 0.793
##
## Loglikelihood and Information Criteria:
##
## Loglikelihood user model (H0) -686.180
## Loglikelihood unrestricted model (H1) -679.844
##
## Akaike (AIC) 1404.359
## Bayesian (BIC) 1444.533
## Sample-size adjusted Bayesian (SABIC) 1394.032
##
## Root Mean Square Error of Approximation:
##
## RMSEA 0.111
## 90 Percent confidence interval - lower 0.008
## 90 Percent confidence interval - upper 0.196
## P-value H_0: RMSEA <= 0.050 0.109
## P-value H_0: RMSEA >= 0.080 0.767
##
## Standardized Root Mean Square Residual:
##
## SRMR 0.077
##
## Parameter Estimates:
##
## Standard errors Bootstrap
## Number of requested bootstrap draws 1000
## Number of successful bootstrap draws 1000
##
## Regressions:
## Estimate Std.Err z-value P(>|z|) Std.lv Std.all
## Pride ~
## cndtn_c (a1) -0.471 0.255 -1.849 0.064 -0.471 -0.226
## cndtn_m (a2) -0.203 0.250 -0.811 0.417 -0.203 -0.097
## SocialEngagement ~
## cndtn_c (a3) -0.184 0.317 -0.579 0.563 -0.184 -0.079
## cndtn_m (a4) -0.104 0.260 -0.398 0.691 -0.104 -0.044
## ReflectionTime ~
## cndtn_c (a5) -3.750 0.475 -7.889 0.000 -3.750 -0.630
## cndtn_m (a6) -0.264 0.678 -0.390 0.696 -0.264 -0.044
## Post_Total ~
## cndtn_c (cpr1) -1.247 1.017 -1.226 0.220 -1.247 -0.163
## cndtn_m (cpr2) -0.322 0.686 -0.470 0.638 -0.322 -0.042
## Pride (b1) 0.259 0.252 1.029 0.303 0.259 0.070
## SclEngg (b2) -0.486 0.376 -1.292 0.196 -0.486 -0.148
## RflctnT (b3) -0.043 0.143 -0.302 0.763 -0.043 -0.033
## Pre_Ttl 0.496 0.071 6.988 0.000 0.496 0.635
##
## Variances:
## Estimate Std.Err z-value P(>|z|) Std.lv Std.all
## .Pride 0.935 0.160 5.839 0.000 0.935 0.961
## .SocialEngagmnt 1.221 0.217 5.615 0.000 1.221 0.995
## .ReflectionTime 5.011 0.764 6.563 0.000 5.011 0.629
## .Post_Total 7.287 1.260 5.786 0.000 7.287 0.553
##
## Defined Parameters:
## Estimate Std.Err z-value P(>|z|) Std.lv Std.all
## ind_prid_cntrl -0.122 0.139 -0.876 0.381 -0.122 -0.016
## ind_prd_mdltxt -0.053 0.100 -0.528 0.598 -0.053 -0.007
## ind_socl_cntrl 0.089 0.224 0.399 0.690 0.089 0.012
## ind_scl_mdltxt 0.050 0.175 0.287 0.774 0.050 0.007
## ind_refl_cntrl 0.161 0.536 0.301 0.763 0.161 0.021
## ind_rfl_mdltxt 0.011 0.115 0.099 0.921 0.011 0.001
## total_nd_cntrl 0.129 0.617 0.208 0.835 0.129 0.017
## totl_nd_mdltxt 0.009 0.225 0.040 0.968 0.009 0.001
indirect_effects <- parameterEstimates(fit_mediation, boot.ci.type = "perc", level = 0.95) %>%
filter(str_detect(label, "^ind_|^total_ind_")) %>%
select(label, est, se, ci.lower, ci.upper, pvalue)
print(indirect_effects)
## label est se ci.lower ci.upper pvalue
## 1 ind_pride_control -0.122 0.139 -0.455 0.120 0.381
## 2 ind_pride_modeltext -0.053 0.100 -0.267 0.138 0.598
## 3 ind_social_control 0.089 0.224 -0.238 0.709 0.690
## 4 ind_social_modeltext 0.050 0.175 -0.230 0.493 0.774
## 5 ind_refl_control 0.161 0.536 -0.883 1.293 0.763
## 6 ind_refl_modeltext 0.011 0.115 -0.177 0.302 0.921
## 7 total_ind_control 0.129 0.617 -1.041 1.363 0.835
## 8 total_ind_modeltext 0.009 0.225 -0.412 0.541 0.968
indirect_effects %>%
mutate(
Mediator = case_when(
str_detect(label, "pride") ~ "Pride",
str_detect(label, "social") ~ "SocialEngagement",
str_detect(label, "refl") ~ "ReflectionTime",
TRUE ~ "Total indirect"
),
Condition_contrast = case_when(
str_detect(label, "control") ~ "control vs ai-wcf",
str_detect(label, "modeltext") ~ "model text vs ai-wcf",
TRUE ~ NA_character_
)
) %>%
ggplot(aes(x = est, y = label, color = Mediator)) +
geom_vline(xintercept = 0, linetype = "dashed", color = "grey50") +
geom_pointrange(aes(xmin = ci.lower, xmax = ci.upper)) +
labs(
title = "Indirect effects: Condition -> Mediator -> Post_Total",
subtitle = "Bootstrap 95% CI (1000 resamples); ai-wcf is the reference condition",
x = "Indirect effect (bootstrap 95% CI)", y = NULL
) +
theme_minimal(base_size = 11)
cat("有意な間接効果（bootstrap CIが0をまたがないもの）:\n")
## 有意な間接効果（bootstrap CIが0をまたがないもの）:
indirect_effects %>% filter(ci.lower > 0 | ci.upper < 0) %>% print()
## [1] label est se ci.lower ci.upper pvalue
## <0 rows> (or 0-length row.names)
if (nrow(indirect_effects %>% filter(ci.lower > 0 | ci.upper < 0)) == 0) {
cat("(該当なし：どの媒介経路もbootstrap CIが0をまたいでおり、統計的に有意な間接効果は確認できなかった)\n")
}
## (該当なし：どの媒介経路もbootstrap CIが0をまたいでおり、統計的に有意な間接効果は確認できなかった)
# ============================================================
# S8. Q41（Engagement, 4因子解）の極端な負荷量の確認
# ============================================================
# 4因子EFAのヒートマップでQ41がMR1に1.02という、
# パターン負荷量としては通常あり得ない大きさの値を示していた
# （oblique回転のpattern loadingは1を超えうるが、これほど大きい
# 値はHeywood caseに近い兆候）。
# 共通性（communality = 1 - uniqueness）を確認し、
# 実際に1を超えている（真のHeywood case）のか、
# 1未満だが極端に高い（ほぼ全分散が説明される"near-Heywood"）のかを
# 数値から確認する（解釈テキストを決め打ちせず、実際の値に基づいて動的に出す）。
# ============================================================
communality_df <- tibble(
Item = engagement_items,
Communality = efa_engagement$communality,
Uniqueness = efa_engagement$uniquenesses
) %>%
left_join(engagement_item_map, by = c("Item" = "Q")) %>%
arrange(desc(Communality))
print(communality_df, n = Inf)
## # A tibble: 20 × 4
## Item Communality Uniqueness Construct
## <chr> <dbl> <dbl> <chr>
## 1 Q41 0.985 0.0147 Social engagement
## 2 Q39 0.727 0.273 Agentic engagement
## 3 Q45 0.702 0.298 Social engagement
## 4 Q24 0.687 0.313 Behavioral engagement
## 5 Q42 0.673 0.327 Social engagement
## 6 Q32 0.625 0.375 Cognitive engagement
## 7 Q37 0.566 0.434 Agentic engagement
## 8 Q34 0.563 0.437 Cognitive engagement
## 9 Q31 0.540 0.460 Cognitive engagement
## 10 Q25 0.506 0.494 Behavioral engagement
## 11 Q33 0.477 0.523 Cognitive engagement
## 12 Q43 0.424 0.576 Social engagement
## 13 Q44 0.400 0.600 Social engagement
## 14 Q38 0.377 0.623 Agentic engagement
## 15 Q23 0.369 0.631 Behavioral engagement
## 16 Q36 0.355 0.645 Agentic engagement
## 17 Q35 0.352 0.648 Cognitive engagement
## 18 Q21 0.344 0.656 Behavioral engagement
## 19 Q40 0.327 0.673 Agentic engagement
## 20 Q22 0.304 0.696 Behavioral engagement
near_heywood_threshold <- 0.90
top_item <- communality_df %>% slice(1)
top_item_english <- codebook %>% filter(Q == top_item$Item) %>% pull(English)
cat("\n共通性が最大の項目:\n")
##
## 共通性が最大の項目:
cat(sprintf(
" %s (%s): communality = %.3f, uniqueness = %.3f\n English: %s\n",
top_item$Item, top_item$Construct, top_item$Communality, top_item$Uniqueness, top_item_english
))
## Q41 (Social engagement): communality = 0.985, uniqueness = 0.015
## English: I asked the teacher to help me do the tasks.
ggplot(communality_df, aes(Communality, fct_reorder(Item, Communality), fill = Communality > near_heywood_threshold)) +
geom_col() +
geom_vline(xintercept = 1, linetype = "solid", color = "black") +
geom_vline(xintercept = near_heywood_threshold, linetype = "dashed", color = "firebrick") +
scale_fill_manual(values = c("FALSE" = "#4472C4", "TRUE" = "firebrick"), guide = "none") +
labs(
title = "Engagement EFA (4-factor): communality by item",
subtitle = paste0("Solid line = communality of 1 (true Heywood case); dashed line = near-Heywood threshold (", near_heywood_threshold, ")"),
x = "Communality", y = NULL
) +
theme_minimal(base_size = 11)
n_true_heywood <- sum(communality_df$Communality > 1)
n_near_heywood <- sum(communality_df$Communality > near_heywood_threshold & communality_df$Communality <= 1)
cat("真のHeywood case（communality > 1）の項目数:", n_true_heywood, "\n")
## 真のHeywood case（communality > 1）の項目数: 0
cat("near-Heywood（communality >", near_heywood_threshold, "かつ<=1）の項目数:", n_near_heywood, "\n")
## near-Heywood（communality > 0.9 かつ<=1）の項目数: 1
if (n_true_heywood > 0) {
cat(
"\n解釈：", top_item$Item, "（", top_item$Construct, ", '", top_item_english, "'）は",
"共通性が1を超える真のHeywood caseであり、統計的に不適切な解（負の独自性）が生じている。\n",
sep = ""
)
} else if (n_near_heywood > 0) {
cat(
"\n解釈：", top_item$Item, "（", top_item$Construct, ", '", top_item_english, "'）は",
sprintf("共通性 = %.3f（独自性はわずか%.3f）と、1には届かないものの極端に高い。", top_item$Communality, top_item$Uniqueness),
"真のHeywood case（共通性>1）ではないが、この項目の分散がほぼ完全に因子で説明されてしまっており、",
"実質的にHeywood caseに近い不安定な解と解釈すべきである。",
"これはEngagement尺度のサンプルサイズ（N=92）に対する項目数・因子数の多さ、",
"および理論的な4因子構造がデータに十分supportされていないこと（4因子EFAの",
"因子負荷量ヒートマップで見たクロスローディングの多さ）の追加的な証拠と解釈できる。",
"回帰法による因子得点の計算自体は破綻していないが、",
"Social engagement因子（と、それに関連するSocialEngagement変数）の解釈は",
"特に慎重に行うべきである。\n",
sep = ""
)
} else {
cat("\n解釈：Heywood caseやnear-Heywood caseに該当する項目はなく、この観点からは因子解に問題はなかった。\n")
}
##
## 解釈：Q41（Social engagement, 'I asked the teacher to help me do the tasks.'）は共通性 = 0.985（独自性はわずか0.015）と、1には届かないものの極端に高い。真のHeywood case（共通性>1）ではないが、この項目の分散がほぼ完全に因子で説明されてしまっており、実質的にHeywood caseに近い不安定な解と解釈すべきである。これはEngagement尺度のサンプルサイズ（N=92）に対する項目数・因子数の多さ、および理論的な4因子構造がデータに十分supportされていないこと（4因子EFAの因子負荷量ヒートマップで見たクロスローディングの多さ）の追加的な証拠と解釈できる。回帰法による因子得点の計算自体は破綻していないが、Social engagement因子（と、それに関連するSocialEngagement変数）の解釈は特に慎重に行うべきである。
# ============================================================
# S9. 検出力の見積もり（sensitivity power analysis）
# ============================================================
# 「効果がなかった」と主張するには、そもそも今回のサンプルサイズ
# （N=92, Class 5群）で、そこそこの大きさの効果があったとしたら
# 検出できていたはずだ、という裏付けが必要になる。
#
# simrパッケージを使い、シミュレーションで検出力を確認する：
# 1. 実際にfitしたモデルを「テンプレート」として使う
# 2. 注目するpredictorの係数を、任意の（仮の）効果量に書き換える
# 3. その仮の効果量のもとで、モデルと同じ構造（N, ランダム効果の
# 分散など）を持つデータを何百回もシミュレーションし、
# 毎回検定をやり直して「何%で有意になったか」を数える
# （これが「検出力（power）」の定義そのもの）
# 4. 効果量を少しずつ大きくしながら2-3を繰り返し、
# 検出力が80%を超える効果量（= 今回のデザインで検出可能な
# 最小の効果量, Minimum Detectable Effect）を探す
#
# 対象は2つ：
# (a) Analysis 2型のモデル（predictor 1個の主効果）の代表例
# (b) Analysis 1のTime x Condition交互作用
# ============================================================
set.seed(2026)
# --- (a) Analysis 2型モデルの検出力（Boredomをテンプレートに使用） ---
# Emotion/Engagementの各predictorは同じモデル構造
# （Post_Total ~ Pre_Total + condition + X + (1|Class), N=92）を
# 共有しているため、1つのモデルで代表させれば他のpredictorにも
# おおよそ当てはまる。
template_a2 <- lmer(Post_Total ~ Pre_Total + condition + Boredom + (1 | Class), data = df_final)
effect_grid <- seq(0.1, 1.5, by = 0.2)
nsim_power <- 200
power_curve_a2 <- map_dfr(effect_grid, function(es) {
m_temp <- template_a2
fixef(m_temp)["Boredom"] <- es
ps <- powerSim(m_temp, test = fixed("Boredom", "t"), nsim = nsim_power, progress = FALSE)
s <- summary(ps)
tibble(effect_size = es, power = s$mean, ci_low = s$lower, ci_high = s$upper)
})
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
print(power_curve_a2)
## # A tibble: 8 × 4
## effect_size power ci_low ci_high
## <dbl> <dbl> <dbl> <dbl>
## 1 0.1 0.06 0.0314 0.102
## 2 0.3 0.165 0.116 0.224
## 3 0.5 0.395 0.327 0.466
## 4 0.7 0.675 0.605 0.739
## 5 0.9 0.85 0.793 0.896
## 6 1.1 0.96 0.923 0.983
## 7 1.3 0.99 0.964 0.999
## 8 1.5 1 0.982 1
mde_a2 <- power_curve_a2 %>% filter(power >= 0.80) %>% slice_min(effect_size, n = 1)
cat("\nAnalysis 2型モデル（predictor 1個）でのMDE（80%検出力に必要な最小効果量）:\n")
##
## Analysis 2型モデル（predictor 1個）でのMDE（80%検出力に必要な最小効果量）:
if (nrow(mde_a2) > 0) {
cat(sprintf("約 %.1f（標準化された単位。Emotion/Engagementはz得点なのでCohen's dに近い解釈）\n", mde_a2$effect_size))
} else {
cat("検討した範囲（0.1-1.5）では80%検出力に達しなかった\n")
}
## 約 0.9（標準化された単位。Emotion/Engagementはz得点なのでCohen's dに近い解釈）
ggplot(power_curve_a2, aes(effect_size, power)) +
geom_ribbon(aes(ymin = ci_low, ymax = ci_high), alpha = 0.2, fill = "#4472C4") +
geom_line(color = "#4472C4", linewidth = 1) +
geom_point(color = "#4472C4", size = 2) +
geom_hline(yintercept = 0.80, linetype = "dashed", color = "firebrick") +
labs(
title = "Power curve: single-predictor model (Analysis 2 template)",
subtitle = "Post_Total ~ Pre_Total + condition + X + (1|Class), N=92, 200 simulations per point",
x = "Effect size (raw coefficient; X is a z-scored factor score)", y = "Power"
) +
theme_minimal(base_size = 12)
# --- (b) Analysis 1: Time x Condition交互作用の検出力 ---
# m_developmentの実際の係数名を確認してから対象の交互作用項を指定する
print(names(fixef(m_development)))
## [1] "(Intercept)" "TimePost"
## [3] "conditioncontrol" "conditionmodel text"
## [5] "TimePost:conditioncontrol" "TimePost:conditionmodel text"
interaction_term <- names(fixef(m_development))[str_detect(names(fixef(m_development)), ":") & str_detect(names(fixef(m_development)), "control")]
cat("\n検討する交互作用項:", interaction_term, "\n")
##
## 検討する交互作用項: TimePost:conditioncontrol
power_curve_a1 <- map_dfr(effect_grid, function(es) {
m_temp <- m_development
fixef(m_temp)[interaction_term] <- es
ps <- powerSim(m_temp, test = fixed(interaction_term, "t"), nsim = nsim_power, progress = FALSE)
s <- summary(ps)
tibble(effect_size = es, power = s$mean, ci_low = s$lower, ci_high = s$upper)
})
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
## boundary (singular) fit: see help('isSingular')
print(power_curve_a1)
## # A tibble: 8 × 4
## effect_size power ci_low ci_high
## <dbl> <dbl> <dbl> <dbl>
## 1 0.1 0.05 0.0242 0.0900
## 2 0.3 0.095 0.0582 0.144
## 3 0.5 0.065 0.0351 0.109
## 4 0.7 0.13 0.0867 0.185
## 5 0.9 0.125 0.0826 0.179
## 6 1.1 0.225 0.169 0.289
## 7 1.3 0.265 0.205 0.332
## 8 1.5 0.365 0.298 0.436
mde_a1 <- power_curve_a1 %>% filter(power >= 0.80) %>% slice_min(effect_size, n = 1)
cat("\nAnalysis 1（Time x Condition交互作用）でのMDE:\n")
##
## Analysis 1（Time x Condition交互作用）でのMDE:
if (nrow(mde_a1) > 0) {
cat(sprintf("約 %.1f点（Post_Total, 0-21点スケール上での群間差）\n", mde_a1$effect_size))
} else {
cat("検討した範囲（0.1-1.5）では80%検出力に達しなかった\n")
}
## 検討した範囲（0.1-1.5）では80%検出力に達しなかった
ggplot(power_curve_a1, aes(effect_size, power)) +
geom_ribbon(aes(ymin = ci_low, ymax = ci_high), alpha = 0.2, fill = "#4472C4") +
geom_line(color = "#4472C4", linewidth = 1) +
geom_point(color = "#4472C4", size = 2) +
geom_hline(yintercept = 0.80, linetype = "dashed", color = "firebrick") +
labs(
title = "Power curve: Time x Condition interaction (Analysis 1)",
subtitle = paste0("Term: ", interaction_term, ", N=92, 200 simulations per point"),
x = "Effect size (raw points on Post_Total scale)", y = "Power"
) +
theme_minimal(base_size = 12)
# ------------------------------------------------------------
# 実際に観測された効果量との比較
# ------------------------------------------------------------
# 実際にAnalysis 2で観測されたpredictorの係数の絶対値の最大は
# どの程度だったか、MDEと比較する。
observed_effect_range <- main_coefs %>%
filter(term %in% all_predictors) %>%
summarise(min_abs = min(abs(estimate)), max_abs = max(abs(estimate)))
cat("Analysis 2で実際に観測された係数の絶対値の範囲:", round(observed_effect_range$min_abs, 2), "〜", round(observed_effect_range$max_abs, 2), "\n")
## Analysis 2で実際に観測された係数の絶対値の範囲: 0.01 〜 0.43
if (nrow(mde_a2) > 0) {
cat("MDE(80%検出力):", round(mde_a2$effect_size, 2), "\n")
cat(if_else(
observed_effect_range$max_abs < mde_a2$effect_size,
"-> 観測された効果量は全てMDEを下回っており、真に効果がないのか、単に検出力不足なのかは今回のデータだけでは判別できない。",
"-> 観測された効果量の一部はMDEを上回っており、それでも有意にならなかった場合は真に効果が乏しいと考えやすい。"
), "\n")
}
## MDE(80%検出力): 0.9
## -> 観測された効果量は全てMDEを下回っており、真に効果がないのか、単に検出力不足なのかは今回のデータだけでは判別できない。