보조 혼합 샘플링을 이용한 베이지안 로짓 모형(Bayesian logit models with auiliary mixture sampling)

오랜만에 돌아온 논문 리뷰.

오늘은 드디어 게재된 나의 소중한 석사 논문을 리뷰해보려고 한다.

(내가 쓰고, 내가 리뷰)

---

Title : Bayesian logit models with auxiliary mixture sampling for analyzing diabetes diagnosis data

Authors : Eun Hee Rhee, Beom Seuk Hwang

---

1. Introduction

 

로짓 모형(Logit model)은 범주형 자료를 분석할 때, 가장 일반적으로 사용되는 모형으로 잘 알려져있다.

이때, 로짓 모형의 모수를 추정하는 방법으로는 Frequentist의 MLE와 Bayesian에서의 MCMC를 고려할 수 있다.

 

베이지안 관점에서 로짓 모형에 대한 추론은 오랜 시간 연구되어 왔으며, 모수에 대한 사후 분포(posterior distribution)이 closed form으로 유도되지 않으므로 주로 Metropolis-Hasting Algorithm을 통해 추론이 수행되었다.

 

물론, Metropplis-Hasting Algorithm(이하 M-H)는 모수에 대한 사후 분포가 closed form이 아니더라도 제안 분포(proposal distribution)을 통해 사후 표본을 샘플링 할 수 있다는 장점이 존재하지만, 생성된 표본들은 독립이 아니며, 모수의 차원이 커질 수록 제안 분포의 적절성을 보장 할 수 없다는 단점이 있다.

 

이에 Frühwirth-Schnatter와 Frühwirth (2007)에서는 Auxiliary latent variable을 도입하여,

모수에 대한 조건부 사후 분포 뿐만 아니라 잠재 변수에 대한 조건부 사후 분포 역시 closed form으로 유도하여 모든 추론이 Gibbs Sampling을 통해 수행되는 Auxiliary Mixture Sampling 방법을 소개하였다.

 

본 논문에서는 제안된 Auxiliary Mixture Sampling 방법의 퍼포먼스를 검증하기 위해 불균형 비율을 조절하여 모의실험을 진행하였고, 실제 데이터로는 국내 젊은 당뇨병 데이터에 적용하여 M-H와 수렴 성능을 비교하고, IRIS 데이터에 적용하여 ML 방법론들과 분류 정확도를 비교하였다.

 

이 리뷰에서는 IRIS 데이터에 적용한 결과만 간단히 리뷰하도록 하겠다.

 

---

2.  Auxiliary mixture sampling for logit model

 

Auxiliary latent variable의 도입과 sampling과정을 살펴 보기 전에, 범주형 자료에 대한 일반적인 로짓 모형의 형태를 살펴보자.

 

$K+1$개의 범주를 갖는 다범주 확률 변수 $y_1,\cdots ,y_N \in \{0,1,\cdots , K\}$가 있을 때, 기준 범주(baseline category)를 '0'으로 설정하고, 공변량 벡터 $\boldsymbol{x}_i$를 사용하여 로짓 모형을 나타내면 다음과 같다.

 

$$y_i \sim Multi(\pi_{0,i},\pi_{1,i},\cdots,\pi_{K,i})$$

$$\log\cfrac{p(y_i = k)}{p(y_i = 0)} = \log\cfrac{\pi_{k,i}}{\pi_{0,i}} = \boldsymbol{x}_i \boldsymbol{\beta}_k$$

$$ \pi_{0,i} = p(y_i=0) = \sum_{\ell=1}^K  \cfrac{1}{1+e^{\boldsymbol{x}_i\boldsymbol{\beta}_\ell}}  ~~~~~~ \pi_{k,i} = p(y_i=k) = \sum_{\ell=1}^K \cfrac{e^{\boldsymbol{x}_i\boldsymbol{\beta}_k}}{1+e^{\boldsymbol{x}_i\boldsymbol{\beta}_\ell}}$$

$$k=1,2,\cdots,K ~~~~ i=1,2,\cdots,N$$

 

베이지안 접근법에선 일반적으로 회귀 계수 $\boldsymbol{\beta}_k$에 대한 사전 분포로 정규분포를 가정하게 되며,

정규 사전 분포를 사용하였을 때, 모수 $\boldsymbol{\beta}_k$에 대한 조건부 사후 분포는 다음과 같다.

 

$$p(\boldsymbol{\beta}_k|\boldsymbol{y},\boldsymbol{X}) \propto p(\boldsymbol{y}|\boldsymbol{\beta}_k,\boldsymbol{X})p(\boldsymbol{\beta}_k)$$

$$~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ \propto \left[\prod_{i=1}^{N}\prod_{k=1}^K \pi_{k,i}^{I(y_i = k)}\right] p(\boldsymbol{\beta}_k)$$

$$~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ \propto \left[\prod_{i=1}^{N}\prod_{k=1}^K \left(\sum_{\ell=1}^{K} \cfrac{e^{\boldsymbol{x}_i\boldsymbol{\beta}_k}}{e^{\boldsymbol{x}_i\boldsymbol{\beta}_\ell}}\right)^{I(y_i=k)}\right]p(\boldsymbol{\beta}_k)$$

 

이때, Likelihood function 부분이 모수 $\boldsymbol{\beta}_k$에 대해 굉장히 복잡한 형태를 띄고 있기 때문에, 모수에 대한 조건부 사후 분포가 closed form으로 나오지 않게된다.

 

이 문제를 해결하기 위해, Frühwirth-Schnatter와 Frühwirth (2007)는 두 단계에 거친 Auxiliary latent variables를 도입하게 된다.

 

• 유틸리티 $y_{k,i}^u$ : 모수 $\boldsymbol{\beta}_k$에 대한 선형성을 만족 시키기 위해 도입됨.
• 성분 지표기 $r_{k,i}$ : 모형의 오차항 $\epsilon_{ki}$이 정규성을 만족하도록 하기 위해 도입됨.

우선, 유틸리티 $y_{k,i}^u$는 McFadden (1973)과 Scott (2011)에서 제안되었으며, 쉽게 말해 '$y_i$의 범주가 그룹 $k$로 선택될 정도' 라고 생각 할 수 있다.

즉, 유틸리티 $y_{k,i}^u$는 관측 값 $y_i$ 마다 총 $K$ 도입되게 되며, $y_{i}$의 그룹이 $k$라는 것은 $max\{ y_{0,i}^u,y_{1,i}^u,\cdots,y_{K,i}^u\}=y_{k,i}^u$ 라는 것과 같다.

도입된 유틸리티는 공변량 $\boldsymbol{x}_i$를 이용해 선형으로 모델링 되게 되며, 기준 그룹 0에 대한 유틸리티 $y_{0,i}^u$는 모든 공변량들과 독립, 즉, $\boldsymbol{\beta}_0=0$으로 가정된다.

(이와 같은 가정은 모수 $\boldsymbol{\beta}_k$에 대한 해석을 일반 로짓 모형과 동일하게 하기 위함이다.)

이때, 유틸리티 $y_{k,i}^u$에 대한 모형식은 다음과 같다.

 

$$y_{k,i}^u = \boldsymbol{x}_i\boldsymbol{\beta}_k+\epsilon_{ki} ,\hspace{0.3cm} k=1,\cdots,K ,~ ~i=1,\cdots,N$$

 

$y_{0,i}^u$와 모든 회귀 오차 $\epsilon_{ki}$의 분포는 standard Gumbel distribution으로 가정되며, 이때 $\pi_{k,i},k=0,\cdots,K$는 위에서 살펴본 일반 로짓 모형에서와 정확히 같은 꼴로 나타나게 된다.

 

두번째로 도입되는 잠재변수는 성분 지표기 $r_{k,i}$ (component indicator)이다.

앞서 언급하였듯이 회귀 오차 $\epsilon_{ki}$는 standard Gumbel distribution을 따른다.

$$p_\epsilon(\epsilon) = exp\{-\epsilon - e^{-\epsilon}\}$$

이때, 우리는 성분 지표기 $r_{k,i}$를 도입하여 $p_\epsilon(\epsilon)$를 $M$개의 성분을 갖는 혼합정규분포 $q_{M,\epsilon}(\epsilon)$으로 근사시켜 회귀 오차의 정규성을 만족시키게 된다.

 

근사를 위한 Objective function은 $p_\epsilon(\epsilon)$와 $q_{M,\epsilon}(\epsilon)$에 대한 쿨백-라이블러 발산에 혼합정규분포의 weight 합이 1이 되도록하는 페널티 항을 추가하여 정의하였다.

 

Frühwirth-Schnatter와 Frühwirth (2007)에서는 근사 결과, 구성성분이 10개($M=10$)인 혼합 정규분포가 가장 적절한 근사 분포로 나타났지만, 본 논문에서는 구성성분이 9개인 혼합정규분포가 쿨백-라이블러 값 0.0012로 가장 적절한 근사 분포로 나타났다.

모든 보조 잠재변수 유틸리티 $y_{k,i}^u$와 성분 지표기 $r_{ki}$가 도입되게 되면 모수 $\boldsymbol{\beta}_k$에 대한 Likelihood function을 Normal quadratic form으로 나타낼 수 있게 되고, Normal prior distribution이 conjugate prior가 되게 된다.

 

최종적으로 로짓 모형에 대한 auxiliary mixture sampling 단계는 다음과 같이 정리 할 수 있다.

---

•  Step 1: Initialize $y_{k,i}^{u(0)}$and $r_{ki}$. $(k=1,\cdots,K,~~i=1,\cdots,N)$
•  Step 2: Sample $\boldsymbol{\beta}_k^{(t)}$ from $p(\boldsymbol{\beta}_k|\boldsymbol{y^u},\boldsymbol{R})$, p(\boldsymbol{\beta}_k)=N_p(\boldsymbol{\mu}_{k0},\boldsymbol{\Sigma}_{k0})$$p(\boldsymbol{\beta}_k|\boldsymbol{y^u},\boldsymbol{R}) = N_p (\boldsymbol{\mu}_N,\boldsymbol{\Sigma}_N)$$$$\boldsymbol{\mu}_N=\boldsymbol{\Sigma}_N\left(\sum_{i=1}^N\frac{\boldsymbol{x}_i^T(y_{k,i}^u-m_{r_{ki}})}{s_{r_{ki}}^2}+\boldsymbol{\Sigma}_{k0}^{-1}\boldsymbol{\mu}_{k0}\right)$$$$\boldsymbol{\Sigma}_N = \left(\sum_{i=1}^N\frac{\boldsymbol{x}_i^T\boldsymbol{x}_i}{s_{r_{ki}}^2}+\boldsymbol{\Sigma}_{k0}^{-1}\right)^{-1}$$
• Step 3: Let $\lambda_{ki}^{(t)}=exp\{\boldsymbol{x}_i\boldsymbol{\beta}_k^{(t)}\},~(k=1,\cdots,K,~~i=1,\cdots,N)$        (3-1) Sample $y_{k,i}^u$$y_{k,i}^{u(0)} = -\log\left(-\frac{\log(U_i)}{1+\sum_{\ell=1}^K\lambda_{\ell i}^{(t)}}-\frac{\log(V_{ki})}{\lambda_{ki}^{(t)}}\boldsymbol{I}(y_i\neq k)\right),~ k=1,\cdots,K$$$$where~~U_i, V_{1i},V_{2i},\cdots, V_{Ki}\sim Uniform\{0,1\}$$        (3-2) Sample $r_{ki}$$$Pr(r_{ki}=j|y_{k,i}^{u(t)},\boldsymbol{\beta}_k)\propto\frac{\omega_j}{s_j}exp\left\{-\frac{1}{2}\left(\frac{y_{k,i}^{u(t)}-\log\lambda_{ki}^{(t)}-m_j}{s_j}\right)^2\right\}$$

---

3. IRIS 데이터에 적용

 

Multinomial logit modeld을 IRIS 데이터에 적용해보았고, 추론 방법으로는 MLE와 auxiliary mixture sampling을 사용하여 각각 방법에 따른 분류 성능을 비교해보았다.

 

모든 추정 과정은 R-software를 통해 수행되었고, MLE는 VGAM (version 1.1-5) 패키지를 사용하여 추정하였다.

 

Method	CV-Accuracy	Species	Correctly classified	Misclassified
Auxiliary mixture sampling	0.9800	Setosa	50	0
		Versicolor	48	2
		Virginica	49	1
MLE	0.9533	Setosa	50	0
		Versicolor	47	3
		Virginica	46	4

 

Stratified 5 - fold cross validation 방법을 이용하여 CV-Accuracy를 계산하였다.

분류 결과를 살펴보면, auxiliary mixture sampling을 추정 방법으로 사용한 multinomial logit model이 MLE를 사용한 것보다 더 좋은 분류 성능을 보인 것을 알 수 있다.

---

4. Conclusion

 

이 논문은 Frühwirth-Schnatter and Frühwirth (2007)에서 제안된 auxiliary mixture sampling의 성능을 모의실험으로 검증해보고 logit model의 두 가지 기능인 추론과 분류의 관점에서 적용하고 비교 분석해보았다.

 

빅데이터가 각광받는 지금 시대에서 베이지안 접근방법은 다소 복잡하고 시간이 오래 걸린다는 이유로 외면받고 있다.

하지만, auxiliary mixture sampling 방법처럼 Gibbs sampling을 활용하여 체인의 수렴 성능을 개선시키고, 축적된 방대한 양의 과거 자료들을 prior distribution으로 활용하여 prior에 대한 객관성을 확보할 수 있게 된다면, 블랙박스 모형들과는 달리 유의미한 해석이 가능한 모형들을 많이 사용할 수 있을 것으로 기대된다.

---

Reference

• Rhee EH and Hwang BS (2022). Bayesian logit models with auxiliary mixture sampling for analyzing diabetes diagnosis data, The Korean Journal of Applied Statistics, 35(1), 129-144
• Frühwirth-Schnatter S and Frühwirth R (2007). Auxiliary mixture sampling with applications to logistic models, Computational Statistics and Data Analysis, 51.7, 3509–3528.
• Held L and Holmes CC (2006). Bayesian auxiliary variable models for binary and multinomial regression, Bayesian Analysis, 1, 145-168
• Scott SL (2011). Data augmentation, frequentist estimation, and the Bayesian analysis of mltinomial logit models, Statistical Papers, 52, 87-109
• McFadden D (1973). Conditional logigt analysis of qualitive choice behavior, Frontiers in Econometrics, Academic press, New York, 105-142
• Shephard N (1994). Partial non-Gaussian state space, Biometrika, 81, 115-131

---

# Define function
## Function for sampling utility
multi.utility = function(y,X,beta.t){ # beta.t : m*d matrix , i : sample , j = category
  beta.t = as.matrix(beta.t)
  n = length(y) ; m = dim(beta.t)[1]
  
  lam.mat = exp(X%*%t(beta.t)); u.lam=rowSums(lam.mat)
  
  U=runif(n);V=matrix(runif(n*m),nrow=n,ncol=m)
  a = -log(U)/(1+u.lam)
  b = -log(V)/lam.mat
  
  I = matrix(NA,nrow=n,ncol=m)
  for( i in 1:n){
    I[i,] = as.numeric(y[i]!=seq(1,m))
  }
  yu = -log(a+b*I)
  return(yu)
}

## Function for computing posterior probabilities of random component
multi.mixprob = function( yu.mat, X, beta.t , mixture.mat){
  K = dim(mixture.mat)[2]
  N = dim(yu.mat)[1] ; M = dim(yu.mat)[2]
  lam.mat = exp(X%*%t(beta.t))
  prob.r = list()
  prob.r.yi = matrix(NA, nrow = M , ncol = K)
  
  for ( i in 1:N){
    yu = yu.mat[i,]
    lam = lam.mat[i,]
    for( m in 1:M){
      yu.m = yu[m]
      lam.m = lam[m]
      pr.r = numeric(K)
      for( k in 1:K){
        w=mixture.mat[1,k];mean=mixture.mat[2,k];s=sqrt(mixture.mat[3,k])
        pr.r[k]=w/s*exp(-((yu.m-log(lam.m)-mean)/s)^2/2)
      }
      prob.r.yi[m,]=pr.r/sum(pr.r)
    }
    prob.r[[i]]=prob.r.yi
  }
  return(prob.r)
}

## Function for sampling random component
multi.mixidx = function(mixprob.mat){
  N=length(mixprob.mat) ; M = dim(mixprob.mat[[1]])[1];K = dim(mixprob.mat[[1]])[2]
  r.idx = matrix(NA,nrow=N,ncol=M)
  
  for( i in 1:N){
    for(m in 1:M){
      r.idx[i,m] = sample(1:K,1,prob = mixprob.mat[[i]][m,])
    }
  }
  return(r.idx)
}


###############################################################################
###############################################################################
#                        Data analysis :IRIS                                  #
###############################################################################
###############################################################################

rm(list=ls())
library(mvtnorm)
library(nnet)
library(dplyr)
library(ggplot2)
library(caret)
library(randomForest)
###############################################################################

### Split data into train and test set.
n=dim(iris)[1]

library(splitstackshape)
library(tibble)

d <- rownames_to_column(DF, var = "ID") %>% mutate_at(vars(ID), as.integer)
flds=list()
set.seed(1234)
for(i in 1:5){
  j=5-(i-1)
  if(j>1){
    a=stratified(d, "Species", size = 1/j)
    flds[[i]]=a$ID
    d=d%>%filter(ID %in% setdiff(d$ID,a$ID))
  } else{
    flds[[i]]=d$ID  
  }
}

#############################################################################################
#                                           MLE                                             #
#############################################################################################
library(foreach)

numCores <- parallel::detectCores() - 1

myCluster <- parallel::makeCluster(numCores)
doParallel::registerDoParallel(myCluster)
idx=1
function_reslut=foreach::foreach(idx=1:5,.combine = 'c',.export = "multinom")%dopar%{
  k.idx=flds[[idx]]
  predict(multinom(Species~.,data=DF[-k.idx,]),DF[k.idx,],type="class")
}

#############################################################################################
#                                   Auxiliary mixture sampling                              #
#############################################################################################
Kfold.Obj=list()
idx.list=list(c(1,2,3),c(4,5,6),c(7,8,9),c(10,11,12),c(13,14,15))
for( k in 1:5){
  # train
  k.idx=flds[[k]]
  y=as.numeric(DF$Species[-k.idx])-1
  X=as.matrix(cbind("Intercept"=rep(1,length(y)),DF[-k.idx,-5]));q=dim(X)[2]
  name<-colnames(X)
  
  Niter=50000;Nburn=Niter/2
  b0=matrix(rep(0,q),ncol=1)
  B0=diag(rep(100,q))
  
  ## Set Initial point ##
  set.seed(1234)
  init.b=matrix(rnorm(q*2),ncol=q) # d*m
  
  yu.mat=multi.utility(y,X,init.b)
  prob.R=multi.mixprob(yu.mat,X,init.b,mix.mat)
  R=multi.mixidx(prob.R)
  
  ## Save Object ##
  iter.beta1 = matrix(NA,nrow=(Niter+1),ncol=q)
  iter.beta2 = matrix(NA,nrow=(Niter+1),ncol=q)
  
  iter.beta1[1,]=init.b[1,]
  iter.beta2[1,]=init.b[2,]
  
  start.t=as.POSIXct(Sys.time())
  
  set.seed(1234)
  for ( it in 2:(Niter+1)){
    
    # (a) Sample beta_k
    
    ## k=1
    b.1N=get.MV(X,R[,1],yu.mat[,1],mix.mat,b0,B0)[[1]]
    B.1N=get.MV(X,R[,1],yu.mat[,1],mix.mat,b0,B0)[[2]]
    
    iter.beta1[it,]=rmvnorm(1,mean=b.1N,sigma=B.1N)
    
    ## k=2
    b.2N=get.MV(X,R[,2],yu.mat[,2],mix.mat,b0,B0)[[1]]
    B.2N=get.MV(X,R[,2],yu.mat[,2],mix.mat,b0,B0)[[2]]
    
    iter.beta2[it,]=rmvnorm(1,mean=b.2N,sigma=B.2N)
    
    ### Make beta.t mat
    beta.t = as.matrix(rbind(iter.beta1[it,],iter.beta2[it,]))
    
    
    # (b1) Sample latent utility
    yu.mat=multi.utility(y,X,beta.t)
    
    # (b2) Sample latent component indicator R
    prob.R=multi.mixprob(yu.mat,X,beta.t,mix.mat)
    R=multi.mixidx(prob.R)
    
    if( it%%1000==0){
      print(paste("iter=",it))
    }
    
  }
  
  R.beta1=Get.estimation(iter.beta1);rownames(R.beta1)=paste("versicolor::",name)
  R.beta2=Get.estimation(iter.beta2);rownames(R.beta2)=paste("virginica::",name)
  
  Kfold.idx=idx.list[[k]]
  
  Kfold.Obj[[Kfold.idx[1]]]=rbind(R.beta1,R.beta2)
  Kfold.Obj[[Kfold.idx[2]]]=iter.beta1
  Kfold.Obj[[Kfold.idx[3]]]=iter.beta2
  print(paste("Finish fold-",k,"Time:",as.POSIXct(Sys.time())-start.t))
}