RGRからゴンペルツ曲線を推定する 2

R

前回はRGRを重量の対数の変化速度として求めたが、そもそもRGRの定義は

{\displaystyle r = \frac{1}{w}\frac{dw}{dt} }

であったので、複数のデータポイントがあれば何らかの方法で平滑化して微分すれば、より精度良くRGRが推定できるだろう。例えば前後2点を含めた3点をラグランジュ補間して得た二次関数を微分するとか。

Rでラグランジュ補間をするパッケージは多分あるが面倒なので二次関数を当てはめてしまう。

r_t <- function(w, t) sum(coef(lm(w ~ t + I(t^2)))[2:3] * c(1, 2*t[2])) / w[2]

r <- numeric(length(g))
r <- NA
for(i in 1:length(g)){
  r[i] <- r_t(g[c(i-1, i, i+1)], t[c(i-1, i, i+1)])
}

そしてあとは前回のヤツを多少いじって実行する。

result <- glm(r ~ t, family = gaussian(link="log"), start = c(1, 0))
result_coef <- coef(result)
plot(t, g)
f <- function(c1) (sum(gompertz(t, c1, result_coef[1], -result_coef[2]) - g))^2
c1 <- optimize(f, interval = c(0, 100))$minimum
points(t, gompertz(t, c1, result_coef[1], -result_coef[2]), col = "red", type = "l")

すると少しはマシな推定値が出る。
f:id:Rion778:20161022170713p:plain

ただ、繰り返しになるがゴンペルツ曲線を求めるなら最初からゴンペルツ曲線をデータに当てはめたほうが良い。例えばgに適当な誤差を与えて推定を繰り返してみると、マズさが分かるだろう。

特定のタイミングにおけるRGRの推定値がほしければ、ゴンペルツ関数を微分したものに推定したパラメータを渡せば良い。微分した結果を関数として使うにはderiv()を使って引数にfunc=TRUEを指定する。

dg <- deriv(~ c1 * exp(-(exp(a-b*t)/b)), "t", func = TRUE)
attr(dg(t), "gradient")
r <- attr(dg(t), "gradient")[,1] / gompertz(t, c1, a, b)
plot(r)

f:id:Rion778:20161022172949p:plain

RGRからゴンペルツ曲線を推定する

R

植物の重量wについて、瞬間的な成長速度は植物の重量に比例すると仮定する。

\frac{dw}{dt} = rw

このときrを相対成長速度(RGR: Relative Growth Rate)と呼ぶ。もしrが変化しないのであれば、植物は指数関数的に増大を続けてしまい、地球は程なくして滅ぶだろう。

w_t = w_0 e^{rt}

しかし、植物は無限に大きくなるわけではないので、一般的には相対成長速度は時間とともに減少する。ここに例えば適当な定数aとbを用意して、rは次のように減少すると仮定してみよう。

r= e^{a-bt}

先の式に代入すれば次のようになる。

\frac{dw}{dt} = we^{a-bt}

そして面倒なのでWolfram Alphaにまかせて微分方程式を問いてもらうとこうなる。

w_t = c_1 e^{-\frac{e^{a-b t}}{b}}

さらに定数項を少し整理して、c_2をこのように定義する。

c_2 = e^{-\frac{e^a}{b}}

そして得られた式をゴンペルツ関数(Gompertz function)と呼ぶ。

w_t = c_1c_2^{e^{-bt}}

瞬間的なRGRの推定

最初の仮定と若干違ってしまうのだが、短い時間\Delta tで植物の体重を測定したためにrが変化しなかったとしよう。すると、最初の式を\Delta tの間で積分することで次の式が得られる。

r = \frac{\Delta \ln w}{\Delta t}

実際の測定で推定できそうな雰囲気がしてきただろう。これで植物体重の変化からRGRの変化が推定できる。そして、RGRの変化についてはr=e^{a-bt}と仮定したので、定数aとbは、対数リンク関数を設定したglmで推定できる。

RGRからゴンペルツ曲線を推定する

まず、aとbが残ったバージョンでゴンペルツ関数を定義し、適当なパラメータでゴンペルツ曲線を描いてみよう。

gompertz <- function(t, c1, a, b){
  c1 * exp(-(exp(a-b*t)/b))
}

t <- seq(1, 20, 0.5)
a <- 0.3
b <- 0.25
c1 <- 10
plot(t, gompertz(t,c1, a, b), type = "l")

f:id:Rion778:20161020222308p:plain

次に重量を対数変換し、その変化量としてRGRの近似値を得てみよう。

g <- gompertz(t, c1, a, b)
r <- diff(log(g))/0.5
plot(t[-1], [r)

f:id:Rion778:20161020224446p:plain]

本当の植物の葉っぱなどを実際に測定をしてみる(重量測定は破壊が伴うので、例えば葉面積で代替する)と、RGRは実際にこのような減衰をすることが多い。

あとはglmでaとbを推定する。c_1についてはoptimize()で頑張って求める。

newt <- t[1:length(r)]+0.5
result <- glm(r ~ newt, family = gaussian(link="log"), start = c(1, 0))
result_coef <- coef(result)
plot(t, g)
f <- function(c1) (sum(gompertz(newt, c1, result_coef[1], -result_coef[2]) - g[-1]))^2
c1 <- optimize(f, interval = c(0, 100))$minimum
c1 <- 10
points(newt, gompertz(newt, c1, result_coef[1], -result_coef[2]), col = "red", type = "l")

するとこうだ。
f:id:Rion778:20161020225004p:plain
測定間隔の間はRGRが変化しないと仮定しているので、少しずれる。測定間隔が広いほどずれは大きくなる。

重量データに直接当てはめる

ちなみに、もし本当にゴンペルツ曲線のパラメータがほしいのであれば、重量の測定データに直接当てはめるべきで、このような回り道をすべきではない。
例えばoptim()を使うとこうなる。

f <- function(c) (sum((gompertz(t, c[1], c[2], c[3]) - g)^2))
res <- optim(c(10, 0.5, 1), f)$par
res # > [1] 9.9993114 0.3007776 0.2500798
plot(t, g)
points(t, gompertz(t, res[1], res[2], res[3]), col = "red", type = "l")

f:id:Rion778:20161020233012p:plain

Excelで表形式に入力されてしまったデータを整形するにはINDEX関数が便利

Excel

f:id:Rion778:20161017191321p:plain
こういう「お前は何がしたかったんだ」みたいなファイルが引き継がれることありますよね？ない？そいつはよかったな！！！！！

クソみたいなデータを整形する

普通データはこう入力するだろう。こうでなくては何も始められない。
f:id:Rion778:20161017191525p:plain
だが不幸なことに現実はこうではないのだ。何かを始めるためには現実を修正しなければならない。

INDEX関数

データの整形にはINDEX関数を使う。

INDEX関数はこのように引数を指定する。

=INDEX(配列, 行番号, 列番号)

すると、配列中の所定の行、列の位置にある値を返す。

Excelに配列などあったかという話だが、ある。例えば、{}の中にコンマ","もしくはセミコロン";"で区切って要素を指定することで直接値を指定して配列を作成できる。コンマで区切ると列が、セミコロンで区切ると行が増える。例えば、次の式は6を返す。

=INDEX({1, 2, 3; 4, 5, 6}, 2, 3)

セル内の数式はデバッグしにくいのでこのようなデータ入力はすべきではないが、もしやるなら改行しておくと分かりやすいだろう。

=INDEX({1,2,3;
        4,5,6},
        2, 3)

ちなみにセル内改行はWindowsだとAlt+Enterだが、MacではCommand+Option+Enterである。Windows版とMac版ではこういう細かい違いがたくさんあって大変不快である。

そして、セルの範囲は配列としても振る舞うので、とりあえずこの場では配列=セル範囲という認識でよい。

作業列を追加する

あとは元の表に対応する行番号と列番号を生成するだけだが、これは作業列を追加すると便利である。作業列などを追加するのはなんだかダサい気がするが、INDEX関数の中でごちゃごちゃやると何かミスが有っても気づきにくいので、デバッグのしやすさを考慮すると作業列(あるいは作業行)は積極的に使っていったほうが良い。

作業列を作成するコツは最初の1ループ分(今回の例では表の1行目分)は手で入力し、それに続く要素は最初のループを参照するように指定することだ。
f:id:Rion778:20161017193915p:plain
そして、必要であれば参照に修正を加える。今回の例では行番号が3ずつ増えるよう修正している。オートフィルは予期せぬ動作をする例が多いので、あまり当てにしてはならない。また、行の追加・削除をする場合には、全体をコピーして値貼り付けして値に変換しておいた方が良い(Excelが空気を読んで調整してくれる場合も多いが、信用してはならない)。

作業列が完成したらINDEX関数を組んでいく。配列にするセル範囲を絶対参照にするのを忘れてはならない。完成したら下に伸ばしていくだけだ。作業列がちゃんと完成していれば、アクティブセルの右下をダブルクリックするだけで適切な範囲までオートフィルしてくれる。ドラッグによるオートフィルは時間がかかる上にミスしやすいので、避けられる場面では避けよう。
f:id:Rion778:20161017195121p:plain
なお、絶対参照にするショートカットはWindowsではF4だが、MacではCommand+tである。いい加減にして欲しい。

Project Euler Problem 1:10

R Project Euler

久しぶりに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))
}

するとこのようなものが出来る。
f:id:Rion778:20161008220158p:plain
空間相関を考慮しないバージョンはこんな感じ。

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]);
}

f:id:Rion778:20161008220315p:plain

参考

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

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

R

前回。
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の扱い

R ggplot2

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

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

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

ggplot() + geom_bar(aes(data))

f:id:Rion778:20160825202906p:plain
ところがfactorにするとこうだ。

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

f:id:Rion778:20160825202941p:plain
どうも集計に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でやってみる。
データ解析のための統計モデリング入門――一般化線形モデル・階層ベイズモデル・MCMC (確率と情報の科学)

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)

f:id:Rion778:20160815014841p:plain

結果のプロット

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

traceplot(fit, inc_warmup = TRUE)

f:id:Rion778:20160815015315p:plain
また、サンプリング結果は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)

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

stan_dens(fit1, separate_chains  = TRUE)

f:id:Rion778:20160817204956p:plain
(※追記ここまで)

テキストの図を再現する

テキストの図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)

f:id:Rion778:20160815020806p:plain

WindowsにRStanをインストール

R Windows

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

環境

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章)

R

決定木

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

決定木を書く

テキストでは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)

f:id:Rion778:20160807214654p:plain

量的変数を用いる

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)

f:id:Rion778:20160807215145p:plain
また、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)

f:id:Rion778:20160807215932p:plain

練習問題解答

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