東京大学・株式会社Nospareの菅澤です．
前回の記事に続き，今回は順位尤度による正規コピュラモデルのベイズ推定法を実装し，数値実験とデータ解析例を紹介していきます．

設定のおさらい

$y_i=(y_{i1},\ldots,y_{ip}) \ (i=1,\ldots,n)$を$p$次元の連続離散混合データとします．このとき，$y_i$の同時分布をモデル化するため，潜在変数$z_i=(z_{i1},\ldots,z_{ip})\sim N(0, {\bf \Psi})$を導入し，$y_i$の各要素を$y_{ik}=G_k(z_{ik})$の形で表現します．ここで，${\bf \Psi}$は潜在変数の相関行列で，$G_k$は適当な非減少関数です．推定対象は相関行列の${\bf \Psi}$です．

前回の記事で紹介したように，順位尤度(rank likelihood)を使うことで，${\bf \Psi}$の事後分布を直接的に与えることができ，パラメータ拡大法による簡便なMCMCアルゴリズム(ギブスサンプラー)によって${\bf \Psi}$の事後サンプルが生成できます．

MCMCアルゴリズム (Rでの実装)

ギブスサンプリングを実装した関数RL_mcmcを以下のように与えます．

##   MCMC algorithm with rank likelihood   ##
# Y: (n, p)-matrix
# mc: MCMC length
# bn: burn-in length
# nu0, S0: hyperparameters of invserse-Wishart prior 
library(truncnorm)
library(MCMCpack)

RL_mcmc <- function(Y, mc=1000, bn=200, nu0=1, S0=diag(p)){
  n <- dim(Y)[1]
  p <- dim(Y)[2]
  Rank_data <- NULL
  for(j in 1:p){ 
    Rank_data <- cbind(Rank_data, match(Y[,j],sort(unique(Y[,j])))) 
  }  
  Rlevels <- apply(Rank_data, 2, max, na.rm=TRUE)
  Ranks <- apply(Y, 2, rank, ties.method="max", na.last="keep")
  N <- apply(!is.na(Ranks), 2, sum)
  U <- t( t(Ranks)/(N+1) )   # normalized rank
  
  # initial values
  Z <- qnorm(U)    # latent variable
  Zfill <- matrix(rnorm(n*p), n, p)
  Z[is.na(Y)] <- Zfill[is.na(Y)]  # imputation of latent variable for missing Y
  Sig <- cov(Z)
  
  ##  MCMC iterations
  Phi.pos <- array(NA, c(mc, p, p))
  for(k in 1:mc) {
    # update Z
    for(j in 1:p){
      Yj <- Y[,j]
      Zj <- Z[,j]
      vv <- Sig[j,j] - Sig[j,-j]%*%solve(Sig[-j,-j])%*%Sig[-j,j]   # conditional variance
      mm <- as.vector( Sig[j,-j]%*%solve(Sig[-j,-j])%*%t(Z[,-j]) )   # conditional expectation 
      for(i in 1:n){
        y <- Yj[i]
        lb <- suppressWarnings( max(Zj[Yj<y], na.rm=T) )
        ub <- suppressWarnings( min(Zj[Yj>y], na.rm=T) )
        Z[i,j] <- rtruncnorm(1, a=lb, b=ub, mean=mm[i], sd=sqrt(vv))
      }
    }
    # update covariance-matrix
    invSig <- rwish(nu0+n, solve(S0+t(Z)%*%Z))
    Sig <- solve( invSig )
    
    # save correlation matrix
    DD <- sqrt(diag(Sig))
    Phi.pos[k,,] <- Sig/(DD%*%t(DD))
  }
  
  # output
  Phi.pos <- Phi.pos[-(1:bn),,]
  return(Phi.pos)
}

この関数は，データ行列Yを入力として，${\bf \Psi}$の事後サンプルが生成される関数となっています．

また，データ行列Yに欠測値があったとしても，ランダムな欠測(いわゆるmissing at random)の仮定のもとで自動的に補完しながら推定できるアルゴリズムになっています．

数値実験

まずは簡単な数値実験を用いて推定の性能を確かめてみます．以下のようにして擬似データを生成します．

library(MASS)
set.seed(123)

n <- 100  # sample size
p <- 6    # dimension 

# true correlation 
Psi_true <- matrix(NA, p, p)
for(j in 1:p){
  for(k in 1:p){
    Psi_true[j, k] <- (0.7)^(abs(k-j))
  }
}

# data generation 
LY <- mvrnorm(n, rep(0, p), Psi_true)　　   # latent variables
Y <- matrix(NA, n, p)

# Y1 (four ordered categories)
Y[,1] <- 1 
Y[LY[,1]>(-1), 1] <- Y[LY[,1]>(-1), 1] + 1
Y[LY[,1]>0, 1] <- Y[LY[,1]>0, 1] + 1
Y[LY[,1]>1, 1] <- Y[LY[,1]>1, 1] + 1

# Y2 & Y3 & Y5 (continuous)
Y[,c(2,3,5)] <- LY[,c(2,3,5)]

# Y4 & Y6 (binary)
Y[,4] <- ifelse(LY[,4]>0, 1, 0)
Y[,6] <- ifelse(LY[,6]>0, 1, 0)

# pairwise scatter plot
pairs(Y)

この設定では，$p=6$次元のデータを$n=100$個生成しています．6変数のうち，1つは4段階の順序カテゴリー変数，2つは二値変数，残り3つは連続変数のデータになっています．生成データの対散布図は以下のようになります．

作成したRL_mcmc関数を用いて事後サンプルを生成し，${\bf \Psi}$の事後平均と95%信用区間を計算します．

posterior <- RL_mcmc(Y, mc=2000, bn=500)

Psi_PM <- apply(posterior, c(2,3), mean)   # posterior mean
Psi_CI <- apply(posterior, c(2,3), quantile, prob=c(0.025, 0.975))   # 95% credible interval

結果は以下のようになりました．(以下では上三角要素をベクトル化した結果を表示します．)

> Upper_tri <- function(mat){ mat[upper.tri(mat)] }    # 上三角成分を抽出する関数を定義しておく
> round(Upper_tri(Psi_true), 3)   # true value
 [1] 0.700 0.490 0.700 0.343 0.490 0.700 0.240 0.343
 [9] 0.490 0.700 0.168 0.240 0.343 0.490 0.700
> round(Upper_tri(Psi_PM), 3)  
 [1] 0.748 0.423 0.702 0.286 0.228 0.525 0.202 0.216
 [9] 0.403 0.632 0.110 0.182 0.378 0.451 0.735
> round(rbind(Upper_tri(Psi_CI[1,,]), Upper_tri(Psi_CI[2,,])), 3)
      [,1]  [,2]  [,3]  [,4]   [,5]  [,6]  [,7]  [,8]  [,9]
lower 0.633 0.239 0.590 0.033 -0.021 0.317 0.004 0.033 0.221
upper 0.836 0.590 0.793 0.502  0.432 0.693 0.391 0.390 0.561
     [,10]  [,11]  [,12] [,13] [,14] [,15]
lower 0.456 -0.174 -0.076 0.150 0.157 0.583
upper 0.774  0.368  0.420 0.585 0.681 0.849

事後平均の値は比較的真値に近い値になっていることがわかります．また，95%信用区間はほとんどの真値をカバーしていることもわかります．今回の例は$n=100$と比較的サンプルサイズが小さい例ですが，ある程度の精度で相関行列${\bf \Psi}$を推定できていることが確認できます．

データ解析 (社会流動性データ)

次にHoff (2007)で扱われたデータ解析例を再現してみます．このデータは1994年に米国で実施された総合的社会調査に基づくものです．回答者の所得(INC)，学士，修士，博士などの最高学位(DEG)，子供の数(CHILD)，年齢(AGE)とともに，回答者の親に関しても所得(PINC)，最高学位(PDEG)，子供の数(PCHILD)のデータが観測されています．詳細についてはHoff (2007)を参照してください．

このデータから上記7つの変数の相関行列を推定することを考えます．サンプルサイズは$n=1002$で，DEGおよびPDEGは5段階の順序カテゴリー変数，CHILDおよびPCHILDはカウント変数となっており，様々な型の変数が混ざった多次元データになっています．

このデータのファイルsocmob.RDataはHoff氏のサイトからダウンロードできます．

関数RL_mcmcを用いて事後サンプルを生成し，事後平均と95%信用区間を計算します．

load("socmob.RData") 
posterior <- RL_mcmc(Y=socmob, mc=1500, bn=500)

variable_name <- c("INC", "DEG", "CHILD", "PINC", "PDEG", "PCHILD", "AGE")
dimnames(posterior)[[2]] <- dimnames(posterior)[[3]] <- variable_name

Psi_PM <- apply(posterior, c(2,3), mean)   # posterior mean
Psi_CI <- apply(posterior, c(2,3), quantile, prob=c(0.025, 0.975))   # 95% credible interval

結果は以下のようになります．

> round(Psi_PM, 3)      # 事後平均
          INC    DEG  CHILD   PINC   PDEG PCHILD    AGE
INC     1.000  0.496  0.308  0.137  0.185 -0.045  0.337
DEG     0.496  1.000 -0.033  0.203  0.475 -0.200  0.047
CHILD   0.308 -0.033  1.000 -0.151 -0.248  0.230  0.601
PINC    0.137  0.203 -0.151  1.000  0.448 -0.214 -0.133
PDEG    0.185  0.475 -0.248  0.448  1.000 -0.287 -0.232
PCHILD -0.045 -0.200  0.230 -0.214 -0.287  1.000  0.119
AGE     0.337  0.047  0.601 -0.133 -0.232  0.119  1.000

> round(Psi_CI[1,,], 3)    # 信用区間の下限
          INC    DEG  CHILD   PINC   PDEG PCHILD    AGE
INC     1.000  0.437  0.244  0.033  0.115 -0.106  0.278
DEG     0.437  1.000 -0.106  0.112  0.413 -0.262 -0.017
CHILD   0.244 -0.106  1.000 -0.243 -0.313  0.166  0.556
PINC    0.033  0.112 -0.243  1.000  0.356 -0.300 -0.216
PDEG    0.115  0.413 -0.313  0.356  1.000 -0.356 -0.296
PCHILD -0.106 -0.262  0.166 -0.300 -0.356  1.000  0.063
AGE     0.278 -0.017  0.556 -0.216 -0.296  0.063  1.000

> round(Psi_CI[2,,], 3)　    # 信用区間の上限
         INC    DEG  CHILD   PINC   PDEG PCHILD    AGE
INC    1.000  0.551  0.369  0.226  0.253  0.022  0.393
DEG    0.551  1.000  0.034  0.291  0.534 -0.134  0.114
CHILD  0.369  0.034  1.000 -0.051 -0.179  0.292  0.641
PINC   0.226  0.291 -0.051  1.000  0.531 -0.125 -0.039
PDEG   0.253  0.534 -0.179  0.531  1.000 -0.220 -0.164
PCHILD 0.022 -0.134  0.292 -0.125 -0.220  1.000  0.178
AGE    0.393  0.114  0.641 -0.039 -0.164  0.178  1.000

このように，複数の型の変数が混ざっていても相関行列の推測が行えることがわかります．今回の例では比較的サンプルサイズが大きい($n=1002$)ため，各要素の信用区間の幅が適度に短くなっている(効率的な推定が行えている)ことも確認できます．

今回は${\bf \Psi}$に対する推定結果を直接的に見ましたが，${\bf \Psi}$の事後サンプルを用いて条件付き独立性に関係する量の事後サンプルを生成することもでき，変数間の関係性をグラフ的に表現することも可能となります．この点の詳細な議論についてはHoff (2007)の元論文をご覧ください．

おわりに

今回は順位尤度(rank likelihood)を用いて従属構造をベイズ推定する方法を実装し，数値実験とデータ解析例について紹介しました．

株式会社Nospareでは統計学の様々な分野を専門とする研究者が所属しております．統計アドバイザリーやビジネスデータの分析につきましては株式会社Nospareまでお問い合わせください．