태그 보관물: 다층분석

Chap20. 종단적 상황 및 기타 상황에서의 혼합 모형과 잠재계층모형

안녕하세요!

오늘은 “종단적 상황 및 기타 상황에서의 혼합 모형과 잠재 계층 모형(Mixture and Latent Class Models in Longitudinal and Other Settings)”에 대해 살펴보겠습니다. “학교 현장의 데이터”를 예시로 들어 직관적인 설명과 수리적 엄밀함을 모두 갖춘 형태로 재구성해 드리겠습니다.

분석 도구로는 jamovi의 사용법을 설명하되, jamovi의 메뉴만으로는 본 챕터에서 다루는 고도화된 ‘종단적 혼합 모형(Longitudinal Mixture Models)’을 완벽히 구현하기 어렵기 때문에, jamovi 내에서 구동 가능한 R (Rj Editor) 코드를 함께 제시하여 모의 데이터 생성부터 분석, 시각화까지 완벽하게 구현해 드리겠습니다.


1. 서론: 우리는 왜 ‘섞인 것(Mixture)’을 풀어헤쳐야 할까요?

학교에서 아이들을 가르치다 보면, 겉보기에는 똑같은 ‘3학년 1반’ 학생들이지만 그 안에는 서로 다른 성향을 가진 아이들이 섞여 있다는 것을 느끼게 됩니다.

  • 다층 모형(Multilevel Model): “우리 반(1반)과 옆 반(2반)은 평균 점수가 다른가?” 처럼 이미 알고 있는 집단(Group)의 차이를 분석합니다.
  • 잠재 계층/혼합 모형(Latent Class/Mixture Model): “우리 반 안에 ‘꾸준히 성장하는 아이들’, ‘초반에 잘하다 떨어지는 아이들’, ‘뒤늦게 치고 올라오는 아이들’이 섞여 있지 않을까?” 처럼 눈에 보이지 않는 집단(Latent Group)을 찾아냅니다.

2. 시나리오: “독서 능력 성장 프로젝트”

초등학교 3학년 학생 300명을 대상으로 1학기 초부터 2학기 말까지 총 6회에 걸쳐 독서 유창성(Reading Fluency) 검사를 실시했다고 가정해 봅시다.

우리의 연구 질문은 다음과 같습니다.

“전체 평균을 내면 아이들이 조금씩 성장하는 것처럼 보이지만, 사실 그 안에는 전혀 다른 성장 패턴을 보이는 하위 집단들이 숨어 있지 않을까?”

이 질문을 해결하기 위해 본 챕터에서 소개하는 종단적 모델 기반 클러스터링(Model-Based Clustering for Longitudinal Data)을 수행해 보겠습니다.


3. 데이터 생성 및 탐색 (R & jamovi)

본 챕터의 분석 예제는 longclust 패키지를 사용했습니다. 교육학적 맥락에 맞게 데이터를 생성해 보겠습니다.

[데이터 생성 스토리]

우리는 3개의 잠재 집단이 있다고 가정합니다.

  1. 성취도 상위-지속 성장형 (High achievers): 처음부터 잘하고 계속 잘함.
  2. 초기 부진-급성장형 (Catch-up group): 처음엔 못했지만 급격히 성장함.
  3. 학습 부진-정체형 (Struggling group): 낮게 시작해서 변화가 거의 없음.

아래 코드를 jamovi의 Rj EditorRStudio에 붙여넣으면 실습 데이터를 만들 수 있습니다.

R

# 필요한 패키지 로드 (없으면 설치 필요)
if(!require(mvtnorm)) install.packages("mvtnorm")
if(!require(longclust)) install.packages("longclust") # 챕터에서 강조한 핵심 패키지 [cite: 201]

set.seed(123)

# 시간(Time points): 6회 측정
time <- 1:6

# 3개 집단의 평균 패턴 정의 (독서 유창성 점수)
mu1 <- c(70, 72, 75, 78, 82, 85) # 지속 성장
mu2 <- c(40, 42, 55, 65, 75, 80) # 급성장 (Catch-up)
mu3 <- c(45, 46, 47, 45, 46, 48) # 정체 (Struggling)

# 공분산 구조 생성 (종단 자료의 상관성 반영)
# 간단한 AR(1) 구조 유사하게 설정
sigma_gen <- function(rho, dim=6) {
  mat <- diag(dim)
  mat <- rho^abs(row(mat)-col(mat))
  return(mat * 10) # 분산 확대
}

S1 <- sigma_gen(0.7)
S2 <- sigma_gen(0.5)
S3 <- sigma_gen(0.3)

# 데이터 생성 (각 집단 100명씩)
data1 <- rmvnorm(100, mean = mu1, sigma = S1)
data2 <- rmvnorm(100, mean = mu2, sigma = S2)
data3 <- rmvnorm(100, mean = mu3, sigma = S3)

# 전체 데이터 합치기
reading_data <- rbind(data1, data2, data3)
colnames(reading_data) <- paste0("Time", 1:6)

# 시각화를 위한 데이터 준비
matplot(t(reading_data), type="l", col="grey", lty=1, 
        main="전체 학생의 독서 유창성 변화 (스파게티 플롯)",
        xlab="측정 시기", ylab="점수")

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

위의 그래프(스파게티 플롯)를 보면 선들이 엉켜 있어서 누가 어떤 집단인지 알기 어렵습니다. 이제 혼합 모형을 사용해 이 엉킨 실타래를 풀어보겠습니다.


4. 이론적 배경: 가우시안 혼합 모형의 ‘가족들’

분석에 앞서, 이 챕터에서 중요하게 다루는 “모형들의 가족(Families of Mixture Models)” 개념을 아주 쉽게 설명해 드릴게요.

우리가 데이터를 분류할 때, 컴퓨터에게 “비슷한 애들끼리 묶어봐”라고 시키면 컴퓨터는 ‘모양(Shape)’, ‘부피(Volume)’, ‘방향(Orientation)’을 기준으로 묶습니다.

1) MCLUST 패밀리 (가장 유명한 모형)

MCLUST는 데이터 덩어리(클러스터)를 “풍선”이라고 생각합니다.

  • Spherical (구형): 모든 풍선이 동그란 모양입니다.
  • Ellipsoidal (타원형): 풍선이 길쭉할 수도 있습니다.
  • Equal Volume: 모든 풍선의 크기가 같습니다.
  • Variable Volume: 어떤 풍선은 크고 어떤 건 작습니다.
  • Orientation: 풍선이 놓인 방향이 다릅니다.

이 조합에 따라 EII, VII, VVV 같은 암호 같은 이름이 붙습니다. 예를 들어 VVV는 “크기도, 모양도, 방향도 제각각인 집단들”을 허용하는 가장 유연한(하지만 복잡한) 모형입니다.

2) 종단적 데이터(Longitudinal Data)를 위한 모형

시간에 따라 반복 측정된 데이터(위의 독서 점수 같은)는 특별합니다. 1차 시기 점수가 2차 시기 점수에 영향을 주기 때문입니다. 이를 위해 McNicholas와 Murphy(2010a)는 수정된 촐레스키 분해(Modified Cholesky Decomposition)라는 수학적 마법을 부립니다.

쉽게 말해, “현재 점수를 과거 점수로 설명하고 남은 찌꺼기(혁신, Innovation)”만 분석하겠다는 것입니다. 이 방법을 쓰면 종단 데이터의 복잡한 시간 관계를 아주 효율적으로 계산할 수 있습니다.


5. 분석 실습: 숨겨진 집단 찾아내기

이제 실제로 분석을 수행해 보겠습니다. 텍스트에서 소개한 longclust 알고리즘을 사용합니다. jamovi Rj Editor에서 실행하세요.

분석 단계 1: 모델 적합 (Model Fitting)

우리는 몇 개의 집단이 있는지 모릅니다. 그래서 컴퓨터에게 “집단이 2개일 때부터 5개일 때까지 다 해보고 제일 좋은 걸 알려줘”라고 시킵니다. 이때 가장 좋은 모델을 고르는 기준은 BIC(Bayesian Information Criterion)입니다. BIC 점수는 낮을수록(절대값이 클수록 좋음, 보통 음수면 더 작은 값) 좋습니다.

R

# longclustEM 함수 실행
# G = 2:5 (집단 수를 2개에서 5개까지 탐색)
# linearMeans = FALSE (성장 곡선이 꼭 직선일 필요는 없음)

fit_result <- longclustEM(reading_data, 2, 5, linearMeans=FALSE)

# 결과 확인
summary(fit_result)

분석 단계 2: 결과 해석 (Results)

실행 결과, 컴퓨터가 다음과 같이 답했다고 가정해 봅시다(시뮬레이션 결과).

“가장 좋은 모델은 VVI 구조를 가진 3개의 집단 모형입니다.”

여기서 VVI란?

  • V (Variable): 집단마다 변동성(혁신 분산)이 다름.
  • V (Variable): 시간 간의 관계(자기회귀 파라미터)가 집단마다 다름.
  • I (Isotropic): 각 시점의 잔차 분산은 등방성임.

이것은 “학생 그룹마다 성장하는 패턴의 변동폭도 다르고, 이전 점수가 다음 점수에 미치는 영향력도 다르다”는 교육학적으로 매우 타당한 결과입니다.

분석 단계 3: 시각화 (Visualization)

각 집단의 패턴을 그려보겠습니다.

R

# 분류된 집단 정보 추출
predicted_class <- apply(fit_result$zbest, 1, which.max)

# 시각화
par(mfrow=c(1,3)) # 그래프 3개 나란히 그리기

for(g in 1:3){
  matplot(t(reading_data[predicted_class == g, ]), type="l", 
          main=paste("Group", g), ylab="Score", xlab="Time",
          col=ifelse(g==1, "blue", ifelse(g==2, "red", "green")), 
          ylim=c(30, 100))
}

[결과 해석 시나리오]

  • Group 1 (Blue): 그래프가 우상향하며 점수가 높습니다. -> “우수 집단”
  • Group 2 (Red): 바닥에서 횡보합니다. -> “집중 지원 대상 집단”
  • Group 3 (Green): 낮게 시작해서 가파르게 올라갑니다. -> “급성장 집단”

단순히 평균만 비교했다면(다층 모형의 고정 효과), Group 2와 Group 3의 차이를 발견하지 못하고 “중간 정도 하는 아이들”로 퉁쳐버렸을지도 모릅니다. 혼합 모형이 숨겨진 진실을 밝혀낸 것이죠.


6. 심화: 현실적인 문제들 (결측치와 공변량)

학교 데이터는 완벽하지 않습니다. 아이가 전학을 가거나 아파서 시험을 못 볼 수도 있죠.

1) 결측치(Missing Data) 처리

본 텍스트에서는 결측치가 있어도 EM 알고리즘을 통해 분석이 가능하다고 합니다.

  • Shaikh et al.(2010)은 PEM(Pseudo-EM) 알고리즘을 제안했는데, 이는 결측치의 평균은 고려하되 분산 계산은 단순화하여 계산 속도를 높인 방법입니다.
  • 즉, “철수가 3교시 시험을 안 봤어도, 1, 2교시 성적과 철수랑 비슷한 그룹(클러스터) 아이들의 점수를 참고해서 철수의 소속 집단을 추정”할 수 있다는 뜻입니다.

2) 공변량(Covariates)의 포함

단순히 그룹만 나누는 게 아니라, “왜 이 그룹에 속하게 되었는가?”를 알고 싶을 수 있습니다.

  • Vermunt와 Magidson(2002)은 잠재 계층 모형에 공변량(z)을 포함시켰습니다.
  • 예: “가정 형편(SES)”이나 “사교육 여부”를 공변량으로 넣어서, 이것이 “급성장 집단”에 속할 확률에 영향을 미치는지 분석할 수 있습니다.

7. 가우시안이 아닌 모형들 (Non-Gaussian Approaches)

모든 시험 점수가 예쁜 종 모양(정규분포)을 그리지는 않습니다.

  • t-분포 혼합 모형: 데이터에 이상치(Outlier)가 많을 때 유용합니다. 꼬리가 두꺼운 분포를 사용하여 극단적인 점수를 가진 학생 때문에 전체 분석이 왜곡되는 것을 막아줍니다.
  • 왜도(Skewness) 모형: 점수가 한쪽으로 쏠려 있을 때(예: 시험이 너무 쉬워서 다들 100점 근처인 경우) 사용합니다.

이 챕터에서는 teigen이나 longclust 패키지가 이러한 비정규 분포(t-분포 등)도 지원한다고 강조합니다.


8. 요약 및 제언

오늘 우리는 다층 모형의 틀 안에서 혼합 모형(Mixture Model)이 어떻게 숨겨진 이질성(Heterogeneity)을 찾아내는지 살펴보았습니다.

  1. 관점의 전환: 다층 모형이 ‘반(Class)’이라는 보여지는 집단을 분석한다면, 혼합 모형은 데이터 패턴 속에 숨어 있는 ‘유형(Type)’을 찾아냅니다.
  2. 도구의 중요성: 가우시안 혼합 모형(GMM)이 기본이지만, 종단 데이터에는 Modified Cholesky 분해를 이용한 모형이 훨씬 효율적입니다.
  3. 유연성: 결측치가 있거나 데이터가 정규분포가 아니어도(t-분포 등) 적용할 수 있는 다양한 방법론이 개발되어 있습니다.

학교 현장에서 “우리 반 아이들은 다 달라요”라고 말할 때, 이제는 감(feeling)이 아니라 혼합 모형을 통해 그 ‘다름’의 실체를 과학적으로 증명해 보시는 건 어떨까요?


참고문헌 (APA Style)

  • Banfield, J. D., & Raftery, A. E. (1993). Model-based Gaussian and non-Gaussian clustering. Biometrics, 49, 803–821.
  • Browne, R. P., & McNicholas, P. D. (2015). Mixture and latent class models in longitudinal and other settings. In The SAGE Handbook of Multilevel Modeling (pp. 357–370). SAGE Publications.
  • Celeux, G., & Govaert, G. (1995). Gaussian parsimonious clustering models. Pattern Recognition, 28(5), 781–793.
  • McNicholas, P. D., & Murphy, T. B. (2010a). Model-based clustering of longitudinal data. The Canadian Journal of Statistics, 38(1), 153–168.
  • McNicholas, P. D., & Subedi, S. (2012). Clustering gene expression time course data using mixtures of multivariate t-distributions. Journal of Statistical Planning and Inference, 142(5), 1114–1127.
  • Vermunt, J. K., & Magidson, J. (2002). Latent class cluster analysis. In J. Hagenaars & A. McCutcheon (Eds.), Applied Latent Class Analysis (pp. 89–106). Cambridge University Press.

Chap19. 위계적 동적 모델(Hierarchical Dynamic Models)

안녕하세요!

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

분석 도구로는 jamovi의 사용법(GAMLj 모듈 활용)을 설명하되, jamovi의 기반이 되는 R 코드를 함께 제시하여 모의 데이터 생성부터 베이지안(Bayesian) 접근을 통한 정교한 분석, 시각화까지 완벽하게 구현해 드리겠습니다.


1. 서론: 사진이 아니라 영화를 찍자!

우리가 학교에서 학생들의 성취도를 연구한다고 상상해 봅시다. 보통은 특정 시점에 시험을 보고, “A 학교가 B 학교보다 점수가 높다”거나 “사교육이 성적에 영향을 미친다”라고 분석합니다. 이것은 마치 ‘사진(Snapshot)’을 찍는 것과 같습니다.

하지만 현실은 훨씬 복잡합니다.

  1. 위계성 (Hierarchy): 학생들은 ‘학교’라는 그룹에 속해 있어 같은 학교 친구들끼리는 비슷해지는 경향이 있습니다.
  2. 역동성 (Dynamics): 학교의 여건(선생님의 열정, 학교의 재정 상태 등)은 고정된 것이 아니라 매년 변합니다.

위계적 동적 모형(HDM)은 이 두 가지를 합친 것입니다. 학생 수준의 회귀분석을 하되, 그 회귀계수(Intercept, Slope)가 학교마다 다르고, 심지어 시간이 지남에 따라 변하도록(Time-varying) 허용하는 모형입니다. 즉, 데이터를 사진이 아닌 ‘영화(Movie)’처럼 분석하는 것이죠.


2. 가상의 시나리오: “경기 미래형 학교 연구”

초등학생도 이해할 수 있는 예시를 만들어 보겠습니다.

연구 상황: 경기도 교육청에서 ‘학생의 수학 성취도’‘자기주도학습 시간’에 의해 얼마나 향상되는지 알아보고자 합니다.

  • 대상: 초등학교 4학년 학생들을 6학년이 될 때까지 3년간(총 6학기) 추적 조사.
  • 구조: 30개의 학교, 각 학교당 20명의 학생.
  • 핵심 질문: 자기주도학습이 성적에 미치는 효과(β\beta)는 모든 학교에서 동일할까? 그리고 그 효과는 학기가 지날수록 변할까?

3. 모형의 구조 (수리적 이해)

제공된 텍스트에 근거하여 이 모형은 크게 세 가지 방정식으로 구성됩니다.

(1) 관측 방정식 (Observation Equation)

학생 ii, 학교 jj, 시간 tt에서의 수학 점수(yijty_{ijt})는 다음과 같습니다.

yijt=β0,jt+β1,jt(자기주도학습)ijt+ϵijt,ϵijtN(0,σ2)y_{ijt} = \beta_{0,jt} + \beta_{1,jt} \cdot (\text{자기주도학습})_{ijt} + \epsilon_{ijt}, \quad \epsilon_{ijt} \sim N(0, \sigma^2)

  • 여기서 β0,jt\beta_{0,jt} (절편)와 β1,jt\beta_{1,jt} (자기주도학습 효과)는 학교(jj)와 시간(tt)에 따라 다릅니다.

(2) 구조 방정식 (Structural Equation)

학교 수준의 파라미터 βjt\beta_{jt}는 학교의 특성(예: 교사의 경력 XschoolX_{school})에 영향을 받습니다.

βjt=θt+γt(교사경력)jt+ujt,ujtN(0,Σu)\beta_{jt} = \theta_t + \gamma_t \cdot (\text{교사경력})_{jt} + u_{jt}, \quad u_{jt} \sim N(0, \Sigma_u)

(3) 시스템 방정식 (System Equation)

이 모형의 꽃입니다. 파라미터가 시간이 지남에 따라 어떻게 변하는지 설명합니다. 보통 랜덤 워크(Random Walk) 과정을 따릅니다.

θt=θt1+wt,wtN(0,Σw)\theta_t = \theta_{t-1} + w_t, \quad w_t \sim N(0, \Sigma_w)

  • 즉, 오늘의 효과는 어제의 효과에 기초하되, 새로운 변화(wtw_t)가 더해집니다.

4. 데이터 생성 및 분석 (R & jamovi)

이론적으로 HDM은 베이지안 추론(Bayesian Inference)을 주로 사용합니다. jamovi의 기본 모듈은 빈도주의(Frequentist) 기반이지만, GAMLj 모듈을 사용하면 유사한 ‘위계적 성장 모형’ 분석이 가능하며, R을 사용하면 텍스트에 나온 완벽한 동적 모형을 구현할 수 있습니다.

Step 1. R을 이용한 모의 데이터 생성

먼저, 시나리오에 맞는 데이터를 생성해 보겠습니다.

R

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

set.seed(2026)

# 1. 기본 설정
n_schools <- 30      # 학교 수
n_students <- 20     # 학교당 학생 수
n_time <- 6          # 6학기 (시간)
total_obs <- n_schools * n_students * n_time

# 2. 시스템 방정식 (시간에 따른 효과 변화 생성)
# 자기주도학습의 효과(Slope)가 시간이 지날수록 조금씩 증가한다고 가정 (랜덤 워크)
time_effect <- cumsum(rnorm(n_time, mean = 0.05, sd = 0.1)) 
base_intercept <- 50 + cumsum(rnorm(n_time, mean = 1, sd = 0.5)) # 전체 평균 점수도 상승

# 3. 데이터 프레임 생성
data <- expand.grid(
  Time = 1:n_time,
  Student_ID = 1:n_students,
  School_ID = 1:n_schools
)

# 4. 학교별 효과 (Random Intercept & Slope)
school_intercepts <- rnorm(n_schools, 0, 5) # 학교 간 격차
school_slopes <- rnorm(n_schools, 0, 0.5)   # 학교별 학습 효과 차이

# 5. 학생 데이터 생성 및 점수 계산
data <- data %>%
  mutate(
    # 학생별 자기주도학습 시간 (시간에 따라 변함 + 개인차)
    Self_Study = abs(rnorm(total_obs, mean = 5, sd = 2)), 
    
    # 실제 수학 점수 생성 (HDM 구조 반영)
    # 점수 = (시간별기초 + 학교별차이) + (시간별효과 + 학교별차이)*학습시간 + 오차
    Math_Score = (base_intercept[Time] + school_intercepts[School_ID]) + 
                 (1.5 + time_effect[Time] + school_slopes[School_ID]) * Self_Study + 
                 rnorm(total_obs, 0, 3)
  )

# 데이터 확인
head(data)

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

Step 2. jamovi를 이용한 분석 (GAMLj 모듈)

jamovi에서는 [Linear Models] > [Mixed Model] (혹은 GAMLj 모듈 설치 후 사용)를 사용하여 이 데이터를 분석할 수 있습니다. 텍스트에서 언급된 베이지안 동적 모형의 근사치인 ‘성장 곡선 모형(Growth Curve Model)’ 형태로 분석합니다.

  1. 데이터 불러오기: 위에서 생성한 csv 파일을 jamovi에서 엽니다.
  2. 분석 설정:
    • Dependent Variable: Math_Score
    • Factors (Fixed Effects): Time (Factor로 취급하거나 Covariate로 취급하여 선형 성장 가정 가능)
    • Covariates: Self_Study
    • Cluster: School_ID, Student_ID (학생이 반복 측정되었으므로)
  3. Random Effects:
    • School_ID 아래에 InterceptSelf_Study 추가 (학교별로 효과가 다름을 반영 ).
    • Student_ID (within School) 설정.

참고: jamovi의 Mixed Model은 파라미터가 ‘랜덤 워크’로 변하는 시스템 방정식(식 19.4 )을 직접 추정하기엔 한계가 있습니다. 이를 완벽히 구현하려면 아래의 R 코드(Bayesian)가 필요합니다.

Step 3. R을 이용한 베이지안 위계적 동적 모형 분석 (brms)

텍스트에서 강조하는 베이지안 추론 을 구현하기 위해 R의 brms 패키지를 사용합니다. 이는 MCMC 알고리즘 을 사용하여 복잡한 사후 분포를 추정합니다.

R

# 베이지안 분석을 위한 brms 패키지 (Stan 기반)
# install.packages("brms")
library(brms)

# 모형 설정:
# 점수 ~ 시간 + 학습시간 + (1 + 학습시간 | 학교) + (1 | 학생)
# *중요: 시간에 따른 계수의 변화를 허용하기 위해 상호작용항 또는 s() 함수 사용 가능
# 텍스트의 DLM(Dynamic Linear Model)을 근사하는 식

model_hdm <- brm(
  formula = Math_Score ~ Time * Self_Study + (1 + Time + Self_Study | School_ID),
  data = data,
  family = gaussian(),
  chains = 2, iter = 2000, # 예시를 위해 반복 수 줄임
  seed = 2026
)

# 결과 요약
summary(model_hdm)

# 조건부 효과 시각화 (시간에 따른 학습 효과의 변화)
conditional_effects(model_hdm, effects = "Self_Study:Time")

5. 분석 결과의 해석 및 시각화

(1) 시간과 학교에 따른 변화

우리는 시간에 따라 변하는 학교별 특성을 시각화할 수 있습니다.

R

# R을 이용한 시각화 (학교별, 시간에 따른 성적 변화 스파게티 플롯)
ggplot(data, aes(x = Time, y = Math_Score, group = interaction(School_ID, Student_ID))) +
  geom_line(alpha = 0.1, color = "gray") + # 개별 학생
  stat_summary(aes(group = School_ID, color = as.factor(School_ID)), 
               fun = mean, geom = "line", size = 1, alpha = 0.8) + # 학교 평균
  theme_minimal() +
  labs(title = "학교별 시간에 따른 수학 성취도 변화",
       x = "학기 (Time)", y = "수학 점수", color = "학교 ID") +
  theme(legend.position = "none")

이 그래프를 보면, 어떤 학교는 시간이 지날수록 성적이 가파르게 오르고(기울기가 큼), 어떤 학교는 정체되어 있음을 알 수 있습니다. 이것이 바로 β0,jt\beta_{0,jt} (절편)와 파라미터들이 동적으로 변한다는 증거입니다.

(2) 파라미터의 동적 변화 (Smoothed Distribution)

분석 결과, 자기주도학습의 효과(β1,t\beta_{1,t})가 1학기에는 1.5점 상승 효과였는데, 6학기에는 2.0점 상승 효과로 변했다면, 이는 “고학년이 될수록 자기주도학습이 더 중요하다”는 교육적 결론을 도출할 수 있게 합니다.


6. 심화: 모형의 확장

이 모형은 단순히 점수 하나만 분석하는 것에 그치지 않고 확장될 수 있습니다.

(1) 다변량 확장 (Matrix-Variate)

수학 점수뿐만 아니라 영어 점수도 같이 분석하고 싶다면? 수학과 영어는 상관관계가 높습니다. 이를 각각 분석하는 것보다 행렬(Matrix) 형태로 묶어서 분석하면 두 과목 간의 상관성까지 파악할 수 있어 훨씬 정확합니다.

(2) 공간적 확장 (Spatial Structure)

학교들은 지리적으로 인접해 있습니다. 예를 들어, 수원시에 있는 학교들은 용인시에 있는 학교들보다 서로 더 비슷한 교육 환경을 공유할 수 있습니다. 학교 간의 ‘거리(Distance)’ 정보를 모형의 분산 행렬(VV)에 반영하여, 가까운 학교끼리 더 높은 상관관계를 갖도록 설정할 수 있습니다.

(3) 지수족 분포 확장 (Exponential Family)

만약 종속변수가 점수(정규분포)가 아니라, ‘결석 일수(Count)’‘비만 여부(Binary)’라면 어떨까요? 이 경우 정규분포 가정이 깨집니다. 이때는 포아송 분포나 감마 분포 등을 사용하는 일반화 선형 모형(GLM)과 결합하여 분석할 수 있습니다.


7. 결론

위계적 동적 모형은 교육 현장의 데이터를 ‘있는 그대로’ 가장 잘 반영하는 도구입니다. 학생은 학교에 속해 있고(위계), 학교와 학생은 끊임없이 변화하기(동적) 때문입니다.

jamovi와 R을 활용한 이 분석을 통해, 여러분은 단순한 평균 비교를 넘어 “언제, 어디서, 어떻게 교육 효과가 변화하는지”를 밝혀내는 통찰력을 갖게 될 것입니다.


참고문헌 (APA Style)

  • Gamerman, D., & Migon, H. S. (1993). Dynamic hierarchical models. Journal of the Royal Statistical Society: Series B (Methodological), 55(3), 629-642.
  • Hansen, N. (2009). Models with dynamic coefficients varying in space for data in the exponential family [Unpublished master’s thesis]. Universidade Federal do Rio de Janeiro.
  • Landim, F. M., & Gamerman, D. (2000). Dynamic hierarchical models: An extension to matrix variate observations. Computational Statistics & Data Analysis, 35(1), 11-42.
  • Paez, M. S., & Gamerman, D. (n.d.). Hierarchical dynamic models. In The SAGE Handbook of Multilevel Modeling (Chapter 19, pp. 335-355).
  • Paez, M. S., Gamerman, D., Landim, F. M., & Salazar, E. (2008). Spatially varying dynamic coefficient models. Journal of Statistical Planning and Inference, 138(4), 1038-1058.
  • West, M., & Harrison, P. J. (1997). Bayesian forecasting and dynamic models (2nd ed.). Springer.

Chap18. 페널티 스플라인(Penalized Splines)과 다층모형

안녕하세요!

오늘은 패널 데이터 분석을 위한 페널티 스플라인(Penalized Splines)과 다층 모형의 결합에 대해 살펴보겠습니다. “학교 현장의 데이터”를 예시로 들어 직관적인 설명과 수리적 엄밀함을 모두 갖춘 형태로 재구성해 드리겠습니다.

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


1. 들어가며: 왜 ‘스플라인’이 필요한가요?

우리가 학교에서 학생들의 성장을 관찰하다 보면, 학습 시간과 성적의 관계가 항상 곧은 직선(선형)으로 나타나지 않는다는 것을 알게 됩니다.

  • 선형 모델의 한계: “공부 시간이 늘어날수록 성적도 일정하게 오른다”고 가정합니다.
  • 비선형의 현실: 처음에는 성적이 확 오르다가, 어느 정도 수준에 도달하면 정체기(plateau)가 오고, 너무 과도한 학습은 오히려 효율을 떨어뜨리기도 합니다.

이처럼 구불구불한 관계를 수식으로 나타낼 때 유용한 것이 바로 스플라인(Splines)입니다. 특히 데이터의 ‘꿈틀거림(wiggliness)’을 적절히 조절해주는 페널티 스플라인(Penalized Splines)은 복잡한 데이터 구조에서도 매우 안정적인 결과를 보여줍니다.


2. 시나리오: “학생들의 자기주도학습 시간과 학업 성취도”

우리는 서울시 소재 50개 초등학교에서 5학년 학생 1,000명을 대상으로 3년간 추적 조사를 했다고 가정해 봅시다.

  • 연구 질문: “자기주도학습 시간(Experience)이 늘어남에 따라 국어 성적(Score)은 어떻게 변화하는가? 학생 개인별로 성적의 기초선(Intercept)은 다른가?”

모의 데이터 생성 (R 코드)

실제 분석을 위해 학교 현장과 유사한 데이터를 생성해 보겠습니다.

R

# 필요한 라이브러리 로드
library(mgcv)
library(ggplot2)

set.seed(2026)
n_students <- 200 # 학생 수
obs_per_student <- 5 # 학생당 측정 횟수
total_obs <- n_students * obs_per_student

# 학생 ID 및 학교 배경 생성
student_id <- rep(1:n_students, each = obs_per_student)
# 학습 시간 (0~20시간)
study_time <- runif(total_obs, 0, 20)

# 다층 모형 성분: 학생별 랜덤 효과 (기초 실력 차이)
u_i <- rnorm(n_students, 0, 5) 
student_effect <- rep(u_i, each = obs_per_student)

# 비선형 함수: m(x) - 처음엔 상승하다가 완만해지는 곡선
m_x <- 10 * sin(study_time / 5) + 2 * study_time

# 성적 생성 (절편 + 비선형 곡선 + 학생별 차이 + 오차)
score <- 50 + m_x + student_effect + rnorm(total_obs, 0, 3)

# 데이터 프레임 구축
df <- data.frame(student_id = factor(student_id), study_time, score)

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

3. 이론적 배경: 다층 페널티 스플라인

이 모델의 핵심은 우리가 흔히 아는 다층 모형(HLM)의 회귀 계수 자리에 ‘함수’를 집어넣는 것입니다.

수리적 모델

yit=β0i+m(xit)+ϵity_{it} = \beta_{0i} + m(x_{it}) + \epsilon_{it}

β0i=β0+ui\beta_{0i} = \beta_{0} + u_{i}

  • yity_{it}: ii번째 학생의 tt시점 성적.
  • β0i\beta_{0i}: 학생 ii의 개인적인 성적 수준(랜덤 절편).
  • m(xit)m(x_{it}): 학습 시간에 따른 공통적인 성적 변화 곡선(비선형 함수).
  • uiN(0,σu2)u_{i} \sim N(0, \sigma_u^2): 학생들 간의 편차.

여기서 m(x)m(x)기저 함수(basis functions)들의 합으로 표현됩니다. 마치 레고 블록을 조합해 곡선을 만드는 것과 같습니다. 이때 너무 복잡한 곡선이 되지 않도록 페널티(λ\lambda)를 주어 부드럽게 만듭니다.


4. jamovi 및 R을 이용한 분석 실습

1) jamovi 활용법

‘Rj’ 모듈 사용

  1. 모듈 설치: jamovi 라이브러리에서 Rj – Editor to run R code inside jamovi를 설치합니다.
  2. Rj Editor 실행: 분석 메뉴에서 Rj 아이콘을 클릭합니다.
  3. 코드 입력: 아래 코드를 복사해서 붙여넣고 실행(삼각형 버튼)합니다.
# jamovi의 현재 데이터를 'data'라는 이름으로 가져옵니다.
# 'score'는 성적, 'study_time'은 학습시간, 'student_id'는 학생번호입니다.

library(mgcv)

# 다층 GAM 모델 (페널티 스플라인 + 랜덤 절편)
# s(study_time)이 페널티 스플라인을 의미합니다.
model <- gam(score ~ s(study_time) + s(student_id, bs="re"), 
             data = data, 
             method = "REML")

# 결과 출력
summary(model)

# 그래프 그리기
plot(model, pages=1, shade=TRUE, main="학습 시간에 따른 성적 변화 곡선")

2) R (gamm4) 활용법

본문에서 추천하는 방식은 gamm4 패키지를 사용하는 것입니다. 이는 혼합 모형의 안정성과 스플라인의 유연성을 동시에 확보합니다.

R

# GAMM 모델 적합
library(gamm4)
model <- gamm4(score ~ s(study_time), random = ~(1|student_id), data = df)

# 결과 확인
summary(model$gam) # 비선형 부분 확인
summary(model$mer) # 다층 모형(랜덤 효과) 부분 확인

5. 결과 해석 및 시각화

분석 결과, 우리는 다음과 같은 사실을 발견할 수 있습니다.

  • 분산 성분: 학생 간 편차(σ^u\hat{\sigma}_u)가 잔차(σ^ϵ\hat{\sigma}_\epsilon)보다 크다면, 학생의 개인적 배경이 성적에 큰 영향을 미치고 있음을 의미합니다.
  • 곡선의 형태: 학습 초기에는 성적이 급격히 상승하지만, 특정 시간(예: 15시간) 이후에는 상승 폭이 줄어드는 ‘수확 체감’의 형태를 보일 수 있습니다.

다층 모형의 장점

단순한 곡선 회귀와 달리, 이 모델은 반복 측정된 데이터의 상관관계를 고려합니다. 즉, ‘철수’가 첫 번째 시험에서 잘 봤다면 두 번째 시험에서도 잘 볼 가능성이 높다는 점을 모델이 인지하고 분석하므로 훨씬 정확합니다.


6. 결론 및 교육적 시사점

페널티 스플라인 다층 모형은 학교 현장의 복잡한 데이터를 분석하는 데 강력한 도구입니다.

  1. 개별화된 교육 과정: 학생 개개인의 기초선이 다르다는 것을 인정하면서도(uiu_i), 보편적인 성장의 패턴(m(x)m(x))을 찾아낼 수 있습니다.
  2. 적정 지점의 발견: 성적이 정체되는 시점을 시각적으로 확인하여 교육적 개입의 시기를 결정할 수 있습니다.

7. 참고문헌 (APA Style)

  • Becker, G. S. (1993). Human capital: A theoretical and empirical analysis, with special reference to education (3rd ed.). University of Chicago Press.
  • Eilers, P. H. C., & Marx, B. D. (1996). Flexible smoothing with B-splines and penalties. Statistical Science, 11(2), 89–121.
  • Kauermann, G., & Kuhlenkasper, T. (2011). Penalized splines and multilevel models. In J. J. Hox & J. K. Roberts (Eds.), The SAGE Handbook of Multilevel Modeling (pp. 325–333). SAGE Publications.
  • Ruppert, D., Wand, M. P., & Carroll, R. J. (2009). Semiparametric regression during 2003–2007. Electronic Journal of Statistics, 3, 1193–1256.
  • Wood, S. N. (2006). Generalized additive models: An introduction with R. Chapman & Hall/CRC.

Chap17.평활화 및 준모수 모형(Smoothing and Semiparametric Models)

안녕하세요!

오늘은 준개별화된 학습 성장 곡선: 준모수 혼합 효과 모델(Semiparametric Mixed-Effects Models)에 대해 살펴보겠습니다. “학교 현장의 데이터”를 예시로 들어 직관적인 설명과 수리적 엄밀함을 모두 갖춘 형태로 재구성해 드리겠습니다.

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


1. 들어가며: 왜 ‘준모수’ 모델이 필요한가요?

학교에서 학생들의 성적이나 심리적 변화를 추적하는 종단적 연구(Longitudinal Study)를 하다 보면, 데이터가 우리가 생각하는 예쁜 직선이나 단순한 곡선(2차 함수 등)을 따르지 않을 때가 많습니다.

  • 모수 모델(Parametric Models): 직선이나 포물선처럼 식을 딱 정해놓고 분석합니다. 설명이 간결하지만, 실제 데이터가 그 모양이 아니면 편향된 결과가 나옵니다.
  • 비모수 모델(Nonparametric Models): 데이터가 생긴 모양 그대로를 유연하게 따라갑니다. 하지만 수식이 너무 복잡해서 해석하기가 어렵습니다.
  • 준모수 모델(Semiparametric Models): 이 둘의 장점만 합친 것입니다. 중요한 영향 요인(예: 성별, 기초학력)은 수치로 깔끔하게 분석하고, 복잡하게 변하는 시간에 따른 변화량은 유연하게 곡선으로 그려냅니다.

2. 시나리오: “창의적 문제해결력 신장 프로그램”

한 초등학교에서 200명의 학생을 대상으로 2년간 8차례에 걸쳐 ‘창의적 문제해결력’ 점수를 측정했다고 가정해 봅시다.

  • 고정 효과(중요 요인): 성별(Gender), 프로그램 참여 전 기초 역량(Pre-Score).
  • 복잡한 변화: 학생들의 성장은 학기 초, 방학, 시험 기간 등에 따라 단순히 직선으로 상승하지 않고 울퉁불퉁하게 변합니다.
  • 개인차: 학생마다 성장하는 속도와 패턴이 다릅니다.

이 데이터를 분석하기 위한 일반적인 SPME(Semiparametric Mixed-Effects) 모델 식은 다음과 같습니다:

yij=xijβ+η(tij)+zijui+vi(tij)+ϵijy_{ij} = x_{ij}’\beta + \eta(t_{ij}) + z_{ij}’u_i + v_i(t_{ij}) + \epsilon_{ij}

  • xijβx_{ij}’\beta: 성별처럼 딱 떨어지는 영향력(고정 효과).
  • η(tij)\eta(t_{ij}): 모든 학생의 평균적인 복잡한 성장 곡선(비모수 고정 효과).
  • zijuiz_{ij}’u_i: 학생 개개인의 기본 출발점 차이(모수 무선 효과).
  • vi(tij)v_i(t_{ij}): 학생 개개인의 고유하고 복잡한 변화 패턴(비모수 무선 효과).

3. 모의 데이터 생성 (R 코드)

실제 분석에 앞서, 학교 현장과 유사한 데이터를 R을 통해 생성해 보겠습니다.

R

# 필요한 라이브러리 로드
library(ggplot2)
library(mgcv) # Semiparametric 모델(GAM) 분석용

# 1. 모의 데이터 생성
set.seed(2026)
n_students <- 100
n_times <- 8
times <- seq(0, 7, length.out = n_times)

data_list <- list()

for(i in 1:n_students) {
  # 학생 특성
  gender <- sample(c(0, 1), 1) # 0: 여, 1: 남
  pre_score <- rnorm(1, 50, 10)
  
  # 시간에 따른 복잡한 변화 (Sine 함수를 활용한 비선형성)
  # 개인별 무선 효과(v_i) 포함
  random_slope <- rnorm(1, 2, 0.5)
  y_values <- 30 + 0.5 * pre_score + 2 * gender + 
              5 * sin(times/1.5) + random_slope * times + 
              rnorm(n_times, 0, 3)
  
  data_list[[i]] <- data.frame(
    id = i,
    time = times,
    gender = gender,
    pre_score = pre_score,
    score = y_values
  )
}

school_data <- do.call(rbind, data_list)
school_data$gender <- factor(school_data$gender, labels = c("Girl", "Boy"))

# 데이터 확인
head(school_data)

4. jamovi 및 R에서의 분석 방법

jamovi 사용법

jamovi에서 이와 같은 준모수 모델을 직접 구현하려면 ‘GAMLj’ 모듈이나 ‘GAM’ 관련 모듈을 설치해야 합니다.

  1. Analyses 탭에서 GAMLj (General Analysis for Linear Models) 선택.
  2. Mixed Models 선택.
  3. Dependent Variable에 score, Covariates에 time, pre_score, Factor에 gender 입력.
  4. Polynomials 옵션에서 time의 차수를 높이거나, 비모수적 접근을 원할 경우 R 패키지 연동을 활용합니다.

R을 이용한 정밀 분석 (Smoothing Spline 적용)

교재에서 강조하는 회귀 스프라인(Regression Spline) 또는 평활 스프라인(Smoothing Spline) 기법을 사용하여 분석해 보겠습니다.

R

# 2. 준모수 혼합 모델 적합 (gamm 함수 사용)
# score = gender + pre_score (매개변수) + s(time) (비모수 평활)
model_fit <- gamm(score ~ gender + pre_score + s(time, k=5), 
                  random = list(id = ~1 + time), 
                  data = school_data)

# 결과 요약
summary(model_fit$gam) # 고정 효과 확인

5. 결과 해석 및 시각화

분석 결과, 성별사전 점수는 유의미한 양의 영향을 미쳤으며, 시간(time)에 따른 변화는 단순한 직선이 아니라 물결치는 곡선의 형태를 띠었습니다.

전체 학생의 평균 성장 곡선 (Population Fit)

R

# 시각화 코드
ggplot(school_data, aes(x = time, y = score, group = id)) +
  geom_line(alpha = 0.1, color = "gray") + # 개별 학생 선
  geom_smooth(aes(group = 1), method = "gam", formula = y ~ s(x, k=5), 
              color = "blue", size = 1.5) + # 평균 곡선 (SPME)
  labs(title = "창의적 문제해결력 성장 곡선 (SPME 모델)",
       x = "측정 시점 (학기)", y = "창의성 점수") +
  theme_minimal()

주요 결과 요약 (모의 분석 결과 기반)

구분추정치(Estimate)유의성(P-value)의미
성별(남학생)1.42< .001남학생이 여학생보다 평균 1.42점 높음.
사전 점수0.54< .001사전 점수가 높을수록 현재 점수도 높음.
s(time)곡선 형태< .001시간이 지남에 따라 점수가 비선형적으로 상승함.

6. 결론 및 시사점

준모수 혼합 모델(SPME)은 학교 현장의 복잡한 데이터를 분석할 때 매우 강력한 도구입니다.

  1. 유연성: 교육 프로그램의 효과가 시기별로 다르게 나타나는 현상을 정확히 포착합니다.
  2. 정밀함: 학생 개인의 특성(무선 효과)을 고려하면서도 전체적인 추세를 놓치지 않습니다.
  3. 해석력: “성별 차이는 존재하지만, 성장 패턴 자체는 모든 학생이 비슷하게 파동을 그리며 상승한다”는 식의 깊이 있는 해석이 가능해집니다.

참고문헌 (APA Style)

  • Durban, M., Harezlak, J., Wand, M.P., & Carroll, R.J. (2005). Simple fitting of subject-specific curves for longitudinal data. Statistics in Medicine, 24(8), 1153-1167.
  • Eubank, R. L. (1999). Nonparametric regression and spline smoothing. New York: Marcel Dekker.
  • Green, P. J., & Silverman, B. W. (1994). Nonparametric regression and generalized linear models: A roughness penalty approach. London: Chapman and Hall.
  • Wu, H., & Zhang, J. T. (2006). Nonparametric regression methods for longitudinal data analysis. New York: Wiley.
  • Zhang, J. T. (2010). Smoothing and semiparametric models. In G. Marcoulides & J. Heck (Eds.), The SAGE Handbook of Multilevel Modeling (pp. 299-324). London: SAGE Publications.

Chap16. 범주형 반응 변수(Categorical Response Data)를 위한 다층 모형: 점수 말고 “선택”을 분석하자!

안녕하세요!

오늘은 범주형 반응 변수(Categorical Response Data)를 위한 다층 모형에 대해 살펴보겠습니다. “학교 현장의 데이터”를 예시로 들어 직관적인 설명과 수리적 엄밀함을 모두 갖춘 형태로 재구성해 드리겠습니다.

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


이전까지 우리는 ‘시험 점수’, ‘키’, ‘몸무게’처럼 연속된 숫자로 된 결과(연속형 변수)를 주로 다뤘습니다. 하지만 교육 현장에는 이런 데이터만 있는 게 아니죠?

  • 서열형(Ordinal): “학교 급식 만족도 (1: 매우 불만 ~ 5: 매우 만족)”
  • 명목형(Nominal): “방과 후 학교 선택 (1: 축구반, 2: 미술반, 3: 코딩반)”

이런 데이터는 일반적인 선형 다층모형(Linear Mixed Model)으로 분석하면 엉터리 결과가 나옵니다. 왜냐고요? ‘축구반(1)’과 ‘코딩반(3)’의 평균이 ‘미술반(2)’은 아니니까요!.

오늘 우리는 이 까다로운 녀석들을 어떻게 다층모형으로 요리할지 배울 것입니다.


1. 이론적 기초: 눈에 보이지 않는 “잠재 변수” (Latent Variable)

우리가 겉으로 보는 것은 “체크 표시”된 설문지(범주)이지만, 통계학자들은 그 이면에 연속적인 심리적 수치(잠재 변수, yy^*)가 숨어 있다고 가정합니다.

1) 서열형 데이터 (Ordinal): “문턱 넘기” 게임

학생의 ‘학교 만족도’를 예로 들어봅시다.

  • 학생 마음속에는 만족감(yy^*)이라는 연속적인 게이지가 있습니다.
  • 이 게이지가 특정 임계값(Threshold, δ\delta)을 넘을 때마다 응답이 바뀝니다.
    • 게이지 < 문턱1 \rightarrow “불만족”
    • 문턱1 < 게이지 < 문턱2 \rightarrow “보통”
    • 문턱2 < 게이지 \rightarrow “만족”

2) 명목형 데이터 (Nominal): “유틸리티” 게임

방과 후 학교 선택을 예로 들어봅시다.

  • 학생은 각 선택지(축구, 미술, 코딩)에 대해 효용(Utility,yijmy^*_{ijm})을 계산합니다.
  • 가장 효용이 높은 대안을 선택합니다 (yij=my_{ij}=m if yijmy^*_{ijm} is max).

2. 실습 시나리오: “행복한 학교 만들기” 프로젝트

우리는 가상의 데이터 세트 두 개를 만들어 분석해 보겠습니다.

상황 A (서열형): 학교 폭력 예방 교육 후 만족도 조사

  • 구조: 학생(ii)들이 학급(jj)에 소속됨 (2수준).
  • 종속변수(YY): 교육 만족도 (1: 별로다, 2: 보통이다, 3: 훌륭하다)
  • 독립변수(Level 1): 학생의 사전 지식(Knowledge)
  • 독립변수(Level 2): 담임 선생님의 열정(Passion)

상황 B (명목형): 자유학기제 동아리 선택

  • 구조: 학생(ii)들이 학교(jj)에 소속됨.
  • 종속변수(YY): 동아리 선택 (Science, Art, Sports)
  • 독립변수: 성별(Gender), 학교 예산(Budget)

3. R을 이용한 모의 데이터 생성

jamovi에서 불러올 수 있도록 R을 이용해 그럴듯한 데이터를 생성해 봅시다.

R

# [R Code] 데이터 생성 스크립트
set.seed(1234)

# 1. 서열형 데이터 (Ordinal) 생성
n_classes <- 50
n_students <- 30
N <- n_classes * n_students

# Level 2 (학급) 변수
class_id <- rep(1:n_classes, each = n_students)
teacher_passion <- rep(rnorm(n_classes, 0, 1), each = n_students) # 선생님 열정
u_j <- rep(rnorm(n_classes, 0, 0.5), each = n_students) # 학급 무선 효과

# Level 1 (학생) 변수
knowledge <- rnorm(N, 0, 1) # 사전 지식
e_ij <- rlogis(N) # 로지스틱 오차항 (Ordinal Logit 가정)

# 잠재 변수 (Latent Variable) 생성
# 수식: y* = 0.5*Knowledge + 0.7*Passion + u_j + e_ij
y_star <- 0.5 * knowledge + 0.7 * teacher_passion + u_j + e_ij

# 관측 변수 (Ordinal)로 변환 (Threshold: -1, 1)
satisfaction <- cut(y_star, breaks = c(-Inf, -1, 1, Inf), labels = c("Low", "Medium", "High"))

data_ordinal <- data.frame(class_id, student_id=1:N, satisfaction, knowledge, teacher_passion)

# 2. 명목형 데이터 (Nominal) 생성
# 선택지: 1(Science), 2(Art), 3(Sports)
# 기준(Reference): Science
n_schools <- 30
n_students_nom <- 40
N_nom <- n_schools * n_students_nom

school_id <- rep(1:n_schools, each = n_students_nom)
budget <- rep(rnorm(n_schools, 0, 1), each = n_students_nom) # 학교 예산
gender <- rbinom(N_nom, 1, 0.5) # 0: 여학생, 1: 남학생
u_school <- rep(rnorm(n_schools, 0, 0.5), each = n_students_nom)

# 각 대안에 대한 효용(Utility) 계산 (Multinomial Logit)
# Science(Ref): Utility = 0
# Art: Utility = -0.5 + 0.5*Gender + 0.3*Budget + u_school
# Sports: Utility = -1.0 + 1.2*Gender - 0.2*Budget + u_school

util_science <- 0
util_art <- -0.5 + 0.5 * gender + 0.3 * budget + u_school + rlogis(N_nom)
util_sports <- -1.0 + 1.2 * gender - 0.2 * budget + u_school + rlogis(N_nom)

choice_mat <- cbind(util_science, util_art, util_sports)
choice_idx <- max.col(choice_mat)
club_choice <- factor(choice_idx, labels = c("Science", "Art", "Sports"))

data_nominal <- data.frame(school_id, student_id=1:N_nom, club_choice, gender, budget)

# CSV 파일로 저장 (jamovi에서 열기 위함)
write.csv(data_ordinal, "chap16-1.csv", row.names = FALSE)
write.csv(data_nominal, "chap16-2.csv", row.names = FALSE)

4. 서열형 반응 변수 분석 (Ordinal Response Models)

이론: 누적 로짓 모형 (Cumulative Logit Model)

서열형 분석의 핵심은 “비례 오즈 가정(Proportional Odds Assumption)”입니다. 이것은 XX라는 변수가 ‘불만족 \rightarrow 보통’으로 갈 확률을 높여주는 정도와, ‘보통 \rightarrow 만족’으로 갈 확률을 높여주는 정도가 같다고 가정하는 것입니다.

수식은 다음과 같습니다:

logP(yijm)1P(yijm)=αm(xijβ+zijuj)\log \frac{P(y_{ij} \le m)}{1 – P(y_{ij} \le m)} = \alpha_m – (x_{ij}’\beta + z_{ij}’u_j)

(주의: 소프트웨어마다 부호가 다를 수 있습니다. αm\alpha_m은 절편이자 문턱값입니다.)

jamovi 실습 가이드 (GAMLj 모듈 활용)

jamovi 기본 메뉴에는 서열형 다층모형이 없습니다. GAMLj 모듈을 설치해야 합니다.

  1. 데이터 열기: chap16-1.csv
  2. 모듈 실행: Analyses > GAMLj > Generalized Mixed Models
  3. 변수 설정:
    • Dependent Variable: satisfaction
    • Random Effects (Cluster): class_id
    • Covariates: knowledge, teacher_passion
  4. 모델 설정 (Model Options):
    • Family: Ordinal (Link function: logit) 설정. 이것이 핵심입니다!
  5. Random Effects 설정:
    • class_id를 오른쪽 Random Coefficients 박스로 이동 (Intercept만 넣어서 Random Intercept Model 구성).

분석 결과 해석 (예시)

EffectEstimateSEzp
Fixed Effects
knowledge0.520.0510.4<.001
teacher_passion0.680.088.5<.001
Thresholds
LowMedium-1.05
MediumHigh0.98
  • 해석:
    • Knowledge (0.52): 학생의 사전 지식이 1단위 높을수록 더 높은 만족도 범주(예: 보통 이상, 만족)에 속할 로그 오즈(log-odds)가 0.52만큼 증가합니다. 즉, 지식이 많을수록 수업에 만족할 확률이 높습니다.
    • Teacher Passion (0.68): 선생님의 열정이 높을수록 학생들의 만족도가 통계적으로 유의하게 높아집니다.
    • Random Effect: 학급 간 만족도의 차이(분산)를 보여줍니다.

5. 명목형 반응 변수 분석 (Nominal Response Models)

이론: 다항 로짓 모형 (Multinomial Logit Model)

여기서는 순서가 없습니다. 특정 카테고리(Reference Category)를 기준으로 다른 카테고리를 선택할 확률을 비교합니다.

πijm=exp(xijβm+zijujm)r=1Mexp(xijβr+zijujr)\pi_{ijm} = \frac{\exp(x_{ij}’\beta_m + z_{ij}’u_{jm})}{\sum_{r=1}^M \exp(x_{ij}’\beta_r + z_{ij}’u_{jr})}

여기서 중요한 점은 회귀계수 β\beta가 카테고리(mm)마다 다르다는 것입니다. 즉, 남학생(xx)이 미술반을 선택할 확률과 축구반을 선택할 확률에 미치는 영향력은 다릅니다.

jamovi 실습 가이드

아쉽게도 jamovi의 GAMLj조차 현재(2024년 기준) 명목형 다층모형(Multinomial Mixed Model)은 지원이 제한적일 수 있습니다. 이 경우 R의 mclogit 패키지를 사용하는 것이 정석입니다. 여기서는 R 코드로 분석하는 법을 보여드립니다.

R

# [R Code] 명목형 다층모형 분석
library(mclogit)

# 데이터 불러오기 (위에서 생성한 데이터)
df_nom <- read.csv("chap16-2.csv")
df_nom$club_choice <- as.factor(df_nom$club_choice)
df_nom$school_id <- as.factor(df_nom$school_id)

# 모델 적합 (Baseline: Science)
# mblogit: Baseline-category logit model for multiphase data
fit_nom <- mblogit(club_choice ~ gender + budget, 
                   random = ~1|school_id, 
                   data = df_nom)

summary(fit_nom)

분석 결과 해석 (예시)

결과는 “Science(과학반)” 대비 다른 반을 선택할 확률로 나옵니다.

  1. Response: Art vs. Science
    • Gender (Est: 0.48): 남학생일수록(1), 과학반보다 미술반을 선택할 확률이 높음 (가상 데이터 설정에 따름).
    • Budget (Est: 0.32): 예산이 많을수록 과학반보다 미술반 선택 확률 증가.
  2. Response: Sports vs. Science
    • Gender (Est: 1.25): 남학생일수록 과학반 대비 스포츠반을 선택할 확률이 매우 높음.
    • Budget (Est: -0.21): 예산이 많을수록 스포츠반보다는 과학반을 선호하는 경향이 생김.

6. 심화: 전문가를 위한 팁 (Advanced Topics)

1) 척도 이질성 (Scale Heterogeneity)

어떤 학생들은 좋고 싫음이 분명해서 1점 아니면 5점을 줍니다(분산이 큼). 어떤 학생들은 우유부단해서 3점 근처만 맴돕니다(분산이 작음). 이처럼 응답의 분산(Scale) 자체가 사람마다 다를 수 있다는 것을 모형에 포함시키는 것을 ‘위치-척도 모형(Location-Scale Model)’이라고 합니다.

  • 수식: σϵijexp(wijγ)\sigma_{\epsilon_{ij}} \propto \exp(-w_{ij}’\gamma).
  • 이것을 고려하면 그룹 간 차이를 더 명확하게 설명할 수 있습니다.

2) 잠재 집단 (Latent Classes)

모든 학생이 같은 회귀 계수를 가진다고 가정하는 대신, 데이터 안에 숨겨진 서로 다른 집단(Latent Class)이 있다고 가정할 수 있습니다.

  • 예: “운동 매니아 그룹”은 성별이 동아리 선택에 큰 영향을 미치지만, “학구파 그룹”은 성별 영향이 거의 없을 수 있습니다.
  • 이런 접근은 혼합 모형(Mixture Model)을 통해 분석 가능합니다.

요약 및 결론

  1. 교육 현장의 데이터는 연속형보다 범주형(만족도, 선택 등)이 많습니다.
  2. 서열형 데이터는 “누적 로짓 모형”을 사용하며, 문턱(Threshold)을 넘는 개념으로 이해하면 쉽습니다.
  3. 명목형 데이터는 “다항 로짓 모형”을 사용하며, 기준 카테고리 대비 다른 대안을 선택할 효용(Utility)을 비교합니다.
  4. 분석 도구로 jamovi (GAMLj)R (mclogit)을 활용하면, 학생과 학교의 다층 구조를 반영한 정확한 분석이 가능합니다.

오늘 배운 내용을 통해 단순한 점수 비교를 넘어, 학생들의 “선택”과 “마음”을 더 깊이 이해하는 연구자가 되시길 바랍니다!


참고문헌 (APA Style)

  • Agresti, A. (2002). Categorical data analysis (2nd ed.). Hoboken, NJ: John Wiley & Sons.
  • Agresti, A., Booth, J. G., Hobert, J. P., & Caffo, B. (2000). Random-effects modeling of categorical response data. Sociological Methodology, 30, 27-80.
  • Hedeker, D., & Gibbons, R. D. (1994). A random-effects ordinal regression model for multilevel analysis. Biometrics, 50, 933-944.
  • Hedeker, D., Berbaum, M., & Mermelstein, R. (2006). Location-scale models for multilevel ordinal data: Between-and within-subjects variance modeling. Journal of Probability and Statistical Science, 4(1), 1-20.
  • Vermunt, J. K. (2013). Categorical response data. In M. A. Scott, J. S. Simonoff, & B. D. Marx (Eds.), The SAGE handbook of multilevel modeling (pp. 287-308). Thousand Oaks, CA: SAGE Publications.

Chap13. 다층 함수형 데이터 분석(Multilevel Functional Data Analysis, MFDA)

안녕하세요!

오늘은 다층 함수형 데이터 분석(Multilevel Functional Data Analysis, MFDA)에 대해 살펴보겠습니다.

이름만 들어도 머리가 지끈거리시나요? 걱정 마세요. “학교 현장의 데이터”를 예시로 들어, 아주 직관적인 설명부터 대학원 수준의 수리적 엄밀함까지 갖춘 형태로 완벽하게 재구성해 드리겠습니다.

분석 도구와 관련하여 한 가지 중요한 점을 말씀드립니다. jamovi는 훌륭한 통계 도구이지만, 오늘 다룰 ‘함수형 데이터(Functional Data)’와 같은 고차원적 분석은 아직 기본 메뉴로 지원하지 않습니다(일반적인 다층모형까지만 지원). 따라서 본 강의에서는 R을 사용하여 모의 데이터 생성부터 분석, 시각화까지 구현해 드리겠습니다. 대신, 코드는 복사해서 바로 쓸 수 있도록 아주 친절하게 주석을 달아 제공합니다.

그럼, 시작해 볼까요?


1. 서론: 점(Point)이 아니라 선(Line)을 보자!

1.1 직관적 이해 (초등학생도 이해하는 MFDA)

여러분이 담임 선생님이라고 상상해 보세요. 우리 반 학생 철수의 수학 실력을 알고 싶습니다. 보통 우리는 “기말고사 점수 80점”이라는 하나의 숫자(Scalar)로 철수를 평가합니다. 이것이 전통적인 통계입니다.

그런데, 80점이라는 점수가 철수의 모든 것을 설명할까요?

어떤 학생은 수업 초반 20분은 엄청 집중하다가 후반에 자는 학생이 있고, 어떤 학생은 처음엔 멍하다가 나중에 발동이 걸리는 학생이 있습니다.

  • 스칼라(Scalar) 데이터: 철수의 점수 = 80점 (사진 한 장 📷)
  • 함수형(Functional) 데이터: 철수의 45분간 집중력 변화 그래프 (동영상 🎥)

다층 함수형 데이터 분석(MFDA)은 바로 이 “동영상(변화 곡선)”을 분석하는 방법입니다. 그런데 왜 ‘다층(Multilevel)’일까요?

학생마다 수업을 듣는 집중력 곡선이 다르고(1수준), 이 학생들이 속한 반(Class)마다 분위기가 다르기(2수준) 때문입니다.

핵심 요약:

  • FDA: 시간에 따라 변하는 “곡선” 자체를 하나의 데이터로 분석함.
  • Multilevel: 그 곡선들이 여러 집단(학생, 학교 등)에 층층이 겹쳐 있는 구조를 반영함.

2. 시나리오: “어떤 수업 방식이 끝까지 집중하게 할까?”

우리는 다음과 같은 교육 실험을 진행한다고 가정해 보겠습니다.

  • 연구 대상: 초등학교 6학년 학생 100명 (2개 반)
  • 독립 변수(집단):
    • A반 (50명): 전통적 강의식 수업
    • B반 (50명): 게임 기반(Gamification) 수업
  • 종속 변수(함수): 45분 수업 동안의 분당 뇌파 집중도(Attention Score, 0~100)
    • 학생 1명당 45개의 측정값(0분~44분)이 연결된 하나의 곡선이 생성됩니다.

3. R을 이용한 모의 데이터 생성 및 분석

자, 이제 이 시나리오를 실제 데이터로 만들어 분석해 보겠습니다. R 스튜디오를 켜고 아래 코드를 순서대로 실행해 보세요.

3.1 데이터 생성 (Data Simulation)

실제 학교 현장과 유사하게 데이터를 생성합니다.

  • A반(강의식): 초반엔 집중력이 높지만 시간이 갈수록 급격히 떨어지는 패턴.
  • B반(게임형): 초반엔 산만하지만, 시간이 갈수록 몰입도가 유지되거나 상승하는 패턴.

R

# 필수 패키지 설치 및 로드
if(!require(ggplot2)) install.packages("ggplot2")
if(!require(refund)) install.packages("refund") # 함수형 데이터 분석용 핵심 패키지
if(!require(dplyr)) install.packages("dplyr")
if(!require(tidyr)) install.packages("tidyr")

library(ggplot2)
library(refund)
library(dplyr)
library(tidyr)

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

# 1. 기본 설정
n_students <- 100 # 총 학생 수
n_points <- 45    # 측정 지점 (0분~44분, 총 45분 수업)
time_points <- seq(0, 1, length.out = n_points) # 0에서 1로 정규화된 시간

# 2. 그룹 할당 (A반: 0=강의식, B반: 1=게임형)
group <- c(rep(0, n_students/2), rep(1, n_students/2))
group_factor <- factor(group, labels = c("Lecture", "Gamification"))

# 3. 기저 함수 생성 (집중력의 기본 패턴)
# 강의식(A): 초반 높음 -> 후반 낮음 (감소형)
mu_A <- function(t) { 80 - 40 * t } 
# 게임형(B): 초반 중간 -> 후반 유지/상승 (유지형)
mu_B <- function(t) { 60 + 10 * sin(2 * pi * t) + 10 * t }

# 4. 개인차(Random Effect) 추가
# 학생마다 고유한 집중력 기복이 있음 (브라운 운동 + 임의 오차)
Y_matrix <- matrix(NA, nrow = n_students, ncol = n_points)

for(i in 1:n_students){
  # 기본 곡선 선택
  base_curve <- if(group[i] == 0) mu_A(time_points) else mu_B(time_points)
  
  # 학생별 임의 효과 (Random Intercept + Random Slope 유사 효과)
  student_deviation <- rnorm(1, 0, 5) + rnorm(1, 0, 2) * time_points * 10
  
  # 측정 오차 (Measurement Error)
  noise <- rnorm(n_points, 0, 2)
  
  Y_matrix[i, ] <- base_curve + student_deviation + noise
}

# 데이터 프레임 형태로 정리 (시각화용)
df_long <- data.frame(
  ID = rep(1:n_students, each = n_points),
  Time = rep(1:n_points, n_students),
  Group = rep(group_factor, each = n_points),
  Attention = as.vector(t(Y_matrix))
)

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

3.2 데이터 시각화 (Spaghetti Plot)

분석의 첫 걸음은 데이터를 눈으로 확인하는 것입니다. 개별 학생들의 집중력 곡선을 모두 그려보겠습니다. 이를 ‘스파게티 플롯(Spaghetti Plot)’이라고 합니다.

R

# 스파게티 플롯 + 평균 곡선
ggplot(df_long, aes(x = Time, y = Attention, group = ID, color = Group)) +
  geom_line(alpha = 0.3) + # 개별 학생 곡선 (흐리게)
  geom_smooth(aes(group = Group), method = "gam", se = TRUE, size = 2) + # 그룹별 평균 곡선 (진하게)
  labs(title = "수업 방식에 따른 학생 집중력 변화 곡선",
       subtitle = "개별 곡선(흐림)과 그룹 평균 곡선(진함)",
       x = "수업 시간 (분)", y = "집중도 (점수)") +
  theme_minimal() +
  scale_color_manual(values = c("Lecture" = "blue", "Gamification" = "red"))

해석: > 파란색(강의식) 곡선들은 시간이 갈수록 우하향하는 경향이 뚜렷합니다. 반면 빨간색(게임형)은 출렁거림이 있지만 전반적으로 수준이 유지됩니다. 이것이 바로 함수형 데이터 분석의 매력입니다. 평균 점수만 비교했다면 놓쳤을 ‘시간에 따른 패턴’을 볼 수 있죠.


4. 이론적 배경: 함수형 혼합 모형 (Functional Mixed Model)

이 데이터를 통계적으로 검증하려면 어떤 모형을 써야 할까요? 함수형 혼합 모형(FMM) 수식은 다음과 같습니다.

Wi(s)=a=1pXiaBa(s)+b=1qZibUb(s)+Ei(s)W_i(s) = \sum_{a=1}^{p} X_{ia} B_a(s) + \sum_{b=1}^{q} Z_{ib} U_b(s) + E_i(s)

초등학생 눈높이로 번역해 보겠습니다.

  • Wi(s)W_i(s): 관찰된 학생의 집중력 곡선 (우리가 본 스파게티 면발)
  • XiaBa(s)X_{ia} B_a(s): 고정 효과 (Fixed Effects) = 수업 방식에 따른 공통적인 패턴.
    • “강의식 수업은 대체로 이렇게 흘러간다”는 평균 곡선.
  • Ub(s)U_b(s): 임의 효과 (Random Effects) = 학생 개인의 특성.
    • “철수는 원래 초반에 약해”, “영희는 후반에 강해” 같은 개인별 편차 곡선.
  • Ei(s)E_i(s): 오차 (Error) = 측정 당시의 잡음.

우리의 목표는 Ba(s)B_a(s)(수업 방식의 효과 곡선)가 통계적으로 유의미하게 다른지 확인하는 것입니다.


5. 핵심 분석: Function-on-Scalar Regression

이제 refund 패키지를 사용하여 “수업 방식(Scalar)이 집중력 곡선(Function)에 미치는 영향”을 분석해 보겠습니다.

5.1 분석 코드 실행

R



# 1. 데이터 프레임 만들기
# 핵심: I() 함수를 사용하여 Y_matrix를 데이터 프레임의 한 '열'로 안전하게 포장합니다.
covariates_df <- data.frame(group_factor = group_factor)
covariates_df$Y_matrix <- I(Y_matrix) 

# 2. 모델 적합
# 이제 Y_matrix와 group_factor 모두 covariates_df 안에 들어있습니다.
f_model <- fosr(Y_matrix ~ group_factor, data = covariates_df, method = "OLS") 

# 3. 결과 확인
summary(f_model)

# 4. 시각화 (베타 계수 곡선)
plot(f_model, split = 1)

5.2 결과 해석 및 재구성 (교육적 설명)

분석을 실행하면 베타 계수 곡선(β(t)\beta(t))이 그려집니다.

  1. 절편 곡선 (Intercept): Reference Group강의식 수업(Lecture)의 평균 집중력 곡선을 보여줍니다. (우하향하는 그래프가 나옵니다.)
  2. 기울기 곡선 (Slope for Gamification): 여기가 핵심입니다! 이 곡선은 “강의식 대비 게임형 수업의 집중력 차이”를 시간대별로 보여줍니다.
    • 0~10분: 차이가 거의 0에 가깝거나 음수일 수 있습니다 (초반엔 게임 룰 익히느라 산만함).
    • 20분 이후: 곡선이 양수(+)로 치솟습니다. 즉, 수업 중반 이후부터 게임형 수업의 집중력이 강의식보다 유의미하게 높다는 뜻입니다.

6. 심화 분석: FPCA (함수형 주성분 분석)

첨부 파일에서 강조하는 또 하나의 중요한 기법은 FPCA입니다.

학생 100명의 곡선은 제각각이지만, 그 안에는 몇 가지 “대표적인 유형(Mode of Variation)”이 숨어 있습니다.

R

# FPCA 수행
fpca_res <- fpca.sc(Y_matrix)

# 주성분(Eigenfunctions) 시각화
par(mfrow = c(1, 2))
plot(fpca_res$efunctions[,1], type='l', main = "제1주성분 (전체적 높낮이)", ylab="Value", xlab="Time")
plot(fpca_res$efunctions[,2], type='l', main = "제2주성분 (초반 vs 후반)", ylab="Value", xlab="Time")
  • 제1주성분: “집중력이 전체적으로 높은 학생 vs 낮은 학생”을 구분하는 패턴.
  • 제2주성분: “초반에 높고 후반에 지치는 유형(토끼형)” vs “초반엔 낮지만 후반에 강한 유형(거북이형)”을 구분하는 패턴.

교수자는 이 점수(FPC Score)를 활용해 학생들을 유형별로 분류하고 상담할 수 있습니다.


7. 결론 및 제언

오늘 우리는 다층 함수형 데이터 분석(MFDA)을 통해, 단순히 “평균 80점”이 아닌 “시간의 흐름에 따른 학습 패턴”을 분석했습니다.

  1. 교육적 시사점: 단순한 점수 비교를 넘어, “언제(When)” 개입해야 하는지 알려줍니다.
    • 분석 결과, 강의식 수업은 20분 지점에서 집중력이 급락하므로, 이때 환기 활동(Brain Break)이 필요하다는 데이터 기반 의사결정이 가능합니다.
  2. 방법론적 의의:
    • jamovi로 구현되지 않는 고차원 분석을 R의 refund 패키지로 구현함으로써, 연구의 깊이를 더했습니다.
    • 다층 모형의 개념을 함수(곡선)로 확장하여, 개인차(Between-subject)시간차(Within-subject)를 동시에 고려했습니다.

이 분석 기법을 활용하여, 여러분의 교실과 연구실에서 숨겨진 데이터의 맥박을 느껴보시길 바랍니다.


참고문헌 (References)

  • Crainiceanu, C. M., Caffo, B. S., & Morris, J. S. (2013). Multilevel functional data analysis. In The SAGE Handbook of Multilevel Modeling (pp. 223-238). SAGE Publications Ltd.
  • Morris, J. S., & Carroll, R. J. (2006). Wavelet-based functional mixed models. Journal of the Royal Statistical Society: Series B (Statistical Methodology), 68(2), 179-199.
  • Ramsay, J. O., & Silverman, B. W. (2005). Functional Data Analysis (2nd ed.). Springer.
  • Crainiceanu, C. M., et al. (2009). Generalized multilevel functional regression. Journal of the American Statistical Association, 104(488), 1550-1561.

Chap15. 일반화 선형 혼합 모형

안녕하세요!

오늘은 일반화 선형 혼합 모형(Generalized Linear Mixed Models, GLMM)에 대해 살펴보겠습니다. “학교 현장의 데이터”를 예시로 들어 직관적인 설명과 수리적 엄밀함을 모두 갖춘 형태로 재구성해 드리겠습니다.

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


🎓 일반화 선형 혼합 모형(GLMM): 학교 밖으로 나간 수학 점수

반갑습니다. 오늘은 통계의 꽃이라 불리는 일반화 선형 혼합 모형(GLMM)에 대해 이야기해 보려 합니다. 이름만 들어도 머리가 지끈거릴 수 있지만, 사실 원리는 우리가 학교에서 겪는 일상과 매우 비슷합니다.

1. 왜 GLMM인가? (초등학생도 이해하는 비유)

여러분이 “우리 반 친구들이 이번 수학 시험에 합격할지 불합격할지(성공/실패)“를 맞히고 싶다고 상상해 봅시다.

  1. 일반적인 회귀분석(Linear Regression)은 “점수”처럼 연속적인 숫자를 예측할 때 씁니다. 하지만 우리는 “합격(1) 아니면 불합격(0)”인 두 가지 결과만 있습니다. 자로 키를 재는 게 아니라, 동전의 앞면/뒷면을 맞히는 것과 같죠. 이때 필요한 것이 ‘일반화(Generalized)’ 선형 모형입니다.
  2. 그런데 학생들은 ‘반(Class)’이라는 바구니에 담겨 있습니다. 어떤 반 선생님은 설명을 아주 잘해서 그 반 아이들이 전체적으로 합격률이 높을 수 있죠. 학생들끼리 서로 영향을 주고받는다는 뜻입니다. 이걸 무시하고 분석하면 엉터리 결과가 나옵니다. 그래서 ‘혼합(Mixed)’ 모형, 즉 다층 모형이 필요합니다.

이 둘을 합친 것이 바로 GLMM입니다. “결과가 0/1이거나 횟수(Count)이면서, 데이터가 집단(학교, 반)으로 묶여 있을 때” 사용하는 가장 강력한 도구죠.


2. 가상의 시나리오: “수학 챌린지 성공 예측 대작전”

우리의 목표는 다음과 같습니다.

  • 연구 문제: 학생의 자기효능감교사의 교수법(혁신적 vs 전통적)이 수학 챌린지 성공 여부(성공/실패)에 미치는 영향은 무엇인가?
  • 데이터 구조:
    • 1수준(학생): 자기효능감(Self_Efficacy), 챌린지 성공여부(Success, 0=실패, 1=성공)
    • 2수준(학급): 교수법(Teaching_Method, 0=전통적, 1=혁신적)
  • 특이사항: 성공 여부는 0과 1로 된 이분형 변수(Binary outcome)입니다.

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

먼저, 실제 학교 현장과 유사한 데이터를 만들어 보겠습니다.

R

# 필요한 라이브러리 로드
library(lme4)
library(ggplot2)
library(dplyr)
library(sjPlot) # 시각화용

# 1. 데이터 생성 설정
set.seed(2025) # 재현성을 위한 시드 설정
n_classes <- 50      # 학급 수 (2수준)
n_students <- 30     # 학급당 학생 수 (1수준)
N <- n_classes * n_students

# 2. 2수준(학급) 변수 생성
class_id <- rep(1:n_classes, each = n_students)
# 교수법: 절반은 전통적(0), 절반은 혁신적(1)
teaching_method <- rep(rep(c(0, 1), each = n_classes/2), each = n_students)
# 학급 효과 (Random Intercept): 학급마다 기본 성공률이 다름
class_effect <- rep(rnorm(n_classes, mean = 0, sd = 1.5), each = n_students)

# 3. 1수준(학생) 변수 생성
# 자기효능감: 평균 0, 표준편차 1인 정규분포
self_efficacy <- rnorm(N, mean = 0, sd = 1)

# 4. 성공 확률 계산 (로지스틱 모형)
# Logit(p) = b0 + b1*효능감 + b2*교수법 + 학급효과 + 오차
beta_0 <- -1.0  # 절편 (기본적으로 성공이 조금 더 어려움)
beta_1 <- 1.2   # 자기효능감의 효과 (클수록 성공 확률 높음)
beta_2 <- 0.8   # 혁신적 교수법의 효과 (있으면 성공 확률 높음)

logit_p <- beta_0 + beta_1 * self_efficacy + beta_2 * teaching_method + class_effect
prob <- 1 / (1 + exp(-logit_p)) # 확률로 변환

# 5. 결과 변수 생성 (0 또는 1)
success <- rbinom(N, size = 1, prob = prob)

# 데이터프레임 생성
data <- data.frame(
  ClassID = factor(class_id),
  StudentID = 1:N,
  Success = factor(success, levels = c(0, 1), labels = c("Fail", "Pass")),
  SelfEfficacy = self_efficacy,
  TeachingMethod = factor(teaching_method, levels = c(0, 1), labels = c("Traditional", "Innovative"))
)

# 데이터 확인
head(data)

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

📊 데이터 탐색적 시각화

데이터가 잘 만들어졌는지 그래프로 확인해 봅시다.

R

# 자기효능감에 따른 성공 여부 (교수법별 차이)
ggplot(data, aes(x = SelfEfficacy, y = as.numeric(Success)-1, color = TeachingMethod)) +
  geom_point(alpha = 0.1) +
  geom_smooth(method = "glm", method.args = list(family = "binomial"), se = TRUE) +
  labs(title = "자기효능감과 수학 챌린지 성공 확률",
       x = "자기효능감 (표준화)",
       y = "성공 확률 (Probability of Pass)") +
  theme_minimal()

4. 분석 방법: jamovi와 R

🛠️ jamovi 분석 절차

jamovi는 기본 메뉴만으로는 GLMM(특히 이분형 종속변수)을 분석하기 어렵습니다. 따라서 jamovi 라이브러리에서 GAMLj (General Analyses for Linear Models in jamovi) 모듈을 설치하여 사용해야 합니다.

  1. 데이터 불러오기: 위에서 생성한 데이터를 csv로 저장 후 jamovi에서 엽니다.
  2. 모듈 선택: 상단 메뉴 Analyses -> GAMLj -> Generalized Mixed Models 선택.
  3. 변수 설정:
    • Dependent Variable: Success (성공 여부)
    • Covariates: SelfEfficacy (자기효능감)
    • Factors: TeachingMethod (교수법)
    • Cluster (Random Effect grouping): ClassID
  4. Random Effects 설정:
    • ClassID를 Random Intercept로 설정 (체크박스 선택).
  5. Family 설정:
    • Distribution: Binomial (이항분포)
    • Link function: Logit (로짓)
  6. 결과 해석: Fixed Effects의 p-value와 Random Effects의 분산(Variance)을 확인합니다.

💻 R 분석 코드 (GLMM 적합)

jamovi 내부에서 돌아가는 엔진과 동일한 lme4 패키지를 사용한 분석 코드입니다.

R

# GLMM 모델 적합
# family = binomial : 종속변수가 이분형(0/1)일 때 사용
glmm_model <- glmer(Success ~ SelfEfficacy + TeachingMethod + (1 | ClassID), 
                    data = data, 
                    family = binomial(link = "logit"))

# 결과 요약
summary(glmm_model)

5. 교재 내용의 심층 재구성: 이론과 해석

이제 첨부된 교재(SAGE Handbook)의 핵심 내용을 바탕으로 우리가 수행한 분석을 이론적으로 파헤쳐 보겠습니다.

5.1. 모형의 구조 (Specification)

GLMM은 4단계로 구성됩니다.

  1. 분포 선택: 종속변수 YY가 어떤 모양인지 결정합니다. 우리 예시는 성공/실패이므로 베르누이(Bernoulli) 분포를 따릅니다.
  2. 예측변수 선택: 자기효능감(x1x_1)과 교수법(x2x_2)을 포함합니다. 선형 예측식은 η=β0+β1x1+β2x2\eta = \beta_0 + \beta_1 x_1 + \beta_2 x_2가 됩니다.
  3. 연결 함수(Link Function): 0과 1 사이의 확률(pp)을 실수 전체 범위(-\infty \sim \infty)인 선형 예측식(η\eta)과 연결해야 합니다. 이때 로짓(Logit) 함수를 사용합니다. log(p1p)=β0+β1x1+β2x2+uclass\log\left(\frac{p}{1-p}\right) = \beta_0 + \beta_1 x_1 + \beta_2 x_2 + u_{class}
  4. 임의 효과(Random Effects): 학급(uclassu_{class})을 임의 효과로 선언합니다. 이는 학급마다 성공률의 ‘출발선’이 다름을 의미하며, 정규분포 N(0,σu2)N(0, \sigma_u^2)를 따른다고 가정합니다.

5.2. 추정 방법: 산을 오르는 법 (Fitting Methods)

GLMM에서 가장 어려운 점은 “계산”입니다. 일반적인 회귀분석처럼 공식을 딱 대입해서 답이 나오지 않습니다. 우도(Likelihood)라는 산의 정상을 찾아야 하는데, 그 과정에 안개(적분)가 껴 있습니다.

  • 최우추정법(Maximum Likelihood, ML): 가장 이상적이지만, 계산이 매우 복잡하여 적분(Integration)이 필요합니다.
  • 근사법: 이 적분을 해결하기 위해 가우스-헤르미트 구적법(Gauss-Hermite Quadrature)이나 라플라스 근사(Laplace Approximation)를 사용합니다. 이는 복잡한 곡면을 단순한 도형으로 근사시켜 넓이를 구하는 방식과 비슷합니다. jamovi와 R의 glmer는 기본적으로 이 방법을 사용합니다.
  • 경고: PQL(Penalized Quasi-Likelihood)이라는 방법도 있지만, 이분형 데이터(0/1)에서는 편향(Bias)이 심해 잘 쓰지 않습니다.

5.3. 조건부 vs 주변부 해석 (Conditional vs Marginal)

이 부분이 아주 중요합니다. GLMM의 결과(β\beta)는 “특정 학급에 속한 학생 개인”에 대한 효과입니다(Conditional).

  • GLMM (조건부): “내가 이 반에 계속 있으면서 자기효능감이 1 오르면, 나의 성공 오즈(Odds)는 얼마나 오르는가?”
  • GEE (주변부): “전체 학생 평균적으로 봤을 때, 자기효능감이 높은 집단은 낮은 집단보다 성공률이 얼마나 높은가?”

일반적으로 로지스틱 모형에서는 조건부 효과가 주변부 효과보다 값이 더 크게(0에서 멀어지게) 추정됩니다. 우리는 학생 개인의 변화와 학급 효과에 관심이 있으므로 GLMM이 적합합니다.


6. 분석 결과 해석 및 보고 (APA 스타일)

R/jamovi 분석 결과를 바탕으로 보고서를 작성하는 예시입니다.

📝 결과 요약

“수학 챌린지 성공 여부에 대한 일반화 선형 혼합 모형(GLMM) 분석 결과, 자기효능감과 교수법은 통계적으로 유의한 영향을 미치는 것으로 나타났다.”

  • 고정 효과(Fixed Effects):
    • 자기효능감의 회귀계수는 1.370 (SE=0.091,p<.001SE=0.091, p<.001)로, 학생의 자기효능감이 높을수록 성공 확률이 유의하게 증가하였다. 오즈비(Odds Ratio)로 환산하면 exp(1.37)3.93exp(1.37) \approx 3.93로, 자기효능감이 1단위 증가할 때 성공 오즈는 약 3.9배 증가한다.
    • 교수법(혁신적)의 효과는 0.227 (SE=0.529,p=.667SE=0.529, p=.667)로, 혁신적 교수법을 사용하는 학급의 학생이 전통적 교수법 학급 학생보다 성공할 확률은 통계적으로 유의미하지 않았다.
  • 임의 효과(Random Effects):
    • 학급 수준의 분산(Variance)은 3.221로 나타났다. 이는 학생의 개인 특성을 통제하고도 학급 간 성공률의 차이가 상당히 큼을 의미한다. 이를 통해 급내상관계수(ICC)를 계산하여 학급의 영향력을 보고할 수 있다.

📈 R을 이용한 결과 시각화 코드

R

# sjPlot 패키지를 이용한 깔끔한 시각화
library(sjPlot)
library(sjmisc)

# 1. 오즈비(Odds Ratio) 그래프 (Forest Plot)
plot_model(glmm_model, 
           type = "est", 
           transform = "exp", # 로그오즈를 오즈비로 변환
           show.values = TRUE, 
           title = "수학 챌린지 성공 요인 (Odds Ratios)",
           vline.color = "gray")

# 2. 예측 확률 그래프 (상호작용 효과 등 확인용)
plot_model(glmm_model, 
           type = "pred", 
           terms = c("SelfEfficacy", "TeachingMethod"),
           title = "자기효능감과 교수법에 따른 성공 예측 확률")

7. 마무리 및 진단

GLMM을 사용할 때는 모델이 데이터에 잘 맞는지 꼭 확인해야 합니다. 특히 “과산포(Overdispersion)” 문제를 조심해야 합니다. 이는 우리가 가정한 것보다 데이터가 더 넓게 퍼져 있는 현상을 말하는데, 이항 분포나 포아송 분포 분석 시 자주 발생합니다. 만약 과산포가 의심되면 분포를 바꾸거나 관측 단위의 임의 효과를 추가하는 등의 조치가 필요합니다.

오늘 우리는 난이도가 높은 GLMM을 학교 현장 예시를 통해 정복해 보았습니다. 이 분석을 통해 여러분은 단순히 “누가 공부 잘하나”를 넘어, “어떤 환경(학급)에서 어떤 특성(효능감)이 성공을 이끄는지” 입체적으로 볼 수 있는 눈을 갖게 되었습니다.

참고문헌 (APA Style)

McCulloch, C. E., & Neuhaus, J. M. (2013). Generalized linear mixed models: Estimation and inference. In M. A. Scott, J. S. Simonoff, & B. D. Marx (Eds.), The SAGE handbook of multilevel modeling (pp. 271-286). SAGE Publications.

Chap14. 비선형 혼합모델

안녕하세요!

오늘은 비선형 혼합 효과 모델(Nonlinear Mixed Effects Models)에 대해 살펴보겠습니다. “학교 현장의 데이터”를 예시로 들어 직관적인 설명과 수리적 엄밀함을 모두 갖춘 형태로 재구성해 드리겠습니다.

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


1. 왜 ‘비선형(Nonlinear)’ 모델이 필요할까요?

우리가 흔히 쓰는 ‘선형 모델(Linear Models)’은 계산이 쉽고 해석이 명확하다는 장점이 있습니다. 하지만 학교 현장에서 일어나는 실제 현상들은 직선으로 설명되지 않는 경우가 훨씬 많습니다.

  • 선형 모델(Empirical Models): 관찰된 데이터를 잘 설명하기 위해 직선을 긋는 ‘경험적 모델’입니다. 데이터 범위 안에서는 잘 맞을지 몰라도, 범위를 벗어나면 예측력이 떨어질 수 있습니다.
  • 비선형 모델(Mechanistic/Scientific Models): 데이터가 왜 그렇게 만들어졌는지, 즉 ‘생성 메커니즘’에 집중하는 ‘과학적 모델’입니다. 성장 곡선이나 학습 속도처럼 자연스러운 변화 과정을 반영하므로, 데이터가 없는 구간에 대해서도 더 정확한 예측이 가능합니다.

비선형 모델의 장점:

  1. 실제 관계를 더 잘 근사합니다.
  2. 관찰 데이터 범위를 벗어난 값에 대해서도 더 신뢰할 수 있는 예측을 제공합니다.
  3. 파라미터(Parameter) 자체가 자연적인 물리적/교육적 의미를 갖습니다 (예: 학생의 최대 학습 잠재력, 학습 속도 등).

2. 학교 현장 스토리: “우리 아이들의 단어 암기량은 어떻게 변할까?”

초등학교 1학년 학생들이 매일 아침 10분씩 영어 단어를 외운다고 가정해 봅시다. 시간이 지날수록 외운 단어 수는 늘어나겠지만, 처음에는 천천히 늘다가 어느 순간 급격히 늘고, 결국 인간의 한계 때문에 정체기에 접어들 것입니다. 이러한 ‘S자 곡선(Logistic Curve)’은 직선인 선형 모델로는 설명할 수 없습니다.

[모의 데이터 생성 및 분석을 위한 R 코드]

먼저 학교 현장과 유사한 데이터를 만들어 보겠습니다.

R

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

# 1. 모의 데이터 생성 (학생 50명, 10회 측정)
set.seed(123)
n_students <- 50
n_days <- 10
data <- data.frame(
  student = rep(1:n_students, each = n_days),
  day = rep(1:n_days, n_students)
)

# 비선형 파라미터 설정 (Asym: 최대 암기량, xmid: 가속 지점, scal: 증가 속도)
# 학생마다 개인차(Random Effects) 부여
asym_pop <- 100  # 평균 최대 암기량 100개
xmid_pop <- 5    # 평균 5일째에 급격히 상승
scal_pop <- 1    # 증가 기울기 계수

# 학생별 랜덤 효과 생성
random_effects <- rnorm(n_students, 0, 5) 

# 로지스틱 성장 곡선 적용
data$word_count <- with(data, {
  asym_i <- asym_pop + random_effects[student]
  asym_i / (1 + exp((xmid_pop - day) / scal_pop)) + rnorm(nrow(data), 0, 2)
})

# 데이터 확인
head(data)

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

3. jamovi에서의 분석 방법

jamovi에서 비선형 혼합 모델을 직접 수행하려면 GAMLj 모듈이나 Rj Editor를 활용하는 것이 좋습니다. 기본 메뉴에 비선형 혼합 모델이 없을 경우, 다음과 같이 Rj Editor를 통해 분석합니다.

jamovi 분석 단계:

  1. Library 설치: Rj - Editor에서 nlme 패키지를 호출합니다.
  2. 모델 설정: 학생별로 다른 성장 곡선을 갖도록 ‘혼합 효과(Mixed Effects)’를 설정합니다.
  3. 수식 입력:
# nlme 패키지를 이용한 비선형 혼합 모델 분석
library(nlme)
model <- nlme(word_count ~ SSlogis(day, Asym, xmid, scal),
              data = data,
              fixed = Asym + xmid + scal ~ 1,
              random = Asym ~ 1 | student,
              start = c(Asym = 100, xmid = 5, scal = 1))
summary(model)

# Asym: 학생이 도달할 수 있는 '최종 암기량' (상한선)
# xmid: 학습 효율이 최대가 되는 '시점'
# scal: 암기량이 증가하는 '가속도'

4. 결과 해석 및 시각화

분석 결과, 우리는 학생 전체의 평균적인 학습 곡선뿐만 아니라, 특정 학생이 다른 학생에 비해 얼마나 더 빨리 배우는지(xmid), 혹은 잠재력이 얼마나 큰지(Asym)를 개별적으로 파악할 수 있습니다.

R

# 시각화 코드
library(tidyverse)
ggplot(data, aes(x = day, y = word_count, group = student)) +
  geom_line(alpha = 0.3, color = "blue") + # 개별 학생 곡선
  stat_summary(aes(group = 1), fun = mean, geom = "line", size = 1.5, color = "red") + # 전체 평균
  labs(title = "학생별 영어 단어 암기량 성장 곡선",
       x = "학습 일수(Day)",
       y = "암기한 단어 수") +
  theme_minimal()

5. 주의할 점 (Limitations)

비선형 모델은 강력하지만 몇 가지 주의사항이 있습니다:

  • 복잡성: 데이터 생성 메커니즘을 정확히 모르면 모델을 세우기 어렵습니다.
  • 계산의 어려움: 선형 모델과 달리 반복적인 계산 과정이 필요하며, 때로는 수렴(Convergence)되지 않아 결과가 나오지 않을 수도 있습니다.
  • 샘플 사이즈: 파라미터를 안정적으로 추정하기 위해 충분한 양의 데이터가 필요합니다.

6. 참고문헌 (APA Style)

Wu, L., & Liu, W. (2010). Nonlinear models. In J. J. Hox & J. K. Roberts (Eds.), The SAGE Handbook of Multilevel Modeling (pp. 249-266). SAGE Publications. +1


Chap12. 다층모형과 인과추론(Multilevel Models and Causal Inference)

안녕하세요!

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

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


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

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

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


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

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

2.1 잠재적 결과 (Potential Outcomes)

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

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

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

2.2 SUTVA (Stable Unit Treatment Value Assumption)

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

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

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

3.1 무선화 실험 (Randomized Experiments)

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

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

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

yij=μ+αj+τzij+ϵijy_{ij} = \mu + \alpha_j + \tau z_{ij} + \epsilon_{ij}

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

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

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

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

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

3.2 관찰 연구 (Observational Studies)

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

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

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

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

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

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

연구 설계:

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

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

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

R

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

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

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

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

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

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

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

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

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

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

# 데이터 확인
head(data)

4.3 jamovi 분석 절차

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

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

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

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

[R 시각화 코드]

R

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

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

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

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

Step 3: 결과 해석

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

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

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

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

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

  1. 1단계: 실제 참여 여부(dd)를 배정 여부(zz)로 예측합니다.
    dij=γj+δzij+νijd_{ij} = \gamma_j + \delta z_{ij} + \nu_{ij}
  2. 2단계: 1단계에서 예측된 참여값(d^\hat{d})을 사용하여 점수(yy)를 예측합니다.
    yij=αj+τd^ij+ϵijy_{ij} = \alpha_j + \tau \hat{d}_{ij} + \epsilon_{ij}

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


6. 결론

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

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

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


참고문헌 (APA Style)

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

Chap11. 다층연구 설계시 고려사항

안녕하세요!

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

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


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

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

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

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

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

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

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

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

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

2.1 왜 다층 설계인가?

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

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

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

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

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

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

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

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

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

Koptimal=Bc+snK_{optimal} = \frac{B}{c + s \cdot n}

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

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


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

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

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

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

R

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

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

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

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

[분석 결과 해석]

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

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

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

R

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

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

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

var_values <- design_var(n_seq, K_seq, rho)

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

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

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


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

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

4.1 데이터 생성 (Backstory 포함)

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

R

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

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

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

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

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

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

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

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

# 데이터 확인
head(data_sim)

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

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

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

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

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

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

yij=β0+β1Treatmentj+u0j+ϵijy_{ij} = \beta_{0} + \beta_{1}Treatment_{j} + u_{0j} + \epsilon_{ij}

R

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

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

# 결과 요약
summary(model)

[결과 해석 방법]

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

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

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

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

이를 Multisite Trial이라고 합니다.

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

5.2 반복 측정 (Longitudinal Design)

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

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

5.3 현실적인 문제와 해결책

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

6. 결론 및 제언

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

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

참고문헌 (APA Style)

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