トップ > R

R でアソシエーション分析 - arules

R

R言語の arules を使ってアソシエーション分析を試してみました。

ソースは http://github.com/fits/try_samples/tree/master/blog/20180108/

はじめに

データセット

今回は、適当に作った下記データセット （100行） を使います。

1行が 1つの取引で , で区切られたアルファベット群が同時に購入された商品とみなします。（例えば、T,Y,C なら T と Y と C 商品を同時購入）

data.basket
C,S,M,R
T,Y,C
P,Y,C,M
O,W,L
R
O,U,R,L
P
W,C
T,O,W,B,C
C,T,W,B,Y,F,D
・・・
R,P,S
B,
B,S
B,F
C,F,N

インストール

install.packages で arules をインストールしておきます。

インストール例
> install.packages("arules")

実装

arules の apriori 関数で Apriori アルゴリズムを使ったアソシエーションルールの抽出を行えます。

read.transactions でデータファイルを transactions オブジェクト化して、apriori の data 引数として使います。

parameter 引数を使って support（支持度）が 0.05 以上 ※、confidence（確信度）が 0.7 以上のものを抽出するようにしました。

 ※ 今回使用するデータセットの件数が 100件なので 0.05 は 5件になる
sample.R
library(arules)

args <- commandArgs(TRUE)

# データファイルを transactions オブジェクト化
tr <- read.transactions(args[1], format = "basket", sep = ",")

# アソシエーションルールの抽出
tr.ap <- apriori(tr, parameter = list(support = 0.05, confidence = 0.7))

# lift 値でソート
inspect(sort(tr.ap, by = "lift"))

実行

Rscript で実行してみます。

実行結果（support = 0.05, confidence = 0.7）
> Rscript sample.R data.basket

・・・
    lhs      rhs support confidence lift     count
[1] {B,O} => {W} 0.05    0.8333333  3.333333  5
[2] {B,T} => {C} 0.05    1.0000000  2.439024  5
[3] {T}   => {C} 0.08    0.8000000  1.951220  8
[4] {N}   => {C} 0.10    0.7142857  1.742160 10

補足

なお、当然ながら support や confidence の条件を緩和すると結果数が大きく変わります。

実行結果2（support = 0.03, confidence = 0.5）
> Rscript sample2.R data.basket

・・・
      lhs        rhs support confidence lift     count
[1]   {B,C,W} => {T} 0.04    0.8000000  8.000000  4
[2]   {B,C,F} => {T} 0.03    0.7500000  7.500000  3
[3]   {B,C,R} => {T} 0.03    0.7500000  7.500000  3
[4]   {F,W}   => {T} 0.03    0.6000000  6.000000  3
[5]   {B,R,S} => {Y} 0.03    0.6000000  6.000000  3
・・・
[128] {O,W}   => {B} 0.05    0.5000000  1.190476  5
[129] {D,W}   => {B} 0.03    0.5000000  1.190476  3
[130] {D,S}   => {B} 0.03    0.5000000  1.190476  3

R の MXNet で iris を分類

R DeepLearning

「MXNet で iris を分類」 と同様の処理を R言語で実装してみました。

ソースは http://github.com/fits/try_samples/tree/master/blog/20171212/

準備

今回は下記サイトの手順に従って MXNet R パッケージの CPU 版を Windows へインストールしました。

MXNet R パッケージの CPU 版を Windows へインストール
cran <- getOption("repos")
cran["dmlc"] <- "https://apache-mxnet.s3-accelerate.dualstack.amazonaws.com/R/CRAN/"
options(repos = cran)
install.packages("mxnet")

インストールした mxnet のバージョンが少し古いようですが（現時点の MXNet 最新バージョンは 1.0）、今回はこれを使います。

バージョン確認
> packageVersion('mxnet')
[1] ‘0.10.1’

学習と評価

MXNet には、階層型ニューラルネットワークの学習処理を簡単に実行するための関数 mx.mlp が用意されているので、今回はこれを使います。

引数 備考
hidden_node 隠れ層のノード（ニューロン）数（デフォルトは 1）
out_node 出力ノード数（今回は分類数）
num.round 繰り返し回数（デフォルトは 10）
array.batch.size バッチサイズ（デフォルトは 128）
learning.rate 学習係数
activation 活性化関数（デフォルトは tanh）

hidden_node にベクトル（例. c(6, 4)）を設定すれば隠れ層が複数階層化されるようです。

iris のデータセットは R に用意されているものを使います。

mx.mlp の入力データには mx.io.DataIter か R の配列 / 行列を使う必要があるようなので（ラベルデータは配列のみ）、data.matrix で行列化しています。

ラベルデータとする iris の種別 iris$Species は因子型ですが、mxnet では因子型を扱えないようなので as.numeric で数値化しています。

ここで as.numeric の結果は 1 ～ 3 の数値になりますが、mxnet で 3種類の分類を行うには 0 ～ 2 でなければならないようなので -1 しています。

一方、predict の結果を max.col(t(<predictの結果>)) で処理すると 1 ～ 3 の数値になるため、評価用のラベルデータは -1 せずに正解率の算出に使っています。

また、array.layout = 'rowmajor' は Warning message 抑制のために設定しています。

iris_hnn.R
library(mxnet)

train_size = 0.7

n = nrow(iris)
# 1 ～ n から無作為に n * train_size 個を抽出
perm = sample(n, size = round(n * train_size))

# 学習用データ
train <- iris[perm, ]
# 評価用データ
test <-iris[-perm, ]

# 学習用入力データ
train.x <- data.matrix(train[1:4])
# 学習用ラベルデータ（0 ～ 2）
train.y <- as.numeric(train$Species) - 1

# 評価用入力データ
test.x <- data.matrix(test[1:4])
# 評価用ラベルデータ（1 ～ 3）
test.y <- as.numeric(test$Species)

mx.set.seed(0)

# 学習
model <- mx.mlp(train.x, train.y, 
                hidden_node = 5, 
                out_node = 3,
                num.round = 100,
                learning.rate = 0.1,
                array.batch.size = 10,
                activation = 'relu',
                array.layout = 'rowmajor',
                eval.metric = mx.metric.accuracy)

# 評価
pred <- predict(model, test.x, array.layout = 'rowmajor')

# 評価用データの分類結果（1 ～ 3）
pred.y <- max.col(t(pred))

# 評価データの正解率を算出
acc <- sum(pred.y == test.y) / length(pred.y)

print(acc)

実行結果は以下の通り。

実行結果
・・・
> model <- mx.mlp(train.x, train.y, 
+                 hidden_node = 5, 
+                 out_node = 3,
+                 num.round = 100,
+                 learning.rate = 0.1,
+                 array.batch.size = 10,
+                 activation = 'relu',
+                 array.layout = 'rowmajor',
+                 eval.metric = mx.metric.accuracy)
Start training with 1 devices
[1] Train-accuracy=0.32
[2] Train-accuracy=0.281818181818182
・・・
[99] Train-accuracy=0.954545454545455
[100] Train-accuracy=0.954545454545455
・・・
> print(acc)
[1] 0.9555556

備考

predict の実行結果は以下のような内容となっています。

> pred

            [,1]        [,2]        [,3]        [,4]
[1,] 0.968931615 0.968931615 0.968931615 0.968931615
[2,] 0.029328469 0.029328469 0.029328469 0.029328469
[3,] 0.001739914 0.001739914 0.001739914 0.001739914
            [,5]        [,6]        [,7]        [,8]
[1,] 0.968931615 0.968931615 0.968931615 0.968931615
[2,] 0.029328469 0.029328469 0.029328469 0.029328469
[3,] 0.001739914 0.001739914 0.001739914 0.001739914
・・・
            [,41]        [,42]        [,43]        [,44]
[1,] 1.762393e-08 7.670556e-06 5.799695e-06 9.349569e-12
[2,] 3.053433e-02 2.679898e-01 1.714197e-01 7.250102e-05
[3,] 9.694657e-01 7.320026e-01 8.285745e-01 9.999275e-01
            [,45]
[1,] 4.420018e-08
[2,] 8.569881e-03
[3,] 9.914301e-01

1 ～ 3 の中で最も数値の高いものが分類結果となりますので、上記を t で転置して max.col すると以下のようになります。

> max.col(t(pred))

 [1] 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 2 2 2 2 3 2 2 2 2 2
[30] 3 2 2 2 2 3 3 3 3 3 3 3 3 3 3 3

ジニ不純度の算出3 - Python, R, CoffeeScript

python R coffeescript nodejs

前々回 と前回 に続き、下記のようなプログラム言語でジニ不純度（ジニ係数）の算出処理を同様に実装してみました。

今回のソースは http://github.com/fits/try_samples/tree/master/blog/20140622/

Python で実装

Python では itertools モジュールの groupby や combinations 関数が使えます。

groupby は Haskell と同様に隣り合う同じ値をグルーピングできます。 （今回のケースでは sorted でソートが必要）

groupby 結果の値部分（グルーピング部分）には直接 len 関数を使えないようなので list 関数でリスト化してから len を適用します。

また、combinations 関数を使用すると Scala の combinations と同様に要素の組み合わせを取得できます。（下記では AB、AC、BC の 3種類）

gini.py
from itertools import *

def size(xs):
    return float(len(xs))

# (a) 1 - (AA + BB + CC)
def giniA(xs):
    return 1 - sum(map(lambda (k, g): (size(list(g)) / size(xs)) ** 2, groupby(sorted(xs))))

def countby(xs):
    return map(lambda (k, v): (k, size(list(v))), groupby(sorted(xs)))

# (b) AB * 2 + AC * 2 + BC * 2
def giniB(xs):
    return sum(map(
        lambda ((xk, xv), (yk, yv)): xv / size(xs) * yv / size(xs) * 2, 
        combinations(countby(xs), 2)
    ))

vlist = ["A", "B", "B", "C", "B", "A"]

print giniA(vlist)
print giniB(vlist)
実行結果
> python gini.py

0.611111111111
0.611111111111

Python 3 で実行するには

Python 3.4 で実行するにはラムダ式と print のところを書き換える必要があります。

Python 3 のラムダ式ではタプル内の複数の要素を個別に引数として取得できないようなので、Python 2 のように lambda (k, v): ・・・ とは書けず、lambda x: ・・・ として個々の要素をインデックスで参照 （x[0] 等） する事になります。

gini3.py （Python 3.4）
・・・
# (a) 1 - (AA + BB + CC)
def giniA(xs):
    return 1 - sum(map(lambda x: (size(list(x[1])) / size(xs)) ** 2, groupby(sorted(xs))))
・・・
# (b) AB * 2 + AC * 2 + BC * 2
def giniB(xs):
    return sum(map(
        lambda x: x[0][1] / size(xs) * x[1][1] / size(xs) * 2, 
        combinations(countby(xs), 2)
    ))

vlist = ["A", "B", "B", "C", "B", "A"]

print(giniA(vlist))
print(giniB(vlist))
実行結果 （Python 3.4）
> python gini3.py

0.6111111111111112
0.611111111111111

R で実装

R では table 関数で要素毎のカウント値を取得でき、combn 関数で Scala や Python の combinations と同様の組み合わせを行列（matrix）として取得できます。

lapply の結果（リスト）には sum 関数を適用できないようなので Reduce を使って合計しています。

また、apply の第2引数を 2 とすれば列単位にデータを処理できます。

gini.R
# (a) 1 - (AA + BB + CC)
giniA <- function(xs) {
  1 - Reduce("+", lapply(table(xs), function(x) (x / length(xs)) ^ 2))
}

# (b) AB * 2 + AC * 2 + BC * 2
giniB <- function(xs) {
  sum(apply(combn(table(xs), 2), 2, function(x) (x[1] / length(xs)) * (x[2] / length(xs)) * 2))
}

list <- c("A", "B", "B", "C", "B", "A")

giniA(list)
giniB(list)
実行結果
・・・
> giniA(list)
[1] 0.6111111

> giniB(list)
[1] 0.6111111

備考

各処理の結果は下記のようになります。

table(list) の結果
> table(list)

list
A B C 
2 3 1
combn(table(list), 2) の結果
> combn(table(list), 2)

     [,1] [,2] [,3]
[1,]    2    2    3
[2,]    3    1    1

ちなみに、上記は以下のような組み合わせのカウント値です。

     [,1] [,2] [,3]
[1,] "A"  "A"  "B" 
[2,] "B"  "C"  "C"

CoffeeScript で実装

CoffeeScript では、Underscore.js 等のライブラリを使用しない限り、グルーピングや組み合わせの処理を自前で実装する事になると思います。

gini.coffee
countBy = (xs) -> xs.reduce (acc, x) ->
    acc[x] ?= 0
    acc[x]++
    acc
, {}

sum = (xs) -> xs.reduce (acc, x) -> acc + x

# (a) 1 - (AA + BB + CC)
giniA = (xs) -> 1 - sum( (v / xs.length) ** 2 for k, v of countBy(xs) )

flatten = (xs) -> xs.reduce (x, y) -> x.concat y

combination = (xs) -> flatten( [x, y] for x in xs when x isnt y for y in xs )

# (b) BA + CA + AB + CB + AC + BC
giniB = (xs) -> sum( (x[1] / xs.length) * (y[1] / xs.length) for [x, y] in combination([k, v] for k, v of countBy(xs)) )

list = ["A", "B", "B", "C", "B", "A"]

console.log giniA(list)
console.log giniB(list)
実行結果
> coffee gini.coffee

0.6111111111111112
0.611111111111111

R でロジスティック回帰とオッズ比の算出 - glm, MCMClogit

R GLM Bayes

以前、glm・MCMCmetrop1R 関数でロジスティック回帰を試みましたが、今回はその時に利用を断念した MCMCpack の MCMClogit 関数を使ってロジスティック回帰を行います。

題材は、書籍 「 データサイエンティスト養成読本 [ビッグデータ時代のビジネスを支えるデータ分析力が身につく! ] (Software Design plus) 」の p.38 と同様の Titanic データセットを使ったロジスティック回帰とオッズ比の算出です。

以前使ったデータは主に数値でしたが、今回は主に因子（factor）データを扱っている点が異なっています。

今回のソースは http://github.com/fits/try_samples/tree/master/blog/20140302/

はじめに

MCMCpack パッケージを R へインストールしておきます。

install.packages("MCMCpack")

Titanic データセット

Titanic のデータセットは R にデフォルトで用意されており、data.frame(Titanic) すると下記のようになります。

data.frame(Titanic) 結果
> data.frame(Titanic)

   Class    Sex   Age Survived Freq
1    1st   Male Child       No    0
2    2nd   Male Child       No    0
3    3rd   Male Child       No   35
4   Crew   Male Child       No    0
5    1st Female Child       No    0
6    2nd Female Child       No    0
7    3rd Female Child       No   17
8   Crew Female Child       No    0
9    1st   Male Adult       No  118
10   2nd   Male Adult       No  154
11   3rd   Male Adult       No  387
12  Crew   Male Adult       No  670
13   1st Female Adult       No    4
14   2nd Female Adult       No   13
15   3rd Female Adult       No   89
16  Crew Female Adult       No    3
17   1st   Male Child      Yes    5
18   2nd   Male Child      Yes   11
19   3rd   Male Child      Yes   13
20  Crew   Male Child      Yes    0
・・・

内容は以下の通りです。

Class（船室等級） Sex（性別） Age（年齢層） Survived（生存可否） Freq（人数）
1st, 2nd, 3rd, Crew Male, Female Child, Adult No, Yes 数値

今回は、ロジスティック回帰の結果から船室等級・性別・年齢層毎の生存率に対するオッズ比を算出します。

(1) glm を使ったロジスティック回帰とオッズ比

書籍の内容とほとんど同じですが、 まずは glm 関数を使ったロジスティック回帰です。

Survived~. としているので Class・Sex・Age を説明変数としたロジスティック回帰を実施します。

オッズ比は exp(<推定値>) で算出できるので、書籍のような epicalc パッケージは使わず、glm 結果の $coefficients を exp 関数へ渡してオッズ比を算出しました。

logiMcmcglmm.R
d <- data.frame(Titanic)

d.data <- data.frame(
  Class = rep(d$Class, d$Freq),
  Sex = rep(d$Sex, d$Freq),
  Age = rep(d$Age, d$Freq),
  Survived = rep(d$Survived, d$Freq)
)

d.res <- glm(Survived~., data = d.data, family = binomial)
summary(d.res)

# オッズ比
exp(d.res$coefficients)

なお、上記では data.frame(Titanic) の Freq（人数） を展開したデータフレーム（下記）を使ってロジスティック回帰を実施しています。

d.data の内容
   Class    Sex   Age Survived
1    3rd   Male Child       No
2    3rd   Male Child       No
・・・
35   3rd   Male Child       No
36   3rd Female Child       No
・・・

実行結果は下記の通りです。

実行結果
・・・
> d.res <- glm(Survived~., data = d.data, family = binomial)
> summary(d.res)

Call:
glm(formula = Survived ~ ., family = binomial, data = d.data)

Deviance Residuals: 
    Min       1Q   Median       3Q      Max  
-2.0812  -0.7149  -0.6656   0.6858   2.1278  

Coefficients:
            Estimate Std. Error z value Pr(>|z|)    
(Intercept)   0.6853     0.2730   2.510   0.0121 *  
Class2nd     -1.0181     0.1960  -5.194 2.05e-07 ***
Class3rd     -1.7778     0.1716 -10.362  < 2e-16 ***
ClassCrew    -0.8577     0.1573  -5.451 5.00e-08 ***
SexFemale     2.4201     0.1404  17.236  < 2e-16 ***
AgeAdult     -1.0615     0.2440  -4.350 1.36e-05 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

(Dispersion parameter for binomial family taken to be 1)

    Null deviance: 2769.5  on 2200  degrees of freedom
Residual deviance: 2210.1  on 2195  degrees of freedom
AIC: 2222.1

Number of Fisher Scoring iterations: 4

> 
> # オッズ比
> exp(d.res$coefficients)
(Intercept)    Class2nd    Class3rd   ClassCrew   SexFemale    AgeAdult 
  1.9844057   0.3612825   0.1690159   0.4241466  11.2465380   0.3459219

オッズ比から、女性（Female）は男性（Male）の 11倍（11.2465380）、大人（Adult）は子供（Child）の 3分の1 （0.3459219）の生存率という結果になりました。

(2) MCMClogit を使ったロジスティック回帰とオッズ比

それでは、MCMCpack の MCMClogit 関数を使って同様の処理を実施してみます。

glm 関数の場合とほとんど同じですが、MCMClogit 関数は応答変数（下記の Survived）に因子型（factor）を扱えないようなので、Survived の No と Yes を 0（= No） と 1（= Yes）へ変換しています。 （as.numeric(d$Survived) - 1 の箇所）

また、MCMClogit 関数の結果も glm とは違って値の分布となるので summary した結果の平均値（期待値） Mean を使ってオッズ比を算出しています。

logistic_odds_mcmclogit.R
library(MCMCpack)

d <- data.frame(Titanic)

d.data <- data.frame(
  Class = rep(d$Class, d$Freq),
  Sex = rep(d$Sex, d$Freq),
  Age = rep(d$Age, d$Freq),
  Survived = rep(as.numeric(d$Survived) - 1, d$Freq)
)

d.res <- MCMClogit(Survived~., data = d.data)

d.summ <- summary(d.res)

d.summ

#オッズ比
exp(d.summ$statistics[, "Mean"])

実行結果は下記の通りです。

実行結果
・・・
> d.res <- MCMClogit(Survived~., data = d.data)
> 
> d.summ <- summary(d.res)
> 
> d.summ

Iterations = 1001:11000
Thinning interval = 1 
Number of chains = 1 
Sample size per chain = 10000 

1. Empirical mean and standard deviation for each variable,
   plus standard error of the mean:

               Mean     SD Naive SE Time-series SE
(Intercept)  0.6775 0.2682 0.002682       0.012059
Class2nd    -1.0367 0.1961 0.001961       0.008828
Class3rd    -1.7935 0.1768 0.001768       0.008054
ClassCrew   -0.8607 0.1596 0.001596       0.006971
SexFemale    2.4486 0.1384 0.001384       0.006279
AgeAdult    -1.0540 0.2392 0.002392       0.010882

2. Quantiles for each variable:

               2.5%     25%     50%     75%   97.5%
(Intercept)  0.1455  0.4992  0.6831  0.8630  1.2028
Class2nd    -1.4185 -1.1683 -1.0331 -0.9027 -0.6626
Class3rd    -2.1481 -1.9130 -1.7926 -1.6760 -1.4476
ClassCrew   -1.1764 -0.9661 -0.8651 -0.7535 -0.5488
SexFemale    2.1543  2.3600  2.4525  2.5442  2.7110
AgeAdult    -1.5275 -1.2144 -1.0573 -0.8904 -0.5907

> 
> #オッズ比
> exp(d.summ$statistics[, "Mean"])
(Intercept)    Class2nd    Class3rd   ClassCrew   SexFemale    AgeAdult 
  1.9690188   0.3546137   0.1663746   0.4228566  11.5715841   0.3485381

glm と同じような結果となりました。

最後に、それぞれの値の分布は下記のようになりました。

f:id:fits:20140302174258p:plain

f:id:fits:20140302174314p:plain

R で個体差のあるロジスティック回帰2 - MCMCglmm

R GLM Bayes

前回 に続き、今回は個体差を考慮したロジスティック回帰を MCMCglmm で試してみます。

実は MCMCglmm の他にも MCMCpack の MCMChlogit や bayesm の rhierBinLogit 等といろいろ試そうとしてみたのですが、イマイチ使い方が分からなかったので今回は断念しました。

今回使用したサンプルソースは http://github.com/fits/try_samples/tree/master/blog/20140211/

はじめに

今回使用するパッケージを R へインストールしておきます。

install.packages("MCMCglmm")

(2) MCMCglmm を使った階層ベイズモデルによるロジスティック回帰（MCMCglmm 関数）

それでは、前回 と同じデータ （data7.csv） を使って推定と予測線のグラフ化を実施してみます。

MCMCglmm() による推定

MCMCglmm() 関数は glmmML() 関数とほぼ同じ使い方ができますが、 今回のケースでは下記の点が異なります。

  • family に multinomial2 を指定する
  • デフォルトで個体差 （units に該当） が考慮されるようなので特に id 列を指定する必要はない

デフォルトで詳細（処理状況）が出力されるので、下記では出力しないよう verbose = FALSE としています。

logiMcmcglmm.R
d <- read.csv('data7.csv')

library(MCMCglmm)

d.res <- MCMCglmm(cbind(y, N - y) ~ x, data = d, family = "multinomial2", verbose = FALSE)

summary(d.res)
・・・

実行結果は下記のようになります。

実行結果
> d.res <- MCMCglmm(cbind(y, N - y) ~ x, data = d, family = "multinomial2", verbose = FALSE)
> 
> summary(d.res)

 Iterations = 3001:12991
 Thinning interval  = 10
 Sample size  = 1000 

 DIC: 667.8652 

 R-structure:  ~units

      post.mean l-95% CI u-95% CI eff.samp
units     6.988    3.898    10.58    517.4

 Location effects: cbind(y, N - y) ~ x 

            post.mean l-95% CI u-95% CI eff.samp  pMCMC    
(Intercept)   -4.2271  -6.3586  -2.5264    814.3 <0.001 ***
x              1.0144   0.6117   1.5237    800.3 <0.001 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

各パラメータの平均値 （post.mean） が { \bar \beta_1 = -4.2271, \bar \beta_2 = 1.0144, \bar s = \sqrt{6.988} = 2.643 } で、前回の結果（glmmML）と近い値になりました。

ここで units は前回の個体差パラメータ { r_i } の分散 { s ^2 } に該当するようですので、標準偏差とするために平方根を取りました。

プログラム上で上記の結果を使う場合は、MCMCglmm() の結果から { \beta_1, \beta_2 } （Intercept, x）の分布を $Sol で、{ s ^2 } （units）の分布を $VCV で取得できます。

なお、実行する度に若干異なる値になるようなので、何らかのオプションを指定して調整する必要があるのかもしれません。

葉数と生存種子数

MCMCglmm() の結果を使って 前回 と同様にグラフ化してみます。

今回はパラメータの平均値（post.mean の値）を使わずに、事後最頻値 （posterior mode） というものを使って予測線を描画してみました。

事後最頻値は posterior.mode() 関数で取得できます。

logiMcmcglmm.R
・・・
# 生存確率の算出
calcProb <- function(x, b, r)
    1.0 / (1.0 + exp(-1 * (b[1] + b[2] * x + r)))

png("logiMcmcglmm_1.png")

plot(d$x, d$y)

xx <- seq(min(d$x), max(d$x), length = 100)

# beta1 と beta2 の事後最頻値
beta <- posterior.mode(d.res$Sol)
# s の事後最頻値
sd <- sqrt(posterior.mode(d.res$VCV))

lines(xx, max(d$N) * calcProb(xx, beta, 0), col="green")
lines(xx, max(d$N) * calcProb(xx, beta, -1 * sd), col="blue")
lines(xx, max(d$N) * calcProb(xx, beta, sd), col="blue")

dev.off()
・・・

f:id:fits:20140211122143p:plain

ついでに、前々回 のようにpredict() を使って予測線を描画しようとしてみましたが、今のところ predict.MCMCglmm() は新しいデータの指定に対応していないようです。

> predict(d.res, data.frame(x = xx))

Error in predict.MCMCglmm(d.res, data.frame(x = xx)) : 
  sorry newdata not implemented yet

葉数 { x_i = 4 } の種子数分布

こちらも事後最頻値を使って、前回 と同様にグラフ化してみます。

logiMcmcglmm.R
・・・
png("logiMcmcglmm_2.png")

yy <- 0:max(d$N)

plot(yy, table(d[d$x == 4,]$y), xlab="y", ylab="num")

# 葉数 x を固定した場合の生存種子数 y の確率分布を算出
calcL <- function(ylist, xfix, n, b, s)
  sapply(ylist, function(y) integrate(
      f = function(r) dbinom(y, n, calcProb(xfix, b, r)) * dnorm(r, 0, s),
      lower = s * -10,
      upper = s * 10
    )$value
  )

lines(yy, calcL(yy, 4, max(d$N), beta, sd) * length(d[d$x == 4,]$y), col="red", type="b")

dev.off()

f:id:fits:20140211122201p:plain

R で個体差のあるロジスティック回帰1 - glmmML

R GLM Bayes

前回 のロジスティック回帰に続き、書籍 「 データ解析のための統計モデリング入門――一般化線形モデル・階層ベイズモデル・MCMC (確率と情報の科学) 」のサンプルを使って個体差を考慮したロジスティック回帰を GLMM と階層ベイズモデルで試してみます。

  • (1) GLMM によるロジスティック回帰 （glmmML 関数）
  • (2) MCMCglmm を使った階層ベイズモデルによるロジスティック回帰（MCMCglmm 関数）

ただし、今回は (1) だけで (2) に関しては次回に書く予定です。

サンプルソースは http://github.com/fits/try_samples/tree/master/blog/20140209/

はじめに

パッケージのインストール

今回使用するパッケージを R へインストールしておきます。

install.packages("glmmML")

データ

データは書籍のサポート web サイトから取得します。 今回は 7章「一般化線形混合モデル（GLMM）- 個体差のモデリング -」のサンプルデータを使います。

データ data7.csv
"N","y","x","id"
8,0,2,1
8,1,2,2
8,2,2,3
・・・

データ内容は下記の通りです。

項目 内容
N 調査種子数
y 生存種子数
x 植物の葉数
id 個体のID

葉数 { x_i } の値が 2 ～ 6 までとなっており、 葉数毎に 20 個体ずつの生存種子数 { y_i } をサンプリングしたデータとなっています。 （合計 100 個体）

二項分布の性質と過分散

ところで、二項分布の性質として下記があります。
今回のケースでは n が調査種子数（= 8）、p が生存確率です。

  • (1) 平均（期待値とするのが適切かも）が { np }
  • (2) 分散が { np(1 - p) }

ここで、今回のデータは下記のようになっています。

葉数 4（x = 4）場合の平均と分散
> mean(d[d$x == 4,]$y)
[1] 4.05

> var(d[d$x == 4,]$y)
[1] 8.365789

np = 4.05 とすると p = 4.05 / 8 = 約 0.5 となり、二項分布として期待される分散は np(1 - p) = 8 × 0.5 × (1 - 0.5) = 2 で、実際の分散 8.36 の方が明らかに大きい値となっています。 （過分散）

二項分布として考えると、原因不明な個体差等による効果を考慮した統計モデルを使ったロジスティック回帰が必要となり、GLMM や階層ベイズを使用するという流れとなるようです。

(1) GLMM によるロジスティック回帰 （glmmML 関数）

それでは glmmML() 関数を使った推定値の算出とグラフ化を実施してみます。

推定方法

個体差を考慮するため、GLMM では個体差のパラメータ { r_i } を追加した線形予測子 { logit(q_i) = \beta_1 + \beta_2 x_i + r_i } を使う事になるようです。
ここで、 { q_i } は生存確率です。

そして、{ r_i } が正規分布の確率分布に従っていると仮定し、
確率密度関数 { p(r_i | s) = \frac{1}{\sqrt{ 2 \pi s_^2} } \exp(- \frac{r_i ^2}{ 2s ^2 }) } を加味して { r_i } を積分した
尤度 { L_i = \int p(y_i | \beta_1, \beta_2, r_i) p(r_i | s) dr_i } に対して
全データの積 { L(\beta_1, \beta_2, s) = \prod_i L_i } の対数尤度 { \log L(\beta_1, \beta_2, s) } が最大になるようなパラメータ { \beta_1, \beta_2, s } の最尤推定値を求めます。

なお、{ p(y_i | \beta_1, \beta_2, r_i) } は二項分布の確率密度関数で、 { s } は個体差 { r_i } のばらつき（標準偏差）です。

glmmML() による推定

glmmML() 関数へ指定するオプションは、前回 の glm() とほとんど同じですが 、下記の点が異なります。

  • 個体差のパラメータ { \r_i } 部分を cluster オプションで指定

今回は個体毎に異なる値が設定されている id 列を cluster オプションへ指定します。

logiGlmmML.R
d <- read.csv('data7.csv')

library(glmmML)

d.res <- glmmML(cbind(y, N - y) ~ x, data = d, family = binomial, cluster = id)

summary(d.res)
・・・

glmmML() による推定結果は下記の通りです。

実行結果
Call:  glmmML(formula = cbind(y, N - y) ~ x, family = binomial, data = d, cluster = id) 


              coef se(coef)      z Pr(>|z|)
(Intercept) -4.190   0.8777 -4.774 1.81e-06
x            1.005   0.2075  4.843 1.28e-06

Scale parameter in mixing distribution:  2.408 gaussian 
Std. Error:                              0.2202 

        LR p-value for H_0: sigma = 0:  2.136e-55 

Residual deviance: 269.4 on 97 degrees of freedom   AIC: 275.4

{ \hat \beta_1 = -4.190, \hat \beta_2 = 1.005, \hat s = 2.408 } と推定されました。

なお、glmmML() の結果から、
{ \hat \beta_1, \hat \beta_2 } の値は $coefficients で、{ \hat s } の値は $sigma で取得できます。

葉数と生存種子数

推定したパラメータを使って葉数 x と生存種子数 y の予測線をグラフ化する処理は下記のようになります。

生存種子数の予測値は 調査種子数（= max(d$N) = 8）× 生存確率 で算出し、
生存確率 { q_i } は { logit(q_i) = \beta_1 + \beta_2 x_i + r_i } の式を元に算出します。

logiGlmmML.R
・・・
# 生存確率の算出
calcProb <- function(x, b, r)
    1.0 / (1.0 + exp(-1 * (b[1] + b[2] * x + r)))

png("logiGlmmML_1.png")

# x と y をプロット
plot(d$x, d$y)

# x の範囲（2 ～ 6）を 100分割
xx <- seq(min(d$x), max(d$x), length = 100)

beta <- d.res$coefficients

# 個体差 r = 0 の葉数と生存種子数との予測線
lines(xx, max(d$N) * calcProb(xx, beta, 0), col="green")
# 個体差 r = -2.408 の葉数と生存種子数との予測線
lines(xx, max(d$N) * calcProb(xx, beta, -1 * d.res$sigma), col="blue")
# 個体差 r = 2.408 の葉数と生存種子数との予測線
lines(xx, max(d$N) * calcProb(xx, beta, d.res$sigma), col="blue")

dev.off()
・・・

個体差 r = 0 （個体差なし）とした場合の予測線を緑色で、 r = ±2.408 の予測線を青色で描画しています。

f:id:fits:20140209172056p:plain

葉数 { x_i = 4 } の生存種子数分布

次に、葉数 x が 4 の場合の生存種子数 y と個体数の分布と予測線をグラフ化する処理は下記のようになります。

glmmML() の結果から、x = 4 の場合の y の確率分布を算出し、x = 4 のサンプル数（= 20）を乗算して予測線を描画しています。

y の確率分布は、0 ～ 8 のそれぞれの値に対して { \int p(y_i | \beta_1, \beta_2, r_i) p(r_i | s) dr_i } で算出します。 （二項分布 × 正規分布 を個体差 { r_i } で積分）

下記では sapply() と integrate() 関数を使って上記の算出処理を実装しています。

なお、lower と upper オプションを使って積分範囲を { -10 \times \hat s } から { 10 \times \hat s } としています。

logiGlmmML.R
・・・
png("logiGlmmML_2.png")

# 0 ～ 8
yy <- 0:max(d$N)

# 生存種子数 y と個体数の分布
plot(yy, table(d[d$x == 4,]$y), xlab="y", ylab="num")

# 葉数 x を固定した場合の生存種子数 y の確率分布を算出
calcL <- function(ylist, xfix, n, b, s)
  sapply(ylist, function(y) integrate(
      f = function(r) dbinom(y, n, calcProb(xfix, b, r)) * dnorm(r, 0, s),
      lower = s * -10,
      upper = s * 10
    )$value
  )

# 予測線
lines(yy, calcL(yy, 4, max(d$N), beta, d.res$sigma) * length(d[d$x == 4,]$y), col="red", type="b")

dev.off()

f:id:fits:20140209172110p:plain

今回はここまで。 MCMCglmm を使った階層ベイズは次回。

R でロジスティック回帰 - glm, MCMCpack

R GLM Bayes

前回 に続き、今回も書籍 「 データ解析のための統計モデリング入門――一般化線形モデル・階層ベイズモデル・MCMC (確率と情報の科学) 」のサンプルを使って GLM とベイズ統計を試してみます。

題材は、6章「GLM の応用範囲をひろげる -ロジスティック回帰など-」のサンプルデータを使ったロジスティック回帰です。

  • (1) GLM によるロジスティック回帰 （glm 関数）
  • (2) MCMCpack を使ったベイズ統計によるロジスティック回帰（MCMCmetrop1R 関数）

ここで、MCMCpack のロジスティック回帰に MCMClogit 関数を使えると思ったのですが、使い方がよく分からなかったので MCMCmetrop1R 関数のみ利用しました。

サンプルソースは http://github.com/fits/try_samples/tree/master/blog/20131229/

はじめに

MCMCpack パッケージを install.packages("MCMCpack") で R の実行環境へインストールしておきます。

ロジスティック回帰を試すためのデータは書籍のサポート web サイトから取得できます。

データ data4a.csv
N,y,x,f
8,1,9.76,C
8,6,10.48,C
8,5,10.83,C
・・・

データ内容は下記のようになっており、個体 i それぞれにおいて 「 { N_i } 個の観察種子のうち生きていて発芽能力があるものは { y_i } 個、死んだ種子は { N_i - y_i } 個」 となっています。

項目 内容
N 観察種子数
y 生存種子数
x 植物の体サイズ
f 施肥処理 （C: 肥料なし, T: 肥料あり）

体サイズ x と肥料による施肥処理 f が種子の生存する確率（ある個体 i から得られた 1個の種子が生存している確率）にどのように影響しているかをロジスティック回帰で解析します。

(1) GLM によるロジスティック回帰 （glm 関数）

まず、glm 関数を使ったロジスティック回帰を行って予測値の曲線を描画します。

特徴は以下のようになります。

  • family に binomial （二項分布）を指定 （デフォルトのリンク関数が logit となる）
  • family が binomial で応答変数 （cbind(y, N - y) の箇所） が 2値（0, 1）以外なら cbind する必要あり （2値なら y だけでよい）

下記では、生存する確率 { q_i } の関数が { q_i = \frac{1}{1 + \exp(-z_i)} } 、
線形予測子 { z_i = \beta_1 + \beta_2 x_i + \beta_3 f_i } で解析する事になります。

logisticGlm.R
d <- read.csv('data4a.csv')

d.res <- glm(cbind(y, N - y) ~ x + f, data = d, family = binomial)

summary(d.res)

png("logisticGlm.png")
# C: 肥料なしを 赤、T: 肥料ありを 青 で描画
plot(d$x, d$y, col = c("red", "blue")[d$f])

xx <- seq(min(d$x), max(d$x), length = 1000)
ft <- factor("T", levels = levels(d$f))
fc <- factor("C", levels = levels(d$f))
# 下記でも可
#ft <- factor("T", levels = c("C", "T"))
#fc <- factor("C", levels = c("C", "T"))

# T: 肥料ありの予測値を算出
qq.t <- predict(d.res, data.frame(x = xx, f = ft), type="response")
# C: 肥料なしの予測値を算出
qq.c <- predict(d.res, data.frame(x = xx, f = fc), type="response")

# T: 肥料ありの予測値の曲線を 緑 で描画
lines(xx, max(d$N) * qq.t, col = "green")
# C: 肥料なしの予測値の曲線を 黄 で描画
lines(xx, max(d$N) * qq.c, col = "yellow")

dev.off()

predict 関数を使えば glm 結果の回帰モデルから新しいデータによる予測値を求めることができますので、"T:肥料あり" と "C:肥料なし" のそれぞれの予測値を求めて曲線を描画しています。

predict の結果は、生存する確率 { q_i } の予測値ですので、縦軸が生存種子数 { y_i } のグラフへ描画するには観察種子数 N （ここでは 8）を乗算する事になります。

ここで、"T:肥料あり" と "C:肥料なし" は factor 関数を使ってファクタ（因子）として作成しています。 d$f と同じ水準を持たせる必要がありますので levels 指定しています。

また、predict には glm で使った説明変数と同じ名前（上記の x と f）を使ったデータを渡す点に注意が必要です。

実行

実行すると下記のような結果になりました。

> R CMD BATCH logisticGlm.R
実行結果 logisticGlm.Rout
・・・
> d.res <- glm(cbind(y, N - y) ~ x + f, data = d, family = binomial)
> 
> summary(d.res)

Call:
glm(formula = cbind(y, N - y) ~ x + f, family = binomial, data = d)

Deviance Residuals: 
    Min       1Q   Median       3Q      Max  
-3.2584  -0.7948   0.1753   0.8757   3.1589  

Coefficients:
            Estimate Std. Error z value Pr(>|z|)    
(Intercept) -19.5361     1.4138  -13.82   <2e-16 ***
x             1.9524     0.1389   14.06   <2e-16 ***
fT            2.0215     0.2313    8.74   <2e-16 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

(Dispersion parameter for binomial family taken to be 1)

    Null deviance: 499.23  on 99  degrees of freedom
Residual deviance: 123.03  on 97  degrees of freedom
AIC: 272.21

Number of Fisher Scoring iterations: 5
・・・

この結果より、線形予測子 { z_i = \beta_1 + \beta_2 x_i + \beta_3 f_i } のパラメータ推定値 { { \hat\beta_1, \hat\beta_2, \hat\beta_3 } } が -19.5361, 1.9524, 2.0215 となります。

f:id:fits:20131229184250p:plain

(2) MCMCpack を使ったベイズ統計によるロジスティック回帰（MCMCmetrop1R 関数）

次は、MCMCpack の MCMCmetrop1R 関数を使ったロジスティック回帰です。

今回のモデルの対数尤度は { \log L ( \beta_j ) = \sum_i \log( \small N_i \large C \small y_i \normalsize ) + y_i \log(q_i) + (N_i - y_i) \log(1 - q_i) } で、
{ q_i = \frac{1}{1 + \exp(-z_i)} } 、
{ z_i = \beta_1 + \beta_2 x_i + \beta_3 f_i } なので、
これらを関数化 （下記の func） して MCMCmetrop1R 関数へ渡しています。

ちなみに、{ \small N_i \large C \small y_i \normalsize = \frac{N_i !}{y_i ! (N_i - y_i) ! } } です。（combination）

前回は、項目毎のデータ（d$x や d$y）を func 関数へ渡すようにしていましたが、今回はデータ d を丸ごと渡すようにしました。 MCMCmetrop1R(・・・, data = d, ・・・)

ファクタ（因子）のデータ （d$f） を直接計算に使えないようなので as.numeric で数値化して使っています。

なお、{ \log( \small N_i \large C \small y_i \normalsize ) + y_i \log(q_i) + (N_i - y_i) \log(1 - q_i) } の部分は dbinom と log 関数を使うとシンプルになります。

パラメータ { \beta_1 } { \beta_2 } { \beta_3 } の初期値は theta.init で c(0, 0, 0) と指定しました。

logisticMcmcMetrop.R
library(MCMCpack)

d <- read.csv('data4a.csv')

func <- function(beta, data) {
  z <- beta[1] + beta[2] * data$x + beta[3] * as.numeric(data$f)
  q <- 1.0 / (1.0 + exp(-z))
  sum(log(choose(data$N, data$y)) + data$y * log(q) + (data$N - data$y) * log(1 - q))

  # 下記でも可
  # sum(log(dbinom(data$y, data$N, q)))
}

d.res <- MCMCmetrop1R(func, theta.init = c(0, 0, 0), data = d, logfun = TRUE)

summary(d.res)

png("logisticMcmcMetrop.png")

plot(d$x, d$y, col = c("red", "blue")[d$f])

xx <- seq(min(d$x), max(d$x), length = 1000)
ft <- factor("T", levels = levels(d$f))
fc <- factor("C", levels = levels(d$f))

zt <- mean(d.res[,1]) + mean(d.res[,2]) * xx + mean(d.res[,3]) * as.numeric(ft)
zc <- mean(d.res[,1]) + mean(d.res[,2]) * xx + mean(d.res[,3]) * as.numeric(fc)

lines(xx, max(d$N) * 1.0 / (1.0 + exp(-zt)), col="green")
lines(xx, max(d$N) * 1.0 / (1.0 + exp(-zc)), col="yellow")

dev.off()

MCMCmetrop1R 結果の各パラメータの mean の値を使って、
肥料の有無（T と C）の生存確率の予測値をそれぞれ算出して、
{ \frac{1}{1 + \exp(-(\hat\beta_1 + \hat\beta_2 x_i + \hat\beta_3 f_i))} }
観察種子数 N （ここでは 8）を乗算し曲線を描画しています。

実行

実行すると下記のような結果になります。

> R CMD BATCH logisticMcmcMetrop.R
実行結果 logisticMcmcMetrop.Rout
・・・
@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
The Metropolis acceptance rate was 0.44590
@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
> 
> summary(d.res)

Iterations = 501:20500
Thinning interval = 1 
Number of chains = 1 
Sample size per chain = 20000 

1. Empirical mean and standard deviation for each variable,
   plus standard error of the mean:

        Mean     SD  Naive SE Time-series SE
[1,] -21.694 1.5583 0.0110185       0.037186
[2,]   1.966 0.1386 0.0009798       0.003322
[3,]   2.027 0.2268 0.0016035       0.005383

2. Quantiles for each variable:

        2.5%     25%     50%     75%   97.5%
var1 -24.836 -22.733 -21.667 -20.596 -18.766
var2   1.703   1.870   1.963   2.059   2.248
var3   1.594   1.873   2.022   2.172   2.494
・・・

Mean の値が glm の結果と似たような値になりました。 グラフにすると違いが分からないくらいです。

f:id:fits:20131229184326p:plain

{ \beta_1 } の分布

f:id:fits:20131229184348p:plain

{ \beta_2 } の分布

f:id:fits:20131229184713p:plain

{ \beta_3 } の分布

f:id:fits:20131229184417p:plain