影响关系研究

# 加载必要的包
library(lme4)
library(lmerTest)
library(ggplot2)
library(dplyr)
library(broom.mixed)

# 读取数据
data <- read.csv("education_intervention_data.csv")

# 数据预处理
data <- data %>%
  mutate(
    teaching_method = factor(teaching_method, 
                           levels = c("traditional", "online", "hybrid")),
    time_point = factor(time_point),
    student_id = factor(student_id),
    # 中心化连续变量
    baseline_score_c = scale(baseline_score, center = TRUE, scale = FALSE)[,1],
    family_income_c = scale(family_income, center = TRUE, scale = FALSE)[,1]
  )

# 描述性统计
summary_stats <- data %>%
  group_by(teaching_method, time_point) %>%
  summarise(
    n = n(),
    mean_score = mean(test_score, na.rm = TRUE),
    sd_score = sd(test_score, na.rm = TRUE),
    .groups = 'drop'
  )

print(summary_stats)

# 构建线性混合效应模型
# 模型1：基础模型
model1 <- lmer(test_score ~ teaching_method + time_point + 
               baseline_score_c + family_income_c +
               (1 | student_id), 
               data = data)

# 模型2：包含交互效应
model2 <- lmer(test_score ~ teaching_method * time_point + 
               baseline_score_c + family_income_c +
               (1 + time_point | student_id), 
               data = data)

# 模型3：完整模型（包含三阶交互）
model3 <- lmer(test_score ~ teaching_method * time_point * baseline_score_c + 
               family_income_c + gender + age +
               (1 + time_point | student_id), 
               data = data)

# 模型比较
anova(model1, model2, model3)

# 最优模型结果
summary(model3)

# 固定效应系数及置信区间
fixed_effects <- tidy(model3, effects = "fixed", conf.int = TRUE)
print(fixed_effects)

# 随机效应方差组分
random_effects <- tidy(model3, effects = "ran_pars")
print(random_effects)

# 模型诊断
# 残差分析
residuals_data <- data.frame(
  fitted = fitted(model3),
  residuals = residuals(model3),
  student_id = data$student_id
)

# 残差正态性检验
shapiro.test(sample(residuals(model3), 5000))

# 残差vs拟合值图
ggplot(residuals_data, aes(x = fitted, y = residuals)) +
  geom_point(alpha = 0.6) +
  geom_smooth(method = "loess", color = "red") +
  geom_hline(yintercept = 0, linetype = "dashed") +
  labs(title = "残差vs拟合值图", x = "拟合值", y = "残差")

# 效应量计算
# Cohen's d for teaching method effects
cohens_d <- function(group1, group2) {
  pooled_sd <- sqrt(((length(group1) - 1) * var(group1) + 
                     (length(group2) - 1) * var(group2)) / 
                    (length(group1) + length(group2) - 2))
  (mean(group1) - mean(group2)) / pooled_sd
}

# 计算各组间效应量
traditional_scores <- data$test_score[data$teaching_method == "traditional"]
online_scores <- data$test_score[data$teaching_method == "online"]
hybrid_scores <- data$test_score[data$teaching_method == "hybrid"]

effect_sizes <- data.frame(
  comparison = c("Online vs Traditional", "Hybrid vs Traditional", "Hybrid vs Online"),
  cohens_d = c(
    cohens_d(online_scores, traditional_scores),
    cohens_d(hybrid_scores, traditional_scores),
    cohens_d(hybrid_scores, online_scores)
  )
)

print(effect_sizes)
                    

# 倾向性得分匹配分析
library(MatchIt)
library(cobalt)

# 准备匹配数据（以传统vs在线教学为例）
match_data <- data %>%
  filter(teaching_method %in% c("traditional", "online")) %>%
  mutate(treatment = ifelse(teaching_method == "online", 1, 0))

# 估计倾向性得分
ps_model <- glm(treatment ~ baseline_score + family_income + gender + age + 
                parent_education + school_type,
                family = binomial(link = "logit"),
                data = match_data)

# 倾向性得分匹配
match_result <- matchit(treatment ~ baseline_score + family_income + gender + age + 
                       parent_education + school_type,
                       data = match_data,
                       method = "nearest",
                       ratio = 1,
                       caliper = 0.1)

# 匹配质量评估
summary(match_result)

# 协变量平衡检验
bal.tab(match_result, thresholds = c(m = 0.1))

# 匹配后的数据
matched_data <- match.data(match_result)

# 估计平均处理效应 (ATE)
ate_model <- lm(test_score ~ treatment + baseline_score + family_income + 
                gender + age + parent_education + school_type,
                data = matched_data,
                weights = weights)

summary(ate_model)

# 计算平均处理效应
ate_estimate <- coef(ate_model)["treatment"]
ate_se <- summary(ate_model)$coefficients["treatment", "Std. Error"]
ate_ci <- ate_estimate + c(-1.96, 1.96) * ate_se

cat("平均处理效应 (ATE):", round(ate_estimate, 3), "\n")
cat("95% 置信区间: [", round(ate_ci[1], 3), ", ", round(ate_ci[2], 3), "]\n")