アールメカブ

• バックアップ一覧
• 差分 を表示
• 現在との差分 を表示
• ソース を表示
• R_multinom へ行く。
  • 1 (2008-12-04 (木) 09:57:50)
  • 2 (2008-12-04 (木) 10:59:58)

Rの備忘録

# 多重ロジットモデルを適用する

library(nnet)

# これは，Q（反応，A〜F）を性別（SEX），質問相手（Boy,Senior），国籍（NAT）に，さらに交互作用をすべて入れて説明しようとするモデル．

dat.m <- multinom(Q ~ SEX * BS * NAT, data = dat)
dat.m2 <- step(dat.m  )

# 結果を見る

summary( dat.m2)

coef(dat.m2) # 特に係数部分を確認
##  (Intercept)         SEX2       NATK
## B    3.446150 -0.495100658  -3.121002
## C   -1.422745 -0.172240284   1.319633
## D    3.002553  0.007826699  -3.894165
## E    1.040829 -0.498292340 -12.054348
## F   -7.508327 -1.267368205   6.150192

## この数値は対数オッズ比といわれる指標です．
## 下で説明しますが，最初の２行

##   (Intercept)         SEX2       NATK
## B    3.446150 -0.495100658  -3.121002

## というのは，基準を A，男性（１？），日本人
## として，この場合に，ABCDEFを選択する確率を
## 次のように求めることができます

(x <- predict(dat.m2, data.frame(SEX = "1", BS = "S", NAT= "J"),
   type = "probs"))
##            A            B                     C            D            
## 1.798904e-02 5.644843e-01 4.336283e-03 3.622434e-01
 ##  E                      F
 5.093713e-02 9.866955e-06
## e-02というのは，小数点を二つ前にずらすと言うことです．

## これは

exp ( c(0, coef(dat.m2)[,1]))  / sum ( exp ( c(0, coef(dat.m2)[,1]) ))

# でも求まります．

# この基準に対して，男性（１？），の韓国人が選択する確率は

(y <- predict(dat.m2, data.frame(SEX = "1", BS = "S", NAT= "K"),
        type = "probs"))
##            A            B                     C            D     
## 2.529461e-01 3.501368e-01 2.281639e-01 1.037065e-01 
##            E              F 
 ##  4.167902e-06 6.504252e-02

# となり，引き算すると

y - x

##           A           B           C           D          
##  0.23495703 -0.21434745  0.22382762 -0.25853689 
    ##     E             F
-0.05093297  0.06503266

# 韓国人男性の場合 B D E を選択する確率が日本人男性より減ります．

# これは以下の出力で，NATK 列の B D E 行の係数がマイナスであることによります

coef( dat.m2 )
##   (Intercept)         SEX2       NATK
## B    3.446150 -0.495100658  -3.121002
## C   -1.422745 -0.172240284   1.319633
## D    3.002553  0.007826699  -3.894165
## E    1.040829 -0.498292340 -12.054348
## F   -7.508327 -1.267368205   6.150192

# この -3.121002 がどうやって計算されるかは，

# まず韓国と日本（そして男）の場合をまとめて計算して

(x2 <- predict(dat.m2, data.frame(SEX = "1", BS = "S",
          NAT= c("J","K")), type = "probs"))

# オッズを求めます．オッズとは医学などで使われるものです．

log(x2[1,1] * x2[2,2] / (x2[1,2] * x2[2,1] ) )

# この意味は，日本人で A を選択する確率に 韓国人で B を選択する確率をかけ，これを日本人で B を選択する確率と韓国人でAを選択する確率を乗じたものを求め，この結果の対数をとったものです．