廿TT

譬如水怙牛過窓櫺 頭角四蹄都過了 因甚麼尾巴過不得

順序のある2次元分割表の比例ハザードモデル

モデルと尤度関数

2次元分割表 n_{ij} (i=1,\ldots,I, j=1,\ldots,J) が得られたとする.

たとえばこんなふうだ.

強く反対 反対 中立 賛成 強く賛成
皆無 15 5 6 0 1
まれ 8 17 13 7 2
ときどき 3 4 4 7 6
しばしば 2 0 3 5 3
常時 0 2 2 9 16

このデータは柳本・清水(1983)からの孫引きで, 行が悪夢を見る頻度, 列が睡眠薬に対する賛否を表している.

原典は西里『質的データの数量化』とのこと.

まず行 i を固定して考える.

「強く反対」から「強く賛成」の背後には見えない連続型の確率変数 X があり, それが適当な間隔で区切られることで, この結果が確認されたものと考える.

その適当な区切りを, C_1=[c_0,c_1), C_2=[c_1,c_2), \ldots , C_J = [c_{J-1},c_J) とする. 等間隔である必要はない.

説明変数を z とし, ハザード関数 \lambda_j は生存関数 S(x|z) を用いて,

 \displaystyle \lambda_j = \lambda(C_j|z) \\
\displaystyle \quad = \Pr(X \in C_j| X > c_{j-1}) \\
\displaystyle \quad = \frac{S(c_{j-1}|z) - S(c_{j}|z)}{S(c_{j-1}|z)}

と定義される.

離散型の比例ハザードモデルは, オッズみたいなものを考え,

 \displaystyle \frac{\lambda_j}{1-\lambda_j} = \exp(\alpha_j+\beta z)

という形で回帰構造を入れる.

このとき, \alpha\beta を推定するための尤度関数は,


 \displaystyle L \propto \prod_{j=1}^{J-1}\lambda_{j}^{n_{ij}} (1-\lambda_{j})^{\sum_{l=j}^{J}n_{il}-n_{ij}}

である.

以上を全部の行について考えるので, 行についての添え字 i を導入し,


 \displaystyle L \propto \prod_{i=1}^{I} \prod_{j=1}^{J-1}\lambda_{ij}^{n_{ij}} (1-\lambda_{ij})^{\sum_{l=j}^{J}n_{il}-n_{ij}}

が最大化すべき尤度である.

対数を取ると,  \displaystyle \lambda_j = \exp(\alpha_j+\beta z)/(1+\exp(\alpha_j+\beta z)) より,


 \displaystyle \log L \propto \sum_{i=1}^{I} \sum_{j=1}^{J-1}\left[ n_{ij} (\alpha_j+\beta z_{i}) - \left(\sum_{l=j}^{J}n_{il}\right)\log \left(1+\exp(\alpha_j+ \beta z_{i})\right) \right].

さて, ハザードではなく分割表の各セルの確率 p_{ij}を評価したい場合がある.

その場合,


 \displaystyle \frac{\lambda_{ij}}{1-\lambda_{ij}} = p_{ij}/\left(1-\sum_{l=1}^{j-1}p_{il}\right)  = \exp(\alpha_j + \beta z_i)

より

\displaystyle p_{ij}=\frac{\exp(\alpha_j + \beta z_i)}{1+ \exp(\alpha_j + \beta z_i)}\left(1 - \sum_{l=1}^{j-1} p_{il}\right)

を見ればよい.

ただし \displaystyle p_{i1}=\frac{\exp(\alpha_1 + \beta z_i)}{1+ \exp(\alpha_1 + \beta z_i)}.

R によるケーススタディ

冒頭の分割表に比例ハザードモデルを当てはめる.

対数尤度を直接 optim にいれてやれば,推定値が得られる.

まずは「皆無」から「常時」を z=(1,2,3,4,5)' とコード化してモデルのパラメータを推定する(model 1).

観測値と予測値を並べてプロットしてやるとこんな感じだ.

f:id:abrahamcow:20170902192017p:plain

library(tidyverse)

n <- matrix(c(15,5,6,0,1,
              8,17,13,7,2,
              3,4,4,7,6,
              2,0,3,5,3,
              0,2,2,9,16),5,5,byrow = TRUE)

LL <- function(n,m,alpha,beta){
  J <- ncol(n)
  I <- nrow(n)
  out <- 0
  for(j in 1:(J-1)){
    for(i in 1:I){
      out <- out +
        n[i,j]*(alpha[j]+beta*i)-m[i,j]*log(1+exp(alpha[j]+beta*i))
      }
  }
  out
}
wrapped_LL <- function(n,m){
  J <- ncol(n)
  function(par){LL(n,m,par[1:(J-1)],par[-c(1:(J-1))])}
}

J <- ncol(n)
I <- nrow(n)
m <- cbind(sapply(1:(J-1),function(j)apply(n[,j:J],1,sum)),n[,J])
ndot <- apply(n,2,sum)
mdot <- apply(m,2,sum)
opt1 <-optim(c(log(ndot/(mdot-ndot))[-J],0),wrapped_LL(n,m),control = list(fnscale=-1,maxit=5000))

alphahat <-opt1$par[1:(J-1)]
betahat <-opt1$par[-c(1:(J-1))]

pmat <- matrix(NA,I,J)
tmp <- t(sapply(1:I,function(i)plogis(alphahat+betahat*i)))
pmat[,1] <- tmp[,1]
for(i in 1:I){
  for(j in 2:(J-1)){
    pmat[i,j] <- tmp[i,j] - sum(pmat[i,1:(j-1)])*tmp[i,j]
  }
}
pmat[,J] <-1-apply(pmat,1,sum,na.rm=TRUE)

emat <-t(sapply(1:J,function(j)pmat[j,]*ndot[j]))

tidy_e <-as_data_frame(emat) %>% 
  mutate(row=row_number()) %>% 
  gather(col,value,-row) %>% 
  mutate(type="pred")

tidy_n <-as_data_frame(n) %>% 
  mutate(row=row_number()) %>% 
  gather(col,value,-row) %>% 
  mutate(type="obs")

fitdf <- bind_rows(tidy_e,tidy_n)
p1 <-ggplot(fitdf,aes(x=type,y=value,fill=type))+
  geom_col()+
  facet_grid(row~col,labeller = labeller(col=function(x)sub("V","",x)))+
  ggtitle("model1")
print(p1)

ところで, 順序尺度のデータを単純に1, 2, 3, 4, 5と置いてしまうことに抵抗のある方も多かろう.

ダミー変数を使って「皆無」から「常時」をそれぞれ別の係数として推定してやるとこうなる(model 2).

LL2 <- function(n,m,alpha,beta){
  J <- ncol(n)
  I <- nrow(n)
  out <- 0
  beta2 <- c(0,beta)
  for(j in 1:(J-1)){
    for(i in 1:I){
      out <- out +
        n[i,j]*(alpha[j]+beta2[i])-m[i,j]*log(1+exp(alpha[j]+beta2[i]))
    }
  }
  out
}
wrapped_LL2 <- function(n,m){
  J <- ncol(n)
  function(par){LL2(n,m,par[1:(J-1)],par[-c(1:(J-1))])}
}
par <- c(log(ndot/(mdot-ndot))[-J],numeric(I-1))

opt2 <-optim(par,wrapped_LL2(n,m),control = list(fnscale=-1,maxit=5000))

alphahat2 <-opt2$par[1:(J-1)]
betahat2 <-c(0,opt2$par[-c(1:(J-1))])

pmat2 <- matrix(NA,I,J)
tmp <- t(sapply(1:I,function(i)plogis(alphahat2+betahat2[i])))
pmat2[,1] <- tmp[,1]
for(i in 1:I){
  for(j in 2:(J-1)){
    pmat2[i,j] <- tmp[i,j] - sum(pmat2[i,1:(j-1)])*tmp[i,j]
  }
}
pmat2[,J] <-1-apply(pmat2,1,sum,na.rm=TRUE)

emat2 <-t(sapply(1:J,function(j)pmat2[j,]*ndot[j]))

tidy_e2 <-as_data_frame(emat2) %>% 
  mutate(row=row_number()) %>% 
  gather(col,value,-row) %>% 
  mutate(type="pred")

fitdf2 <- bind_rows(tidy_e2,tidy_n)
p2 <-ggplot(fitdf2,aes(x=type,y=value,fill=type))+
  geom_col(position = "dodge")+
  facet_grid(row~col)+
  facet_grid(row~col,labeller = labeller(col=function(x)sub("V","",x)))+
  ggtitle("model2")
print(p2)

f:id:abrahamcow:20170902192925p:plain

当てはまりには大差なく, このデータに関してはいえば, あえてパラメータを増やすことにさして意味はなさそうである.

model 1, model 2, それぞれの AIC は以下のようになった.

> #model 1
> -2*opt1$value + 2*length(opt1$par)
[1] 396.6306
> #model 2
> -2*opt2$value + 2*length(opt2$par)
[1] 402.1568

付録

柳本・清水(1983)では対数尤度を偏微分して, ニュートン法で =0 となる解を求めている.

model 1 のパラメータを求めるコードを R で書くとこんな感じだ.

fitph <-function(n,maxit=1000,tol=1e-8){

dLLalpha <- function(j,n,m,beta,x){
  function(alpha_j){
    sum(n[,j]-m[,j]*exp(alpha_j+beta*x)/(1+exp(alpha_j+beta*x)))
  }
}

dLLbeta <- function(n,m,alpha,x){
  function(beta){
  out <- 0
  for(j in 1:(J-1)){
    for(i in 1:I){
      out <- out +
        n[i,j]*x[i]-m[i,j]*x[i]*exp(alpha[j]+beta*x[i])/(1+exp(alpha[j]+beta*x[i]))
    }
  }
  out
  }
}
J <- ncol(n)
I <- nrow(n)
x=1:I
m <- cbind(sapply(1:(J-1),function(j)apply(n[,j:J],1,sum)),n[,J])
ndot <- apply(n,2,sum)
mdot <- apply(m,2,sum)
alpha = log(ndot/(mdot-ndot))[-J]
beta = 0
for(i in 1:maxit){
  lalpha <-lapply(1:(J-1),function(j)uniroot(dLLalpha(j,n,m,beta,x),c(-50,50)))
  alphanew <- sapply(lalpha,function(x)x$root)
  lbeta <-uniroot(dLLbeta(n,m,alphanew,x),c(-50,50))
  betanew <- lbeta$root
#  print(i)
  if(all(abs(alpha-alphanew)<tol & abs(beta-betanew)<tol)){
    break
    }
  alpha <-alphanew
  beta <-betanew
}

return(list(alpha=lalpha,beta=lbeta,iter=i))
}

fit1 <-fitph(n)

alphahat1_2 <-sapply(fit1$alpha,function(x)x$root)
betahat1_2 <-fit1$beta$root

実行してみると optim でやったのとほぼ同じ結果が得られることがわかるだろう.

参考文献

柳本武美・清水央子. (1983). 2次元分割表における比例ハザードモデルの適用. 応用統計学, 12(1), 17-29.

abrahamcow.hatenablog.com

Cox比例ハザードモデル (医学統計学シリーズ)

Cox比例ハザードモデル (医学統計学シリーズ)

中村剛『Cox比例ハザードモデル 』は, 今回のエントリで直接参照したわけではないが, 数学的なところも丁寧に書いてあるので, 比例ハザードモデルそのものについて学びたい人にはおすすめ。

質的データの数量化―双対尺度法とその応用 (統計ライブラリー)

質的データの数量化―双対尺度法とその応用 (統計ライブラリー)

西里静彦『質的データの数量化』は読んでません。

むすび

ディズニーランド一緒に行ってくれる人募集中です。

abrahamcow.hatenablog.com