[R] multinomial glmnet K모수화 문제

본 게시글은 glmnet의 family='multinomial' 옵션을 사용했을 때, 출력되는 모수들을 어떻게 해석할 수 있을까 고민하며 작성했다.

 

0. 왜 문제인가?

 

하나의 예시를 들어보자.

현재 데이터에 정상 그룹 (0), 입원 그룹 (1), 사망 그룹(2)이 존재하며, 

일반적으로 multinomial regression을 해석할 때 정상그룹을 Reference 삼아 나이와 같은 특정 변수의 계수를 해석한다.

 

(ex) 나이가 많을수록 정상 그룹보다 입원 그룹에 속할 확률이 높아진다. or 나이가 많을수록 정상 그룹보다 사망 그룹에 속할 확률이 높아진다.

 

하지만, multinomial glmnet은 예제와 같이 해석할 수 없다.

왜냐하면, 첫 번째 그룹 혹은 마지막 그룹을 reference 삼은 k-1개의 그룹에 대한 계수가 아닌, 각 k개의 그룹에 대한 계수를 출력하기 때문이다. (계산상 문제 같음)

아래 glmnet의 공식 문서를 참고하면, 회귀계수가 K개의 쌍으로 존재함을 알 수 있다.

 

그럼 이를 어떻게 해석하는가?? -> 모르겠다

각 계수를 이용해 k-1 모수화 시켰을 때의 값을 얻을 수 있는가?? -> 없는 것 같다 (찾아본 바로)

 

각 그룹에 대한 계수를 모두 합하면 0이 되도록 모수화가 되어 있으며, 이는 우리의 직관적인 해석과는 거리가 멀다.

 

1. K-1 모수화 코드 짜기

 

아래의 코드는 해외의 질문글을 바탕으로 재구성됐다.

테스트를 위해 세 가지 라이브러리를 로드하고, iris 데이터를 사용하였다.

iris 데이터는 각각의 꽃을 setosa, versicolor, virginica 세 그룹 중 하나로 분류하는 문제이다.

library(glmnet)
library(pemultinom)
library(nnet)

iris.x <- as.matrix(iris[1:2])
iris.y <- as.matrix(iris[5])

# split
idx <- sample(150)
tr_idx <- idx[1:120]; test_idx <- setdiff(idx, tr_idx)

 

1-1. nnet

만약, penalized multinomial regression이 아니라면 (lambda=0), nnet 패키지의 multinom 함수를 사용하여 k-1 모수화된 계수를 얻을 수 있다.

이를 multinomial glmnet의 결과와 비교해 보자. 

tr_x <- iris.x[tr_idx, ]; tr_y <- iris.y[tr_idx, ]
test_x <- iris.x[test_idx, ]; test_y <- iris.y[test_idx, ]

# fitting via glmnet
mod.glmnet <- glmnet::glmnet(
  x = tr_x, y = tr_y,
  family = "multinomial",
  lambda = 0,
  type.multinomial = "grouped"
)

tr_data <- iris[tr_idx, ]; test_data <- iris[test_idx, ]
# fitting via nnet
mod.nnet <- nnet::multinom(
  Species ~ Sepal.Length + Sepal.Width,
  data = tr_data
)

 

먼저 multinomial glmnet의 경우 setosa, versicolor, virginica 세 그룹에 대한 회귀계수 값을 출력하는 것을 알 수 있다.

다음으로, nnet의 경우 setosa 그룹 대비 versicolor의 회귀계수와, setosa 그룹 대비 virginica 그룹의 회귀계수를 출력한다.

 

1-2. pemultinom (2023)

 

penalty항이 없다면, nnet의 multinom을 사용하면 된다. 하지만 penalty 항이 있다면 어떻게 해야 할까? 

찾아보니, 정말 최근에 penalized multinom을 지원하는 패키지가 나왔다.

다만, cv.glmnet과 달리 lambda를 직접 지정은 못 하고 lambda 개수만 지정이 가능하다.

 

적절한 비교를 위해, pemultinom과 glmnet의 lambda를 통일시켜 주었다

set.seed(1)
# default ref 는 마지막 값으로 되어 있으므로, 첫 번째 값을 ref삼고 싶다면 직접 지정해야 함.
mod.pemulti <- cv.pemultinom(x=iris.x, y=iris.y, ref='setosa', nlambda=100, nfold=10)

mod.glmnet <- cv.glmnet(
  x = iris.x, y = iris.y,
  family = "multinomial",
  lambda = mod.pemulti$lambda,
  type.multinomial = "grouped",
  nfold=10
)

 

최신에 나온 패키지라 얼마나 믿을만한지 모르겠어서 테스트를 진행해 보았다.

pemultinom으로 적합된 모형은 glmnet과 같이 'lambda.1se'와 'lambda.min' 두 가지 옵션을 제공한다.

우리가 원하는 'setosa' 그룹을 reference삼은 회귀계수가 출력됨을 확인할 수 있다.

 

해당 패키지가 최신의 패키지인데, glmnet과 완전히 같은 알고리즘으로 작동이 될까?

회귀계수 값은 모수화로 인한 문제이니 다르더라도, 같은 알고리즘이라면 최소한 선택되는 Lambda 값이나 예측값은 같아야 할 것 같다.

 

 (1) 먼저, 선택된 lambda 들을 비슷하나 조금은 다른 것을 확인할 수 있다. 

 

 (2) 다른 lambda를 선택한 1se 예측 결과를 비교해 보자. 간단한 예제여서 그럴 수 있지만, 예측 결과는 상당히 유사한 것을 알 수 있다. 

 

그리고 "prob" 옵션으로 사용했을 때, 그룹의 순서가 섞이는 것을 확인할 수 있다. prob으로 예측할 때 조심해야 한다.

추측하건대, 디폴트 옵션이 마지막 그룹을 레퍼런스로 잡는 거기 때문에... 순서가 섞이는 것 같다.