Project Euler Problem 1:10

久しぶりにProject Eulerをやってみようと思ったが、せっかくなので進捗をリセットして最初からにした。
とりあえず10問やってみたが大分脳がヤバくなってる感あったのでちゃんと継続してボケ防止に努めたい。
解答コードはgithubで管理することにして、ここには要点をメモっていく。

Problem 1: Multiples of 3 and 5

簡易版FizzBuzz。添字操作と条件式の組み合わせの基本確認。

Problem 2: Even Fibonacci numbers

RでProject Eulerをやる場合、gmpパッケージを縛るか縛らないかで難易度に雲泥の差が生ずるが、基本的には使っていく方針で進めていきたい。便利そうなパッケージがあれば積極的に利用していく。
フィボナッチ数はgmp::fibnum()で求められる。

Problem 3: Largest prime factor

素因数分解はgmp::factorize()で実行できる。

Problem 4: Largest palindrome product

文字列の反転はstringi::stri_reverse()で実行できる。
組み合わせ毎に調べたほうが効率が良いが、探索範囲が狭いので999*999からスタートして下っていく方針でもすぐに終わる。
クリアすると読めるPDFファイルで勉強するための問題だろう。

Problem 5: Smallest multiple

大きな階乗はgmp::factorialZ()を使うと良い。
Big Integerに演算を施す場合は専用の関数が必要になる場合がしばしばある。剰余を求めるにはgmp::mod.bigz()を使う。

Problem 6: Sum square difference

ベクトル演算が出来るRでは基本操作レベル。

Problem 7: 10001st prime

gmpにはgmp::nextprime()という便利な関数がある。

Problem 8: Largest product in a series

文字列の操作方法が問題となる。

• 置換: sub(最初の1つだけ)、gsub(全て)
• 切り出し: substring
• 分割: strsplit(結果はlistなのでunlistする)

Problem 9: Special Pythagorean triplet

多重ループを抜けるにはlabeledLoopパッケージによるラベル付きループが便利だ。
探索範囲が狭いのでbrute forceでもすぐに終わってしまうのだが、これもProblem 4と同様、クリア後にPDFファイルで学習する用に簡単になっているのだろう。

Problem 10: Summation of primes

エラトステネスの篩は覚えるまで何度か自分で実装してみたほうが良いが、numbersというパッケージにエラトステネスの篩が入っている。numbers::Primes()は2つの引数を与えるとその数以下の、2つの値を与えると2数間の素数を列挙する。2つの値を与える場合に2数が20万くらい以上離れているとエラーが出るようだが、一つだけなら2^53-1まで指定できるようだ(あまり大きいと終わらなくなるが)

「データ解析のための統計モデリング入門」第11章メモ

p.254の図11.7をstanを使って再現したかった。

stanコードはこんな感じ。

data {
   int<lower=1> N;     // 欠測値を含むデータ数
   int<lower=1> N_obs; // 欠測値以外のデータの数
   int<lower=1, upper=N> Idx_obs[N_obs]; // 欠測値以外のデータのインデックス
   int<lower=0> Y_obs[N_obs]; // 欠測値を除いたデータベクトル
}

parameters {
   real r[N]; 
   real<lower=0> s_r; 
}

model {
   for (j in 2:N)
      r[j] ~ normal(r[j-1], s_r); // 1次のCAR model
   for (j in 1:N_obs) 
      Y_obs[j] ~ poisson_log(r[Idx_obs[j]]);
}

generated quantities { 
   real<lower=0> Y_mean[N];
   for (j in 1:N)
      Y_mean[j] = exp(r[j]);
}

R側のコードはこんな感じ。

# データを欠測させる
no <- c(6, 9, 12, 13, 26:30)
Idx_obs <- (1:50)[-no]
Y_obs <- Y[Idx_obs]
N <- length(Y)
N_obs <- length(Y_obs)

# サンプリング
dlist <- list(N = N, N_obs = N_obs, Idx_obs = Idx_obs, Y_obs = Y_obs)
fit2 <- stan("sec11_2.stan", data = dlist)

# 結果から中央値と信頼区間を算出
pred <- extract(fit2)
ci <- apply(pred$Y_mean, 2, quantile, probs = c(0.5, 0.025, 0.975))

# プロット
plot(Y, type = "n", ylim = c(0, 25))
polygon(c(1:N, N:1), 
        c(ci[2, ], rev(ci[3, ])),
        border = NA,
        col = "gray90")
points(Y)
points(no, Y[no], pch = 19)
lines(m, lty = 2)
lines(ci[1, ])
for(n in no){
  polygon(c(n-0.5, n+0.5, n+0.5, n-0.5),
          c(0, 0, 100, 100), 
          border = NA, 
          col = rgb(.3, .3, .3, .3))
}

するとこのようなものが出来る。

空間相関を考慮しないバージョンはこんな感じ。

data {
   int<lower=1> N;     // 欠測値を含むデータ数
   int<lower=1> N_obs; // 欠測値以外のデータの数
   int<lower=1, upper=N> Idx_obs[N_obs]; // 欠測値以外のデータのインデックス
   int<lower=0> Y_obs[N_obs]; // 欠測値を除いたデータベクトル
}

parameters {
   real r[N]; 
   real<lower=0> s_r; 
   real mu;
}

model {
   for (j in 1:N)
      r[j] ~ normal(0, s_r);
   for (j in 1:N_obs) 
      Y_obs[j] ~ poisson_log(mu + r[Idx_obs[j]]);
}

generated quantities {
   real<lower=0> Y_mean[N];
   for (j in 1:N)
      Y_mean[j] = exp(mu + r[j]);
}

参考

stanコードの記述にあたっては下記サイトが参考になったというかほとんどそのまんまです。

• 久保緑本11章のCAR model（空間構造のあるベイズモデル） - StatModeling Memorandum
• 不等間隔の状態空間モデル - StatModeling Memorandum

「データ解析のための統計モデリング入門」第9章メモつづき

前回。
rion778.hatenablog.com
これをrstanarmでやってみた。

準備

install.packages("rstanarm")

割と時間かかった。

library(rstan)
library(rstanarm)
rstan_options(auto_write = TRUE)
options(mc.cores = parallel::detectCores())

# データダウンロード
download.file(destfile = "ch9.Rdata", "http://hosho.ees.hokudai.ac.jp/~kubo/stat/iwanamibook/fig/gibbs/d.RData")
load("ch9.Rdata")

実行

fit <- stan_glm(y ~ x, family = poisson(), data = d)

細かい設定は省略したけど1行で済むのは楽。あと早い。

プロット

rstanに入っている作図関数がそのまま使える。

plot(fit)
stan_dens(fit, separate_chains = TRUE) + geom_rug()
stan_trace(fit)

出力は特に違いがないので略。
stan_traceにinclude_warmup=TRUEを指定した時の挙動が若干おかしい気はする(thinが無視されている様子？)。

出力

デフォルトの出力はかなり要約されている。summaryで少し詳しい物が見られるが、$stanfitの中にrstanの時と同じ形式のやつが入っていたので、例えばextract()を使いたい時なんかはこれを使えば良い様だ。

fit                # デフォルト
summary(fit)       # 少し詳しく
fit$stanfit        # rstanと同じで

テキストの図を再現する

上述のように、$stanfitをextract()してやれば良い。変数名が何故かalphaとbetaだ。

post <- extract(fit$stanfit)
add.mean <- function(bb1, bb2, lty = 2, lwd = 1, ...){
  lines(d$x, exp(bb1 + bb2 * d$x), lty = lty, lwd = lwd, ...)
}
plot(d$x, d$y, type = "n", xlab = "x_i", ylab = "y_i")
for (i in 1:length(post$alpha)) {
  add.mean(post$alpha[i], post$beta[i], lty = 1, col = "#00000004")
}
points(d$x, d$y)
add.mean(median(post$alpha), median(post$beta), lty = 1, lwd = 2)

これも出力に差は無いので略。

geom_barにおけるNAの扱い

geom_barというかstat_countの方だと思うけど、例えばこんなデータがあったとして

data <- c(1, 2, 2, 3, 3, 3, NA, NA, NA, NA)

普通に棒グラフを描くとこうだ。

ggplot() + geom_bar(aes(data))

ところがfactorにするとこうだ。

ggplot() + geom_bar(aes(factor(data)))

どうも集計にnaを含むかどうかはfactorであるかどうかに依存しているらしい(正しい動作かどうかはよく分からない)。ちなみにgeom_barとstat_countにはna.rm引数が用意されているが、これは警告の有無を制御するだけでNAが除外されることには変わりがない。
これはLearnBayesパッケージに含まれるstudentdataのDrinkをプロットしようとして気付いたのだが、もしこのようなデータをstat_countに渡す必要が出てきたら、factorで無くするか、NAを除外する必要がある。逆にNAを含めたかったらfactorにしてやる。

install.packages("LearnBayes")
library(LearnBayes)
data(studentdata)
## NAを除外する方法
ggplot(subset(studentdata, !is.na(Drink)), aes(x = Drink)) + geom_bar()
## characterに変換する方法
ggplot(studentdata, aes(x = as.character(Drink))) + geom_bar()

「データ解析のための統計モデリング入門」第9章メモ

緑本9章の内容について、RStanでやってみる。

Stanコード

次の内容で作成し、sample.stanという名前で保存した。

data {
  int<lower=1> N;
  vector[N] x;
  int<lower=0> y[N];
  real mean_x;
}
parameters {
  real beta1;
  real beta2;
}
model {
  for(i in 1:N){
    y[i] ~ poisson(exp(beta1 + beta2 * (x[i] - mean_x)));
  }
}

コードはPyStanによるはじめてのマルコフ連鎖モンテカルロ法 - Qiitaを参考に少しだけいじった。やはり平均値は引いた方が少し早いようだ。Stanでは変数名にドットを使用できないので注意が必要。

Rコード

library(rstan)

rstan_options(auto_write = TRUE)
options(mc.cores = parallel::detectCores())

# データのダウンロード
download.file(destfile = "ch9.Rdata", "http://hosho.ees.hokudai.ac.jp/~kubo/stat/iwanamibook/fig/gibbs/d.RData")
load("ch9.Rdata")

# データlistの準備
list.data <- list(
  x = d$x,
  y = d$y,
  N = nrow(d),
  mean_x = mean(d$x)
)

# stanで計算
fit <- stan(file = "sample.stan", data = list.data,
            iter = 1500, chains = 3, thin = 3, warmup = 100)

緑本メモ - 盆栽日記を参考にした。

結果の確認

> print(fit)
Inference for Stan model: sample.
3 chains, each with iter=1500; warmup=100; thin=3; 
post-warmup draws per chain=467, total post-warmup draws=1401.

        mean se_mean   sd   2.5%    25%    50%    75%  97.5%
beta1   1.97    0.00 0.08   1.81   1.92   1.98   2.03   2.13
beta2   0.08    0.00 0.07  -0.06   0.03   0.08   0.13   0.22
lp__  143.99    0.03 1.02 141.56 143.62 144.31 144.70 144.95
      n_eff Rhat
beta1  1287    1
beta2  1263    1
lp__   1315    1

Samples were drawn using NUTS(diag_e) at Mon Aug 15 00:45:10 2016.
For each parameter, n_eff is a crude measure of effective sample size,
and Rhat is the potential scale reduction factor on split chains (at 
convergence, Rhat=1).

plot(fit)

結果のプロット

トレースプロットはrstanパッケージのtraceplot()で描ける。バーンイン期間の表示の有無等も設定できるようだ。

traceplot(fit, inc_warmup = TRUE)

また、サンプリング結果はextract()で数値として取り出せるので、これを用いることで事後分布のプロット等も描ける。

post <- extract(fit)

old <- par(mfrow = c(1, 2))
plot(density(post$beta1), main = "beta 1")
rug(post$beta1)
plot(density(post$beta2), main = "beta 2")
rug(post$beta2)
par(old)

(※2016/08/17追記)
単に事後分布の密度推定をプロットしたいだけであればstan_dens()を使用するのが簡単だった。chain毎の推定結果を重ねて表示するオプションなんかもある。可視化のためのツールは一通り用意されている(?stan_plotを確認するとよい)。

stan_dens(fit1, separate_chains  = TRUE)

(※追記ここまで)

テキストの図を再現する

テキストの図9.6については、生態学データ解析 - 本/データ解析のための統計モデリング入門のfig09_06.Rを少し修正すれば描くことができる。

# テキストの図を再現
add.mean <- function(bb1, bb2, lty = 2, lwd = 1, ...){
  lines(d$x, exp(bb1 + bb2 * (d$x - mean(d$x))), lty = lty, lwd = lwd, ...)
}

plot(d$x, d$y, type = "n", xlab = "x_i", ylab = "y_i")
for (i in 1:length(post$beta1)) {
  add.mean(post$beta1[i], post$beta2[i], lty = 1, col = "#00000004")
}
points(d$x, d$y)
add.mean(median(post$beta1), median(post$beta2), lty = 1, lwd = 2)

WindowsにRStanをインストール

基本的にはRStan Getting Started (Japanese) · stan-dev/rstan Wiki · GitHubに従っていけた。

環境

• Windows 10 Pro
• R: 3.3.1
• RStudio: 0.99.903

Rtoolsのインストール

• Building R for WindowsよりRtools34.exeをダウンロード。
• セットアップ中、"Edit the system PATH"にチェック。

Makevarsファイルを作成するため、次のコードをRから実行。

dotR <- file.path(Sys.getenv("HOME"), ".R")
if (!file.exists(dotR)) dir.create(dotR)
M <- file.path(dotR, "Makevars")
if (!file.exists(M)) file.create(M)
cat("\nCXXFLAGS=-O3 -mtune=native -march=native -Wno-unused-variable -Wno-unused-function",
    file = M, sep = "\n", append = TRUE)

RStanのインストール

以下のコードをRから実行。

Sys.setenv(MAKEFLAGS = "-j4")
install.packages('rstan', repos = 'https://cloud.r-project.org/', dependencies = TRUE)

インストール完了後、Rを再起動。

C++が動くかどうかの確認

fx <- inline::cxxfunction( signature(x = "integer", y = "numeric" ) , '
  return ScalarReal( INTEGER(x)[0] * REAL(y)[0] ) ;
                           ')
fx( 2L, 5 )
# => 10

最初"cygheap base mismatch detected"が出たので、cygwin1.dllで検索してRtoolsに入っていたもの以外をリネームした(gnupackに入っていた)。
また、「もしg++のバージョンが4.9かそれ以降ならば...」に続くコードを実行してしまうと"Cannot export __gnu_lto_v1: symbol not defined"というようなエラーが出た。Makevarsファイルから該当箇所を削除したら直った(cf.Google グループ)

カテゴリカルデータ解析読書メモ(第10章)

決定木

データ全体を説明変数を用いて段階的にグループ分けする手法。

決定木を書く

テキストではmvpartパッケージを使用しているが、現在は公開が停止されているので代わりにpartykitパッケージを使用する。
rpart()の代わりに使う関数はctree()であるが、引数の入れ方はほとんど変わらない。結果はpartyオブジェクトになっているので、そのままプロットすれば決定木が表示される。

library(partykit)
library(vcd)
result <- ctree(survival ~ stage + operation + xray, weight = Freq,
                data = OvaryCancer)
result
# Model formula:
# survival ~ stage + operation + xray
#
# Fitted party:
# [1] root
# |   [2] stage in early: yes (n = 158, err = 19.6%)
# |   [3] stage in advanced: no (n = 141, err = 16.3%)
# 
# Number of inner nodes:    1
# Number of terminal nodes: 2
plot(result)

量的変数を用いる

ctree()を使った結果はテキストと若干異なる。

# 量的変数を用いたグループ分け
data(iris)
result <- ctree(Species ~ ., data = iris)
result
# Model formula:
#   Species ~ Sepal.Length + Sepal.Width + Petal.Length + Petal.Width
# 
# Fitted party:
#   [1] root
# |   [2] Petal.Length <= 1.9: setosa (n = 50, err = 0.0%)
# |   [3] Petal.Length > 1.9
# |   |   [4] Petal.Width <= 1.7
# |   |   |   [5] Petal.Length <= 4.8: versicolor (n = 46, err = 2.2%)
# |   |   |   [6] Petal.Length > 4.8: versicolor (n = 8, err = 50.0%)
# |   |   [7] Petal.Width > 1.7: virginica (n = 46, err = 2.2%)
# 
# Number of inner nodes:    3
# Number of terminal nodes: 4
plot(result)

また、predict()にtypeオプションの指定は必要ない。

iris$est <- predict(result)
xtabs( ~ Species + est, data = iris)
#             est
# Species      setosa versicolor virginica
#   setosa         50          0         0
#   versicolor      0         49         1
#   virginica       0          5        45

順序カテゴリカルデータ

パッケージrpartOrdinalも既に公開されていないが、これもctree()で解析できる。目的変数が順序カテゴリカルデータになっていれば特にオプションの指定は必要ない様だ。
データセットはパッケージRSADBEに同様の物が入っているのでこれを利用する。

library(RSADBE)
data("lowbwt")
lowbwt$Category <- factor(ifelse(lowbwt$BWT <= 2500, 3,
                                 ifelse(lowbwt$BWT <= 3000, 2,
                                        ifelse(lowbwt$BWT <= 3500, 1, 0))),
                          ordered = TRUE)
ord.rpart <- ctree(Category ~ AGE + LWT + RACE + SMOKE + PTL + HT + UI + FTV,
                   data = lowbwt)
ord.rpart
 
# Model formula:
#   Category ~ AGE + LWT + RACE + SMOKE + PTL + HT + UI + FTV
# 
# Fitted party:
#   [1] root
# |   [2] LWT <= 109: 3 (n = 42, err = 50.0%)
# |   [3] LWT > 109
# |   |   [4] UI <= 0
# |   |   |   [5] SMOKE <= 0
# |   |   |   |   [6] RACE <= 1: 0 (n = 39, err = 43.6%)
# |   |   |   |   [7] RACE > 1: 1 (n = 44, err = 63.6%)
# |   |   |   [8] SMOKE > 0: 3 (n = 47, err = 66.0%)
# |   |   [9] UI > 0: 3 (n = 17, err = 41.2%)
# 
# Number of inner nodes:    4
# Number of terminal nodes: 5

plot(ord.rpart)

練習問題解答

カテゴリカルデータ解析第10章練習問題解答 · GitHub

Rで逆推定(2)

はてなから1年前のブログとかいうメールが来てこんな記事を書いていたのを思い出した。
rion778.hatenablog.com
再び少し調べたところ、これはinvestrというパッケージで実行できるらしい。

データとモデルオブジェクトの作成

データは前回と同じものを使うが、データフレームにしておく。これは、モデルオブジェクト中の$dataがNULLだと逆推定を行う際にデータの指定が必要となるためだ。

temp <- rep(c(15, 20, 25), c(5, 5, 5))
day <- c(34, 33, 36, 37, 35,
         30, 28, 29, 25, 28,
         23, 25, 22, 24, 26)
d <- data.frame(temp, day)
result <- lm(day ~ temp, data = d)

逆推定

逆推定はinvest()関数で行う。さらに、plotFit()関数で推定区間と共にプロットできる。

library(investr)
(inv <- invest(result, y0 = 30, interval = "inversion", mean.response = TRUE))
# estimate    lower    upper 
# 19.09091 18.09070 19.99693 
plotFit(result, interval = "confidence", shade = TRUE)
abline(h = 30, v = c(inv$lower, inv$estimate, inv$upper), lty = 2)

ちなみに、JMPで実行した場合の推定値は

• 予測値: 19.09091
• 下側95%:18.09069
• 上側95%:19.99694

であるので、ほぼ一致している。

参考資料

• CRAN - Package investr
• investr: An R Package for Inverse Estimation

カテゴリカルデータ解析読書メモ(第9章)

対応分析

カテゴリカルデータにおいて各カテゴリーの変数に適当な数字を割り当てると、2つのカテゴリ間で散布図を書いたり相関係数を計算することができる。このとき、相関係数が最大となるような数値をカテゴリー変数に割り当てることで、カテゴリー間の関係を表現、探索する方法が対応分析である。
まず、カテゴリー変数に適当に数字を割り当ててプロットしてみる。テキストでは円の大きさを度数の対数に対応させているが、分かりにくいので普通に度数に比例させた。

# カテゴリーに数値を割り当て、度数に円の大きさを対応させてプロットする
plot(c(0,6), c(0,4), type = "n")
for(j in 1:5){
  for(i in 1:3){
    points(j,i,cex = tabledata[i,j]/mean(tabledata))
  }
}

この状態では行と列の間に相関はほとんど見られない。
相関係数の最大化は、MASSパッケージのcorresp()関数で実行できる。

# corresp()を使って相関係数を最大化する
result <- corresp(tabledata, nf = 1)
result
plot(-3:2, -3:2, type = "n", xlab = "Urban", ylab = "Status")
for(i in 1:3){
  for(j in 1:5){
    points(result$cscore[j], result$rscore[i],
           cex = (tabledata[i, j])/mean(tabledata))
  }
}

plotしてみると、度数の多いセルが対角線に近い位置になるよう並んでいることが分かる。

同様の図は単にcorresp()の結果をplotするだけでも得られる。

plot(result)

nf = 2とすることで、各カテゴリーを二次元データで表現できる。これをplotすると自動的にバイプロットが作成され、カテゴリ間の関係を視覚的に確認することができる。

result <- corresp(tabledata, nf = 2)
plot(result)

多重対応分析

テキストのデータだと少し分かりにくいので、MASSパッケージのfarmsデータを例にする。

> summary(farms)
 Mois   Manag  Use    Manure
 M1:7   BF:3   U1:7   C0:6  
 M2:4   HF:5   U2:8   C1:3  
 M4:2   NM:6   U3:5   C2:4  
 M5:7   SF:6          C3:4  
                      C4:3  

farmsは牧草地の環境や管理、利用方法についてのデータで、4つのカテゴリカル変数を持つ。
テキストではMASSパッケージのmca()関数を使用しているが、ここではFactoMineRパッケージのMCA()関数を使用してみる。

library(FactoMineR)
result <- MCA(farms)

MCA()ではplot = FALSEと明示的にOFFにしなければ実行時にプロットが複数表示される。カテゴリーによるバイプロットの結果の部分を示そう。軸名の部分には固有値の割合も同時に示されている。

この例では、第2軸までではデータの変動の40%程度までしか説明できていない。固有値と累積割合はresult$eigで確認できる。

> result$eig
         eigenvalue percentage of variance
dim 1  6.499174e-01           2.166391e+01
dim 2  5.551954e-01           1.850651e+01
dim 3  5.169428e-01           1.723143e+01
dim 4  3.819977e-01           1.273326e+01
dim 5  3.102940e-01           1.034313e+01
dim 6  2.208944e-01           7.363148e+00
dim 7  1.332712e-01           4.442372e+00
dim 8  8.908661e-02           2.969554e+00
dim 9  7.744688e-02           2.581563e+00
dim 10 4.752489e-02           1.584163e+00
dim 11 1.742866e-02           5.809553e-01
dim 12 2.442300e-32           8.140999e-31
       cumulative percentage of variance
dim 1                           21.66391
dim 2                           40.17043
dim 3                           57.40185
dim 4                           70.13511
dim 5                           80.47825
dim 6                           87.84139
dim 7                           92.28377
dim 8                           95.25332
dim 9                           97.83488
dim 10                          99.41904
dim 11                         100.00000
dim 12                         100.00000

また、各カテゴリーがどの程度軸に寄与しているかはresult$var$contribで確認できる。

> result$var$contrib
          Dim 1      Dim 2        Dim 3      Dim 4      Dim 5
M1  1.997032543  6.3987899 6.653719e-04  2.8555287 21.7714144
M2  1.332351513  6.1585438 9.578125e+00 11.0329611  0.1482187
M4  1.942796388  2.6231723 1.344942e+00  4.5402859 46.1550354
M5  9.185471009  2.3084618 8.911297e+00  0.1010620  0.5528613
BF  1.253592590  8.4488996 5.772589e+00 12.9978140 11.0988568
HF  0.474068031 12.2938426 3.762828e+00  3.7546643  2.2497386
NM 20.591151562  3.0315099 1.913005e+00  1.6407989  0.6060488
SF  9.718852258 12.3551134 1.719340e+00  0.2505076  3.1151790
U1  6.648880395  1.2997710 1.072406e+01  0.2151174  2.9499596
U2  9.699621774  6.0048359 2.352817e-02  4.1868505  0.9439940
U3  0.789421444  3.0648808 1.655492e+01  9.8408916  0.6452055
C0 20.591151562  3.0315099 1.913005e+00  1.6407989  0.6060488
C1  0.237784097 12.3464560 7.134416e+00 17.8537482  5.6832374
C2  5.424873542  4.5789897 4.094891e-01 24.3956217  2.2873335
C3  0.002754719  0.1617089 2.943082e+01  3.1784495  0.8359484
C4 10.110196574 15.8935143 8.069635e-01  1.5148998  0.3509199

これを見ると、第1軸にはNM、C0、C4などの影響が大きく、第2軸にはHF、SF、C1、C4などの影響が大きい。

カテゴリカルデータ解析読書メモ(第8章)

対数線形モデル

対数線形モデルにおいては、分割表における各セルの期待度数を予測する。すなわち、複数のカテゴリカル変数の中から特定の目的変数を設定するのでなく、全ての変数について相互の関係を調べることを目的としている。

lonlin()とloglm()

対数線形モデルによる解析を行う関数としてloglin()とloglm()が用意されている。
loglin()は以下のように使用する。

> library(vcd)
> tabledata <- xtabs(Freq ~ Status + Urban, data = DanishWelfare)
> result <- loglin(tabledata, margin = list(1, 2), param = TRUE, fit = TRUE)
2 iterations: deviation 2.273737e-13 
> result$para
$`(Intercept)`
[1] 5.383959

$Status
     Widow    Married  Unmarried 
-0.7693390  1.0654469 -0.2961079 

$Urban
   Copenhagen SubCopenhagen     LargeCity          City       Country 
   -0.4836304    -0.3771835    -0.4102991     0.6787275     0.5923855 

対してloglm()は次のようにモデル式で入力する。

> result <- loglm(~ Status + Urban, data = tabledata)
> result
Call:
loglm(formula = ~Status + Urban, data = tabledata)

Statistics:
                      X^2 df P(> X^2)
Likelihood Ratio 153.9505  8        0
Pearson          158.1145  8        0
> coef(result)
$`(Intercept)`
[1] 5.383959

$Status
     Widow    Married  Unmarried 
-0.7693390  1.0654469 -0.2961079 

$Urban
   Copenhagen SubCopenhagen     LargeCity          City       Country 
   -0.4836304    -0.3771835    -0.4102991     0.6787275     0.5923855 

計算される値はどちらも同じだが、入出力の理解のしやすさはloglm()の方が上の様に思う。以降、loglm()を使用する。

期待度数と残差を求める

期待度数はfitted()、残差はresiduals()を用いて計算できる。

> round(fitted(result), 1)
Re-fitting to get fitted values
           Urban
Status      Copenhagen SubCopenhagen LargeCity   City Country
  Widow           62.2          69.2      67.0  199.0   182.5
  Married        389.9         433.6     419.5 1246.5  1143.4
  Unmarried       99.9         111.1     107.5  319.4   293.0
> round(residuals(result), 1)
Re-fitting to get frequencies and fitted values
           Urban
Status      Copenhagen SubCopenhagen LargeCity City Country
  Widow            6.1           2.5      -1.5  1.0    -6.3
  Married         -5.4          -0.2       0.0  0.0     3.1
  Unmarried        4.5          -1.8       1.2 -0.8    -1.7

三元表での対数線形モデル

変数が3つ以上になると作成できるモデルの数がかなり多くなる。
交互作用項について細かく検討するのであればupdate()を利用すると入力の手間を節約できる。

> Modelindependent <- loglm(~ Gender + Admit + Dept, data = UCBAdmissions)
> Modelindependent
Call:
loglm(formula = ~Gender + Admit + Dept, data = UCBAdmissions)

Statistics:
                      X^2 df P(> X^2)
Likelihood Ratio 2097.671 16        0
Pearson          2000.328 16        0
> # GenderとAdmitの交互作用を加える
> ModelGA <- update(Modelindependent, .~. + Gender:Admit)
> ModelGA
Call:
loglm(formula = . ~ Gender + Admit + Dept + Gender:Admit, data = UCBAdmissions)

Statistics:
                      X^2 df P(> X^2)
Likelihood Ratio 2004.222 15        0
Pearson          1748.160 15        0

テキストでは全てのモデルを作成して比較しているが、面倒なのでstep()を使用してモデル選択をしてみよう。
まずフルモデルを作成する。

> ModelFull <- loglm(~ Gender * Admit * Dept, data = UCBAdmissions)
> ModelFull
Call:
loglm(formula = ~Gender * Admit * Dept, data = UCBAdmissions)

Statistics:
                 X^2 df P(> X^2)
Likelihood Ratio   0  0        1
Pearson            0  0        1

フルモデルでは期待度数と観測度数が一致する。
では、AICによるモデル選択を試みる。

> step(ModelFull)
Start:  AIC=48
~Gender * Admit * Dept

                    Df    AIC
<none>                 48.000
- Gender:Admit:Dept  5 58.204
Call:
loglm(formula = ~Gender * Admit * Dept, data = UCBAdmissions)

Statistics:
                 X^2 df P(> X^2)
Likelihood Ratio   0  0        1
Pearson            0  0        1

AICではフルモデルが選択された。
次にBICによるモデル選択を試みる。

> step(ModelFull, k = log(sum(UCBAdmissions)))
Start:  AIC=202.02
~Gender * Admit * Dept

                    Df    AIC
- Gender:Admit:Dept  5 180.14
<none>                 202.02

Step:  AIC=180.14
~Gender + Admit + Dept + Gender:Admit + Gender:Dept + Admit:Dept

               Df     AIC
- Gender:Admit  1  173.25
<none>             180.14
- Admit:Dept    5  901.45
- Gender:Dept   5 1266.75

Step:  AIC=173.25
~Gender + Admit + Dept + Gender:Dept + Admit:Dept

              Df     AIC
<none>            173.25
- Admit:Dept   5  986.49
- Gender:Dept  5 1351.78
Call:
loglm(formula = ~Gender + Admit + Dept + Gender:Dept + Admit:Dept, 
    data = UCBAdmissions, evaluate = FALSE)

Statistics:
                      X^2 df    P(> X^2)
Likelihood Ratio 21.73551  6 0.001351993
Pearson          19.93841  6 0.002840164

BICを基準にした場合は、Admitとgenderの相互作用がモデルに含まれていない。これは第5章での解析結果とも一致する。
AICとBICで結果が異なったが、データ数が多い場合にAICを基準とすると、フルモデルを選択しやすい傾向があり、BICをモデル選択に使用する場合が多い、とのことである。

練習問題解答

カテゴリカルデータ解析第8章練習問題解答 · GitHub