종속변수가 이항변수인 집단중심추세모델

1 개요

지난 포스트에서는 종단궤적모형¹의 기초에 대해 다루었다. 기초라 함은 종속변수가 연속형 변수인 경우를 의미한다. 그러나 쑬모있는 변수는 종종 이항변수인 경우가 많다. 해석도 편하다.

이 포스트에서는 종속변수가 이항변수인 경우의 GBTM에 대해 다룬다.

2 연결함수

다양한 형태의 종속변수를 모델에 포괄하기 위한 대표적인 방법은 연결함수(link function)를 사용하는 것이다.

평범한 회귀분석(OLS)은 사실, 종속변수가 연속변수인 경우에만 적용되는 특수한 회귀분석이다. 이 특수한 케이스에서 벗어나, 보다 광범위하고 일반적인 모델로 활용될 수 있도록 확장한 개념이 일반화 회귀분석(GLM)이다. 핵심은 연결함수.

사회과학연구에서 흔히 사용되는 로지스틱 회귀분석도 일반화회귀분석 중 하나다. 연결함수를 로지스틱 함수로 설정한 회귀분석에 해당한다.

예를 들어보자. 평범한(?) 회귀분석은 다음과 같은 관계식을 추정한다.

\[ 소득_{i}=\beta_{0}+\beta_{1}나이_{i}+e_{i} \]

그런데 만약 종속변수가 \(소득_{i}\)이 아닌 \(보험가입여부_{i}\)라면, 추정식의 해석이 모호해진다. 종속변수는 0과 1 사이의 값을 갖게 되므로 일종의 선형확률모형(LPM: linear probability model)으로 여겨진다. 그러나 나이 1살당 가입확률이 20% 오른다고 가정해보자(즉, \(\beta_{1}=0.20\)). 이렇게 되면 어떤 나이에서는 확률이 100%을 넘어간다. 반대로, 0%보다 낮아지기도 한다.

이를 방지키 위해 종속변수를 다시 또 하나의 함수에 넣어서 변환시키는 과정이 연결함수에 해당한다. 로지스틱 함수는 이렇게 생겼다.

\[ g(x)=\frac{e^x}{1+e^x} \]

이 함수 \(g(x)\)의 \(x\)자리에 원래의 추정식 \(소득_{i}\)를 넣으면 된다. 이렇게:

\[ 보험에~가입할~확률=\frac{e^{\beta_{0}+\beta_{1}나이_{i}+e_{i}}}{1+e^{\beta_{0}+\beta_{1}나이_{i}+e_{i}}} \]

즉, 로지스틱 회귀분석은 위 식을 추정하는 것.

GBTM도 마찬가지 개념이 적용된다. 보통 GBTM은 종속변수를 연속형 변수(continuous)로 구성한다.² 그러나 종단추세모형도 여타 회귀분석모델의 종속변수 문제와 동일하게, 연결함수를 사용하여 여러 형태의 종속변수를 활용할 수 있다.

다만, 이 글에서는 간단히 이항변수(binomial)에 사용할 수 있는 thresholds 함수에 대해서 알아볼 것이다. 바로 윗문단에서 언급했듯, 종속변수가 이향변수인 경우 로지스틱 함수를 연결함수로 흔히 사용한다. 깔끔하고 직관적이기 때문. 그러나 이항변수에 활용할 수 있는 연결함수로 thresholds 함수도 존재한다. 로지스틱 함수는 이미 잘 알려져 있으므로, 이 포스트에서는 thresholds 함수를 중심으로 설명한다. 이유는 logit 함수보다 thresholds 함수가 더 다양한 종속변수에 대응할 수 있기 때문이다:

• logit 함수
  확률 \(P(Y=1)\)을 직접 선형예측식으로 모델링.

• thresholds 함수
  실제로는 이항변수인 종속변수 \(Y\)를 가상적(잠재적)으로 연속변수라고 가정하고, \(Y=1\)이 되는 문턱을 설정하여 모델링. 종속변수가 2개 이상의 이산변수(categorical)라면 문턱을 여러단계 설정하여 모델링.

예를 들어보자. 종속변수가 건강보험 가입여부(\(Y_{i}=1~or~2\))인 경우, thresholds 함수를 사용한다는 것은, 개인 \(i\)의 건강보험 가입여부를 직접적으로 결정하는 잠재된 연속변수(\(Y_{i}^*\))가 있음을 가정한다는 것이다. 이 \(Y_{i}^*\)가 어떤 문턱을 넘게되면 \(Y_{i}=1\)로 판정하도록 그 문턱을 정하는 것이 thresholds 함수의 기본적인 원리다. probit 함수와 비슷하다.

3 예시 모델

종속변수가 이항변수인 GBTM을 어떻게 추정하는지 실습해보자.

3.1 개요

보험은 “불확실하고 불특정하게 발생하는 사건으로 인해 발생하는 손실분을 다수가 사전에 분배함으로써 얻어지는 확실성의 프리미엄”을 상품화한 것이다.

따라서 보험 가입자들 사이에는 “위험의 재분배”가 일어난다. 반면 보험 가입자와 미가입자 사이에는 사건발생에 따른 위험의 격차가 더 벌어진다. 한국의 4대보험은 전국민 의무가입이지만, 일부 불안정 노동자는 배제된다. 재분배 정책이 불평등을 오히려 악화시키는 것.

이에, 시간에 따라 사회보험 가입확률이 변화하는 패턴을 유형화하고자 한다. 이 때, 저소득여부(불안정 노동자)가 패턴 변화에 어떤 영향을 미치는지 확인해보기 위해 조절효과모델을 활용하여 함께 살펴보고자 한다. 요컨대, 연결함수방식을 이해한다면 GBTM은 여타 회귀분석과 동일한 바리에이션을 섞을 수 있다.

즉, 다음과 같은 GBTM 기본 모델에서:

\[ Y_{it}^{k} = \beta_{0}^{k}+\beta_{1}^{k}~time_{it}+\epsilon_{it} \]

고임금 여부(\(wage_{it}\))를 추가한 모델을 추정하는 것:³

\[ Y_{it}^{k} = \beta_{0}^{k}+\beta_{1}^{k}~time_{it}+\beta_{2}^{k}~time_{it}\times wage_{it}+\epsilon_{it} \]

3.2 실습

lcmm 패키지를 사용한다. MLE를 수행하여 추정하므로, 분석 재현성을 위해 시드값을 설정해주면 좋다. 먼저 \(k=1\)인 모델, 즉 하위집단이 없는 모델을 추정해준다. 이렇게 추정된 파라미터를 \(k=2\) 모델 추정을 위한 초기값으로 사용하는 것이 기본 골자.

몇몇 디테일은 다음과 같다:

• cat()은 현재 몇개짜리 모델을 추정하고 있는지 콘솔에 출력해준다.
• start_time, end_time은 시작시간, 완료시간을 출력한다. GBTM 연산은 시간이 많이 든다.
• lcmm() 내부의 인자 B는 초기값에 해당한다. 직접입력해도 되고, 아래 코드와 같이 \(k=1\)의 파라미터가 자동으로 입력되도록 해도 무방하다.
• saveRDS() 오랜 연산이 걸리므로, 분석결과를 별도로 저장한뒤 필요할 때 불러온다.

lcmm 패키지를 활용한 실습

set.seed(0509)
for(k in 1:7) {
  start_time <- Sys.time()
  if(k == 1) {
    cat("\n```\n")
    cat("fitting:\n")
    model09[[k]] <- lcmm(fixed = socins ~ year + higwag:year,
                         subject = "pid",
                         ng = k,
                         link = "thresholds",
                         data = df09[[5]],
                         nproc = 6)
  } else {
    model09[[k]] <- lcmm(fixed = socins ~ year + higwag:year,
                         mixture = ~ year + higwag:year,
                         subject = "pid",
                         ng = k,
                         link = "thresholds",
                         data = df09[[5]],
                         B = model09[[1]],
                         nproc = 6)
  }
  names(model09)[k] <- paste0("g", k)
  end_time <- Sys.time()
  cat("G =", k)
  cat(" takes", 
      round(as.numeric(end_time) - as.numeric(start_time, 2)),
      "second(s),\n")
  rm(k, start_time, end_time)
}

saveRDS(model09, "model09.rds")

3.3 최종 모델선정

당연히, \(k\)의 갯수를 늘릴수록 원데이터를 잘 설명하는 모델이 나온다. 100명짜리 데이터를 설명하기 위해 100개의 패턴이라면 완벽하겠다. 그러나 GBTM은 유의미한 자료요약이 목적.

모델의 설명력(goodness of fits)은 통상 AIC, BIC로 판단한다. 이외에도 \(k=n\)모델이 \(k=n-1\) 모델보다 유익한지 검증해주는 방법 들도 있는데, 자세한 내용은 지난 포스트를 참고하자. 여기서는 간단히 AIC, BIC, SABIC만 확인해보자 (Figure 1). 그림에서 보듯 \(k=4\)개 지점에서 설명량 증가폭이 감소한다(팔꿈치). 반면 Entropy는 \(k=5\)일 때 손해가 크다. \(k=4\)로 최종 결정.

Figure 1: Goodness of fits

3.4 분석결과

Figure 2 을 확인해보자. 각 \(k\)를 색으로 표현하엿다. 실선은 추정치를, 점선은 실제값을 의미한다. 그림에서 확인되듯 추정함수는 실제값 추세의 가운데를 적절히 관통하고 있는 것으로 보인다.

4개의 각 조절효과 모델을 시각화하되 좌측에는 저소득 노동자, 우측에는 고소득 노동자(\(wage_{it}=1\))의 조절효과를 시각화하였다. 양쪽 그림 모두 Class 2은 가입확률 1에 매우 가까우며(가입 유지), Class 4는 0에 매우 가깝다(미가입 유지). 소득조절효과가 잘 드러나는 집단은 Class 1이다. 이들 집단은 저소득 노동자일 경우 사회보험에서 배제될 확률이 증가한다. Class 1에 속했더라도 고소득, 즉 안정적인 일자리라면 4대보험 가입확률도 높다.

요컨대 Class 1은 사회보험 역전현상이 나타나는 집단이다. 정책수혜가 절실한 저소득 노동자가 오히려 사회정책에서 배제되고 있는 현상이 드러난다.

Figure 2: 시간에 따른 4대보험 가입확률 변화와 소득수준에 따른 조절효과

4 마치며

위 분석결과를 보며 어떤 생각이 드는가? GBTM을 통해 나타난 패턴으로부터 특정한 현상을 통찰해내는 것은 연구자의 몫이다. 분석결과를 토대로 영향요인을 분석할 수도 있고, 각 집단들의 특성을 분석할 수도 있다. 물론 이 자체로 new finding 이기도 하고.

Footnotes

1. 이하에서는 간략히 GBTM으로 명명하자.

2. 사회과학 분야에서는 연속형 변수를 가정한 모델에 순서형 변수(ordinal)까지도 사용하는 듯 하다.

3. GBTM은 이미 \(t=0\) 시점에서 다양한 \(Y\) 절편을 \(k\)개 추정한다. 따라서 모형의 경제성을 위해 \(\beta_{3}^{k}wage_{it}\)를 별도로 추정하지 않았다.