안녕하세요!

오늘은 사회연결망 및 관계형 데이터의 다층모형(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): 이러한 네트워크 데이터는 정방행렬인 $Y$로 표현할 수 있으며, 여기서 $Y_{ij}$는 행위자 i로부터 행위자 j로 향하는 관계를 나타냅니다. 무방향 네트워크 데이터의 경우 이 행렬은 대칭을 이룹니다.

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

• 개인 중심 연결망(Personal or Egocentric Network): 표본으로 추출된 개인이 하나 이상의 타인(alter)과의 관계에 대해 보고하는 데이터입니다. (예: “지훈아, 네 친한 친구들을 다 적어볼래?”)
• 양자 데이터(Dyadic Data): 두 명의 행위자와 그들 사이의 관계인 양자(dyad), 즉 $(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수준(타인 수준)의 회귀 방정식은 다음과 같습니다.

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

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

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

$y_{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: 양자 데이터의 다층모형 (너와 나의 연결 고리)

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

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

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

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

• 이러한 교차 내재 구조를 다루기 위해 사회적 관계 모형(Social Relations Model, SRM)을 사용합니다.
• 이 모형은 행위자의 두 가지 역할(송신자와 수신자)을 모두 고려하기 위해 $u_{1}$과 $u_{2}$라는 무선 효과 벡터를 추가합니다. $y_{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)인 다층 로지스틱 모형을 사용할 수 있습니다.
• 양자 데이터의 상호성을 반영하기 위해 고정된 송신자/수신자 매개변수와 명시적인 상호성 매개변수를 사용하는 $p_{1}$ 모형이나 , 이를 확장하여 정규 분포를 가정하는 상관된 무선 송신자/수신자 매개변수를 갖는 $p_{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). $p_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.

안녕하세요!

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

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

---

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

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

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

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

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

---

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

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

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

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

$y(s) = x(s)’\beta + u(s) + \epsilon(s)$

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

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

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

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

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

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

$y(s) = X(s)’\beta + u(s) + \epsilon(s)$

이때 구조는 동일하지만, $y(s)$, $u(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.

---

안녕하세요!

오늘은 생존분석과 노쇠모형(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)를 추정할 것입니다. 순간 위험률이 높다는 것은 짧은 시간 구간 안에서 사건(자퇴)이 발생할 확률이 높다는 것을 의미합니다.

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

$h_{ij}(t)=h_{0}(t)\exp[x_{ij}(t)\beta+u_{j}]$

여기서 $h_0(t)$는 기저 위험(baseline hazard)이며, $u_j$는 학교(군집)별 취약성을 나타내는 무작위 효과로 정규분포 $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)

나. 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. 노쇠모형(다층분석) 실행: 다층 구조의 학교 효과($u_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.26	0.17	.12
SES	-0.32	0.07	< .001	SES 점수가 1점 높을수록 자퇴 위험이 유의미하게 감소합니다.
SchoolID (Frailty Variance, $\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.

안녕하세요!

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

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

---

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

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

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

---

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

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

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

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

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

$y_{ijkl}=x_{ijkl}\beta+u_{l}+\delta_{kl}+v_{jkl}+\epsilon_{ijkl}$

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

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

$y_{ijkl}=\beta_{0}+x_{1jkl}\beta_{1}+x_{2ijkl}\beta_{2}+u_{l}+\delta_{0kl}+x_{1jkl}\delta_{1kl}+\nu_{jkl}+\epsilon_{ijkl}$

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

---

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

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

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

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

$y_{it}=x_{it}\beta+u_{i}+\epsilon_{it}$

• 여기서 $u_{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.

안녕하세요!

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

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

---

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

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

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

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

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

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

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

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

---

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

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

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

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

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

고정효과 vs 무선효과 모형

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

---

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

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

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

---

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

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

• 모형: $\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의 $p$-값이 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.

안녕하세요!

오늘은 다층분석 시 이항 내생 변수(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)이란 무엇인가요?

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

쉽게 예를 들어볼까요?

• 연구 목적: 학생의 ‘매일의 또래 갈등(원인)’이 ‘매일의 학업 스트레스(결과)’에 미치는 영향을 알고 싶습니다.
• 숨겨진 진짜 원인(내생성 발생): 사실 이 학생은 원래 ‘예민한 기질(관찰되지 않은 무선 효과, $u_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.

---

안녕하세요!

오늘은 다층모형(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)

---

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$)과 무선 효과($u$)는 정규분포를 따른다.
2. 모형 구조: 이상치(Outlier)가 없으며 모든 관측치가 모형을 따른다.
3. 함수 형태: 평균 반응은 선형적이다($E(y) = X\beta$).

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

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

1. 주변 잔차 (Marginal Residuals, $e^{(1)}$):
   • $y – (고정효과 예측값)$
   • “전체 평균 학생”과 “실제 학생 점수”의 차이입니다. 여기에는 학생 개인의 특성(무선 효과)도 오차에 포함되어 버립니다.
2. 조건부 잔차 (Conditional Residuals, $e^{(2)}$):
   • $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점을 받습니다. 다른 학생들보다 $\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\beta$). 하지만 실제로는 처음에 급격히 늘다가 고학년이 되면 정체될 수도 있습니다(곡선).

🔍 진단 방법

• 잔차 그림: 잔차도(Residual Plot)에서 $U$자 형태나 뒤집힌 $U$자 형태가 나타나면 선형성 위배입니다.
• 해결책: Time의 제곱항($Time^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.

---