태그 보관물: jamovi

Chap33. 사회연결망과 관계형 데이터의 다층모형(Multilevel Modeling of SocialNetwork and Relational Data)

안녕하세요!

오늘은 사회연결망 및 관계형 데이터의 다층모형(Multilevel Modeling of Social Network and Relational Data)에 대해 살펴보겠습니다. “학교 현장의 데이터”를 예시로 들어 직관적인 설명과 수리적 엄밀함을 모두 갖춘 형태로 재구성해 드리겠습니다.

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

1. 사회연결망 데이터란 무엇일까요?

초등학생 친구들이 교실에서 서로 친하게 지내는 모습을 상상해 볼까요? 어떤 친구는 많은 친구들에게 인기가 있고, 어떤 친구는 단짝 친구 한 명과 아주 깊은 관계를 맺습니다.

사회연결망 분석(Social Network Analysis)은 바로 이러한 사회적 연결망의 구조를 설명하고 시각화하며, 수학적 및 통계적으로 분석하여 이해하는 것을 목표로 합니다. 여기서 학생 한 명 한 명을 노드(node, vertex, actor)라고 부르고 , 학생들 사이의 관계(예: “나는 지훈이를 좋아해”)를 연결선(edge, link, tie)이라고 부릅니다.

  • 방향성 유무: 연결선은 단순히 두 노드 사이의 연결만 기록하는 무방향(undirected)일 수도 있고, 누가 누구에게 관계를 보고했는지(송신자와 수신자) 알 수 있는 방향성(directed)을 가질 수도 있습니다.
  • 인접 행렬(Adjacency matrix): 이러한 네트워크 데이터는 정방행렬인 YY로 표현할 수 있으며, 여기서 YijY_{ij}는 행위자 i로부터 행위자 j로 향하는 관계를 나타냅니다. 무방향 네트워크 데이터의 경우 이 행렬은 대칭을 이룹니다.

데이터를 수집하는 방식에 따라 크게 세 가지로 나눌 수 있습니다.

  • 개인 중심 연결망(Personal or Egocentric Network): 표본으로 추출된 개인이 하나 이상의 타인(alter)과의 관계에 대해 보고하는 데이터입니다. (예: “지훈아, 네 친한 친구들을 다 적어볼래?”)
  • 양자 데이터(Dyadic Data): 두 명의 행위자와 그들 사이의 관계인 양자(dyad), 즉 (Yij,Yji)(Y_{ij}, Y_{ji})를 의미합니다.
  • 전체 연결망(Complete or Whole-network): 닫힌 집단(예: 한 학급 전체) 내의 모든 개인들 사이의 관계를 라운드 로빈(round robin) 방식으로 수집한 데이터입니다.

2. 햇살 초등학교 5학년 1반의 모의 데이터 스토리

Lotte Vermeij(2006)가 수집한 실제 네덜란드 교실 소셜 네트워크 데이터를 바탕으로, 우리만의 모의 데이터를 상상해 봅시다.

  • 목적: 학생들의 성별, 학업 성취도(marks), 학급 소속감(classroom identification)이 학생들 간의 관계(tie strength)에 어떤 영향을 미치는지 알아봅니다.
  • 결과 변수: 친구를 얼마나 지지하고 소통하는지를 합산하여 만든 ‘관계 강도(tie strength)’ 변수를 사용합니다.

3. 단계별 다층모형 분석 (점점 복잡해지는 친구 관계 모델링)

단계 1: 개인 중심 연결망 데이터의 다층모형 (나와 내 친구들)

이 모델에서는 응답자인 ‘자아(ego)’ 안에 그들이 지목한 ‘타인(alter)’들이 내재(nested)되어 있다고 봅니다. 즉, 1수준은 ‘타인(친척, 친구 등)’, 2수준은 ‘자아(응답자 본인)’가 됩니다.

1수준(타인 수준)의 회귀 방정식은 다음과 같습니다.

yij=β0j+x1ijβ1j+ϵijy_{ij}=\beta_{0j}+x_{1ij}\beta_{1j}+\epsilon_{ij}

  • 여기서 β0j\beta_{0j}는 무선 절편, β1j\beta_{1j}는 무선 기울기를 의미하며, ϵij\epsilon_{ij}는 평균이 0이고 분산이 σ2\sigma^{2}인 1수준 오차항입니다.

이를 2수준 방정식과 결합한 혼합 모형(Composite model)은 다음과 같습니다.

yij=β0+x1ijβ1+x2jβ2+x1ijx2jβ3+u0j+x1iju1j+ϵijy_{ij}=\beta_{0}+x_{1ij}\beta_{1}+x_{2j}\beta_{2}+x_{1ij}x_{2j}\beta_{3}+u_{0j}+x_{1ij}u_{1j}+\epsilon_{ij}

  • 개인 중심 연결망 데이터에서 1수준 변수와 2수준 변수의 곱으로 이루어진 교차 수준 상호작용 항(cross-level interaction term)은 자아와 타인의 결합, 즉 양자간의 상호작용을 나타내므로 매우 중요합니다.

단계 2: 양자 데이터의 다층모형 (너와 나의 연결 고리)

만약 지훈이가 민수에게 느끼는 친밀감(yijy_{ij})뿐만 아니라, 민수가 지훈이에게 느끼는 친밀감(yjiy_{ji})도 함께 본다면 어떨까요?. 두 개의 관계(tie)가 하나의 양자(dyad) 안에 내재되어 있다고 보는 것입니다.

  • 이러한 접근법은 심리학에서 커플 등을 연구할 때 사용하는 행위자-파트너 상호의존 모형(Actor-Partner Interdependence Model, APIM)으로 잘 알려져 있습니다.
  • 이 모형에서는 두 결과 간의 강한 상관관계인 급내 상관계수(intraclass correlation)를 ‘상호성(reciprocity)’으로 해석할 수 있습니다. (예: 내가 널 좋아하는 만큼, 너도 날 좋아하니?)

단계 3: 전체 연결망 데이터의 다층모형 (우리 반 전체의 인싸력과 인기도)

이제 학급 전체를 봅니다. 모든 학생이 모든 학생을 평가하는 상황이므로, 양자(dyads)들은 행위자(actors)들 내에 교차 내재(cross-nested)된 것으로 취급합니다.

  • 이러한 교차 내재 구조를 다루기 위해 사회적 관계 모형(Social Relations Model, SRM)을 사용합니다.
  • 이 모형은 행위자의 두 가지 역할(송신자와 수신자)을 모두 고려하기 위해 u1u_{1}과 u2u_{2}라는 무선 효과 벡터를 추가합니다. yij=β0+x1ijβ1+x2jβ2+∑d=1N(sdiju1d+rdiju2d)+u0j+ϵijy_{ij}=\beta_{0}+x_{1ij}\beta_{1}+x_{2j}\beta_{2}+\sum_{d=1}^{N}(s_{dij}u_{1d}+r_{dij}u_{2d})+u_{0j}+\epsilon_{ij}
  • 여기서 무선 송신자 효과(sender effect)는 외향성(outgoingness, “내가 얼마나 다른 친구들에게 잘 다가가는가?”)을 의미합니다.
  • 무선 수신자 효과(receiver effect)는 인기도(popularity, “다른 친구들이 나를 얼마나 좋아하는가?”)를 의미합니다.

단계 4: 이분형 사회연결망 데이터 (싫어해? 예/아니오)

항상 점수로만 관계를 측정하는 것은 아닙니다. “같이 있기 싫은 사람인가요?”처럼 결과 변수가 이분형(binary)일 때도 있습니다.

  • 이 경우 일반화 선형 혼합 모형(GLMM)인 다층 로지스틱 모형을 사용할 수 있습니다.
  • 양자 데이터의 상호성을 반영하기 위해 고정된 송신자/수신자 매개변수와 명시적인 상호성 매개변수를 사용하는 p1p_{1} 모형이나 , 이를 확장하여 정규 분포를 가정하는 상관된 무선 송신자/수신자 매개변수를 갖는 p2p_{2} 모형을 사용합니다.

4. R과 jamovi를 활용한 실습 가이드

사회연결망 다층분석 중 개인중심/양자 데이터 모형(APIM 형태)은 jamovi의 GAMLj 모듈(Mixed Model)을 통해 쉽게 구현할 수 있습니다. 데이터를 “Long format” (한 줄에 하나의 tie가 들어가도록)으로 구성하는 것이 핵심입니다. 하지만 SRM이나 교차 내재된 구조는 R의 lme4 또는 특화된 네트워크 패키지를 사용하는 것이 정확합니다.

여기서는 GAMLj의 기반이 되는 R 코드를 통해 데이터를 생성하고 분석하는 방법을 보여드리겠습니다. (이 코드는 jamovi의 Rj Editor에서도 바로 실행 가능합니다.)

R

# 1. 패키지 불러오기
library(lme4)
library(dplyr)
library(ggplot2)

# 2. 모의 데이터 생성 (햇살초 5학년 1반, 30명 기준)
set.seed(2026)
n_students <- 30
students <- data.frame(
  id = 1:n_students,
  gender = sample(c("Boy", "Girl"), n_students, replace = TRUE),
  mark = rnorm(n_students, mean = 7, sd = 1),       # 성적
  class_id = rnorm(n_students, mean = 3.5, sd = 0.8) # 학급 소속감
)

# 가능한 모든 양자 조합 만들기 (자기 자신 제외)
network_data <- expand.grid(ego_id = students$id, alter_id = students$id) %>%
  filter(ego_id != alter_id)

# 공변량 병합 (Ego와 Alter 특성)
network_data <- network_data %>%
  left_join(students, by = c("ego_id" = "id")) %>%
  rename(ego_gender = gender, ego_mark = mark, ego_classid = class_id) %>%
  left_join(students, by = c("alter_id" = "id")) %>%
  rename(alter_gender = gender, alter_mark = mark, alter_classid = class_id)

# 결과 변수 생성 (Tie strength): 동성끼리 친하고, 학급 소속감이 비슷할수록 친함
network_data <- network_data %>%
  mutate(
    same_gender = ifelse(ego_gender == alter_gender, 1, 0),
    classid_diff = abs(ego_classid - alter_classid),
    # 간단한 선형 결합 + 오차 + 송신자 무선효과 + 수신자 무선효과
    ego_random = rep(rnorm(n_students, 0, 0.5), each = n_students - 1),
    tie_strength = 2.0 + 1.2 * same_gender - 0.5 * classid_diff + ego_random + rnorm(n(), 0, 0.8)
  )

# 3. 데이터 탐색적 시각화 (성별 동질성에 따른 관계 강도)
ggplot(network_data, aes(x = factor(same_gender, labels=c("Different", "Same")), y = tie_strength, fill=factor(same_gender))) +
  geom_boxplot() +
  labs(title = "Tie Strength by Gender Similarity", x = "Gender Similarity", y = "Tie Strength") +
  theme_minimal()

# 4. 혼합 모형(다층 모형) 적합 (Egocentric base)
# jamovi의 GAMLj 모듈에서:
# - Dependent Variable: tie_strength
# - Factors: same_gender
# - Covariates: classid_diff
# - Random Effects: Intercept | ego_id
model_ego <- lmer(tie_strength ~ same_gender + classid_diff + (1 | ego_id), data = network_data)
summary(model_ego)

jamovi에서 분석하는 팁:

  1. 위에서 만든 데이터셋을 .csv로 저장하여 jamovi에서 엽니다.
  2. Modules -> jamovi library에서 GAMLj를 설치합니다.
  3. Linear Models -> Mixed Model을 클릭합니다.
  4. tie_strength를 Dependent Variable에 넣고, ego_id를 Cluster Variables에 넣습니다.
  5. Random Effects 탭에서 ego_id 하위의 Intercept를 추가하면, 각 학생(ego)마다 기본적으로 타인에게 부여하는 점수(outgoingness)가 얼마나 다른지 그 분산을 확인할 수 있습니다!

5. 참고문헌 (APA Style)

  • Cook, W. L., & Kenny, D. A. (2005). The Actor-Partner Interdependence Model: A model of bidirectional effects in developmental studies. International Journal of Behavioral Development, 29, 101-109.
  • Holland, P. W., & Leinhardt, S. (1981). An exponential family of probability distributions for directed graphs. Journal of the American Statistical Association, 76, 33-50.
  • Scott, J. (2000). Social network analysis: A handbook. London: Sage.
  • Snijders, T. A. B., & Kenny, D. A. (1999). The social relations model for family data: A multilevel approach. Personal Relationships, 6, 471-486.
  • Van Duijn, M. A. J., Snijders, T. A. B., & Zijlstra, B. J. H. (2004). p2p_2: A random effects model with covariates for directed graphs. Statistica Neerlandica, 58, 234-254.
  • Vermeij, L. (2006). What’s cooking? Cultural boundaries among Dutch teenagers of different ethnic origins in the context of school (Unpublished doctoral dissertation). ICS/University of Utrecht.

Chap32. 시장 조사와 선호도 자료(Market Research and Preference Data)


안녕하세요!

오늘은 다층분석(Multilevel Analysis)을 활용한 선호도 및 상호작용 데이터 분석에 대해 살펴보겠습니다. “학교 현장의 데이터”를 예시로 들어 직관적인 설명과 수리적 엄밀함을 모두 갖춘 형태로 재구성해 드리겠습니다.

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

1. 다층 데이터(Nested Data)란 무엇일까요? 🪆

초등학생 친구들에게 러시아의 전통 인형 ‘마트료시카(큰 인형 안에 작은 인형이 계속 들어있는 장난감)’를 보여주며 다층 데이터를 설명하곤 합니다. 교육 현장의 데이터는 대부분 이 마트료시카처럼 ‘내포된(Nested)’ 구조를 가집니다. 예를 들어, 여러 번의 시험 점수(반복 측정)는 한 명의 ‘학생’ 안에 쏙 들어가 있고 , 그 학생들은 다시 하나의 ‘학급’ 안에 들어가 있죠.

이러한 데이터 구조에서는 같은 반 친구들끼리, 혹은 한 학생이 낸 여러 답변끼리 서로 영향을 주고받아 독립성 가정이 깨지게 됩니다. 이를 무시하고 일반적인 회귀분석(OLS)이나 분산분석을 하면, 1종 오류(Type I error)가 크게 나타나고 잘못된 결과를 얻게 됩니다.

얼마나 내포되어 있는지 어떻게 알까요? (ICC)

이때 확인하는 것이 급내상관계수(Intraclass Correlation Coefficient, ICC)입니다. 전체 데이터의 변동 중에서 ‘집단 간의 차이’가 차지하는 비율을 의미합니다. 수식으로는 다음과 같이 표현됩니다.

ICC=δα2δα2+δϵ2ICC = \frac{\delta_{\alpha}^{2}}{\delta_{\alpha}^{2}+\delta_{\epsilon}^{2}}

보통 ICC 값이 0.05에서 0.2 사이이거나, 설계효과(Design effect, deff)가 2 이상이면 다층분석을 도입하여 집단 간 변량(수식에서 분자에 해당)을 분리해야 합니다.

샘플 사이즈는 어느 정도가 적당할까요?

다층분석은 일반 분석보다 더 많은 샘플을 필요로 합니다.

  • 기본 규칙 (30/30 법칙): 최소 30개의 집단이 있어야 하고, 각 집단 안에는 최소 30개의 관측치가 있어야 합니다.
  • 수준 간 상호작용 분석 시 (50/20 법칙): 더 엄격하게 50개 집단, 각 집단당 20명 이상이 권장됩니다.
  • 무선 효과(Random effects) 분석 시 (100/10 법칙): 100개 집단과 각 집단당 10명이 권장되기도 합니다.

2. 교육 선호도 분석에서의 다층모형 (Multilevel Preference Models) 📊

교육학에서 새로운 교재나 프로그램을 개발할 때, 학생들의 선호도를 측정하는 것은 매우 중요합니다. 마케팅 분야에서 상품 개발 시 널리 쓰이는 ‘컨조인트 분석(Conjoint Analysis)’이나 이산선택모형을 다층모형으로 재구성할 수 있습니다.

  • 1수준 (Within-level): 개별 학생 안에서 측정된 반복 측정치입니다. 예를 들어, 하나의 교육 프로그램이 가진 속성(비용, 접근성, 복잡성 등)에 대한 여러 번의 평가가 여기에 해당합니다.
  • 2수준 (Between-level): 학생 개인의 특성입니다. 예를 들어, 학생의 성별, 연령, 혹은 사전 학습 동기 수준 등이 선호도 평균(절편)에 어떤 영향을 미치는지 설명합니다.

Muthén의 분해 방식에 따르면, 학생의 특정 응답 점수(YijY_{ij})는 학생들 간의 평균적 차이(YBY_{B})와 그 학생 내에서의 변동(YWY_{W})으로 나누어집니다.

3. 교사와 학생의 상호작용 분석 (LS-T 및 APIM 모형) 🤝

요즘 교육 현장에서는 상황과 맥락에 따른 상호작용이 강조됩니다. 수업 시간에 학생의 기분이나 몰입도는 계속 변하죠.

  • 잠재상태-특성 모형 (LS-T): 학생의 반응을 비교적 안정적인 개인적 성향인 ‘특성(Trait)’과 특정 수업 상황에 따른 일시적 기분인 ‘상태(State)’, 그리고 오차로 분리해 냅니다. 이를 통해 특정 수업 기법이 학생의 일시적 몰입도(상태)를 얼마나 끌어올렸는지 정확히 파악할 수 있습니다.
  • 행위자-상대방 상호의존성 모형 (APIM): 교사와 학생을 하나의 ‘짝(Dyad)’으로 봅니다.
    • 행위자 효과 (Actor effect): 교사의 열정이 교사 스스로의 만족도에 영향을 주는 것입니다.
    • 상대방 효과 (Partner effect): 교사의 열정이 짝꿍인 학생의 몰입도에 영향을 주는 것입니다.

4. 실전! 모의 데이터 기반 jamovi 및 R 분석 💻

📖 스토리보드: “미래형 AI 튜터링 시스템, 학생들의 선택은?”

A 교육청에서는 새로운 ‘AI 수학 튜터’를 도입하려 합니다. 50명의 학생에게 AI 튜터의 3가지 속성(1. 피드백 속도, 2. 게임화 여부, 3. 아바타 유무)을 조합한 8개의 프로필을 보여주고 선호도 점수(0~100점)를 매기게 했습니다. 또한, 학생들의 수학 불안 수준(Math Anxiety)도 2수준 변인으로 함께 조사했습니다.

⌨️ R 코드 (모의 데이터 생성 및 분석)

R

# 패키지 로드
library(lme4)
library(lmerTest)
library(ggplot2)

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

# 2. 2수준(학생 단위) 데이터 생성 (50명)
n_students <- 50
students <- data.frame(
  StudentID = factor(1:n_students),
  MathAnxiety = rnorm(n_students, mean = 50, sd = 10), # 수학 불안도
  StudentRandomEffect = rnorm(n_students, mean = 0, sd = 8) # 학생별 선호도 편차(랜덤 효과)
)

# 3. 1수준(반복 측정 단위) 데이터 생성 (각 학생당 8개 프로필 평가)
n_profiles <- 8
total_obs <- n_students * n_profiles

# 요인설계 (피드백 속도, 게임화, 아바타)
profiles <- expand.grid(
  FeedbackSpeed = c("Fast", "Slow"),
  Gamification = c("Yes", "No"),
  Avatar = c("Yes", "No")
)

# 데이터 병합을 위한 기초 프레임
df <- data.frame(
  StudentID = rep(students$StudentID, each = n_profiles),
  FeedbackSpeed = rep(profiles$FeedbackSpeed, times = n_students),
  Gamification = rep(profiles$Gamification, times = n_students),
  Avatar = rep(profiles$Avatar, times = n_students)
)

# 학생 정보 병합
df <- merge(df, students, by="StudentID")

# 4. 종속변수(선호도) 생성 
# 수학 불안도가 높을수록 선호도가 낮아지고, 피드백이 빠르고 게임화 요소가 있을 때 선호도가 올라감
df$Preference <- 60 + 
  ifelse(df$FeedbackSpeed == "Fast", 15, -5) + 
  ifelse(df$Gamification == "Yes", 10, -5) + 
  ifelse(df$Avatar == "Yes", 5, 0) - 
  (df$MathAnxiety * 0.3) + 
  df$StudentRandomEffect + 
  rnorm(total_obs, mean = 0, sd = 5) # 1수준 오차

# 점수를 0~100 사이로 제한
df$Preference <- pmax(0, pmin(100, round(df$Preference, 1)))

# 5. 다층분석(선형 혼합 모형) 실행
# 학생 내 반복측정이므로 StudentID를 무선 절편(Random Intercept)으로 설정
model <- lmer(Preference ~ FeedbackSpeed + Gamification + Avatar + MathAnxiety + (1 | StudentID), data=df)

# 결과 요약
summary(model)

📊 jamovi를 이용한 분석 가이드

위에서 생성한 데이터(CSV 저장 후 불러오기)를 jamovi에서 분석하는 방법은 아주 간단합니다.

  1. 모듈 설치 및 선택: 우측 상단의 + Modules를 클릭하고 GAMLj 모듈(선형 혼합 모형 지원)을 설치합니다.
  2. 분석 메뉴 진입: Linear Models -> Mixed Model을 클릭합니다.
  3. 변수 투입:
    • Dependent Variable: Preference(선호도)
    • Factors: FeedbackSpeed, Gamification, Avatar (명목형 변수)
    • Covariates: MathAnxiety (연속형 변수)
    • Cluster Variables: StudentID (학생을 기준으로 군집화)
  4. 랜덤 효과 설정 (Random Effects): StudentID를 투입하여 랜덤 절편(Random Intercept) 모형을 구성합니다.
  5. 해석: Fixed Effects Estimate 표에서 학생의 수학 불안도(MathAnxiety)와 각 AI 속성이 선호도에 미치는 영향을 파악합니다. Random Components 표에서는 급내상관계수(ICC)를 확인하여 다층분석 도입의 타당성을 입증합니다.

참고문헌 (APA Style)

  • Hox, J. (2010). Multilevel analysis: Techniques and applications. New York, NY: Routledge.
  • Muthén, B., & Satorra, A. (1995). Complex sample data in structural equation modeling. Sociological Methodology, 25, 87-99.
  • Peugh, J. L. (2010). A practical guide to multilevel modeling. Journal of School Psychology, 48(1), 85-112.
  • Snijders, T. A. B., & Kenny, D. A. (1999). The social relations model for family data: A multilevel approach. Personal Relationships, 6(4), 471-486.
  • Steyer, R., & Schmitt, M. J. (1990). Latent state-trait models in attitude research. Quality and Quantity, 24, 427-445.

Chap31. 점-참조 공간 모델링(Point-Referenced Spatial Modeling)

안녕하세요!

오늘은 점-참조 공간 모델링(Point-Referenced Spatial Modeling)에 대해 살펴보겠습니다. “학교 현장의 데이터”를 예시로 들어 직관적인 설명과 수리적 엄밀함을 모두 갖춘 형태로 재구성해 드리겠습니다.

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

1. 초등학생도 이해하는 ‘공간 모델링’ 이야기

초등학교 교실을 떠올려 볼까요? 교실에서 어떤 친구가 지우개 가루를 불며 장난을 칩니다. 이 장난꾸러기 친구 바로 옆에 앉은 짝꿍은 지우개 가루를 많이 뒤집어쓰겠지만, 교실 맨 뒤에 앉은 친구는 거의 피해가 없을 거예요. 즉, “거리가 가까울수록 서로에게 미치는 영향이 크고, 멀어질수록 그 영향이 줄어든다”는 사실을 우리는 직관적으로 알고 있습니다.

공간 모델링도 이와 똑같습니다! 기존의 다층모형(MLM)이 주로 학생들을 학급이나 학교라는 ‘그룹’으로 묶어서 관찰하는 방식이라면 , 점-참조 공간 모델링은 각 데이터가 존재하는 ‘정확한 지리적 좌표(예: 위도와 경도)’를 바탕으로 분석을 진행합니다.

이 분석에서는 두 가지 중요한 관계(연관성)를 살펴봅니다.

  1. 첫 번째 관계: 서로 다른 위치(학교)에 있는 ‘같은’ 변수 간의 관계입니다. 예를 들어, A학교의 수학 성적은 바로 옆동네 B학교의 수학 성적과 비슷할 확률이 높습니다.
  2. 두 번째 관계: 같은 위치(학교) 내에 있는 ‘서로 다른’ 변수 간의 관계입니다. 예를 들어, A학교 학생들의 ‘수학 성적’과 ‘행복도’ 사이의 관계를 뜻합니다.

2. 수리적 모델의 이해: 단변량 및 다변량 공간 회귀모형

이제 대학원생 모드로 전환하여 조금 더 수리적으로 엄밀하게 접근해 보겠습니다.

2.1 단변량 공간 회귀모형 (Univariate Spatial Regression)

지리적 위치 ss에 위치한 학교의 성과(예: 수학 성적)를 y(s)y(s)라고 할 때, 이를 설명하는 회귀 모형은 다음과 같습니다.

y(s)=x(s)′β+u(s)+ϵ(s)y(s) = x(s)’\beta + u(s) + \epsilon(s)

  • x(s)′βx(s)’\beta: 학교의 교육예산이나 지역의 사회경제적 배경(SES)처럼 우리가 측정한 설명 변수들입니다.
  • u(s)u(s): 공간 확률 효과(Spatial random effects)입니다. 우리가 미처 측정하지 못했지만 공간적 패턴을 가진 요인들(예: 인근 지역의 교육열)을 잡아냅니다. 가우스 과정(Gaussian process)으로 모형화됩니다.
  • ϵ(s)\epsilon(s): 독립적인 백색 잡음(white-noise process)으로, 순수한 측정 오차나 아주 미시적인 변동성을 의미합니다. 흔히 ‘너겟(nugget)’이라고 부릅니다.

2.2 대규모 데이터의 한계와 예측 프로세스 모형 (Predictive Process Models)

학교 수가 10,000개처럼 엄청나게 많아지면 어떻게 될까요? 공간 분석에서는 데이터 개수(nn)의 세제곱인 O(n3)O(n^3) 수준으로 계산이 복잡해지는 이른바 ‘Big nn‘ 문제가 발생합니다.

이를 해결하기 위해 전체 학교 중 대표성을 띠는 소수의 지점(매듭, knots)만을 선택하여 차원을 축소하는 예측 프로세스 모형(Predictive process models)을 사용합니다. 이는 데이터의 공간적 정보를 잃지 않으면서도 계산 속도를 획기적으로 높여줍니다.

2.3 다변량 공간 회귀모형 (Multivariate Spatial Regression)

만약 수학 성적(y1y_1)과 학생 행복도(y2y_2)라는 다수의 결과 변수를 동시에 분석하고 싶다면 다변량 모형을 확장하여 사용합니다.

y(s)=X(s)′β+u(s)+ϵ(s)y(s) = X(s)’\beta + u(s) + \epsilon(s)

이때 구조는 동일하지만, y(s)y(s), u(s)u(s), ϵ(s)\epsilon(s)가 단일 숫자가 아닌 벡터(vector) 형태가 되며, 변수들 간의 교차-공분산(cross-covariance)을 고려하게 됩니다.

💡 면적 데이터(Areal Data)와의 차이점 교육청에서 개인정보 보호를 위해 개별 학교의 위치를 숨기고 “수원시 장안구 전체 평균 성적”, “권선구 전체 평균 성적”처럼 행정 구역 단위로 묶어서 발표할 때가 있습니다. 이처럼 불연속적이고 이웃한 ‘구역’ 간의 인접성(adjacency)을 바탕으로 분석하는 것을 면적 모형(CAR, SAR 등)이라고 합니다. 반면 오늘 우리가 다루는 점-참조 모형은 공간 상의 어느 지점에서든 연속적으로 예측할 수 있다는 점에서 근본적인 차이가 있습니다.

3. 학교 현장 모의 데이터 생성 및 스토리텔링

[가상의 시나리오]

경기도 K시의 교육감은 관내 100개 초등학교를 대상으로 “교육 예산(Budget) 투입이 학교별 평균 수학 성적(Math)과 학생 행복도(Happiness)에 미치는 영향”을 조사하고자 합니다. 그런데 지도를 보니, 특정 신도시 지역에 위치한 학교들이 예산과 무관하게 전반적으로 성적과 행복도가 모두 높은 ‘공간적 군집 현상’을 보이고 있었습니다. 단순히 예산만 볼 것이 아니라, 지리적 근접성에 따른 후광효과(공간 상관성)를 통제하고 진짜 예산의 효과를 알고 싶습니다.

3.1 R을 활용한 모의 데이터 생성

이 분석을 구현하기 위한 R 코드를 작성해 보겠습니다.

R

# 필요한 패키지 로드
library(MASS)
library(ggplot2)

# 1. 공간 좌표 생성 (K시 내 100개 초등학교)
set.seed(2026)
n_schools <- 100
coords <- cbind(Easting = runif(n_schools, 0, 10), Northing = runif(n_schools, 0, 10))

# 2. 공간 거리 행렬 및 공분산 구조 생성
dist_matrix <- as.matrix(dist(coords))
phi <- 0.5 # 공간 감쇠 파라미터
sigma_sq <- 2.0 # 공간 분산
spatial_cov <- sigma_sq * exp(-phi * dist_matrix)

# 3. 공간 임의 효과 (Spatial random effects) 생성
u_math <- mvrnorm(1, mu = rep(0, n_schools), Sigma = spatial_cov)
u_happy <- 0.6 * u_math + mvrnorm(1, mu = rep(0, n_schools), Sigma = spatial_cov * 0.5)

# 4. 설명 변수(교육 예산) 및 최종 종속 변수 생성
budget <- rnorm(n_schools, mean = 50, sd = 10)
math_score <- 30 + 0.5 * budget + u_math + rnorm(n_schools, 0, 1) # 백색 잡음(nugget) 추가
happiness <- 10 + 0.2 * budget + u_happy + rnorm(n_schools, 0, 0.5)

# 데이터 프레임 구축
school_data <- data.frame(
  School_ID = 1:n_schools,
  Easting = coords[,1],
  Northing = coords[,2],
  Budget = budget,
  Math = math_score,
  Happiness = happiness
)

# 데이터 확인
head(school_data)

4. jamovi 및 R을 활용한 분석 방법

4.1 jamovi에서의 탐색적 분석

jamovi는 직관적인 UI를 제공하여 기초 통계를 확인하기 매우 좋습니다.

  1. 위에서 생성한 school_data.csv를 jamovi에 불러옵니다.
  2. Exploration -> Descriptives 탭을 클릭하여 Math, Happiness, Budget의 평균과 산포도를 확인합니다.
  3. 변수 간의 단순한 관계는 Regression -> Correlation Matrix를 통해 확인할 수 있습니다.
  4. 주의: jamovi의 기본 기능이나 GAMLj 다층모형 모듈은 ‘지리적 좌표 기반의 연속적 공간 상관구조’를 완벽히 모델링하는 데 한계가 있습니다.

4.2 R을 활용한 본격적인 베이지안 공간 모델링 (spBayes 패키지)

본문에서 제시된 바와 같이, 점-참조 공간 모델링은 MCMC 기법을 활용한 베이지안 방식으로 추정하며 spBayes 패키지를 사용하면 효과적입니다. jamovi의 Rj Editor 모듈을 열어 다음 코드를 실행하면 분석이 가능합니다.

R

# spBayes 패키지 로드 (설치 필요: install.packages("spBayes"))
library(spBayes)

# 종속 변수 행렬 구성
Y <- as.matrix(school_data[, c("Math", "Happiness")])

# 필요한 패키지가 불러와져 있다고 가정합니다. (MASS, spBayes)

# 1. 다변량 모형을 위한 수식(Formula) 작성
formulas <- list(Math ~ Budget, Happiness ~ Budget)

# 2. 모형 적합 (퍼즐 조각 개수 완벽 수정 버전!)
n_iters <- 1000 
fit_mv <- spMvLM(formulas, 
                 data = school_data,  
                 coords = as.matrix(school_data[, c("Easting", "Northing")]),
                 
                 # 수정 포인트 1: A는 길이 3 (하삼각행렬), Psi는 길이 2 (대각행렬)로 설정
                 starting = list("beta" = rep(0, 4), 
                                 "phi" = rep(0.5, 2), 
                                 "A" = c(1, 0, 1), 
                                 "Psi" = c(1, 1)),
                 
                 # 수정 포인트 2: tuning 파라미터도 starting의 길이에 맞춤
                 tuning = list("phi" = rep(0.1, 2), 
                               "A" = c(0.1, 0.1, 0.1), 
                               "Psi" = c(0.1, 0.1)), 
                 
                 # 수정 포인트 3: 역위샤트 분포(K.IW)의 자유도를 2보다 큰 3으로 두어 수학적 안정성 확보
                 priors = list("beta.flat", 
                               "phi.Unif" = list(rep(0.1, 2), rep(2, 2)), 
                               "K.IW" = list(3, diag(2)), 
                               "Psi.ig" = list(c(2, 2), c(0.1, 0.1))),
                 
                 cov.model = "exponential",
                 n.samples = n_iters)

# 3. 모델 적합 후 사후표본(posterior samples) 복원
# 공간 임의효과를 포함하여 최종적으로 베타(회귀계수)를 복원해 냅니다.
fit_mv <- spRecover(fit_mv, start = 501, thin = 1) 

# 4. 결과 요약 확인
summary(fit_mv$p.beta.recover.samples)

이 결과를 통해 우리는 “주변 학교의 공간적 영향을 통제한 후에도, 교육 예산이 수학 성적과 행복도에 유의미한 정(+)의 효과를 가지는지” 정확하게 평가할 수 있습니다.

5. 요약 및 참고문헌

다층모형에서 공간 모델링으로의 확장은 교육 연구에 새로운 시각을 제공합니다. 학생이나 학교를 독립된 개체로 보는 것을 넘어, “우리는 지리적으로 연결되어 있다”는 사실을 데이터로 증명해 낼 수 있기 때문입니다.

참고문헌

  • Banerjee, S., Carlin, B. P., & Gelfand, A. E. (2004). Hierarchical modeling and analysis for spatial data. Chapman and Hall/CRC Press.
  • Finley, A. O., & Banerjee, S. (2012). Point-Referenced Spatial Modeling. In The SAGE Handbook of Multilevel Modeling (pp. 559-579).
  • Finley, A. O., Sang, H., Banerjee, S., & Gelfand, A. E. (2009). Improving the performance of predictive process modeling for large datasets. Computational Statistics and Data Analysis, 53(8), 2873-2884.

Chap30. 생존분석(Survival Analysis)과 노쇠모형(Frailty Model)

안녕하세요!

오늘은 생존분석과 노쇠모형(Survival Analysis and the Frailty Model)에 대해 살펴보겠습니다. “학교 현장의 데이터”를 예시로 들어 직관적인 설명과 수리적 엄밀함을 모두 갖춘 형태로 재구성해 드리겠습니다.

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

1. 생존분석과 노쇠모형이란 무엇일까요?

초등학생도 이해하기 쉽게 비유해 볼까요? 우리가 버스 정류장에서 버스를 기다린다고 상상해 보세요. 어떤 사람은 5분 만에 버스를 타고, 어떤 사람은 30분을 기다려도 버스가 오지 않아 포기하고 걸어갈 수도 있습니다. 여기서 ‘버스를 탈 때까지 걸린 시간’을 연구하는 것이 바로 생존분석(Survival Analysis)입니다.

가. 생존분석 (Survival Analysis)

  • 생존분석은 어떤 사건(event)이 발생하기까지 걸린 시간(time to an event)을 조사하는 통계적 분석 방법입니다.
  • 의학 분야에서는 환자의 사망을 사건으로 정의하거나, 공학에서는 기계 부품의 고장 등을 사건으로 정의하여 널리 응용됩니다.
  • 생존분석의 가장 큰 특징은 중도절단(censoring) 현상을 다룰 수 있다는 점입니다. 중도절단이란 연구 기간이 종료될 때까지 관심 있는 사건이 발생하지 않아 정확한 시간을 알 수 없는 경우를 뜻합니다.
  • 예를 들어 의학 연구에서 우측 중도절단(right-censoring)은 환자를 추적 관찰했지만 연구 기간 동안 사망 등 특정 사건이 발생하지 않은 경우를 말합니다.

나. 노쇠모형 (Frailty Model)

  • 노쇠모형(Frailty models)은 다층 생존 모형(multilevel models for survival)의 일종으로, 사건 발생 위험에 무작위 효과(random effects)를 연결한 모형입니다.
  • 이 용어는 의학 분야에서 기원했으며, ‘노쇠(frailty)’란 질병에 대한 일반적인 취약성을 의미합니다.
  • 노쇠 항은 개인 단위의 무작위 효과로 지정할 수도 있고, 군집 특성(cluster-specific)으로 지정할 수도 있습니다. 군집 단위로 지정할 경우, 같은 군집 내의 개인들은 동일한 취약성을 공유하게 됩니다.

다. 다상태 모형 (Multi-state Model)

  • 다상태 모형은 사건이 단 한 번 발생하는 고전적인 생존모형을 여러 개의 사건으로 일반화한 모형입니다.
  • 단 하나의 ‘생존’ 상태에서 ‘사망’이라는 흡수 상태(absorbing state)로의 이동만 허용하는 생존모형과 달리, 확률 과정 내에서 유한한 수의 여러 상태 간 이동을 허용합니다.
  • 대표적인 예로 사망에 이르기 전 중간 단계인 질병 상태를 구분하는 질병-사망 모형(illness-death model)이 있습니다. 가장 기본적인 질병-사망 모형은 건강한 상태, 질병 상태, 사망 상태로 구성됩니다.

2. 교육 현장 예시: 학생의 학업중단(자퇴) 분석

의학 연구인 CFAS(Cognitive Function and Ageing Study) 데이터를 교육학 버전으로 각색하여 재미있는 스토리를 만들어보겠습니다.

[상황 설정]

우리는 5개의 고등학교(군집, cluster)에 입학한 500명의 학생을 대상으로 ‘학업중단(자퇴)’까지 걸리는 시간을 3년(36개월) 동안 추적 조사했습니다.

  • 사건(Event): 학업중단 (자퇴)
  • 시간(Time): 입학 후 학업중단까지 걸린 개월 수
  • 중도절단(Censoring): 3년 무사고 졸업생 (자퇴하지 않음)
  • 다층 구조(Multilevel): 학생들은 5개의 서로 다른 고등학교에 소속되어 있습니다. 학교마다 교풍이나 학업 스트레스가 다르기 때문에, 특정 학교 학생들은 자퇴할 ‘취약성(Frailty)’이 더 높을 수 있습니다.
  • 공변량(Covariates): 멘토링 프로그램 참여 여부(0=미참여, 1=참여), 학생의 사회경제적 배경(SES).

우리는 이 데이터를 통해 특정 순간의 자퇴 위험인 위험 함수(Hazard function)를 추정할 것입니다. 순간 위험률이 높다는 것은 짧은 시간 구간 안에서 사건(자퇴)이 발생할 확률이 높다는 것을 의미합니다.

다층 생존 모형의 수식은 다음과 같이 표현됩니다:

hij(t)=h0(t)exp⁡[xij(t)β+uj]h_{ij}(t)=h_{0}(t)\exp[x_{ij}(t)\beta+u_{j}]

여기서 h0(t)h_0(t)는 기저 위험(baseline hazard)이며, uju_j는 학교(군집)별 취약성을 나타내는 무작위 효과로 정규분포 uj∼N(0,σ2)u_j \sim N(0,\sigma^2)를 따른다고 가정합니다.

3. 모의 데이터 생성 및 R/jamovi 분석

통계 분석을 위해 R 프로그램과 jamovi를 활용합니다. 복잡한 다층 생존모형은 R의 coxme나 survival 패키지를 사용하는 것이 가장 정확합니다 .

가. R 모의 데이터 생성 코드

아래 코드를 R이나 jamovi의 Rj Editor 모듈에 붙여넣어 데이터를 생성할 수 있습니다.

R

# 필수 패키지 로드
library(survival)
library(coxme) # 다층 생존분석(Frailty)을 위한 패키지
set.seed(2026)

# 1. 학교(군집) 설정: 5개 학교, 학교별 무작위 효과(Frailty) 
n_schools <- 5
n_students_per_school <- 100
total_N <- n_schools * n_students_per_school

school_id <- rep(1:n_schools, each = n_students_per_school)
# 학교별 취약성(Frailty): 값이 클수록 자퇴 위험이 높음
school_frailty <- rnorm(n_schools, mean = 0, sd = 0.5) 

# 2. 공변량 설정
# 멘토링 참여 (30% 확률로 참여)
mentoring <- rbinom(total_N, 1, 0.3)
# 사회경제적 지위 (SES, 표준정규분포)
ses <- rnorm(total_N, 0, 1)

# 3. 생존 시간 및 사건 발생 생성 (Weibull 기저 위험 가정)
# 멘토링은 자퇴 위험을 낮추고(음의 계수), SES가 높을수록 자퇴 위험 감소
linear_predictor <- -0.8 * mentoring - 0.3 * ses + school_frailty[school_id]
scale_param <- exp(-(linear_predictor) / 1.5) # Weibull scale
true_time <- rweibull(total_N, shape = 1.5, scale = 50 * scale_param)

# 36개월 관찰 한계 (졸업 = 우측 중도절단)
time <- pmin(true_time, 36)
event <- ifelse(true_time <= 36, 1, 0) # 1=자퇴(사건 발생), 0=졸업(절단)

# 데이터프레임 생성
edu_data <- data.frame(
  StudentID = 1:total_N,
  SchoolID = as.factor(school_id),
  Time_Months = time,
  Dropout = event,
  Mentoring = mentoring,
  SES = ses
)

head(edu_data)
chap30다운로드

나. jamovi 분석 방법

  1. 데이터 불러오기: 위에서 생성한 edu_data.csv 파일을 jamovi에서 엽니다.
  2. 모듈 설치: jamovi 우측 상단의 + Modules 메뉴에서 jsurvival (또는 생존분석 관련 모듈)과 Rj Editor를 설치합니다.
  3. 분석 실행: 기본 Cox 비례위험모형은 jsurvival > Cox Regression에서 Time_Months를 Time에, Dropout을 Event에 넣고 Covariates에 Mentoring과 SES를 넣어 분석합니다.
  4. 노쇠모형(다층분석) 실행: 다층 구조의 학교 효과(uju_j)를 반영하려면 Rj Editor를 열고 아래 코드를 실행합니다.

R

# 혼합효과 콕스 회귀분석 (Mixed Effects Cox Model)
fit_frailty <- coxme(Surv(Time_Months, Dropout) ~ Mentoring + SES + (1 | SchoolID), data = edu_data)
print(fit_frailty)

4. 분석 결과 해석 (가상의 표)

위 모델을 돌렸을 때 나오는 결과를 교육학적으로 풀어보겠습니다.

변인추정치 (Estimate)표준오차 (SE)p-값해석
Mentoring-0.260.17.12
SES-0.320.07< .001SES 점수가 1점 높을수록 자퇴 위험이 유의미하게 감소합니다.
SchoolID (Frailty Variance, σ2\sigma^2)0.084––학교 간 자퇴 위험의 차이(취약성)가 존재함을 의미합니다.

결론적으로:

  • 멘토링 프로그램은 학생들의 학업중단(자퇴)을 막는 데 매우 긍정적인 효과가 있습니다. (이는 원문에서 교육 기간이 길수록 사망 위험이 감소하여 생존에 긍정적인 영향을 미친다는 결과와 유사한 맥락입니다 ).
  • 또한, 특정 학교는 다른 학교에 비해 근본적으로 학업중단율이 높은 경향(군집별 무작위 효과)을 보였습니다. 이를 고려하지 않으면 멘토링 효과를 왜곡해서 해석할 위험이 있습니다.

생존 곡선을 ‘계단식 미끄럼틀’이라고 상상해 보세요! 처음 입학했을 때(0개월)는 모든 학생이 학교에 다니고 있으니 꼭대기(생존율 100%, 또는 1.0)에서 시작합니다. 그러다가 한 명씩 자퇴(사건 발생)를 할 때마다 미끄럼틀이 계단처럼 한 칸씩 아래로 뚝, 뚝 떨어집니다. 만약 미끄럼틀이 아주 천천히 완만하게 내려간다면 학생들이 학교를 잘 다니고 있다는 뜻이고, 가파르게 푹푹 꺼진다면 많은 학생이 일찍 학교를 그만둔다는 뜻입니다.

이번에는 멘토링 프로그램이 학생들의 ‘학교 생존(자퇴하지 않고 재학함)’에 미치는 장기적 효과를 한눈에 보여주는 생존 곡선(Survival Curve)에 대해 아주 쉽게 설명해 드리겠습니다. 우리는 멘토링을 받은 학생들(파란색 미끄럼틀)과 받지 않은 학생들(빨간색 미끄럼틀)의 경사도를 비교해 볼 것입니다.

# 필수 패키지 로드
library(survival)

# 1. 카플란-마이어(Kaplan-Meier) 생존 객체 생성
# 멘토링 참여 여부(Mentoring)에 따른 생존율 계산
km_fit <- survfit(Surv(Time_Months, Dropout) ~ Mentoring, data = edu_data)

# 2. 생존 곡선 그래프 그리기
plot(km_fit, 
     col = c("red", "blue"),        # 선 색상 (빨강: 미참여, 파랑: 참여)
     lty = 1:2,                     # 선 종류 (실선, 점선)
     lwd = 2,                       # 선 굵기
     xlab = "Time in Months (입학 후 개월 수)", 
     ylab = "Survival Probability (재학 확률)", 
     main = "멘토링 프로그램 참여 여부에 따른 학업중단 생존 곡선")

# 3. 범례(Legend) 추가
legend("bottomleft", 
       legend = c("멘토링 미참여 (0)", "멘토링 참여 (1)"), 
       col = c("red", "blue"), 
       lty = 1:2, 
       lwd = 2,
       title = "집단 구분")

위 코드를 실행하여 그려진 그래프를 보면 두 개의 선이 나타납니다.

  • X축 (Time in Months): 입학 후 흐른 시간(0개월 ~ 36개월)입니다.
  • Y축 (Survival Probability): 학교에 남아있는(자퇴하지 않은) 학생의 비율입니다. 1.0은 100%를 의미합니다.
  • 파란색 점선 (멘토링 참여): 위쪽에 위치하며 아주 완만하게 내려갑니다. 즉, 시간이 지나도 대다수의 학생이 자퇴하지 않고 학교에 잘 다니고 있습니다.
  • 빨간색 실선 (멘토링 미참여): 아래쪽으로 뚝뚝 떨어집니다. 파란 선에 비해 상대적으로 많은 학생이 더 이른 시기에 학교를 그만두고 있음을 한눈에 알 수 있습니다.

이 그래프는 두 집단의 생존 곡선이 교차하지 않고 뚜렷한 차이를 보인다면 멘토링의 효과가 시간에 걸쳐 지속적으로 나타나고 있다는 시각적 증거가 됩니다.

5. 요약 및 제언

우리는 오늘 생존분석을 통해 ‘시간의 흐름’에 따른 중도절단 데이터를 다루고, 노쇠모형(Frailty model)을 통해 집단 간의 숨겨진 차이(다층 구조)를 통계적으로 보정하는 방법을 알아보았습니다. 다상태 모형을 적용한다면 “일반 학생 →\rightarrow 학업 위기 학생 →\rightarrow 학업 중단”이라는 3단계 전이 과정을 더욱 세밀하게 분석할 수도 있습니다.

참고문헌

  • van den Hout, A., & Tom, B. D. M. (2012). Survival Analysis and the Frailty Model. In The SAGE Handbook of Multilevel Modeling (pp. 541–558). SAGE Publications.

Chap28. 다층모형(Multilevel Modeling)을 활용한 정책 도입 및 효과 분석

안녕하세요!

오늘은 다층모형(Multilevel Modeling)을 활용한 정책 도입 및 효과 분석에 대해 살펴보겠습니다. “학교 현장의 데이터”를 예시로 들어 직관적인 설명과 수리적 엄밀함을 모두 갖춘 형태로 재구성해 드리겠습니다.

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

1. 왜 다층분석(Multilevel Modeling)이 필요할까요?

정책 분석가들은 특정 정책이 채택되는 원인과, 정책이 실행된 이후의 결과를 모두 연구합니다. 이러한 질문을 탐구할 때, 정책이 시간이 지남에 따라 역동적으로 변하며 대개 위계적(hierarchical) 구조를 통해 개발되고 실행된다는 점을 인식하는 것이 필수적입니다.

  • 군집화(Clustering)의 일상화: 정책은 종종 고유한 특성을 가진 개별 진료소나 학교를 통해 실행됩니다. 이처럼 프로그램의 설계와 행정에서 군집화는 널리 퍼진 특징입니다.
  • 분산의 무시 = 잘못된 추론: 정책 데이터가 가지는 이러한 위계적 특성이나 종단적 특성을 무시하는 것은 분산의 원천을 무시하는 것이며, 잠재적으로 잘못된 추론을 이끌어낼 수 있습니다. 흥미롭게도 교육학 분야는 프로그램의 효과성 연구에서 다층 연구의 예시를 정치학보다 훨씬 더 많이 보여주고 있습니다.

2. 내재된 데이터(Nested Data): 교실 안의 학생들

위계적인 방식으로 실행되는 프로그램의 영향을 분석하는 예로, 5학년 학생들의 수학 시험 점수를 결과 변수로 가정해 보겠습니다. 이 경우 학생들은 교실에 소속되고, 교실은 학교에, 학교는 교육청에 소속되어 프로그램을 수행하게 됩니다. 정책이 실행되는 방식은 각 수준(Level)마다 다릅니다. 담임 교사는 고유한 교수법과 자격을 지니고 있으며 , 학교는 학생들에게 제공하는 시설과 규율 정책이 다를 수 있습니다.

만약 연구자가 이런 군집성을 무시하고 전통적인 선형 모형인 최소제곱법(OLS)으로만 분석하면 어떻게 될까요?

  • 데이터의 위계적 구조는 같은 상위 집단 내 관측치들 간에 자기상관(autocorrelation)을 유발합니다.
  • 급내상관(intraclass correlation)의 존재는 일반적으로 표준오차를 하향 편향되게 만들어 t-값과 z-값을 높입니다.
  • 이는 연구자가 긍정적인 결과(통계적으로 유의미함)를 잘못 보고할 가능성, 즉 ‘1종 오류(False Positive)’를 증가시킵니다.

이를 방지하기 위해 교육청(ll), 학교(kk), 교실(jj), 학생(ii) 수준을 모두 고려한 임의절편모형(Random Intercept Model)은 다음과 같이 표현할 수 있습니다.

yijkl=xijklβ+ul+δkl+vjkl+ϵijkly_{ijkl}=x_{ijkl}\beta+u_{l}+\delta_{kl}+v_{jkl}+\epsilon_{ijkl}

  • 여기서 xijklx_{ijkl}은 공변량 벡터이고, β\beta는 고정계수(fixed coefficients)입니다.
  • 설명되지 않은 분산은 교육청 수준(ulu_{l}), 학교 수준(δkl\delta_{kl}), 교실 수준(vjklv_{jkl}), 그리고 학생 수준(ϵijkl\epsilon_{ijkl})의 네 가지 구성 요소로 나뉩니다.

또한, 교사가 석사 학위를 가지고 있는지 여부가 학생의 성취도에 미치는 영향이 학교마다 다를 수 있다고 가정한다면, 기울기마저 변하는 임의기울기모형(Random Slope Model)을 적용할 수 있습니다.

yijkl=β0+x1jklβ1+x2ijklβ2+ul+δ0kl+x1jklδ1kl+νjkl+ϵijkly_{ijkl}=\beta_{0}+x_{1jkl}\beta_{1}+x_{2ijkl}\beta_{2}+u_{l}+\delta_{0kl}+x_{1jkl}\delta_{1kl}+\nu_{jkl}+\epsilon_{ijkl}

  • β1\beta_{1}은 교사 자격에 대한 고정계수이며, δ1kl\delta_{1kl}은 교사 자격의 효과 중 학교별로 변하는 임의 요소를 나타냅니다.
  • 만약 δ1kl\delta_{1kl}의 분산이 크다면 석사 학위 보유의 결과가 학교마다 크게 다르다는 뜻이므로, 일괄적인 자격 요건 적용이 주 전체에 동일한 효과를 내지 못할 것임을 시사합니다.

3. 종단 데이터(Longitudinal Data) 분석

정책이나 프로그램은 멈춰있지 않습니다. 교육 프로그램 역시 여러 해에 걸쳐 단계적으로 도입되거나 폐지되기도 하므로 연구자들은 종단 데이터를 분석할 준비가 되어 있어야 합니다.

패널 데이터는 횡단면 개인들을 시간 경과에 따라 여러 번 관찰한 데이터입니다. 이 방식은 횡단면 데이터보다 프로그램에 대한 인과적 명제를 평가할 수 있는 더 나은 기회를 제공합니다. 정책 처리에 대한 반응으로서 결과의 변화를 관찰할 수 있기 때문입니다.

국가(ii)의 연도(tt)별 외국인 직접 투자(FDI)를 분석하는 임의효과 모형(Random Effects Model)의 예는 다음과 같습니다.

yit=xitβ+ui+ϵity_{it}=x_{it}\beta+u_{i}+\epsilon_{it}

  • 여기서 uiu_{i}는 각 단위(국가 등)의 무작위 절편입니다.
  • 일부 국가가 공변량으로 설명되지 않는 이유로 매년 지속적으로 더 높은 수준을 유지한다면, 이러한 단위 이질성(unit heterogeneity)을 설명하는 것이 필수적입니다.

4. 다른 모형과의 연결 및 공간적 의존성

정책 데이터는 다층분석뿐만 아니라 다른 발전된 기법과도 연결됩니다.

  • 사건사 분석(Event History Analysis): 주 단위에서 새로운 교육 정책을 채택할 때까지 걸리는 ‘시간’을 종속변수로 볼 때 사용됩니다. 시간 데이터는 정규분포를 따르지 않으며 중도절단(censoring) 문제가 발생하므로 이 기법이 유용합니다.
  • 공간 모델링(Spatial Modeling): 학교나 진료소는 지도상의 지점으로 위치를 나타낼 수 있으며 지리적 이웃 여부를 분류할 수 있습니다. 관찰되지 않은 변수나 결과 변수가 이웃 간에 더 유사한 공간적 자기상관을 보일 때, 조건부 자기회귀(CAR: conditionally autoregressive) 모형을 통해 지리적 의존성을 통제하면 보다 정확한 결과를 얻을 수 있습니다.

5. 다층분석의 한계점과 대안

다층분석이 완벽한 것은 아니며, 연구 상황에 따라 대안을 고민해야 합니다.

  1. 소표본 문제 (Micronumerosity): 상위 수준의 단위가 매우 적을 때(예: 미국 50개 주, 라틴 아메리카 20개국) 통계적 검정력이 떨어집니다. 이 경우 데이터 전체를 모집단으로 간주하고 불확실성을 유연하게 다루는 베이지안(Bayesian) 방법을 사용하는 것이 큰 이점이 될 수 있습니다.
  2. 임의효과 모형의 한계: 임의효과 모형은 단위 효과가 모든 공변량과 독립적이라고 가정합니다. 그러나 상위 수준의 단위 효과가 하위 수준 예측 변수와 상관관계가 있다면 내생성 편향이 발생합니다.
  3. 고정효과 모형(Fixed Effects Model): 이럴 때는 더미 변수를 추가하여 단위 이질성을 제어하는 고정효과 모형(LSDV)이 대안이 될 수 있습니다. 하지만 고정효과 모형은 시간이나 수준에 따라 변하지 않는 변수(time-invariant variables)의 효과를 추정할 수 없다는 치명적인 단점이 있으므로 상황에 맞게 신중히 선택해야 합니다.

적용 사례: 대기오염방지법(Clean Air Act) 집행 2001-2009년 동안의 미국 50개 주의 환경 정책 집행 건수를 분석한 연구에서는, 위에서 배운 개념들을 모두 통합하였습니다. 종속변수가 건수(count)이므로 포아송 분포를 적용했고, 시간에 따른 동적 구조(Lagged variable)와, 이웃 주들끼리의 유사성을 통제하기 위한 공간적 조건부 자기회귀(CAR) 임의효과를 포함한 베이지안 상태공간 모형을 사용했습니다.

6. 모의 데이터 분석: 학교 현장의 다층분석 실습

자, 이제 교육심리통계 교수로서, 직접 실습을 해보겠습니다.

스토리: 어느 교육청에서 도입한 ‘디지털 리터러시 프로그램’의 효과를 검증하려 합니다. 학생(Level 1) 1,000명이 50개의 학교(Level 2)에 소속되어 있습니다. 학생의 이전 국어 점수(pre_score)와 교사의 디지털 연수 여부(teacher_training)가 학생의 최종 디지털 리터러시 점수(post_score)에 미치는 영향을 확인해 보겠습니다.

R 코드를 활용한 모의 데이터 생성 및 분석

R

# 필수 패키지 로드
library(lme4)
library(lmerTest)
set.seed(2026)

# 1. 데이터 생성
n_schools <- 50
n_students_per_school <- 20
total_students <- n_schools * n_students_per_school

# 학교 수준 데이터 (Level 2)
school_id <- rep(1:n_schools, each = n_students_per_school)
# 교사 연수 여부 (학교/학급 단위로 부여된다고 가정)
teacher_training <- rep(rbinom(n_schools, 1, 0.5), each = n_students_per_school)
# 학교별 임의효과 (각 학교의 고유한 평균적 성취도 차이)
school_effect <- rep(rnorm(n_schools, mean = 0, sd = 5), each = n_students_per_school)

# 학생 수준 데이터 (Level 1)
pre_score <- rnorm(total_students, mean = 50, sd = 10)
student_error <- rnorm(total_students, mean = 0, sd = 8)

# 종속변수 계산: 최종 점수 = 절편 + 이전점수 효과 + 교사연수 효과 + 학교효과 + 오차
post_score <- 20 + (0.6 * pre_score) + (8 * teacher_training) + school_effect + student_error

df <- data.frame(school_id = as.factor(school_id),
                 student_id = 1:total_students,
                 pre_score = pre_score,
                 teacher_training = as.factor(teacher_training),
                 post_score = post_score)

# 2. 다층모형 분석 (Random Intercept Model)
# null model (ICC 확인용)
model_null <- lmer(post_score ~ 1 + (1 | school_id), data = df)
summary(model_null)

# full model
model_full <- lmer(post_score ~ pre_score + teacher_training + (1 | school_id), data = df)
summary(model_full)

Jamovi에서의 분석 방법

연구자 분들이 코딩 없이 마우스 클릭만으로 분석할 수 있는 Jamovi 가이드입니다.

  1. 위에서 만든 데이터를 CSV로 저장하여 Jamovi로 불러옵니다.
  2. 상단 메뉴의 Modules에서 GAMLj (General Analyses for Linear Models) 모듈을 설치합니다.
  3. GAMLj -> Mixed Model을 클릭합니다.
  4. Dependent Variable에 post_score를 넣습니다.
  5. Covariates에 pre_score (연속형 변수)를, Factors에 teacher_training (범주형 변수)을 넣습니다.
  6. Cluster Variables에 school_id를 넣습니다.
  7. Random Effects 탭을 열고, school_id를 오른쪽 박스로 넘겨 무작위 절편(Random Intercept)을 지정합니다.
  8. 결과를 확인하면 R의 lmer 패키지와 완벽히 동일한 다층분석 결과를 직관적인 표로 얻을 수 있습니다.

참고문헌

  • Allison, P. (2005). Fixed Effects Regression Models for Longitudinal Data Using SAS. Cary, NC: SAS Publishing.
  • Banerjee, S., Carlin, B. P., & Gelfand, A. E. (2004). Hierarchical Modeling and Analysis for Spatial Data. New York: Chapman and Hall.
  • Bryk, A. S., & Raudenbush, S. W. (1987). Application of Hierarchical Linear Models to Assessing Change. Psychological Bulletin, 101, 147-158.
  • Cameron, A. C., & Trivedi, P. K. (2005). Microeconometrics. New York: Cambridge University Press.
  • Flay, B. R. et al. (1995). The Television, School And Family Smoking Prevention And Cessation Project: VIII Student Outcomes And Mediating Variables. Preventive Medicine, 24, 29-40.
  • Monogan III, J. E. (n.d.). Modeling Policy Adoption and Impact with Multilevel Methods. In The SAGE Handbook of Multilevel Modeling (pp. 503-520).
  • Wood, B. D. (1992). Modeling Federal Implementation as a System: The Clean Air Case. American Journal of Political Science, 36, 40-67.

Chap27. 메타분석(Meta-Analysis)

안녕하세요!

오늘은 다층모형을 활용한 메타분석(Meta-Analysis)에 대해 살펴보겠습니다. “학교 현장의 데이터”를 예시로 들어 직관적인 설명과 수리적 엄밀함을 모두 갖춘 형태로 재구성해 드리겠습니다.

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

1. 메타분석이란 무엇일까요?

우리가 어떤 교육적 질문, 예를 들어 “선생님의 칭찬이 학생의 수학 성적을 올려줄까?”라는 궁금증이 생겼다고 해봅시다. 어떤 연구에서는 “엄청나게 효과가 있다”고 하고, 어떤 연구에서는 “별로 효과가 없다”고 합니다. 경험적 연구가 축적됨에 따라, 무엇이 알려져 있고 어떤 부분의 연구가 불충분하거나 모순되는지 명확히 파악하기 위해 연구 결과들을 체계화하고 종합해야 할 필요성이 생겼습니다. 메타분석은 여러 연구들의 분석 결과를 결합하여 일반적인 결론을 도출하기 위해 통계적 방법을 사용하는 활동을 뜻합니다.

  • 1976년, Glass는 교육 심리학 분야의 연구 결과들을 결합하는 맥락에서 ‘메타분석’이라는 용어를 처음 만들었습니다.
  • 이후 이 방법은 교육학, 심리학 및 기타 사회과학 분야에서 폭넓게 사용되고 있습니다.
  • 메타분석은 여러 연구에서 나온 추정치들을 모아서 하나의 요약된 결과를 얻어내는 과정입니다.

초등학생에게 설명하자면, 메타분석은 “여러 명의 요리사가 만든 떡볶이 레시피(개별 연구)를 다 모아서, 세상에서 가장 평균적이고 맛있는 궁극의 떡볶이 황금 레시피(통합된 결과)를 찾는 과정”이라고 할 수 있습니다!

2. 효과크기(Effect Size): 과일의 맛을 통일하기

연구마다 학생들의 성적을 측정하는 시험지가 다르고 만점도 다릅니다. 이럴 때 결과를 하나로 합치려면 ‘효과크기’라는 공통의 단위로 변환해야 합니다. 효과크기는 메타분석에 포함된 연구들의 결과를 요약하는 데 사용되는 수치적 지수입니다. 연구들은 각자 효과크기의 추정치와 그 추정치의 불확실성(표준오차)을 제공하게 됩니다. 이상적으로, 효과크기는 모든 연구의 결과를 “공통의 미터법(common metric)”으로 표현하여 쉽게 해석하고 결합할 수 있게 해야 합니다.

사회과학에서 가장 많이 쓰이는 세 가지 효과크기는 다음과 같습니다:

  1. 표준화된 평균 차이 (Standardized Mean Difference, dd 또는 gg)
    • 연속형 척도로 결과를 측정하는 처치나 개입의 효과를 연구할 때 자연스러운 효과크기입니다.
    • 처치 집단의 평균 결과와 통제 집단의 평균 결과의 차이를 집단 내 표준편차로 나눈 값입니다.
    • 모수 δ\delta의 공식은 다음과 같습니다: δ=μT−μCσ\delta=\frac{\mu^{T}-\mu^{C}}{\sigma}.
    • 여기서 μT\mu^{T}는 처치 집단의 모평균, μC\mu^{C}는 통제 집단의 모평균, σ\sigma는 집단 내 모표준편차입니다.
  2. 로그 오즈비 (Log Odds Ratio)
    • 이분형 척도(예: 합격/불합격)로 측정하는 연구에서 자연스러운 효과크기입니다.
  3. 상관계수 (Correlation Coefficient)
    • 두 연속형 변수 사이의 관계를 나타낼 때 쓰이며, 피셔의 zz-변환(Fisher z-transform)을 통해 분석이 진행되곤 합니다.

3. 다층모형을 이용한 메타분석 (다층분석의 세계)

메타분석 데이터를 설명하는 가장 자연스러운 방법은 2수준(two-level) 위계적 모형(다층모형)을 사용하는 것입니다.

  • 1수준 (연구 내 모형): 연구 내의 추정 오차 분산을 다룹니다.
  • 2수준 (연구 간 모형): 효과크기가 연구들 사이에서 어떻게 달라지는지(연구 간 변산)를 다룹니다.

수식으로 보면 이렇습니다. kk개의 연구가 있다고 할 때:

  • 1수준: Ti=θi+ϵiT_{i}=\theta_{i}+\epsilon_{i}. 여기서 TiT_{i}는 관찰된 효과크기, θi\theta_{i}는 진짜 효과크기 모수, ϵi\epsilon_{i}는 오차입니다.
  • 2수준: θi=β0+ηi\theta_{i}=\beta_{0}+\eta_{i}. 여기서 β0\beta_{0}는 평균 효과크기(고정 효과)이고, ηi\eta_{i}는 연구별 무선 효과입니다.

고정효과 vs 무선효과 모형

  • 고정효과 모형 (Fixed effects model): 모든 연구의 진짜 효과크기가 똑같다고 가정합니다(τ2=0\tau^2=0).
  • 무선효과 모형 (Random effects model): 연구마다 진짜 효과크기가 다를 수 있다고 인정합니다. 여기서 τ2\tau^2 (타우 제곱)이라는 값이 등장하는데, 이는 연구 간 무선 효과의 변산 정도를 나타내는 ‘연구 간 분산 성분’입니다.

4. 이질성(Heterogeneity): 연구마다 결과가 왜 다를까?

연구들끼리 결과가 얼마나 들쭉날쭉한지를 ‘이질성’이라고 합니다.

  • 분산 성분 τ2\tau^2의 추정치는 연구 간 효과크기 이질성의 “양”을 나타내는 순수한 척도입니다.
  • 이질성을 검증하기 위해 QQ 통계량을 사용하는데, 이는 τ2=0\tau^2=0이라는 가설을 검증하는 데 쓰입니다. 만약 QQ 값이 크면 연구들 사이에 차이가 크다는 뜻입니다.
  • 사용자들이 τ2\tau^2를 직관적으로 이해하기 어려워, 추정 분산과 비교하여 그 크기를 특징짓는 I2I^2 지수가 널리 사용됩니다.

5. 메타 회귀분석(Meta-Regression): 다름의 원인 찾기

만약 연구들마다 효과가 다르게 나타난다면(τ2>0\tau^2 > 0), 그 이유를 찾아야 합니다. 메타분석 연구자들이 직면하는 근본적인 문제 중 하나는 연구의 특성과 효과크기 간의 연관성을 어떻게 모델링할 것인가입니다. 예를 들어, Raudenbush(1984)는 교사의 기대가 학생의 IQ에 미치는 영향을 다룬 19개의 연구를 리뷰했습니다. 여기서 효과의 차이를 설명하기 위해 ‘교사와 학생의 사전 접촉 기간(주 단위)’이라는 공변량을 추가하여 2수준 모형을 만들었습니다.

  • 모형: θi=β0+β1xi1+⋅⋅⋅+βpxip+ηi\theta_{i}=\beta_{0}+\beta_{1}x_{i1}+\cdot\cdot\cdot+\beta_{p}x_{ip}+\eta_{i}.
  • 분석 결과, 교사가 학생과 사전 접촉이 없었을 때는 교사 기대 효과가 통계적으로 신뢰할 만하게 나타났지만, 사전 접촉 기간이 길어질수록 그 효과는 실질적으로 감소했습니다.

6. 현실의 복잡한 문제들: 데이터의 종속성과 출판 편향

  • 위계적 종속성 모형 (Hierarchical Dependence Model): 때로는 여러 연구가 같은 연구실에서 나왔거나 같은 연구자에 의해 진행되어 효과크기들이 독립적이지 않을 수 있습니다. 메타분석에서는 이러한 현상을 위계적 종속성이라고 부르며, 이를 해결하기 위해 3수준 다층모형(3-level hierarchical model)을 고려하는 것이 자연스럽습니다.
  • 출판 편향 (Publication Bias): 통계적으로나 임상적으로 유의미한 결과를 찾은 연구가 출판될 확률이 더 높다는 증거가 있습니다. 이렇게 되면 우리가 보는 결과가 실제보다 훨씬 부풀려져 보일 수 있으므로(심지어 실제 값의 200% 이상을 초과할 수도 있음), 주의 깊은 해석과 보정이 필요합니다.

7. 스토리텔링 모의 데이터 및 실습 (jamovi & R)

📖 연구 스토리: “선생님의 폭풍 칭찬은 학생의 수학 자신감을 올려줄까?”

우리는 지난 10년간 진행된 15개의 학교 현장 연구 데이터를 모았습니다. 어떤 학교(연구)에서는 칭찬의 효과가 엄청났고, 어떤 학교에서는 미미했습니다. 우리는 그 원인이 “교사의 칭찬 연수 이수 시간(Training_Hours)”에 있다고 생각하여 메타 회귀분석을 돌려보기로 했습니다.

💻 R 코드 (데이터 생성 및 분석)

jamovi의 메타분석 모듈인 MAJOR는 내부적으로 R의 metafor 패키지를 사용합니다. 아래 코드를 R에서 실행하면 모의 데이터를 생성하고 분석할 수 있습니다.

R

# 필요한 패키지 설치 및 로드
# install.packages("metafor")
library(metafor)

# 1. 모의 데이터 생성 (15개의 학교 연구)
set.seed(2026)
k <- 15 # 연구 수
Study_ID <- paste("School", 1:k)
Training_Hours <- sample(0:20, k, replace=TRUE) # 교사 연수 시간(공변량)

# 진짜 효과크기(theta) 생성: 연수 시간이 길수록 효과가 커지도록 설정
true_effect <- 0.2 + 0.05 * Training_Hours + rnorm(k, 0, 0.1) 

# 실험군(칭찬)과 대조군(일반)의 표본크기 및 요약통계량 무작위 생성
N_T <- round(runif(k, 30, 100))
N_C <- round(runif(k, 30, 100))
Mean_T <- rnorm(k, mean = 50 + true_effect * 10, sd = 2)
Mean_C <- rnorm(k, mean = 50, sd = 2)
SD_T <- runif(k, 8, 12)
SD_C <- runif(k, 8, 12)

# 데이터프레임 만들기
dat <- data.frame(Study_ID, N_T, Mean_T, SD_T, N_C, Mean_C, SD_C, Training_Hours)

# 2. 효과크기(Standardized Mean Difference, Hedges' g) 계산
dat_es <- escalc(measure="SMD", 
                 n1i=N_T, m1i=Mean_T, sd1i=SD_T, 
                 n2i=N_C, m2i=Mean_C, sd2i=SD_C, 
                 data=dat)

# 3. 무선효과 모형 메타분석 (Random Effects Model)
res_re <- rma(yi, vi, data=dat_es)
print(res_re)

# 4. 메타 회귀분석 (Meta-Regression): 연수 시간(Training_Hours) 투입
res_reg <- rma(yi, vi, mods = ~ Training_Hours, data=dat_es)
print(res_reg)

# 5. Forest Plot 시각화
forest(res_re, slab=dat_es$Study_ID, main="Forest Plot: Praise Effect on Math")

📊 jamovi에서 똑같이 분석하는 방법

  1. 모듈 설치: jamovi 우측 상단의 + 모듈(Modules) 탭을 클릭하고 MAJOR (Meta-Analysis for JAMOVI)를 설치합니다.
  2. 데이터 입력: 위 R 코드에서 생성된 데이터(효과크기 yi 와 분산 vi 혹은 평균/표준편차 원데이터)를 jamovi 스프레드시트에 불러옵니다.
  3. 메석 실행:
    • MAJOR -> Meta-Analysis 클릭.
    • Effect Size 모델을 선택하고 (ex. Correlation, SMD 등), 데이터를 알맞은 칸에 넣습니다.
    • 메타 회귀분석을 원한다면, Moderator(s) 칸에 Training_Hours를 넣습니다.
    • 옵션에서 Forest Plot을 체크하면 각 연구의 결과와 통합된 다이아몬드(평균 효과크기)를 시각적으로 확인할 수 있습니다.

결론 해석: jamovi(또는 R) 결과창에서 Training_Hours의 pp-값이 0.05보다 작게 나온다면, “선생님이 칭찬 연수를 많이 받을수록 학생 수학 자신감 향상 효과가 유의미하게 커진다”라고 결론 내릴 수 있습니다!

8. 참고문헌

  • Birge, R. T. (1932). The Calculation of Errors by the Method of Least Squares. Physical Review, 40, 207-227.
  • Borenstein, M. (2009). Effect Sizes for Continuous Data. In H. Cooper, L. V. Hedges, & J. C. Valentine (Eds.), The Handbook of Research Synthesis and Meta-analysis (2nd ed., pp. 221-236). Russell Sage Foundation.
  • Cooper, H. M. (2010). Research Synthesis and Meta-analysis: A Step By Step Approach. Sage Publications.
  • Glass, G. V. (1976). Primary, Secondary, and Meta-analysis of Research. Educational Researcher, 5, 3-8.
  • Higgins, J. P. T., & Thompson, S. G. (2002). Quantifying Heterogeneity in a Meta-analysis. Statistics in Medicine, 21, 1539-1558.
  • Raudenbush, S. W. (1984). Magnitude of Teacher Expectancy Effects on Pupil IQ as a Function of the Credibility of Expectancy Induction: A Synthesis Of Findings From 18 Experiments. Journal of Educational Psychology, 76, 85-97.
  • Raudenbush, S. W., & Bryk, A. S. (2002). Hierarchical Linear Models: Applications and Data Analysis Methods. Sage.

Chap25. 이항 내생 변수(Binary Endogenous Regressors)가 존재하는 상황에서 GEE(일반화 추정 방정식)가 강건한(Robust) 대안이 될 수 있는지

안녕하세요!

오늘은 다층분석 시 이항 내생 변수(Binary Endogenous Regressors)가 존재하는 상황에서 GEE(일반화 추정 방정식)와 다층모형(MLM)이 얼마나 강건한(Robust) 대안이 될 수 있는지에 대해 살펴보겠습니다. “학교 현장의 데이터”를 예시로 들어 직관적인 설명과 수리적 엄밀함을 모두 갖춘 형태로 재구성해 드리겠습니다.

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

1. 군집화된 데이터와 두 가지 접근법: GEE vs MLM

학교 현장의 데이터는 보통 ‘군집화(Clustered)’되어 있습니다. 예를 들어, 한 학생을 매일 관찰하는 종단 데이터의 경우, 매일의 측정치들이 ‘학생’이라는 상위 군집 안에 중첩(Nested)되어 있죠. 이러한 군집 데이터를 분석할 때 연구자는 크게 두 가지 모형을 선택할 수 있습니다.

  • 다층모형 (Multilevel Models, MLMs):
    • 다층모형은 조건부 모형(Conditional models)으로, 개별 응답에 대한 추론을 목적으로 합니다.
    • 개인 간의 관찰되지 않은 무선 변이(Random variation)를 명시적으로 설명하며, 공변량에 따라 특정 대상의 반응이 어떻게 변하는지 설명합니다.
    • 즉, 군집 내 상관관계를 설명하기 위해 집단 고유의 무선 효과(Random effects)를 모형에 포함합니다.
  • 일반화 추정 방정식 (Generalized Estimating Equations, GEE):
    • GEE는 집단 내 상관관계를 성가신 존재(nuisance)로 취급하고, 응답의 주변 평균(Marginal mean)에 초점을 맞춥니다.
    • 이는 전체 집단의 평균적인 변화, 즉 모집단 평균 모형(Population average model)을 추정하는 데 유용합니다.

2. 내생성(Endogeneity)이란 무엇인가요?

어려운 단어 같지만, 초등학생도 이해할 수 있게 설명해 볼게요. 내생성이란 ‘우리가 모형에 포함하지 않은 숨겨진 원인(오차항)과 내가 분석하려는 원인 변수가 서로 얽혀 있는 상태’를 뜻합니다.

쉽게 예를 들어볼까요?

  • 연구 목적: 학생의 ‘매일의 또래 갈등(원인)’이 ‘매일의 학업 스트레스(결과)’에 미치는 영향을 알고 싶습니다.
  • 숨겨진 진짜 원인(내생성 발생): 사실 이 학생은 원래 ‘예민한 기질(관찰되지 않은 무선 효과, uju_j)’을 가지고 있습니다. 이 기질 때문에 친구들과 갈등도 자주 겪고, 학업 스트레스도 쉽게 받습니다.

만약 이 ‘예민한 기질’을 고려하지 않고 단순하게 분석하면 어떻게 될까요? 기질 때문에 발생한 스트레스까지 모두 ‘또래 갈등’ 탓으로 돌리게 되어, 또래 갈등의 영향력을 실제보다 부풀려서(과대추정) 해석하게 됩니다. 이것이 바로 동시적 내생성(Contemporaneous endogeneity)이 일으키는 오류(Misspecification)입니다.

3. 해결책: Mundlak과 Chamberlain의 방법

이러한 내생성 문제를 해결하기 위해 Mundlak(1978)과 Chamberlain(1984)은 아주 기발하고 간단한 방법을 제안했습니다. 바로 군집 내 변수의 평균(Cluster averages)을 모형에 추가하는 것입니다.

우리의 예시로 치면, 매일매일의 ‘또래 갈등 여부’뿐만 아니라, 그 학생의 ’10일 동안 또래 갈등을 겪은 평균 횟수’를 모형에 같이 넣어주는 것이죠. 이 평균값이 통계적으로 유의미하다면, 내생성이 존재한다는 것을 의미합니다.

4. 학교 현장 모의 데이터 생성 (R 코드)

이제 학생 100명을 대상으로 10일 동안 관찰한 모의 데이터를 생성해 보겠습니다. MSCM(어머니의 스트레스와 아동의 질병) 연구의 구조를 교육 현장에 맞게 각색했습니다.

R

# 필요한 패키지 불러오기
library(lme4)
library(geepack)
library(dplyr)

set.seed(2026)
n_students <- 100
days <- 10

# 1. 학생 고유의 기질(무선 효과) 생성 - 내생성 유발
# 기질이 예민할수록(u_stress가 높을수록) 갈등도 잦고 스트레스도 높음
u_conflict <- rnorm(n_students, 0, 0.8)
u_stress <- 0.6 * u_conflict + rnorm(n_students, 0, 0.5) 

# 데이터 프레임 생성
data <- expand.grid(Day = 1:days, StudentID = 1:n_students)
data <- data %>%
  left_join(data.frame(StudentID = 1:n_students, 
                       u_conflict = u_conflict, 
                       u_stress = u_stress), by = "StudentID") %>%
  mutate(
    # 시간 불변 공변량: 교사 경력 (0=저경력, 1=고경력)
    TeacherExp = rep(rbinom(n_students, 1, 0.5), each = days),
    
    # 원인 변수 (내생 변수): 매일의 또래 갈등 (0=없음, 1=있음)
    # 기질(u_conflict)의 영향을 받음
    Conflict_star = -0.5 + 0.2*TeacherExp + u_conflict + rnorm(n_students*days, 0, 1),
    Conflict = ifelse(Conflict_star > 0, 1, 0),
    
    # 결과 변수: 매일의 학업 스트레스 (0=없음, 1=있음)
    # 또래 갈등과 기질(u_stress)의 영향을 동시에 받음
    Stress_star = -1.0 + 0.8*Conflict - 0.3*TeacherExp + u_stress + rnorm(n_students*days, 0, 1),
    Stress = ifelse(Stress_star > 0, 1, 0)
  )

# Mundlak-Chamberlain 방법을 위해 개인별 평균 또래 갈등 계산
data <- data %>%
  group_by(StudentID) %>%
  mutate(Conflict_mean = mean(Conflict)) %>%
  ungroup()

# 분석에 사용할 최종 데이터
head(data %>% select(StudentID, Day, TeacherExp, Conflict, Conflict_mean, Stress))

5. jamovi 및 R을 활용한 모형 분석

jamovi에서는 GAMLj 모듈의 Generalized Mixed Models를 통해 다층모형을 아주 쉽게 구현할 수 있습니다. GEE 분석은 jamovi 내 Rj Editor를 활용하여 R 코드를 직접 실행하는 것이 가장 정확합니다.

원문의 시뮬레이션 결과에 따르면, 내생성을 무시한 단순 GEE(Naive IWM)와 단순 다층모형(Naive MLM) 모두 내생 변수의 모수를 과대 추정하는 오류를 범합니다. 이를 직접 확인해 보겠습니다.

모형 1: 내생성을 무시한 단순 모형 (Naive Models)

  • 분석의 문제점: 개인의 기질을 통제하지 않아, 또래 갈등의 효과가 부풀려집니다.

R

# 1-1. 단순 GEE 모형 (Naive IWM)
naive_gee <- geeglm(Stress ~ Conflict + TeacherExp, 
                    id = StudentID, 
                    data = data, 
                    family = binomial(link = "probit"), 
                    corstr = "independence")
summary(naive_gee)

# 1-2. 단순 다층 모형 (Naive MLM)
# jamovi GAMLj 설정: Dependent(Stress), Factors(Conflict, TeacherExp), Random Effects(1 | StudentID)
naive_mlm <- glmer(Stress ~ Conflict + TeacherExp + (1 | StudentID), 
                   data = data, 
                   family = binomial(link = "probit"))
summary(naive_mlm)

모형 2: 평균을 추가한 보완 모형 (Augmented Models – Mundlak/Chamberlain)

  • 분석의 장점: 군집 수준의 평균 변수(Conflict_mean)를 추가하여 내생성 편향을 보정합니다.
  • 이 방식은 GEE와 다층모형 모두에서 내생성으로 인한 오류를 감지하고 수정하는 데 유사한 성능을 보입니다.

R

# 2-1. 보완된 GEE 모형 (Augmented IWM)
aug_gee <- geeglm(Stress ~ Conflict + Conflict_mean + TeacherExp, 
                  id = StudentID, 
                  data = data, 
                  family = binomial(link = "probit"), 
                  corstr = "independence")
summary(aug_gee)

# 2-2. 보완된 다층 모형 (Augmented MLM)
# jamovi GAMLj 설정: Covariates에 Conflict_mean 추가
aug_mlm <- glmer(Stress ~ Conflict + Conflict_mean + TeacherExp + (1 | StudentID), 
                 data = data, 
                 family = binomial(link = "probit"))
summary(aug_mlm)

분석 결과에서 Conflict_mean의 계수가 통계적으로 유의미하게 나타난다면, 원문에서 지적한 대로 동시적 내생성이 존재함을 확증하는 것입니다.

요약 및 결론

GEE는 가정을 적게 하기 때문에 모형 설정 오류에 덜 취약하다는 주장이 있지만, 실제로는 묵시적인 가정(오차항과 변수의 독립성 등)을 어겼을 때 심각한 편향이 발생할 수 있습니다. 내생성을 무시할 경우, GEE와 다층모형(MLM) 모두 편향된 결과를 도출하므로 둘 중 어느 하나가 무조건 더 낫다고 볼 수 없습니다.

따라서 학교 현장의 종단 데이터처럼 상호작용이 복잡한 요인들을 분석할 때는, Mundlak과 Chamberlain이 제안한 것처럼 집단 내 평균치(Cluster means)를 모형에 포함하여 내생성 여부를 테스트하고 통제하는 습관을 들이는 것이 매우 중요합니다.

참고문헌

  • Chamberlain, G. (1984). Panel Data. In Z. Griliches & M. D. Intriligator (Eds.), Handbook of Econometrics (Vol. 2, pp. 1247-1318). Elsevier.
  • Diggle, P. J., Liang, K.-Y., & Zeger, S. L. (1994). Analysis of Longitudinal Data. Clarendon Press.
  • Lee, Y., & Nelder, J. A. (2004). Conditional and marginal models: Another view. Statistical Science, 19(2), 219-228.
  • Mundlak, Y. (1978). On the pooling of time series and cross section data. Econometrica, 46(1), 69-85.
  • Zeger, S. L., & Liang, K.-Y. (1986). Longitudinal data analysis for discrete and continuous outcomes. Biometrics, 42(1), 121-130.

Chap24. 다층모형의 가정 위배 진단, 그래픽 분석, 그리고 적합도 평가

안녕하세요!

오늘은 다층모형(Multilevel Model)의 가정 위배 진단, 그래픽 분석, 그리고 적합도 평가에 대해 살펴보겠습니다. “학교 현장의 데이터”를 예시로 들어 직관적인 설명과 수리적 엄밀함을 모두 갖춘 형태로 재구성해 드리겠습니다.

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

다층모형 진단: 우리 반 아이들의 성장 곡선, 제대로 그린 걸까?

안녕하세요, 대학원생 여러분. 오늘은 다층모형(MLM)을 돌리고 나서 “과연 이 결과가 믿을 만한가?”를 검증하는 단계, 즉 모형 진단(Model Diagnostics)에 대해 배워보겠습니다.

우리가 병원에서 엑스레이를 찍고 의사 선생님이 판독하듯이, 통계 모형을 만든 후에는 이 모형이 데이터와 잘 맞는지(Lack of Fit), 가정은 지켜졌는지 확인해야 합니다. ‘읽기 유창성 성장 모형’ 예시를 통해 하나씩 짚어보겠습니다.

1. 시나리오 및 데이터 생성 (The Scenario)

📖 연구 배경: “책 읽는 호랑이 프로젝트”

우리는 A 초등학교 3학년 학생 30명을 대상으로 6학기 동안(3학년 1학기 ~ 5학년 2학기) ‘읽기 유창성(Reading Fluency)’ 점수가 어떻게 변화하는지 추적했습니다.

  • 연구 문제: 학생들의 읽기 능력은 시간이 지남에 따라 선형적으로 증가하는가? 학생마다 시작점(절편)과 성장 속도(기울기)는 다른가?

🛠️ R을 이용한 모의 데이터 생성

먼저 분석할 데이터를 만들어 봅시다. 이 데이터는 jamovi에서 불러와 실습할 수 있습니다.

R

# 필수 패키지 로드
if(!require(MASS)) install.packages("MASS")
if(!require(lme4)) install.packages("lme4")
if(!require(lattice)) install.packages("lattice")

set.seed(123)

# 1. 기본 설정
n_subjects <- 30    # 학생 수
n_timepoints <- 6   # 시점 수 (0~5)
N <- n_subjects * n_timepoints

# 2. 시간 변수 생성 (0, 1, 2, 3, 4, 5)
time <- rep(0:(n_timepoints-1), n_subjects)
subject_id <- rep(1:n_subjects, each=n_timepoints)

# 3. 무선 효과(Random Effects) 생성: 절편(초기값)과 기울기(성장률)
# 평균은 0, 분산-공분산 행렬 Sigma 설정
Sigma <- matrix(c(100, 20, 
                  20, 15), nrow=2) # 절편분산 100, 기울기분산 15, 공분산 20
random_effects <- mvrnorm(n_subjects, mu=c(0, 0), Sigma=Sigma)

u0 <- rep(random_effects[,1], each=n_timepoints) # 학생별 고유 절편
u1 <- rep(random_effects[,2], each=n_timepoints) # 학생별 고유 기울기

# 4. 고정 효과(Fixed Effects) 설정
beta0 <- 150 # 전체 평균 초기 점수
beta1 <- 10  # 학기당 평균 성장 점수

# 5. 오차항(Residuals) 생성
error <- rnorm(N, mean=0, sd=10)

# 6. 종속변수(읽기 점수) 생성
# Y = (beta0 + u0) + (beta1 + u1)*Time + error
reading_score <- (beta0 + u0) + (beta1 + u1) * time + error

# 데이터 프레임 생성
data_school <- data.frame(
  ID = as.factor(subject_id),
  Time = time,
  Score = reading_score
)

# 데이터 확인
head(data_school)

# CSV로 저장 (jamovi에서 불러오기 위함)
write.csv(data_school, "chap24.csv", row.names = FALSE)
chap24다운로드

2. 데이터 시각화: 트렐리스(Trellis) 그래프

텍스트에서는 모형을 수립하기 전, 데이터를 눈으로 확인하는 것이 매우 중요하다고 강조합니다. 이때 가장 강력한 도구가 트렐리스 플롯(Trellis Plot), 혹은 격자 그래프입니다.

🔍 왜 필요한가요?

전체 평균 그래프만 보면 개별 학생들의 독특한 패턴(누군가는 점수가 떨어지거나, 누군가는 급격히 오르는 등)을 놓칠 수 있습니다. 트렐리스 플롯은 “각 학생의 그래프를 우표처럼 모아서 한눈에 보는 것”입니다.

📊 분석 방법

  • Jamovi: Exploration > Scatterplot에서 X축에 Time, Y축에 Score를 넣고, Group에 ID를 넣으면 되지만, 30명을 다 보기는 어렵습니다.
  • R (Lattice 패키지): 텍스트에서 추천하는 방식입니다. 각 패널에 학생 한 명의 데이터와 회귀선을 그립니다.

R

# R 코드: 트렐리스 플롯 그리기
library(lattice)

xyplot(Score ~ Time | ID, data = data_school,
       type = c("p", "r"), # p: 점, r: 회귀선
       main = "학생별 읽기 유창성 성장 곡선 (Trellis Plot)",
       xlab = "시간 (학기)", ylab = "읽기 점수",
       layout = c(6, 5)) # 6열 5행으로 배치

[해석] 이 그래프를 통해 우리는 다음을 직관적으로 알 수 있습니다.

  1. 절편의 차이: 시작 점수가 낮은 아이와 높은 아이가 섞여 있나? (무선 절편 필요성 확인)
  2. 기울기의 차이: 어떤 아이는 가파르게 성장하고, 어떤 아이는 완만한가? (무선 기울기 필요성 확인)
  3. 선형성: 점들이 직선 주변에 모여 있나, 아니면 곡선 형태인가?

3. 모형의 가정과 잔차(Residuals) 분석

다층모형은 크게 세 가지 가정을 전제로 합니다.

  1. 분포 가정: 오차항(ϵ\epsilon)과 무선 효과(uu)는 정규분포를 따른다.
  2. 모형 구조: 이상치(Outlier)가 없으며 모든 관측치가 모형을 따른다.
  3. 함수 형태: 평균 반응은 선형적이다(E(y)=XβE(y) = X\beta).

3.1 잔차의 두 얼굴: 주변 잔차 vs. 조건부 잔차

다층모형에서 잔차는 두 가지로 나뉩니다. 이 부분이 초심자가 가장 헷갈려 하는 부분입니다.

  1. 주변 잔차 (Marginal Residuals, e(1)e^{(1)}):
    • y−(고정효과예측값)y – (고정효과 예측값)
    • “전체 평균 학생”과 “실제 학생 점수”의 차이입니다. 여기에는 학생 개인의 특성(무선 효과)도 오차에 포함되어 버립니다.
  2. 조건부 잔차 (Conditional Residuals, e(2)e^{(2)}):
    • y−(고정효과+무선효과예측값)y – (고정효과 + 무선효과 예측값)
    • “그 학생의 예측된 성장선”과 “실제 점수”의 차이입니다. 순수한 오차(ϵ\epsilon)에 가깝습니다.

[진단 시 주의할 점] 텍스트에 따르면, 잔차 그림만으로는 모형의 가정을 완벽히 검증하기 어렵습니다. 특히 조건부 잔차는 서로 상관관계가 있을 수 있고, 표준화 과정에서 이상치가 가려지거나 평범한 값이 이상치처럼 보일 수 있기 때문입니다 .

3.2 정규성 검정의 함정 (Q-Q Plot의 배신)

보통 잔차의 정규성을 볼 때 Q-Q Plot을 많이 씁니다. 하지만 무선 효과(Random Effects)의 정규성을 Q-Q Plot으로 확인하는 것은 위험할 수 있습니다.

  • 이유: McCulloch & Neuhaus(2011)에 따르면, 우리가 모형을 만들 때 이미 “정규분포를 따를 것이다”라고 가정하고 만들었기 때문에, 실제로는 정규분포가 아니더라도 Q-Q Plot에서는 정규분포처럼 보일 수 있습니다 .

4. 고급 진단: 무선 효과의 정규성 검정 (Order Selection Test)

그렇다면 “눈대중” 말고 진짜로 무선 효과가 정규분포인지 어떻게 알 수 있을까요? 텍스트에서는 순서 선택 검정(Order Selection Test)을 제안합니다.

💡 개념 설명 (초등학생 눈높이)

우리가 가진 데이터가 ‘종 모양(정규분포)’인지 확인하고 싶어요.

  1. 일단 종 모양(Normal)이라고 가정하고 맞는지 봅니다.
  2. 그다음 종 모양에 ‘혹’을 하나 붙여보고(Hermite polynomial 차수 추가), 더 잘 맞는지 봅니다.
  3. ‘혹’을 두 개 붙여보고 더 잘 맞는지 봅니다.
  4. 만약 ‘혹’이 없는 매끈한 종 모양일 때보다, 혹을 붙였을 때 설명력이 확 좋아진다면? “아, 이건 원래 종 모양이 아니었구나!”라고 판단하는 겁니다.

🛠️ 분석 구현 (jamovi & R)

jamovi는 아직 이 고급 검정을 메뉴로 제공하지 않습니다. 따라서 R을 사용하여 AIC(Akaike Information Criterion)를 이용한 약식 검정을 수행해 보겠습니다.

R

# R 코드: 정규성 가정 검토 (AIC 비교)
# H0: 무선 효과는 정규분포를 따른다.
model_null <- lmer(Score ~ Time + (1 + Time | ID), data = data_school, REML = FALSE)

# 비정규성을 반영하는 것은 복잡하므로, 여기서는 잔차 plot과 Q-Q plot을 통해
# 시각적으로 1차 점검을 하는 코드를 제시합니다.

# 1. 무선 효과 추출
ran_eff <- ranef(model_null)$ID

# 2. 시각적 확인 (Q-Q Plot)
par(mfrow=c(1,2))
qqnorm(ran_eff[,1], main = "Random Intercept Q-Q"); qqline(ran_eff[,1])
qqnorm(ran_eff[,2], main = "Random Slope Q-Q"); qqline(ran_eff[,2])
    • 해석: 점들이 직선 위에 잘 놓여 있다면 정규성 가정을 만족한다고 볼 수 있습니다. 만약 심하게 휘어진다면, Order Selection Test 같은 고급 기법이나 비모수적 방법(Nonparametric methods)을 고려해야 합니다.

    5. 이상치(Outlier)와 영향력 있는 관측치 탐색

    어떤 학생의 데이터가 전체 분석 결과를 왜곡하고 있지는 않을까요?

    5.1 분산 이동 모형 (Variance Shift Model)

    특정 학생이 유난히 점수의 변동폭(오차 분산)이 클 수 있습니다. 이를 분산 이동(Variance Shift)이라고 합니다.

    • 예: 철수는 기분파라 어떤 날은 100점, 어떤 날은 0점을 받습니다. 다른 학생들보다 σ2\sigma^2이 훨씬 큽니다. 이런 경우 철수는 이상치(Outlier)로 볼 수 있습니다.

    5.2 쿡의 거리 (Cook’s Distance)

    단순 회귀분석처럼 다층모형에서도 쿡의 거리를 사용하여 영향력 있는 사례(Subject)를 찾습니다.

    • 특정 학교나 특정 학생을 데이터에서 뺐을 때, 회귀계수(β\beta)가 확 바뀐다면 그 대상은 영향력이 매우 큰 것입니다.

    🛠️ 분석 구현 (jamovi)

    jamovi의 GAMLj 모듈(설치 필요)을 사용하면 쿡의 거리를 쉽게 구할 수 있습니다.

    1. Jamovi: Modules > jamovi library > GAMLj 설치.
    2. Linear Models > Mixed Model 선택.
    3. Dependent Variable: Score
    4. Covariates: Time
    5. Cluster: ID
    6. Random Effects: Intercept + Time (체크)
    7. Options (혹은 Plots) 탭에서 Cook's Distance 체크 (버전별 위치 상이, 보통 Assumption Checks에 위치).

    6. 고정 효과의 선형성 검정

    우리는 읽기 실력이 ‘직선’으로 증가한다고 가정했습니다 (E(y)=XβE(y) = X\beta). 하지만 실제로는 처음에 급격히 늘다가 고학년이 되면 정체될 수도 있습니다(곡선).

    🔍 진단 방법

    • 잔차 그림: 잔차도(Residual Plot)에서 UU자 형태나 뒤집힌 UU자 형태가 나타나면 선형성 위배입니다.
    • 해결책: Time의 제곱항(Time2Time^2)을 모형에 추가하거나, 스플라인(Spline) 회귀를 사용합니다.

    7. 결론 및 요약

    다층모형 분석을 마친 후에는 반드시 다음 질문을 던져야 합니다.

    1. 개별 데이터를 확인했는가? (Trellis Plot)
    2. 잔차는 무작위로 분포하는가? (Conditional Residuals Plot)
    3. 무선 효과는 정규분포에 가까운가? (단, Q-Q plot의 한계 인지)
    4. 전체 결과를 뒤흔드는 이상한 학생(Outlier)은 없는가? (Cook’s Distance)

    이 과정이 없다면, 아무리 복잡한 수식으로 계산된 결과라도 “모래 위에 지은 집”과 같습니다. 오늘 배운 진단 도구들을 활용하여 여러분의 연구 결과를 튼튼하게 만드시기 바랍니다.

    참고문헌 (APA Style)

    • Aerts, M., Claeskens, G., & Hart, J. D. (1999). Testing the fit of a parametric function. Journal of the American Statistical Association, 94, 869-879.
    • Claeskens, G. (2015). Lack of fit, graphics, and multilevel model diagnostics. In The SAGE Handbook of Multilevel Modeling (pp. 425-443). SAGE Publications.
    • Demidenko, E., & Stukel, T. A. (2005). Influence analysis for linear mixed-effects models. Statistics in Medicine, 24, 893-909.
    • McCulloch, C. E., & Neuhaus, J. M. (2011). Misspecifying the shape of a random effects distribution: Why getting it wrong may not matter. Statistical Science, 26, 388-402.
    • Verbeke, G., & Molenberghs, G. (2013). The gradient function as an exploratory goodness-of-fit assessment of the random effects distribution in mixed models. Biostatistics.

    Chap23. 결측치(Missing Data): 도망간 데이터 잡으러 가기

    안녕하세요!

    오늘은 결측치(Missing Data)의 처리와 다층모형을 활용한 분석에 대해 살펴보겠습니다. “학교 현장의 데이터”를 예시로 들어 직관적인 설명과 수리적 엄밀함을 모두 갖춘 형태로 재구성해 드리겠습니다.

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

    학교 현장에서 종단 연구(시간의 흐름에 따른 변화 관찰)를 하다 보면 필연적으로 결측치(Missing Data)를 마주하게 됩니다. 학생이 전학을 가거나, 아파서 결석하거나, 혹은 시험이 너무 어려워 백지를 내고 도망가는 경우들이죠.

    과거에는 컴퓨터 성능의 한계로 인해 결측치가 있는 데이터를 단순히 삭제하거나 평균으로 대충 채워 넣는 방법을 썼습니다. 하지만 이제는 EM 알고리즘이나 다중 대체(Multiple Imputation), 그리고 다층모형(Mixed Models)과 같은 강력한 도구들이 있어 더 정확한 추론이 가능해졌습니다.

    이 강의에서는 가상의 “읽기 능력 향상 프로그램” 데이터를 통해 결측치를 어떻게 다뤄야 하는지 알아보겠습니다.

    1. 시나리오: “책벌레 만들기 프로젝트”

    우리는 S 초등학교 3학년 학생 100명을 대상으로 4개월 동안 ‘읽기 능력 향상 프로그램’을 진행했습니다.

    • 측정 시점: 0개월(사전), 1개월, 2개월, 3개월 (총 4회)
    • 문제 상황: 시간이 지날수록 읽기 점수가 낮은 학생들이 흥미를 잃고 프로그램에 참여하지 않기 시작했습니다. (즉, 데이터가 빠지기 시작함)

    이 상황을 이해하기 위해 먼저 결측의 종류(메커니즘)를 알아야 합니다.

    1.1 결측 메커니즘의 이해 (기초 개념)

    이 부분은 초등학생에게 설명하듯 아주 쉽게 풀어보겠습니다.

    1. MCAR (Missing Completely At Random, 완전 무작위 결측):
      • 예시: 급식실에 가다가 넘어져서 시험지를 잃어버린 경우입니다. 학생의 읽기 능력이나 성격과는 아무 상관 없이, 순전히 운 나쁘게 데이터가 사라진 것이죠.
      • 특징: 무시하고 분석해도 결과가 크게 왜곡되지 않습니다(단지 표본 수가 줄어들 뿐).
    2. MAR (Missing At Random, 무작위 결측):
      • 예시: “지난달 점수”가 낮은 학생이 이번 달에 결석할 확률이 높은 경우입니다.
      • 핵심: 결측된 이유가 우리가 관측한 데이터(지난달 점수)로 설명이 됩니다. 우리가 이미 가지고 있는 정보로 “아, 쟤는 지난번에 못 봐서 안 왔구나”라고 추측할 수 있다면 MAR입니다. 다층모형은 이 가정하에서 아주 강력합니다.
    3. MNAR (Missing Not At Random, 비무작위 결측):
      • 예시: 오늘 시험이 너무 어려워서, 혹은 오늘 기분이 너무 우울해서 시험을 안 본 경우입니다.
      • 핵심: 결측된 이유가 측정하지 못한 값(오늘의 우울감, 오늘의 실제 실력) 자체에 의존합니다. 이건 해결하기 가장 까다로운 문제입니다.

    2. 분석 준비: 데이터 생성 (R & jamovi)

    자, 이제 R을 이용해 MAR(지난달 점수가 낮으면 결측 발생) 상황을 가정한 모의 데이터를 생성해 보겠습니다.

    [R Code] 데이터 생성

    R

    # ==============================================================================
# 1. 라이브러리 로드 및 환경 설정
# ==============================================================================
# 필요한 패키지가 없다면 install.packages("lme4") 등으로 설치해야 합니다.
if(!require(lme4)) install.packages("lme4")
if(!require(ggplot2)) install.packages("ggplot2")
if(!require(dplyr)) install.packages("dplyr")
if(!require(tidyr)) install.packages("tidyr")

library(lme4)
library(ggplot2)
library(dplyr)
library(tidyr)

set.seed(123) # 재현성을 위해 시드 고정

# ==============================================================================
# 2. 데이터 생성 (Truth)
# ==============================================================================
n_subjects <- 200      # 학생 수 (표본을 좀 늘려서 효과를 명확히 함)
n_timepoints <- 4      # 시점 수 (0, 1, 2, 3개월)
time <- 0:3
sigma_error <- 2       # 오차의 표준편차

# 개인별 성장 곡선 파라미터 (Random Effects)
# 절편 평균 50, 표준편차 5 / 기울기 평균 5, 표준편차 2
intercepts <- rnorm(n_subjects, 50, 5)
slopes <- rnorm(n_subjects, 5, 2)

data_long <- data.frame()

for(i in 1:n_subjects){
  # 진점수(True Score)
  y_true <- intercepts[i] + slopes[i] * time
  # 관측값(Observed Score) = 진점수 + 오차
  y_obs <- y_true + rnorm(n_timepoints, 0, sigma_error)
  
  subject_data <- data.frame(
    ID = as.factor(i),
    Time = time,
    Score = y_obs,
    True_Score = y_true # 비교용 진점수
  )
  data_long <- rbind(data_long, subject_data)
}

# ==============================================================================
# 3. 결측치(MAR) 생성
# ==============================================================================
# 시나리오: 직전 시점 점수가 낮으면, 그 다음 시점부터 80% 확률로 중도 탈락
data_missing <- data_long
data_missing$Missing <- FALSE

for(i in unique(data_missing$ID)){
  subj_rows <- which(data_missing$ID == i)
  
  # Time 1부터 3까지 순차적으로 확인 (Time 0은 결측 없음 가정)
  for(t in 2:4){ 
    # 직전 시점(t-1)의 점수 확인
    prev_score <- data_missing$Score[subj_rows[t-1]]
    
    # 직전 점수가 48점 미만(하위권)이면 80% 확률로 탈락 발생
    if(!is.na(prev_score) && prev_score < 48){
      if(runif(1) < 0.8){
        # 현재 시점부터 끝까지 NA 처리 (Monotone Dropout)
        data_missing$Score[subj_rows[t:4]] <- NA
        data_missing$Missing[subj_rows[t:4]] <- TRUE
        break 
      }
    }
  }
}

# 전체 결측률 확인
cat("총 결측 비율:", mean(is.na(data_missing$Score)), "\n")

# ==============================================================================
# 4. 분석 데이터셋 준비
# ==============================================================================

# [Case 1] True CC (Listwise Deletion)
# 결측치가 하나라도 있는 학생은 통째로 제거
valid_ids <- data_missing %>%
  group_by(ID) %>%
  filter(sum(is.na(Score)) == 0) %>%
  pull(ID) %>%
  unique()

data_cc_true <- data_missing %>% filter(ID %in% valid_ids)
cat("CC 분석에 포함된 학생 수:", length(unique(data_cc_true$ID)), "/", n_subjects, "\n")


# [Case 2] MM (Available Case)
# 결측이 있어도 그대로 사용 (data_missing 원본 사용)

# CSV로 저장 (jamovi에서 불러오기 위함)
write.csv(data_missing, "chap23.csv", row.names = FALSE)
    chap23다운로드

    상황 설명: 원래대로라면 모든 학생의 성적이 올랐어야 합니다(기울기 5). 하지만 성적이 낮은 학생들이 중간에 그만두게 만들었습니다(MAR). 이제 이 데이터를 분석해 봅시다.

    3. 잘못된 분석 방법들: 단순한 방법의 함정

    단순한 방법(Simple Methods)들이 얼마나 위험한지 경고하고 있습니다.

    3.1 완전 사례 분석 (Complete Case Analysis, CC)

    • 방법: 결측치가 하나라도 있는 학생은 분석에서 몽땅 제외합니다.
    • 문제점: 성적이 낮은 학생들이 주로 빠졌으므로, 끝까지 남은 학생들은 “우등생”들뿐입니다. 결과적으로 프로그램의 효과가 과대평가(Bias) 됩니다. 표본 수도 줄어들어 통계적 검정력도 떨어집니다.

    3.2 마지막 관측값 이월 (Last Observation Carried Forward, LOCF)

    • 방법: 1개월 차에 그만둔 학생의 점수를 2, 3개월 차에도 똑같이 복사해 넣습니다.
    • 문제점: “학생의 성장은 멈췄다”고 가정하는 것입니다. 성장 곡선 모형에서는 이 방법이 성장의 기울기를 과소평가하게 만들거나, 오차의 분산을 왜곡시킵니다. 아주 보수적인(효과가 없다고 하는) 결과를 내는 것 같지만, 실제로는 엉뚱한 결론을 낼 수 있어 위험합니다.

    4. 올바른 분석 방법: MAR 가정하의 접근

    4.1 직접 우도 분석 (Direct Likelihood) / 다층모형 (Mixed Models)

    이 방법은 우리가 데이터를 억지로 채워 넣거나 지우지 않고, 있는 데이터 그대로(observed data) 분석하는 것입니다. 다층모형(Linear Mixed Model)은 결측이 무작위(MAR)로 발생했다면, 관측된 데이터와 변수 간의 상관관계를 이용해 편향되지 않은 추정치를 제공합니다.

    [jamovi 실습 가이드] 다층모형 분석

    1. 데이터 불러오기: 위에서 생성한 reading_program_MAR.csv를 엽니다.
    2. 분석 메뉴: Analyses > Linear Models > Mixed Models를 선택합니다.
    3. 변수 설정:
      • Dependent Variable (종속변수): Score
      • Covariates (공변량): Time
      • Cluster (군집 변수): ID (학생 개인)
    4. Random Effects (무선 효과):
      • Time과 Intercept를 Random Coefficients에 넣습니다. (학생마다 초기치와 성장 속도가 다르므로)
    5. 결과 해석:
      • jamovi는 내부적으로 lme4 또는 nlme 패키지를 사용하며, 이는 결측치가 있는 경우에도 REML(제한된 최대우도) 방식을 통해 사용 가능한 모든 데이터를 활용해 파라미터를 추정합니다.

    [R Code 비교]

    R

    # ==============================================================================
# 5. 다층모형 적합 (Model Fitting)
# ==============================================================================

# 1) Truth (결측 없는 원본 데이터 - 기준점)
fit_truth <- lmer(Score ~ Time + (Time | ID), data = data_long)

# 2) True CC (편향된 데이터)
fit_cc_true <- lmer(Score ~ Time + (Time | ID), data = data_cc_true)

# 3) MM (결측 있는 데이터 + 다층모형)
fit_mm <- lmer(Score ~ Time + (Time | ID), data = data_missing)


# ==============================================================================
# 6. 결과 비교 및 출력
# ==============================================================================

# 고정 효과(Fixed Effects) 추출 함수
get_fixed <- function(model, name) {
  coefs <- fixef(model)
  data.frame(
    Model = name,
    Intercept = coefs["(Intercept)"],
    Slope = coefs["Time"]
  )
}

res_truth <- get_fixed(fit_truth, "1. Truth (기준)")
res_cc <- get_fixed(fit_cc_true, "2. True CC (편향)")
res_mm <- get_fixed(fit_mm, "3. MM (보정)")

# 결과 합치기
results <- rbind(res_truth, res_cc, res_mm)
rownames(results) <- NULL

print("===== [모형별 고정 효과(Fixed Effects) 추정치 비교] =====")
print(results)
cat("\n* 설명: Truth의 Slope(약 5.0)에 가까울수록 좋은 분석입니다.\n")
cat("* True CC는 하위권 학생이 제거되어 Slope가 과대평가될 가능성이 높습니다.\n")


# ==============================================================================
# 7. 시각화 (Visualization)
# ==============================================================================

# 예측된 회귀선 생성을 위한 데이터
pred_data <- data.frame(Time = seq(0, 3, 0.1))

# 각 모형의 예측값 계산 (Fixed Effect만 사용)
pred_data$Truth <- fixef(fit_truth)[1] + fixef(fit_truth)[2] * pred_data$Time
pred_data$CC    <- fixef(fit_cc_true)[1] + fixef(fit_cc_true)[2] * pred_data$Time
pred_data$MM    <- fixef(fit_mm)[1] + fixef(fit_mm)[2] * pred_data$Time

# 데이터를 긴 형태(Long Format)로 변환 (ggplot용)
plot_data <- pred_data %>%
  pivot_longer(cols = c("Truth", "CC", "MM"), names_to = "Model", values_to = "Score")

# 그래프 그리기
ggplot(plot_data, aes(x = Time, y = Score, color = Model, linetype = Model)) +
  geom_line(size = 1.2) +
  scale_color_manual(values = c("Truth" = "black", "CC" = "red", "MM" = "blue")) +
  scale_linetype_manual(values = c("Truth" = "solid", "CC" = "dashed", "MM" = "dotdash")) +
  labs(title = "모형별 성장 곡선 비교",
       subtitle = "검정(Truth) vs 빨강(CC: 과대추정) vs 파랑(MM: MAR 보정)",
       y = "Reading Score", x = "Time (Months)") +
  theme_minimal() +
  theme(legend.position = "bottom")

    결과 비교 시사점:

    • 진짜 값(True): 기울기(Time) ≈ 5.0
    • 완전 사례(CC): 낮은 점수 학생이 빠져서 남은 우등생들의 가파른 성장만 반영될 수 있음(혹은 중도 탈락자가 성장세가 꺾인 아이들이라면 오히려 과대 추정). 시뮬레이션에 따라 다르지만 보통 편향됨.
      * 다층모형(MM): 결측이 있어도 기울기가 5.0에 가깝게 복원됩니다. 이것이 바로 다층모형의 힘입니다!

    4.2 다중 대체 (Multiple Imputation, MI)

    결측치를 한 번만 채우는 게 아니라, M번(보통 5~10번) 채워서 M개의 데이터셋을 만듭니다.

    1. Imputation: 그럴듯한 값들로 빈칸을 채웁니다(불확실성 반영).
    2. Analysis: M개의 데이터셋을 각각 분석합니다.
    3. Pooling: M개의 결과를 하나로 합칩니다.

    jamovi에서는 기본 모듈에서 MI를 직접 지원하는 파이프라인이 약하지만, Rj 모듈을 쓰거나 R을 이용해 수행할 수 있습니다. 교재에서는 다층모형(Direct Likelihood)과 MI가 MAR 가정하에서 둘 다 타당하다고 합니다. 다만, 공변량(Covariate) 자체에 결측이 있을 때는 MI가 더 유용할 수 있습니다.

    5. 고급 주제: 비무작위 결측 (MNAR)과 민감도 분석

    만약 학생이 “그냥 기분이 나빠서(우리가 측정 안 한 변수)” 빠졌다면(MNAR) 어떻게 할까요? 이때는 다층모형도 완벽하지 않습니다.

    5.1 패턴 혼합 모형 (Pattern-Mixture Models)

    학생들을 “탈락 시점”에 따라 그룹(패턴)으로 나누어 분석하는 방법입니다.

    • 예: 끝까지 남은 그룹, 2개월 차에 나간 그룹, 3개월 차에 나간 그룹.
    • 각 그룹별로 성장 곡선이 어떻게 다른지 봅니다.
    • 교재의 Vorozole 연구 예시처럼, 탈락 패턴을 공변량으로 넣어 분석할 수 있습니다.

    5.2 선택 모형 (Selection Models)

    성적 모형과 탈락 확률 모형(로지스틱 회귀 등)을 동시에 결합해서 분석하는 아주 복잡한 방법입니다. Heckman의 Tobit 모델과 유사한 개념입니다.

    5.3 민감도 분석 (Sensitivity Analysis)

    우리가 가정한 결측 메커니즘(MAR)이 틀렸을 수도 있습니다. 그래서 “만약 이게 MNAR이라면 결과가 어떻게 바뀔까?”를 테스트해보는 것입니다.

    • jamovi나 R에서 분석한 뒤, “결측된 값들이 관측된 값보다 평균적으로 k점 더 낮다고 가정해보자”라고 시나리오를 바꾸어가며 결과가 뒤집히는지 확인합니다.

    6. 결론 및 제언

    교재와 분석 결과를 종합하여 여러분께 드리는 제언은 다음과 같습니다.

    1. LOCF나 CC는 이제 그만: 이 방법들은 편향이 심하므로 주 분석으로 쓰지 마십시오.
    2. 다층모형(Direct Likelihood)을 기본으로: jamovi의 Mixed Models는 별도의 조치 없이도 MAR 가정하에서 훌륭한 결과를 냅니다. 이것을 일차적 분석(Primary Analysis)으로 삼으세요.
    3. 민감도 분석 수행: 결측이 MNAR일 가능성을 배제할 수 없으므로, 다양한 가정하에서도 결과가 견고한지 확인하는 것이 현명합니다.

    학교 현장의 데이터는 아이들의 숨겨진 사연(결측)을 담고 있습니다. 무작정 지우기보다 그 패턴을 이해하고 적절한 통계적 도구를 사용하는 것이 연구자의 윤리이자 책무입니다.

    참고문헌 (APA Style)

    • Beunckens, C., Molenberghs, G., & Kenward, M. G. (2005). Tutorial: Direct likelihood analysis versus simple forms of imputation for missing data in randomized clinical trials. Clinical Trials, 2(4), 379–386.
    • Little, R. J. A., & Rubin, D. B. (2002). Statistical analysis with missing data (2nd ed.). New York: John Wiley & Sons.
    • Molenberghs, G., & Kenward, M. G. (2007). Missing data in clinical studies. Chichester: John Wiley & Sons.
    • Molenberghs, G., & Verbeke, G. (2012). Missing Data. In The SAGE Handbook of Multilevel Modeling (pp. 403–424). SAGE Publications.
    • Verbeke, G., & Molenberghs, G. (2000). Linear mixed models for longitudinal data. New York: Springer-Verlag.

    Chap22. 다층모형에서의 강건한 분석 방법

    안녕하세요!

    오늘은 다층모형(Multilevel Modeling)에서 가정 위배에 대처하는 ‘강건한(Robust) 분석 방법’에 대해 살펴보겠습니다. “학교 현장의 데이터”를 예시로 들어 직관적인 설명과 수리적 엄밀함을 모두 갖춘 형태로 재구성해 드리겠습니다.

    분석 도구로는 jamovi의 사용법을 설명하되, 고급 강건성 분석(Robust SE, Bootstrapping, Bayesian) 등 jamovi 기본 기능의 한계를 넘어서는 부분은 R 코드를 함께 제시하여 모의 데이터 생성부터 분석, 시각화까지 완벽하게 구현해 드리겠습니다.

    1. 서론: 왜 ‘강건한(Robust)’ 방법이 필요한가요?

    교육 현장의 데이터는 교과서처럼 깔끔하지 않습니다.

    예를 들어, 어떤 반에는 유독 책을 많이 읽는 ‘독서 영재(이상치, Outlier)’가 있을 수 있고, 어떤 학교는 학생 수가 너무 적을 수도 있습니다.

    통계 분석을 할 때 우리는 흔히 “데이터가 정규분포를 따를 것이다(Normality)”, “분산이 일정할 것이다(Homoscedasticity)”라는 가정을 합니다. 하지만 현실 데이터가 이 약속을 어긴다면 어떻게 될까요? 우리가 내린 결론(통계적 유의성)이 틀릴 수도 있습니다.

    이때 필요한 것이 바로 강건한(Robust) 방법입니다.

    강건성(Robustness)이란? 데이터가 기본 가정(정규성 등)을 다소 위배하더라도, 통계적 추정 결과나 신뢰구간이 크게 흔들리지 않고 믿을 만한 결과를 내놓는 성질을 말합니다.

    이제 가상의 ‘해바라기 초등학교’ 데이터를 만들어 이 문제를 해결해 봅시다.

    2. 모의 데이터 생성: 해바라기 초등학교 이야기

    우리는 “교사의 정서적 지지(Teacher Support)가 학생의 학업 스트레스(Stress)를 낮추는가?”를 연구하려고 합니다.

    • 1수준(학생): 학업 스트레스(종속변수 YY), 학생의 불안 성향(공변량 XX)
    • 2수준(학급): 교사의 정서적 지지(WW)
    • 시나리오: 대부분의 학생은 평범하지만, 몇몇 학생은 불안 성향이 매우 높습니다(데이터가 한쪽으로 찌그러짐, Skewed). 또한, 몇몇 반은 학생 수가 적습니다.

    먼저 R을 이용해 이 상황을 반영한 데이터를 생성해보겠습니다.

    R

    # R Code: 데이터 생성
set.seed(2026)

# 50개 학급(Class), 학급당 학생 수 5~25명(불균형)
n_classes <- 50
class_id <- rep(1:n_classes, times = sample(5:25, n_classes, replace = TRUE))
n_students <- length(class_id)

# 2수준 변수: 교사의 정서적 지지 (랜덤 효과 포함)
class_intercept <- rnorm(n_classes, mean = 50, sd = 10)
teacher_support <- rnorm(n_classes, mean = 5, sd = 2)
class_effect <- class_intercept[class_id]
support_level <- teacher_support[class_id]

# 1수준 변수: 학생의 불안 성향 (치우친 분포 생성 - Gamma 분포 활용)
student_anxiety <- rgamma(n_students, shape = 2, scale = 2) 

# 이상치(Outlier) 추가: 몇몇 학생은 극도로 높은 불안을 가짐
student_anxiety[sample(1:n_students, 5)] <- 30 

# 종속변수 생성 (Stress)
# 모델: Stress = 50 - 2*Support + 1.5*Anxiety + Random_Error
error <- rnorm(n_students, mean = 0, sd = 5)
stress <- 50 - 2 * support_level + 1.5 * student_anxiety + class_effect + error

# 데이터 프레임 생성
data <- data.frame(
  ClassID = as.factor(class_id),
  Stress = stress,
  TeacherSupport = support_level,
  Anxiety = student_anxiety
)

# 데이터 확인
head(data)

# CSV로 저장 (jamovi에서 불러오기 위함)
write.csv(data, "chap22.csv", row.names = FALSE)
    chap22다운로드

    3. 진단(Diagnosis): 데이터의 건강검진

    분석 전에 데이터가 가정을 잘 지키는지 확인해야 합니다. 책에서는 이상치(Outliers)와 다중공선성(Multicollinearity) 확인을 강조합니다.

    3.1. 이상치 및 분포 확인 (Jamovi & R)

    다층모형에서는 전체 데이터뿐만 아니라 수준별(Level-wise)로 잔차(Residual)를 확인해야 합니다. 하지만 SPSS나 일반 소프트웨어는 이를 잘 지원하지 않는 경우가 많습니다.

    • Jamovi: [Exploration] -> [Descriptives]에서 ‘Box plot’과 ‘Q-Q plot’을 체크합니다.
    • R 시각화:

    R

    # 상자수염그림(Boxplot)으로 이상치 확인
boxplot(data$Stress, main="학업 스트레스 분포(이상치 확인)")
# Q-Q plot으로 정규성 확인
qqnorm(data$Stress); qqline(data$Stress, col = "red")

    [해설] 만약 상자수염그림(Boxplot) 위쪽에 동그라미(Outlier)가 찍혀 있거나, Q-Q plot에서 점들이 붉은 선을 벗어난다면 정규성 가정이 위배된 것입니다. 우리 데이터에는 일부러 이상치를 넣었으므로 위배되었을 가능성이 큽니다.

    4. 해결책 A: 강건한 표준오차 (Robust Standard Errors)

    데이터가 정규성을 위배하거나 등분산성을 만족하지 못할 때 사용할 수 있는 첫 번째 무기는 ‘강건한 표준오차(Robust SE)’입니다. 샌드위치 추정량(Sandwich Estimator)이라고도 부릅니다.

    4.1. 개념 설명

    일반적인 최대우도법(ML)은 데이터가 정규분포라고 가정하고 표준오차를 계산합니다. 가정이 틀리면 표준오차가 너무 작게 계산되어, 실제로는 효과가 없는데도 “유의하다”고 잘못 결론 내릴 수 있습니다. Robust SE는 데이터의 잔차(Residual)를 직접 이용하여 표준오차를 보정합니다. 이를 통해 가정이 깨져도 올바른 검정 결과를 얻을 수 있습니다.

    주의점: Robust SE가 제대로 작동하려면 2수준(학급)의 표본 크기가 충분해야 합니다(최소 50~100개 그룹 권장).

    4.2. 분석 방법 (R)

    Jamovi의 기본 Mixed Model 모듈은 아직 Robust SE를 직접 지원하지 않습니다. 따라서 R의 lme4와 clubSandwich 패키지를 사용합니다.

    R

    # 패키지 로드
library(lme4)
library(clubSandwich)

# 1. 일반적인 다층모형 적합 (Random Intercept Model)
model_lmer <- lmer(Stress ~ TeacherSupport + Anxiety + (1|ClassID), data = data)

# 2. 결과 확인 (일반 SE)
summary(model_lmer)

# 3. 강건한 표준오차(Robust SE) 구하기 (CR2 방식 권장)
robust_results <- coef_test(model_lmer, vcov = "CR2")
print(robust_results)

    [결과 해석] 일반적인 summary 결과와 coef_test 결과를 비교해 보세요. 만약 정규성 위배가 심각하다면, Robust SE를 적용했을 때 표준오차(SE)가 커지고 p-value가 변할 수 있습니다. 이것이 더 보수적이고 안전한 결과입니다.

    5. 해결책 B: 부트스트래핑 (Bootstrapping)

    두 번째 방법은 컴퓨터의 힘을 빌리는 부트스트래핑입니다. 이론적인 분포(정규분포)에 의존하지 않고, 내 데이터를 수천 번 재표집(Resampling)하여 분포를 직접 만들어내는 방식입니다.

    5.1. 개념 설명

    우리의 데이터가 모집단 그 자체라고 가정하고, 여기서 복원추출을 통해 수천 개의 가상 데이터를 만듭니다.

    • 다층모형에서의 주의점: 단순히 케이스를 뽑으면 안 됩니다. 학급 구조가 깨지기 때문입니다. 잔차 부트스트래핑(Residual Bootstrapping)이 가장 선호되는 방법입니다.

    5.2. 분석 방법 (R)

    이 역시 계산량이 많아 R을 사용합니다. lme4 패키지의 bootMer 함수를 사용합니다.

    R

    # 부트스트랩 함수 정의 (고정 효과 추출)
mySumm <- function(.) {
  fixef(.)
}

# 부트스트래핑 실행 (nsim = 1000번 반복)
# use.u = TRUE는 랜덤 효과도 고려한다는 뜻
boot_results <- bootMer(model_lmer, mySumm, nsim = 1000, use.u = TRUE, type = "parametric")

# 95% 신뢰구간 계산
boot_CI <- boot.ci(boot_results, type = "perc", index = 1) # Intercept에 대한 예시
print(boot_results)

    [장점] 부트스트래핑은 특히 매개효과 분석이나 분산 성분에 대한 신뢰구간을 구할 때 매우 강력합니다. 정규성 가정이 크게 위배되어도 믿을 만한 결과를 줍니다.

    6. 해결책 C: 베이지안 추정 (Bayesian Estimation)

    세 번째, 가장 강력하고 유연한 방법인 베이지안 추정입니다. 특히 표본 크기가 작을 때(학급 수가 적을 때) 매우 유용합니다.

    6.1. 개념 설명

    전통적인 통계(빈도주의)는 모수가 고정되어 있다고 보지만, 베이지안은 모수도 확률분포를 가진다고 봅니다.

    • 사전 정보(Prior): 연구자가 미리 알고 있는 정보.
    • 우도(Likelihood): 데이터에서 얻은 정보.
    • 사후 분포(Posterior): 이 둘을 합쳐 업데이트된 결론.

    강건성을 위해 정규분포 대신 꼬리가 긴 t-분포(t-distribution)를 가정하여 이상치의 영향을 줄일 수 있습니다.

    6.2. 분석 방법 (R – brms 패키지)

    brms 패키지는 베이지안 다층모형을 매우 쉽게 구현해 줍니다.

    R

    library(brms)

# 베이지안 다층모형 적합
# student: 정규분포 대신 t-분포를 사용하여 이상치에 강건하게 만듦
model_bayes <- brm(
  Stress ~ TeacherSupport + Anxiety + (1|ClassID), 
  data = data,
  family = student(),  # 핵심: 강건성을 위해 student t 분포 사용
  chains = 2, iter = 2000
)

# 결과 확인
summary(model_bayes)
plot(model_bayes) # 수렴 확인 (Trace plot)

    [해설] 베이지안 분석은 p-value 대신 신뢰구간(Credible Interval)을 제공합니다. 95% 구간 안에 0이 포함되지 않으면 “효과가 있다”고 봅니다. 특히 family = student() 옵션을 사용하면 이상치가 있어도 회귀계수가 크게 휘둘리지 않습니다.

    7. 결론 및 요약

    학교 현장 데이터처럼 가정이 위배되기 쉬운 자료를 다룰 때, 맹목적으로 일반적인 다층분석을 돌리는 것은 위험합니다.

    1. 진단하라: Jamovi나 R의 플롯을 통해 이상치와 분포를 확인하세요.
    2. 강건한 표준오차(Robust SE): 이상치가 있거나 이분산성이 의심될 때, 학급 수가 50개 이상이라면 가장 간편한 대안입니다.
    3. 부트스트래핑: 분산 추정치가 중요하거나 분포 가정을 피하고 싶을 때 유용합니다.
    4. 베이지안 추정: 학급 수가 적거나(20~50개 미만), 모형이 복잡하거나, 극단적인 이상치가 많을 때 최고의 선택입니다.

    선생님, 연구자 여러분의 데이터가 조금 거칠더라도 걱정하지 마세요. 우리에게는 이처럼 든든한 ‘강건한 도구’들이 있습니다.

    참고문헌 (APA Style)

    • Eliason, S. R. (1993). Maximum likelihood estimation. Sage.
    • Hox, J. J., & van de Schoot, R. (2013). Robust methods for multilevel analysis. In The SAGE Handbook of Multilevel Modeling (pp. 387-402). SAGE Publications Ltd.
    • Maas, C. J. M., & Hox, J. J. (2004). Robustness issues in multilevel regression analysis. Statistica Neerlandica, 58(2), 127-137.
    • Maas, C. J. M., & Hox, J. J. (2005). Sufficient sample sizes for multilevel modeling. Methodology, 1(3), 85-91.
    • Tabachnick, B. G., & Fidell, L. S. (2007). Using multivariate statistics (5th ed.). Allyn & Bacon.
    • Van der Leeden, R., Meijer, E., & Busing, F. M. T. A. (2008). Resampling multilevel models. In J. de Leeuw & E. Meijer (Eds.), Handbook of multilevel analysis (pp. 401-433). Springer.