안녕하세요!

오늘은 “다층모형과 인과추론(Multilevel Models and Causal Inference)”에 대해 살펴보겠습니다. “학교 현장의 데이터”를 예시로 들어 직관적인 설명과 수리적 엄밀함을 모두 갖춘 형태로 재구성해 드리겠습니다.

분석 도구로는 jamovi의 사용법을 설명하되, jamovi의 기반이 되는 R 코드를 함께 제시하여 모의 데이터 생성부터 분석, 시각화까지 완벽하게 구현해 드리겠습니다.

---

1. 들어가며: 왜 다층모형인가?

전통적인 회귀분석은 모든 학생이 서로 독립적이라고 가정합니다. 하지만 교육 현장은 그렇지 않습니다. 같은 학교, 같은 반 학생들은 급훈, 담임 선생님, 학교 분위기 등을 공유합니다. 이를 위계적 구조(Hierarchical Structure) 또는 집단 의존성(Group Dependencies)이라고 합니다.

인과추론(Causal Inference)의 관점에서, 데이터가 이러한 계층 구조를 가질 때 다층모형(Multilevel Model)을 사용하는 것은 단순한 통계적 선호가 아니라, 편향(Bias)을 줄이고 정확한 표준오차를 추정하기 위한 필수 전략입니다.

---

2. 인과추론의 기초 개념과 교육적 예시

본격적인 분석에 앞서, 인과추론의 핵심 개념을 학교 상황에 빗대어 정의해 봅시다.

2.1 잠재적 결과 (Potential Outcomes)

어떤 학생 철수($i$)가 있습니다.

• $y_i(1)$: 철수가 ‘방과후 보충수업($z=1$)’을 들었을 때의 성적
• $y_i(0)$: 철수가 ‘방과후 보충수업($z=0$)’을 듣지 않았을 때의 성적

인과 효과(Causal Effect)는 이 둘의 차이 $y_i(1) – y_i(0)$입니다. 하지만 현실에서 우리는 철수가 수업을 듣거나, 듣지 않거나 둘 중 하나의 결과만 볼 수 있습니다. 이를 “인과추론의 근본적인 문제(Missing Data Problem)”라고 합니다.

2.2 SUTVA (Stable Unit Treatment Value Assumption)

이 가정은 “철수가 보충수업을 받았는지 여부가, 옆 짝꿍 영희의 성적에 영향을 주지 않아야 한다(상호간섭 없음)”는 것입니다.

• 문제점: 학교에서는 이 가정이 자주 깨집니다. 철수가 보충수업에서 배운 내용을 영희에게 알려줄 수 있기 때문입니다. 이를 해결하기 위해 집단(학교/학급) 단위 무선화가 권장되기도 합니다.

---

3. 연구 설계에 따른 다층모형 적용

3.1 무선화 실험 (Randomized Experiments)

가장 이상적인 상황입니다. 처치(Treatment)가 무작위로 배정되면, 평균적으로 두 집단은 성향이 비슷해집니다($E[y_i(1)] = E[y_i|z_i=1]$).

A. 개인 단위 무선배정 (학생별 제비뽑기)

학생들에게 무작위로 새로운 ‘독서 프로그램’을 배정했습니다. 하지만 학생들은 학교($j$)라는 집단에 속해 있습니다. 학교마다 평균 독서 능력이 다를 수 있으므로, 이를 반영한 다층모형(Random Intercept Model)이 필요합니다.

$y_{ij} = \mu + \alpha_j + \tau z_{ij} + \epsilon_{ij}$

• $\alpha_j$: $j$번째 학교의 고유한 특성(학교 효과, 랜덤 절편)
• $\tau$: 독서 프로그램의 효과 (우리가 알고 싶은 값)

이 모형을 쓰면 학교 간 차이($\alpha_j$)를 통제하고 순수한 프로그램 효과($\tau$)를 더 정밀하게 추정할 수 있습니다.

B. 집단 단위 무선배정 (학교별 제비뽑기)

교육 정책 연구에서는 흔히 “A학교는 실험군, B학교는 대조군”으로 배정합니다. 이를 군집 무선화(Cluster Randomized Experiments)라고 합니다.

• 이유: ‘학교 폭력 예방 캠페인’처럼 학교 전체 분위기를 바꾸는 처치는 학생 개인별로 쪼개서 적용할 수 없기 때문입니다.
• 분석: 처치 변수($z_j$)가 학생 수준($i$)이 아닌 학교 수준($j$)에 들어갑니다.

3.2 관찰 연구 (Observational Studies)

현실적으로 무선 배정이 불가능할 때(예: 사립학교 진학 효과), 우리는 무시가능성(Ignorability) 가정을 도입합니다. 즉, “부모의 소득, 지능 등 공변량($x$)이 같다면, 사립학교와 공립학교 학생은 비교 가능하다”고 가정하는 것입니다.

• 성향점수(Propensity Score) 활용: 다층 구조에서는 성향점수를 추정할 때도 다층모형을 사용하는 것이 좋습니다.

---

4. [실습] jamovi & R을 활용한 다층 인과 분석

이제 가상의 시나리오를 통해 실제 데이터를 생성하고 분석해 보겠습니다.

4.1 시나리오: “아침 독서 마라톤” 효과 분석

연구 배경: 경기도 교육청은 초등학생의 어휘력 향상을 위해 매일 아침 20분간 책을 읽는 ‘아침 독서 마라톤’ 프로그램을 개발했습니다.

연구 설계:

• 총 20개 학교, 학교당 30명의 학생(총 600명).
• 군집 무선화(Cluster RCT): 학교 단위로 제비뽑기를 하여 10개 학교는 ‘프로그램 시행(Treatment)’, 10개 학교는 ‘기존 자습(Control)’을 하도록 했습니다.
• 데이터 구조:
  • Level 1: 학생 (사후 어휘력 점수 score)
  • Level 2: 학교 (school_id)
  • 처치: program (1=시행, 0=미시행)

4.2 R을 이용한 모의 데이터 생성

jamovi는 R 기반이므로, 아래 코드로 데이터를 생성하여 CSV로 저장한 뒤 jamovi에서 불러오면 됩니다.

R

# 필수 라이브러리 로드
library(lme4)
library(tidyverse)

set.seed(2026) # 재현성을 위한 시드 설정

# 1. 파라미터 설정
n_schools <- 20       # 학교 수
n_students <- 30      # 학교당 학생 수
n_total <- n_schools * n_students

# 2. 학교 수준 효과 (Level 2)
# 학교마다 평균 어휘력이 다름 (표준편차 5)
school_intercept <- rnorm(n_schools, mean = 0, sd = 5)

# 처치 배정 (학교 단위 무선화)
# 1~10번 학교: 통제군(0), 11~20번 학교: 실험군(1)
school_treatment <- c(rep(0, 10), rep(1, 10))

# 학교 데이터 프레임
school_data <- data.frame(
  school_id = 1:n_schools,
  school_eff = school_intercept,
  program = school_treatment
)

# 3. 학생 수준 데이터 생성 (Level 1)
data <- data.frame(
  student_id = 1:n_total,
  school_id = rep(1:n_schools, each = n_students)
)

# 학교 정보 병합
data <- left_join(data, school_data, by = "school_id")

# 4. 결과 변수 생성 (어휘력 점수)
# 기본 점수 70점 + 프로그램 효과 8점 + 학교 효과 + 개인 오차(sd=8)
# y_ij = 70 + 8 * z_j + u_j + e_ij
data <- data %>%
  mutate(
    error = rnorm(n_total, mean = 0, sd = 8),
    score = 70 + 8 * program + school_eff + error
  )

# 팩터 변환
data$school_id <- as.factor(data$school_id)
data$program <- factor(data$program, levels = c(0, 1), labels = c("Control", "Treatment"))

# 데이터 확인
head(data)

4.3 jamovi 분석 절차

이 데이터는 학교 간 차이(School Effect)가 존재하고, 처치가 학교 단위로 부여되었으므로 다층모형(Linear Mixed Model)을 사용해야 정확합니다.

Step 1: 데이터 탐색 및 시각화

분석 전에 데이터의 구조를 눈으로 확인해야 합니다.

1. jamovi 메뉴: Exploration > Descriptives
2. Variables에 score를 넣고, Split by에 program을 넣습니다.
3. Box Plot: 학교별 차이를 보기 위해 Box plot을 체크하고, X축에 program을 둡니다. (※ jamovi 기본 기능으로는 학교별 boxplot을 한 번에 그리기 어려우므로 R 모듈인 seolmatrix나 scatr 모듈을 설치하여 시각화하면 좋습니다.)

[R 시각화 코드]

R

# 학교별 점수 분포 시각화
ggplot(data, aes(x = school_id, y = score, fill = program)) +
  geom_boxplot() +
  theme_minimal() +
  labs(title = "학교별 어휘력 점수 분포", y = "어휘력 점수", x = "학교 ID")

해석: 상자 그림을 보면 같은 처치 집단 내에서도 학교마다 점수의 높낮이가 다름을 알 수 있습니다. 이것이 바로 $\alpha_j$ (학교 효과)입니다.

Step 2: 다층모형 분석 (Linear Mixed Models)

1. 모듈 선택: 상단 메뉴에서 Linear Models > Mixed Model을 클릭합니다. (보이지 않으면 jamovi Library에서 GAMLj 모듈을 설치하는 것을 강력 추천합니다. 여기서는 기본 Mixed Model 기준으로 설명합니다.)
2. 변수 설정:
   • Dependent Variable (종속변수): score
   • Covariates (공변량) 또는 Factors: program (처치 변수)
   • Cluster (군집 변수): school_id
3. Random Effects (랜덤 효과) 설정:
   • 왼쪽의 program을 오른쪽으로 옮기지 않고, Intercept만 Random Coefficients에 둡니다. (기본적으로 (Intercept | school_id)로 설정됨)
   • 이는 학교마다 평균 점수(절편)가 다름을 허용하는 것입니다.
4. Fixed Effects (고정 효과) 설정:
   • program을 Model Terms에 넣습니다. 이것이 우리가 알고 싶은 ‘독서 마라톤 효과’입니다.

Step 3: 결과 해석

jamovi의 결과표(Estimates)는 다음과 유사하게 나옵니다.

Effect	Estimate	SE	t	p
Intercept	72.574	0.972	74.681	< .001
program (Treatment)	10.162	1.944	5.228	< .001

• Fixed Effects: program의 Estimate가 약 10.162입니다. 즉, 독서 마라톤을 한 학교 학생들이 하지 않은 학교보다 평균적으로 약 10.162점 더 높은 어휘력을 보입니다. $p < .05$이므로 통계적으로 유의합니다.
• Random Components (Variance):
  • $\sigma^2_{school}$ (School Intercept): 학교 간 분산. 이 값이 0보다 크다면 학교 효과가 존재한다는 뜻입니다.
  • ICC (Intraclass Correlation Coefficient): 전체 분산 중 학교가 설명하는 비율입니다.

---

5. 심화: 불응(Noncompliance)과 도구변수(IV)

실험을 했는데, 독서 프로그램을 하라고 배정받은 학교의 일부 학생이 땡땡이를 쳤다면(Noncompliance) 어떻게 될까요? 이때는 “배정된 상태($z$)”를 도구변수(Instrument)로 사용하여, 실제 “참여한 상태($d$)”의 효과를 추정해야 합니다.

jamovi/R 구현 (2단계 최소자승법 개념)

1. 1단계: 실제 참여 여부($d$)를 배정 여부($z$)로 예측합니다.
   $d_{ij} = \gamma_j + \delta z_{ij} + \nu_{ij}$
2. 2단계: 1단계에서 예측된 참여값($\hat{d}$)을 사용하여 점수($y$)를 예측합니다.
   $y_{ij} = \alpha_j + \tau \hat{d}_{ij} + \epsilon_{ij}$

이 분석은 jamovi의 sem (구조방정식) 모듈이나 R의 AER 패키지(ivreg)를 통해 수행할 수 있습니다. 중요한 건 배정($z$)은 오직 참여($d$)를 통해서만 결과($y$)에 영향을 미쳐야 한다(배제 제한)는 가정입니다.

---

6. 결론

다층모형을 활용한 인과추론은 교육 현장과 같이 “집단 속에 개인이 속한 데이터”를 분석할 때 가장 강력한 도구입니다.

1. 설계: 가능하다면 학교 단위 무선화(Cluster RCT)가 상호간섭(SUTVA 위배) 문제를 피하는 데 유리합니다.
2. 분석: 단순히 평균을 비교하는 t-test 대신, 학교의 무선 절편(Random Intercept)을 포함한 혼합 모형을 사용해야 표준오차의 과소추정을 막을 수 있습니다.
3. 해석: 결과는 “개인 수준의 효과”인지 “학교 수준의 효과”인지 명확히 구분하여 해석해야 합니다.

이 장의 내용이 여러분의 연구에 튼튼한 방법론적 기초가 되기를 바랍니다.

---

참고문헌 (APA Style)

• Almond, D., Chay, K., & Lee, D. (2005). The costs of low birth weight. The Quarterly Journal of Economics, 120(3), 1031-1083.
• Angrist, J. D., Imbens, G. W., & Rubin, D. B. (1996). Identification of causal effects using instrumental variables. Journal of the American Statistical Association, 91(434), 444-472.
• Cornfield, J. (1978). Randomization by group: A formal analysis. American Journal of Epidemiology, 108(2), 100-102.
• Gelman, A., & Hill, J. (2007). Data analysis using regression and multilevel/hierarchical models. Cambridge University Press.
• Hill, J. (2013). Multilevel models and causal inference. In The SAGE Handbook of Multilevel Modeling (Chapter 12, pp. 201-219).
• Hong, G., & Raudenbush, S. W. (2006). Evaluating kindergarten retention policy: A case study of causal inference for multilevel observational data. Journal of the American Statistical Association, 101(475), 901-910.
• Kim, J., & Seltzer, M. (2007). Causal inference in multilevel settings in which selection processes vary across schools (Tech. Rep.). CRESST, UCLA.
• Rubin, D. B. (1978). Bayesian inference for causal effects: The role of randomization. The Annals of Statistics, 6(1), 34-58.
• Rubin, D. B. (1990). Formal modes of statistical inference for causal effects. Journal of Statistical Planning and Inference, 25(3), 279-292.
• Slavin, R. E., Madden, N. A., Dolan, L. J., & Wasik, B. A. (1996). Every child, every school: Success for all. Corwin Press.

안녕하세요!

오늘은 “다층 연구 설계의 최적화(Sample Size and Power Analysis in Multilevel Designs)”에 대해 살펴보겠습니다. “학교 현장의 데이터”를 예시로 들어 직관적인 설명과 수리적 엄밀함을 모두 갖춘 형태로 재구성해 드리겠습니다.

분석 도구로는 jamovi의 기반이 되는 R 코드를 함께 제시하여 ‘최적 표본 크기 산출’부터 ‘모의 데이터 생성’, 그리고 ‘분석’까지 완벽하게 구현해 드리겠습니다. (참고: jamovi는 데이터 분석에는 강력하지만, 연구 설계 단계의 복잡한 ‘최적 표본 산출’ 기능은 제한적이므로, 이 부분은 R의 수리적 계산 기능을 활용하고 분석은 jamovi의 혼합모형 논리로 설명하겠습니다.)

---

1. 서론: 연구자의 영원한 딜레마, “돈이냐, 정확성이냐?”

교육 연구를 진행할 때 우리는 항상 두 가지 제약 조건 사이에서 고민합니다.

1. 통계적 검증력(Power): 효과가 있다면 있다고 말할 수 있는 힘 (높을수록 좋음).
2. 예산(Budget): 연구비와 시간은 한정되어 있음.

이 장에서는 다층 모형(Multilevel Modeling) 상황, 즉 학생이 학교에 소속되어 있는 구조에서 어떻게 하면 가장 적은 비용으로 가장 높은 정확성을 얻을 수 있는지, 그 최적 설계(Optimal Design) 방법을 알려드리겠습니다.

우리의 가상 시나리오를 소개합니다.

[시나리오: 프로젝트 ‘독서왕’]

교육심리학자 김 교수는 새로운 독서 프로그램이 초등학생의 ‘문해력’을 높이는지 검증하려 합니다.

• 총 예산: 1,000만 원 (가상의 화폐 단위)
• 비용 구조:
  • 학교 하나를 섭외하는 비용(행정 절차, 학교 보상 등): 200만 원 (c)
  • 학생 한 명을 검사하는 비용(검사지, 간식 등): 10만 원 (s)
• 연구 질문: 몇 개의 학교를 섭외하고, 학교당 몇 명의 학생을 뽑아야 내 돈 1,000만 원 안에서 가장 정확한 결과를 얻을까요?

---

2. 군집 무작위 배정(Cluster Randomized Trial)의 설계

2.1 왜 다층 설계인가?

가장 쉬운 방법은 전국의 학생 명부에서 무작위로 학생을 뽑는 것입니다. 하지만 현실적으로 불가능합니다.

• 실행 가능성: 학교 단위로 프로그램을 돌려야 합니다.
• 오염(Contamination): 한 반에서 철수는 실험집단, 영희는 통제집단이면 서로 이야기하며 효과가 섞여버립니다.

그래서 우리는 학교(Cluster)를 통째로 실험군 혹은 대조군으로 배정하는 군집 무작위 배정(Cluster Randomized Trial)을 사용합니다.

2.2 급내상관계수(ICC)와 설계 효과

문제는 같은 학교 아이들끼리는 서로 비슷하다는 점입니다(학교 분위기, 선생님의 영향 등). 이를 급내상관계수(ICC, $\rho$)라고 합니다.

• ICC가 높다 = 학교 간 차이가 크다 = 같은 학교 아이들은 매우 비슷하다.
• ICC가 높으면, 학생을 100명 더 뽑는 것보다 학교를 1개 더 섭외하는 게 훨씬 중요해집니다.

2.3 최적 표본 크기 공식 (The Magic Formula)

주어진 예산($B$) 하에서 처치 효과($\beta_1$)의 분산($Var(\hat{\beta}_1)$)을 최소화하는 최적의 학교 수($K$)와 학교당 학생 수($n$)는 다음과 같습니다.

$n_{optimal} = \sqrt{\frac{c(1-\rho)}{s\rho}}$

$K_{optimal} = \frac{B}{c + s \cdot n}$

• $c$: 학교당 비용 (200)
• $s$: 학생당 비용 (10)
• $\rho$: ICC (가정된 값, 보통 0.05~0.10)

이 공식을 보면, ICC($\rho$)가 커질수록 학교당 학생 수($n$)는 줄여야 합니다. 왜냐하면 같은 학교에서 많이 뽑아봤자 정보가 중복되기 때문입니다.

---

3. R과 jamovi를 활용한 최적 설계 및 데이터 생성

이제 김 교수의 ‘독서왕’ 프로젝트를 위해 R을 사용하여 최적 표본을 계산하고, 이를 분석할 수 있는 모의 데이터를 생성해 보겠습니다.

3.1 최적 표본 크기 계산 (R Code)

김 교수의 상황: 예산 10,000, $c=200, s=10$, 그리고 ICC($\rho$)는 선행연구를 통해 0.05로 가정합니다.

R

# [R Code] 최적 표본 크기 산출
# 파라미터 설정
B <- 10000   # 총 예산
c <- 200     # 학교(Cluster)당 비용
s <- 10      # 학생(Person)당 비용
rho <- 0.05  # 급내상관계수 (ICC)

# 1. 최적의 학생 수 (n) 계산 [cite: 120]
n_opt <- sqrt((c * (1 - rho)) / (s * rho))

# 2. 최적의 학교 수 (K) 계산 [cite: 119]
K_opt <- B / (c + s * n_opt)

# 결과 출력
cat("최적의 학교당 학생 수 (n):", round(n_opt, 2), "명\n")
cat("최적의 학교 수 (K):", round(K_opt, 2), "개교\n")

[분석 결과 해석]

• 계산 결과, $n \approx 19.49$, $K \approx 25.32$가 나옵니다.
• 현실적으로 반올림하여 학교당 19명, 총 26개 학교(실험 13, 통제 13)를 섭외하는 것이 최적입니다.

3.2 비용 효율성 시각화 (Design Efficiency Plot)

아래 그래프는 학교당 학생 수($n$)를 변화시킬 때, 추정의 오차(분산)가 어떻게 변하는지 보여줍니다. 우리는 분산이 가장 낮은 지점을 찾아야 합니다.

R

# [R Code] 시각화 생성
library(ggplot2)

n_seq <- seq(5, 50, by = 1) # 학생 수를 5명에서 50명까지 변화시킴
K_seq <- B / (c + s * n_seq) # 예산 제약에 따른 학교 수

# 분산 계산 공식 (단순화된 형태) [cite: 128]
# g_rho 부분과 예산 부분을 결합
design_var <- function(n, K, rho) {
  # Standard Error calculation based on Eq 11.3 & 11.6 logic relation
  # Here we look at relative variance proportional to the function
  design_effect <- 1 + (n - 1) * rho
  total_N <- n * K
  return(design_effect / total_N) 
}

var_values <- design_var(n_seq, K_seq, rho)

df_plot <- data.frame(n = n_seq, Variance = var_values)

ggplot(df_plot, aes(x = n, y = Variance)) +
  geom_line(color = "blue", linewidth = 1.2) +
  geom_vline(xintercept = 19, linetype = "dashed", color = "red") +
  annotate("text", x = 19, y = max(var_values), label = "Optimal n = 19", vjust = -1) +
  labs(title = "최적 표본 설계: 학교당 학생 수 vs. 추정 오차",
       subtitle = "예산 고정 시, 학생을 19명 뽑을 때 오차가 최소화됨",
       x = "학교당 학생 수 (n)", y = "치료 효과의 분산 (낮을수록 좋음)") +
  theme_minimal()

이 그래프를 통해 김 교수는 무작정 학생을 많이 뽑는다고 좋은 게 아니라, 학교 수와 학생 수의 황금 비율을 맞춰야 함을 알 수 있습니다.

---

4. 모의 데이터 생성 및 분석 (Linear Mixed Model)

이제 최적 설계($K=26, n=19$)에 따라 데이터를 수집했다고 가정하고, 이를 jamovi(또는 R의 lmer)에서 분석하는 방법을 보여드리겠습니다.

4.1 데이터 생성 (Backstory 포함)

• 상황: 26개 학교를 선정, 절반은 ‘독서왕 프로그램(Treatment)’, 절반은 ‘기존 수업(Control)’ 진행.
• 효과: 프로그램은 문해력 점수를 평균 5점 정도 높여줄 것으로 기대.
• 변산: 학교 간 차이(ICC) 반영.

R

# [R Code] 모의 데이터 생성
set.seed(123) # 재현성을 위해 시드 설정

K <- 26       # 학교 수
n <- 19       # 학교당 학생 수
N <- K * n    # 총 학생 수

# 학교 ID 및 치료 집단 배정 (0: 통제, 1: 처치)
school_id <- rep(1:K, each = n)
treatment <- rep(c(rep(0, K/2), rep(1, K/2)), each = n)

# 랜덤 효과 생성 (학교 간 차이)
u0j <- rep(rnorm(K, mean = 0, sd = sqrt(5)), each = n) # 학교 분산 = 5

# 오차항 생성 (학생 간 차이)
# ICC = 0.05 이려면, 학교분산/(학교분산+오차분산) = 0.05
# 5 / (5 + 95) = 0.05 -> 오차 분산은 95, SD는 약 9.75
eij <- rnorm(N, mean = 0, sd = sqrt(95))

# 고정 효과 (진짜 치료 효과 = 5점)
beta0 <- 50  # 평균 점수
beta1 <- 5   # 치료 효과

# 종속 변수 (문해력 점수) 생성 
# y_ij = beta0 + beta1*x_j + u_0j + e_ij
y <- beta0 + beta1 * treatment + u0j + eij

# 데이터 프레임 생성
data_sim <- data.frame(
  SchoolID = factor(school_id),
  StudentID = 1:N,
  Treatment = factor(treatment, labels = c("Control", "Program")),
  Score = y
)

# 데이터 확인
head(data_sim)

# CSV로 저장 (jamovi에서 불러오기 위함)
write.csv(data_sim, "chap11.csv", row.names = FALSE)

4.2 jamovi 식 분석 절차 (R lme4 문법 활용)

jamovi의 Linear Mixed Models (GAMLj 모듈) 메뉴에서 분석하는 것과 동일한 R 코드입니다.

데이터 구조는 다음과 같습니다.

• Level 1 (Person): 학생별 문해력 점수 (Score)
• Level 2 (Cluster): 학교 (SchoolID)
• Predictor: 처치 여부 (Treatment, Level 2 변수)

수식 모형은 다음과 같습니다15:

$y_{ij} = \beta_{0} + \beta_{1}Treatment_{j} + u_{0j} + \epsilon_{ij}$

R

# [R Code] 다층 모형 분석
library(lme4)
library(lmerTest)

# 모형 적합: 절편에 대해서만 학교별 무선 효과(Random Intercept) 허용
model <- lmer(Score ~ Treatment + (1 | SchoolID), data = data_sim)

# 결과 요약
summary(model)

[결과 해석 방법]

1. Fixed Effects: TreatmentProgram의 Estimate가 5에 가까운지, p-value가 0.05보다 작은지 확인합니다. (우리가 5로 설정했으므로 유의하게 나올 것입니다.)
2. Random Effects: SchoolID의 Variance가 설계 시 가정한 분산과 유사한지 확인합니다.

---

5. 심화: 더 복잡한 상황들

5.1 다기관 임상시험 (Multisite Trials)

만약 김 교수가 학교 전체를 배정하는 게 아니라, 각 학교 안에서 철수는 실험반, 영희는 통제반으로 나눌 수 있다면 어떨까요?

이를 Multisite Trial이라고 합니다.

• 장점: 학교 효과($u_{0j}$)가 치료 효과 추정에서 사라지므로 훨씬 강력한 검증력(Power)을 가집니다.
• 설계: 이 경우 $Var(\hat{\beta}_1)$은 학교 분산에 영향을 받지 않으므로, 더 적은 예산으로도 유의한 결과를 얻을 수 있습니다. 하지만 ‘오염’ 문제가 없어야만 가능합니다.

5.2 반복 측정 (Longitudinal Design)

만약 김 교수가 프로그램을 1년 동안 진행하면서 학생들의 변화를 보고 싶다면 몇 번 측정해야 할까요?

• 선형 변화(직선 성장)를 가정할 때: 시작(Baseline)과 끝(End), 딱 2번 측정하거나, 중간 지점 하나를 추가하는 것이 비용 대비 가장 효율적입니다.
• 이차 함수(곡선 성장)를 가정할 때: 최소 3번의 측정이 필요하며, 등간격으로 측정하는 것이 좋습니다.

5.3 현실적인 문제와 해결책

1. ICC를 모를 때: 보통 문헌 연구를 통해 보수적으로(약간 높게) 잡습니다. ICC를 실제보다 절반 정도로 낮게 잘못 예측했더라도, 최적 설계 대비 효율성 손실은 약 10% 내외로 크지 않다는 연구가 있습니다.
2. 학교 크기가 다를 때: 모든 학교가 19명일 수는 없습니다. 학교 간 크기 편차(CV)가 0.5 정도라면, 계산된 학교 수($K$)보다 약 11% 정도 더 많은 학교를 섭외하여 이를 보정해야 합니다.

---

6. 결론 및 제언

오늘 살펴본 내용을 요약하면 다음과 같습니다.

1. 학교 현장 연구에서는 학생 수만큼이나 학교 수(Cluster Number)가 중요하다.
2. 비용 함수와 ICC를 고려하면, 무조건 많은 표본보다 최적의 비율을 찾는 것이 경제적이다.
3. 학교를 통째로 배정(CRT)하는 것보다 학교 내 무선 배정(Multisite)이 통계적으로는 더 유리하나, 오염 가능성을 고려해야 한다.

---

참고문헌 (APA Style)

• Maas, C. J. M., & Hox, J. J. (2005). Sufficient sample sizes for multilevel modeling. Methodology, 1(3), 86-92.
• Moerbeek, M., & Teerenstra, S. (2011). Optimal design in multilevel experiments. In J. J. Hox & J. K. Roberts (Eds.), Handbook of Advanced Multilevel Analysis (pp. 257-281). New York: Routledge.
• Raudenbush, S. W. (1997). Statistical analysis and optimal design for cluster randomized trials. Psychological Methods, 2(2), 173-185.
• Snijders, T. A. B., & Bosker, R. J. (1993). Standard errors and sample sizes for two-level research. Journal of Educational Statistics, 18(3), 237-259.
• Van Breukelen, G. J. P., & Moerbeek, M. (2013). Design considerations in multilevel studies. In M. A. Scott, J. S. Simonoff, & B. D. Marx (Eds.), The SAGE Handbook of Multilevel Modeling (pp. 183-200). SAGE Publications.

안녕하세요!

오늘은 개체 내 오차 구조(Within-Individual Error Structures)의 복잡성과 모델링에 대해 살펴보겠습니다. “학교 현장의 데이터”를 예시로 들어 직관적인 설명과 수리적 엄밀함을 모두 갖춘 형태로 재구성해 드리겠습니다.

분석 도구로는 jamovi의 사용법을 설명하되, jamovi의 기반이 되는 R 코드를 함께 제시하여 모의 데이터 생성부터 분석, 시각화까지 완벽하게 구현해 드리겠습니다.

---

우리가 학교에서 학생들의 성장을 추적할 때(종단 연구), 단순히 “평균 점수가 올랐나?”만 보는 것은 반쪽짜리 분석입니다. 학생 개개인이 어떻게 변화하는지, 그 변화의 폭(분산)은 일정한지, 어제의 성적이 오늘의 성적에 얼마나 영향을 미치는지(상관)를 파악해야 합니다. 제공해주신 텍스트는 이러한 공분산 구조(Covariance Structure)를 어떻게 모델링할 것인가에 대한 깊이 있는 통찰을 제공합니다.

1. 교육 현장의 예시: “읽기 유창성 성장 프로젝트”

이해를 돕기 위해 가상의 시나리오를 설정하겠습니다.

상황: A 초등학교에서 3학년 학생 50명을 대상으로 ‘읽기 유창성(1분당 읽은 단어 수)’을 1년 동안 4회(3월, 6월, 9월, 12월) 측정했습니다.

핵심 질문:

1. 모든 학생의 읽기 실력 격차(분산)는 3월이나 12월이나 똑같을까요? (등분산성)
2. 3월 성적이 좋은 학생은 6월에도 좋을까요? 12월까지 그 영향이 갈까요? (계열 상관)

---

2. 기본 모델과 개념 (The General Model)

우리는 데이터를 다음의 선형 회귀 모델로 표현할 수 있습니다.

$y = X\beta + Zu + \epsilon$

• $y$: 학생들의 읽기 점수 벡터
• $X\beta$ (고정 효과): 전체 학생들의 평균적인 성장 곡선 (예: 시간이 지날수록 점수가 오른다).
• $Zu$ (임의 효과): 학생 개인별 특성 (예: 어떤 학생은 시작부터 잘하고, 어떤 학생은 성장 속도가 빠르다).
• $\epsilon$ (오차 항): 설명되지 않는 나머지 변동. 오늘의 핵심 주제는 바로 이 $\epsilon$의 구조인 $\Sigma$를 파헤치는 것입니다.

---

3. 데이터 시각화와 진단 (Diagnostics)

모델링을 하기 전에 눈으로 확인해야 합니다. 텍스트에서는 프로파일 도표(Profile Plot), OSM, PRISM을 추천합니다.

3.1 프로파일 도표 (Profile Plot)

학생 개개인의 성장 궤적을 그린 그래프입니다.

• 해석: 선들이 서로 꼬이지 않고 나란히 간다면? $\rightarrow$ 학생 내 상관이 높음(잘하던 애가 계속 잘함).
• 분산: 시간이 갈수록 선들의 폭이 넓어진다면? $\rightarrow$ 이분산성(Heterogeneity) 존재.

3.2 OSM (Ordinary Scatterplot Matrix)

모든 시점 간의 산점도 행렬입니다.

• 특징: 대각선에서 멀어질수록(시간 격차가 클수록) 상관관계가 어떻게 변하는지 보여줍니다. 읽기 데이터에서는 보통 시간이 지날수록(대각선에서 멀어질수록) 상관이 낮아지는 패턴을 보입니다.

3.3 PRISM (Partial Regression-on-Intervenors Scatterplot Matrix)

이것은 저자(Zimmerman)가 강조하는 고급 도구입니다.

• 개념: $t_1$과 $t_3$의 관계를 볼 때, 중간 시점인 $t_2$의 영향을 제거하고(통제하고) 봅니다.
• 해석: 만약 $t_2$를 통제했을 때 $t_1$과 $t_3$의 관계가 사라진다면(무작위 산포), 이는 “조건부 독립”을 의미하며, Antedependence 모델이나 AR(1) 모델이 적합하다는 강력한 신호입니다.

---

4. 공분산 구조 모델링: 유연성 vs 간명성

어떤 모델을 선택해야 할까요? 너무 단순하면 현실을 반영 못 하고(과소적합), 너무 복잡하면 추정이 어렵습니다(과대적합). 이 둘 사이의 균형(Balance)이 핵심입니다.

4.1 가장 단순한 모델들 (Parsimonious)

1. 복합 대칭 (Compound Symmetry, CS):
   • 가정: 3월의 분산 = 12월의 분산. 3월-6월 상관 = 3월-12월 상관. 모든 것이 일정함.
   • 교육적 현실: 글쎄요, 보통 고학년이 될수록 격차(분산)가 커지고, 먼 시점 간의 상관은 떨어지기 때문에 현실적이지 않을 수 있습니다.
   • 구조: $\sigma^2 \begin{pmatrix} 1 & \rho & \rho & \rho \\ \rho & 1 & \rho & \rho \\ \rho & \rho & 1 & \rho \\ \rho & \rho & \rho & 1 \end{pmatrix}$
2. 1차 자기회귀 (AR(1)):
   • 가정: 바로 전 시험 성적이 다음 시험에 가장 큰 영향을 줌. 시간이 멀어질수록 상관이 지수적으로 감소 ($\rho, \rho^2, \rho^3…$).
   • 교육적 현실: 성장 데이터에 매우 적합합니다. 단, 분산이 일정하다고 가정하는(Homogeneous) 한계가 있습니다.

4.2 중간 단계의 모델들 (Heterogeneous & Banded)

3. 이분산 자기회귀 (Heterogeneous AR, ARH(1)):
   • AR(1)과 같지만, 시점별로 분산(시험 점수의 퍼짐 정도)이 다를 수 있게 허용합니다. 학년이 올라갈수록 실력 격차가 벌어지는 현상을 설명하기 좋습니다.
4. 토플리츠 (Toeplitz, TOEP):
   • AR(1)처럼 시간이 지날수록 상관이 변하지만, 그 패턴을 $\rho^k$처럼 강제로 수학적 공식에 넣지 않고 자유롭게 둡니다(Band 구조). 더 유연하지만 파라미터가 많아집니다.

4.3 가장 유연한 모델들 (Flexible)

5. 무구조 (Unstructured, UN):
   • 가정: 아무런 제약이 없습니다. 모든 시점의 분산과 상관을 각각 다 추정합니다.
   • 장단점: 가장 정확하지만, 측정 시점이 많아지면 추정해야 할 파라미터가 기하급수적으로 늘어나($n(n+1)/2$) 모델이 터질 수 있습니다.
6. Antedependence (AD) 및 Structured AD (SAD) 모델:
   • 텍스트의 핵심 제안: 1차 자기회귀(AR1)는 너무 딱딱하고, 무구조(UN)는 너무 무겁습니다.
   • AD 모델: 시점 간의 “조건부 독립”을 가정하여 파라미터를 줄이면서도, 비정상성(Non-stationarity: 분산과 상관이 변하는 성질)을 허용하는 똑똑한 모델입니다.
   • SAD 모델은 이를 더 구조화하여 파라미터 수를 $O(1)$ 수준으로 줄여 간명성을 확보합니다.

---

5. 분석 실습: R과 jamovi 활용

이제 가상의 “읽기 유창성 데이터”를 생성하고, 위에서 배운 모델들을 비교 분석해 보겠습니다.

(참고: jamovi의 GAMLj 모듈이나 Linear Mixed Models는 CS, AR1, UN, Toeplitz 등을 지원하지만, PRISM 도표나 고급 AD 모델은 R이 필요합니다. 따라서 R 코드를 중심으로 설명하고 jamovi 활용법을 덧붙이겠습니다.)

5.1 R을 이용한 모의 데이터 생성 및 시각화

R

# 필요한 라이브러리 로드
library(MASS)
library(nlme)
library(lattice)
library(ggplot2)

# --- 1. 가상의 읽기 유창성 데이터 생성 (Backstory) ---
# 50명의 학생, 4회 측정 (t=1, 2, 3, 4)
# 가정: 시간이 갈수록 분산이 커지고(이분산), 상관은 멀어질수록 낮아짐(AR 구조)
set.seed(2024)
N <- 50
Time <- 4
Mean_Vector <- c(40, 55, 65, 72) # 시간별 평균 점수 상승
Sigma <- matrix(0, Time, Time)

# 공분산 행렬 생성 (ARH(1) 구조와 유사하게 설정)
# 분산은 100, 140, 180, 220으로 증가
vars <- c(100, 140, 180, 220)
cor_rho <- 0.7 # 1-lag 상관계수
for(i in 1:Time){
  for(j in 1:Time){
    Sigma[i,j] <- sqrt(vars[i]) * sqrt(vars[j]) * (cor_rho^abs(i-j))
  }
}

# 데이터 생성
data_matrix <- mvrnorm(n = N, mu = Mean_Vector, Sigma = Sigma)
colnames(data_matrix) <- c("T1", "T2", "T3", "T4")
df_wide <- data.frame(ID = 1:N, data_matrix)

# Long format으로 변환 (분석용)
df_long <- reshape(df_wide, direction = "long", varying = list(2:5), 
                   v.names = "Score", timevar = "Time", times = 1:4)
df_long <- df_long[order(df_long$ID, df_long$Time), ]

# --- 2. 진단 시각화 ---

# (1) 프로파일 도표 (Profile Plot) [cite: 65]
ggplot(df_long, aes(x = Time, y = Score, group = ID)) +
  geom_line(alpha = 0.5) +
  geom_smooth(aes(group = 1), method = "loess", color = "red", size = 1.5) +
  labs(title = "학생 읽기 유창성 프로파일 도표", subtitle = "시간에 따른 개별 성장 곡선") +
  theme_minimal()

# (2) OSM (Ordinary Scatterplot Matrix) 
splom(df_wide[,2:5], 
      main = "OSM: 시점 간 산점도 행렬",
      pscales = 0, varname.cex = 0.8)

# (3) PRISM (간략화된 개념적 구현) [cite: 127]
# T1과 T3의 관계에서 T2를 통제한 편상관을 시각적으로 확인
res_T1_T2 <- residuals(lm(T1 ~ T2, data = df_wide))
res_T3_T2 <- residuals(lm(T3 ~ T2, data = df_wide))
plot(res_T1_T2, res_T3_T2, 
     main = "PRISM 예시 (T1 vs T3 | T2)", 
     xlab = "T1 | T2 residuals", ylab = "T3 | T2 residuals")
abline(lm(res_T3_T2 ~ res_T1_T2), col="blue")
# 설명: 이 점들이 무작위로 퍼져있으면(기울기가 0에 가까우면) 조건부 독립 -> AD 모델 적합

# --- 3. 모델 적합 및 비교 (gls 함수 사용) ---

# Model 1: 복합 대칭 (Compound Symmetry) - 너무 단순함 [cite: 232]
fit_cs <- gls(Score ~ factor(Time), data = df_long,
              correlation = corCompSymm(form = ~ 1 | ID))

# Model 2: AR(1) - 등분산 가정 [cite: 250]
fit_ar1 <- gls(Score ~ factor(Time), data = df_long,
               correlation = corAR1(form = ~ 1 | ID))

# Model 3: Unstructured (무구조) - 가장 유연함 [cite: 498]
fit_un <- gls(Score ~ factor(Time), data = df_long,
              correlation = corSymm(form = ~ 1 | ID),
              weights = varIdent(form = ~ 1 | Time)) # 이분산 허용

# Model 4: ARH(1) (이분산 AR1) - 현실적 대안 
fit_arh1 <- gls(Score ~ factor(Time), data = df_long,
                correlation = corAR1(form = ~ 1 | ID),
                weights = varIdent(form = ~ 1 | Time))

# --- 4. 모델 선택 (AIC/BIC 비교) [cite: 47, 55] ---
anova(fit_cs, fit_ar1, fit_un, fit_arh1)

5.2 jamovi에서의 분석 절차

jamovi에서는 [Linear Mixed Models] 메뉴를 사용합니다.

1. 데이터 불러오기: 위에서 생성한 데이터를 csv로 저장하여 jamovi에서 엽니다.
2. 분석 설정:
   • Dependent Variable: Score
   • Covariates: Time (또는 Factors로 설정)
   • Cluster: ID
3. Random Effects & Covariance Structure:
   • jamovi 옵션 중 Residual Covariance (혹은 Random Effect의 상관 구조) 설정 탭을 찾습니다.
   • 여기서 Unstructured, AR(1), Compound Symmetry 등을 선택할 수 있습니다.
4. 모델 비교: 각 모델별로 AIC/BIC 값을 확인하여 가장 작은 값을 가진 모델을 선택합니다. (보통 교육 성장 데이터는 ARH(1)이나 Unstructured가 잘 나옵니다).

---

6. 결론 및 제언 (Conclusions)

이 장의 저자들(Núñez-Antón & Zimmerman)은 다음과 같이 결론짓습니다.

1. 평균과 공분산을 함께 보라: 평균적인 성장만 보지 말고, 오차의 구조(변동성, 상관성)를 함께 모델링해야 정확한 추론이 가능합니다.
2. 균형을 찾아라: 무조건 복잡한 모델(Unstructured)이 좋은 것이 아닙니다. 데이터의 특성(비정상성, 계열 상관 등)을 반영하되, 가능한 한 파라미터 수가 적은 모델(예: Structured Antedependence 등)을 찾는 것이 최선입니다.
3. 시각화를 믿어라: 통계 수치(AIC/BIC)도 중요하지만, 프로파일 도표와 PRISM 등을 통해 데이터의 실제 패턴을 눈으로 확인하는 것이 강력한 진단 도구가 됩니다.

[User를 위한 다음 단계]

위의 R 코드를 복사하여 RStudio에서 실행해 보시겠습니까? 생성된 그래프(프로파일 도표, PRISM)를 통해 방금 설명한 “조건부 독립” 개념이 시각적으로 어떻게 나타나는지 직접 확인해 보시면 이해가 훨씬 빠를 것입니다.

---

참고문헌 (References)

• Akaike, H. (1974). A new look at the statistical model identification. IEEE Transactions on Automatic Control, 19, 716-723.
• Diggle, P. J., Heagerty, P. J., Liang, K.-Y., & Zeger, S. L. (2002). Analysis of longitudinal data (2nd ed.). Oxford University Press.
• Gabriel, K. R. (1962). Ante-dependence analysis of an ordered set of variables. Annals of Mathematical Statistics, 33, 201-212.
• Goldstein, H. (2003). Multilevel statistical model (3rd ed.). Arnold.
• Kenward, M. G. (1987). A method for comparing profiles of repeated measurements. Applied Statistics, 36, 296-308.
• Laird, N. M., & Ware, J. H. (1982). Random-effects models for longitudinal data. Biometrics, 38, 963-974.
• Núñez-Antón, V., & Zimmerman, D. L. (2001). Parametric modeling of growth curve data: An overview (with discussion). Test, 10, 1-73.
• Zimmerman, D. L. (2000). Viewing the correlation structure in longitudinal data through a PRISM. The American Statistician, 54, 310-318.
• Zimmerman, D. L., & Núñez-Antón, V. (2010). Antedependence models for longitudinal data. Chapman & Hall/CRC.