TIM Labs

sato.yuichiroによるエントリー一覧

このブログ記事は『データ解析のための統計モデリング入門』(久保拓弥 著、岩波書店)という、とても分かりやすい統計モデリングの入門書を、さらに分かりやすくするために読書メモをまとめたものです。 今回は第6章、GLMの応用範囲をひろげるについてのまとめの七回目です。

この章では様々なタイプのGLMについて説明がされています。 その中でガンマ回帰について説明がされています。 なので、いろいろなパラメータについてガンマ分布をプロットしてくれるコードを用意しました。 コードはRで書きました。

rs <- seq(0.5, 4.5, by = 1.0)
ys <- seq(0.5, 19.5, by = 1.0)
xl <- "y"
yl <- "Probability"
legends <- paste0("r = ", rs)
for (s in 1:20) {
  title <- paste0("Gamma distribution, p(y|shape = ", s, ", rate = r)")
  plot(0, 0, type = "n", xlim = c(0, 20), ylim = c(0.0, 1.0), main = title, xlab = xl, ylab = yl)
  for (i in 1:5) {
    lines(ys, dgamma(ys, s, rs[i]), type = "l")
    points(ys, dgamma(ys, s, rs[i]), pch = i)
  }
  legend("topright", legend = legends, pch = 1:5)
}

実行するとグラフが20枚プロットされます。 例えば下図のようなグラフがプロットされます。

gamma_4.png

ガンマ分布はパラメータの組によっては、指数分布のような分布や、ポアソン分布のような山なりの分布を取れることが分かります。

このブログ記事は『データ解析のための統計モデリング入門』(久保拓弥 著、岩波書店)という、とても分かりやすい統計モデリングの入門書を、さらに分かりやすくするために読書メモをまとめたものです。 今回は第6章、GLMの応用範囲をひろげるについてのまとめの六回目です。

この章では様々なタイプのGLMについて説明がされています。 その中でオフセット項わざという手法について説明がされているので、サンプルデータに対してこの手法を試してみるコードを用意しました。 コードはRで書きました。

dataSize <- 100
createData <- function(n, createY) {
  createY <- match.fun(createY)
  x <- runif(dataSize, 0, 5)
  A <- runif(dataSize, n, 10 * n)
  y <- createY(x, A)
  data.frame(x, A, y)
}

linkFunction <- function(x, A) rpois(dataSize, exp(0.1 + 0.3 * x + log(A)))
dataGenerator1 <- function() createData(1, linkFunction)
dataGenerator2 <- function() createData(100, linkFunction)

doGlm <- function(dataGen) {
  dataGen <- match.fun(dataGen)
  data <- dataGen()
  fit <- glm(y ~ x, offset = log(A), data = data, family = poisson)
  c(fit$coefficients[1], fit$coefficients[2])
}

showResult <- function(dataGen) {
  dataGen <- match.fun(dataGen)
  m <- replicate(100, doGlm(dataGen))
  meanm <- apply(m, 1, mean)
  sdm <- apply(m, 1, sd)
  cat("mean b1", meanm[1], "\n")
  cat("mean b2", meanm[2], "\n")
  cat("sd b1", sdm[1], "\n")
  cat("sd b2", sdm[2], "\n")
}

cat("case1\n")
showResult(dataGenerator1)
cat("\ncase2\n")
showResult(dataGenerator2)

サンプルデータに含まれる変数は xAy です。 説明変数は xA です。 応答変数は y です。 ポアソン回帰のリンク関数は λ = exp(b1 + b2 * x + log(A)) という形です。 なので、log(A) 部分がオフセット項になります。 本の例で言うと A は面積のような、y との割算で密度を計算できるような値です。

本の中では1000打数300安打の打者と10打数3安打の打者を同じ三割打者と見なしてよいかという問題に触れられています。 サンプルデータの変数名で言うと、A = 1000y = 300 の場合と、A = 10y = 3 の場合とを同じと見なしてよいかという問題です。 感覚的には、10回中の3回打つ打者より、1000回中の300回打つ打者の方がより正確に三割打者であると言えます。 なので、これらの状況には何らかの差が出てきて欲しいと思うところです。

この問題の状況を見るために、ポアソン回帰を2つの場合について行いました。 ケース1においては A1 から 10 までのランダムな値を取ります。 ケース2においては A100 から 1000 までのランダムな値を取ります。 x はどちらの場合も 0 から 5 までのランダムな値です。 yλ = exp(0.1 + 0.3 * x + log(A)) つまり λ = A * exp(0.1 + 0.3 * x) というリンク関数で得られる λ を平均に持つポアソン分布からサンプリングされます。 なので、A が100倍になれば、y の平均も100倍になります。

コードを実行すると次のような結果が得られます。

case1
mean b1 0.09576121 
mean b2 0.3013721 
sd b1 0.05874369 
sd b2 0.01892569 

case2
mean b1 0.1002007 
mean b2 0.299965 
sd b1 0.006423389 
sd b2 0.00181942 

推定されるパラメータの平均値はどちらも大差がありませんが、その標準偏差は10倍近い差があります。 確かに A の大きさの違いが推定値のばらつきに影響を与えているのが分かります。

このブログ記事は『データ解析のための統計モデリング入門』(久保拓弥 著、岩波書店)という、とても分かりやすい統計モデリングの入門書を、さらに分かりやすくするために読書メモをまとめたものです。 今回は第6章、GLMの応用範囲をひろげるについてのまとめの五回目です。

この章では様々なタイプのGLMについて説明がされています。 その中で、ロジスティック回帰についてAICを使ってモデル選択をする様子が説明されています。 なので、サンプルデータを用意して、いろいろなモデルについてAICを計算してくれるコードを用意しました。 コードはRで書きました。

N <- 10
dataSize <- 100
createData <- function(createY) {
  createY <- match.fun(createY)
  x <- runif(dataSize, 0, 20)
  c <- runif(dataSize, 0, 20)
  y <- createY(x)
  data.frame(x, c, y)
}

dataGenerator <- function()
  createData(function(x) rbinom(dataSize, N, 1 / (1 + exp(5.0 - 0.5 * x))))

data <- dataGenerator()

calcAIC <- function(formula) {
  fit <- glm(formula, data = data, family = binomial)
  fit$aic
}

cat("model\t\taic\n")
cat("null\t\t", calcAIC(cbind(y, N - y) ~ 1), "\n")
cat("x\t\t", calcAIC(cbind(y, N - y) ~ x), "\n")
cat("c\t\t", calcAIC(cbind(y, N - y) ~ c), "\n")
cat("x + c\t\t", calcAIC(cbind(y, N - y) ~ x + c), "\n")
cat("xc\t\t", calcAIC(cbind(y, N - y) ~ x : c), "\n")
cat("x + xc\t\t", calcAIC(cbind(y, N - y) ~ x + x : c), "\n")
cat("c + xc\t\t", calcAIC(cbind(y, N - y) ~ c + x : c), "\n")
cat("x + c + xc\t", calcAIC(cbind(y, N - y) ~ x * c), "\n")

サンプルデータは xcy を変数として含んでいます。 説明変数は x で応答変数は y です。 そこに y と相関のない変数として c を用意しました。 c は必要ない変数なので、モデルに組み込むとバイアスが大きくなるはずです。 コードを実行すると以下のような結果が出力されます。

model       aic
null         885.6348 
x            264.6336 
c            887.2804 
x + c        265.7918 
xc           667.2655 
x + xc       265.8728 
c + xc       381.9791 
x + c + xc   267.7918

確かに x だけがモデルに使われている時に一番AICが低くなるのが分かります。

このブログ記事は『データ解析のための統計モデリング入門』(久保拓弥 著、岩波書店)という、とても分かりやすい統計モデリングの入門書を、さらに分かりやすくするために読書メモをまとめたものです。 今回は第6章、GLMの応用範囲をひろげるについてのまとめの四回目です。

この章では様々なタイプのGLMについて説明がされています。 まずはロジスティック回帰について説明されています。 なので、サンプルデータを用意して、ロジスティック回帰をやってみるコードを用意しました。 コードはRで書きました。

N <- 100
dataSize <- 100
createData <- function(createY) {
  createY <- match.fun(createY)
  x <- runif(dataSize, 0, 200)
  y <- createY(x)
  data.frame(x, y)
}

dataGenerator <- function()
  createData(function(x) rbinom(dataSize, N, 1 / (1 + exp(5.0 - 0.05 * x))))

data <- dataGenerator()
fit <- glm(cbind(y, N - y) ~ x, data = data, family = binomial)
b1 <- fit$coefficients[1]
b2 <- fit$coefficients[2]

linkFunction <- function(x) 1 / (1 + exp(-b1 - b2 * x))

xs <- 1:200
ys <- 0:N
title <- "Logistic regression"
xl <- "x"
yl <- "The most probable y"
maxY <- numeric(length = 200)
for (x in xs) {
  q <- linkFunction(x)
  maxY[x] <- which.max(dbinom(ys, N, q)) - 1
}
plot(0, 0, type = "n", xlim = c(0, 200), ylim = c(0, N), main = title, xlab = xl, ylab = yl)
lines(xs, maxY, type="l")

前回もサンプルデータに対してロジスティック回帰をやってみて結果をプロットしました。 今回は説明変数 x に対して、どの応答変数 y が一番高い確率で得られるかをプロットしました。 コードを実行すると下図のような図がプロットされます。

logistic_regresson_2.png

ロジスティック関数のような曲線が得られました。 もちろん y の値はこの曲線の周囲に分布しているわけですが、一番取りやすい値をつなげると、このような曲線になるということです。

このブログ記事は『データ解析のための統計モデリング入門』(久保拓弥 著、岩波書店)という、とても分かりやすい統計モデリングの入門書を、さらに分かりやすくするために読書メモをまとめたものです。 今回は第6章、GLMの応用範囲をひろげるについてのまとめの三回目です。

この章では様々なタイプのGLMについて説明がされています。 まずはロジスティック回帰について説明されています。 なので、サンプルデータを用意して、ロジスティック回帰をやってみるコードを用意しました。 コードはRで書きました。

N <- 10
dataSize <- 100

createData <- function(createY) {
  createY <- match.fun(createY)
  x <- runif(dataSize, 0, 20)
  y <- createY(x)
  data.frame(x, y)
}

dataGenerator <- function()
  createData(function(x) rbinom(dataSize, N, 1 / (1 + exp(5.0 - 0.5 * x))))

data <- dataGenerator()

fit <- glm(cbind(y, N - y) ~ x, data = data, family = binomial)
b1 <- fit$coefficients[1]
b2 <- fit$coefficients[2]

linkFunction <- function(x) 1 / (1 + exp(-b1 - b2 * x))

xs <- seq(5, 20, by = 3)
ys <- 0:N
title <- "Logistic regression"
xl <- "y"
yl <- "Probability"
legends <- paste0("x = ", xs)
plot(0, 0, type = "n", xlim = c(0, N), ylim = c(0.0, 1.0), main = title, xlab = xl, ylab = yl)
for (x in xs) {
  q <- linkFunction(x)
  y <- dbinom(ys, N, q)
  lines(ys, y, type="l")
  points(ys, y, pch = x)
}
legend("topleft", legend = legends, pch = xs)

作ったサンプルデータにロジスティック回帰をやって、リンク関数のパラメータを推定しているのは以下の部分です。

fit <- glm(cbind(y, N - y) ~ x, data = data, family = binomial)
b1 <- fit$coefficients[1]
b2 <- fit$coefficients[2]

こうして得られた b1b2 を使って以下のリンク関数を作れば、説明変数 x に対して、事象の発生確率 q が得られます。

linkFunction <- function(x) 1 / (1 + exp(-b1 - b2 * x))

これで事象を N 回繰り返した時の応答変数 y の分布が求まるので、あとはそれをプロットします。 コードを実行すると、下図のようなグラフがプロットされます。

logistic-regression.png

推定される b2 は正の数なので、x が大きくなるにしたがって y の分布も大きな値を取るようになっています。

このブログ記事は『データ解析のための統計モデリング入門』(久保拓弥 著、岩波書店)という、とても分かりやすい統計モデリングの入門書を、さらに分かりやすくするために読書メモをまとめたものです。 今回は第6章、GLMの応用範囲をひろげるについてのまとめの二回目です。

この章では様々なタイプのGLMについて説明がされています。 まずはロジスティック回帰について説明されています。 なので、いろいろなパラメータについてロジスティック回帰のモデルをプロットしてくれるコードを用意しました。 コードはRで書きました。

N <- 20

linkFunction <- function(b1, b2) {
  function(x) {
    1 / (1 + exp(-b1 - b2 * x))
  }
}

distribution <- function(x, linkFunction) {
  l <- match.fun(linkFunction)
  function(y) {
    dbinom(y, N, l(x))
  }
}

xs <- seq(0, 10, by = 2)
ys <- 0:N
xl <- "y"
yl <- "Probability"
legends <- paste0("x = ", xs)
for (b1 in seq(-0.5, 0.5, by = 0.2)) {
  for (b2 in seq(-0.5, 0.5, by = 0.2)) {
    l <- linkFunction(b1, b2)
    title <- paste0("q = 1 / (1 + exp(-(", b1, " + ", b2, "x)))")
    plot(0, 0, type = "n", xlim = c(0, N), ylim = c(0.0, 1.0), main = title, xlab = xl, ylab = yl)
    for (x in xs) {
      d <- distribution(x, l)
      lines(ys, d(ys), type = "l")
      points(ys, d(ys), pch = x)
    }
    legend("topright", legend = legends, pch = xs)
  }
}

リンク関数はロジスティック関数です。 コードでは以下の部分です。

linkFunction <- function(b1, b2) {
  function(x) {
    1 / (1 + exp(-b1 - b2 * x))
  }
}

これが二項分布のパラメータの一つである事象の発生確率 q を与えます。 もう一つのパラメータは試行回数 N ですが、これは 20 にしました。 ということで、リンク関数のパラメータ b1b2 さえ与えれば、説明変数 x に対して応答変数 y の分布がどのように変化するかが分かります。

コードを実行すると様々な b1b2 の値に対してモデルをプロットします。 例えば下図のような図がプロットされます。

binomial_model.png

これは b1 = -0.3b2 = 0.3 の場合です。 x0 ならq = 1 / (1 + exp(0.3)) = 0.425... です。 x10 ならq = 1 / (1 + exp(-2.7)) = 0.937... です。 b2 が正なので x の増加に従って exp の引数部分がどんどん小さくなって、q 全体としては 1 に近づくわけです。 事象の発生確率が高くなれば、N 回の試行でその事象が起きる回数は大きくなります。 図を見ると、確かに x の増加に従って y の分布が大きな値を取るようになっているのが分かります。

このブログ記事は『データ解析のための統計モデリング入門』(久保拓弥 著、岩波書店)という、とても分かりやすい統計モデリングの入門書を、さらに分かりやすくするために読書メモをまとめたものです。 今回は第6章、GLMの応用範囲をひろげるについてのまとめの一回目です。

この章では様々なタイプのGLMについて説明がされています。 まずはロジスティック回帰について説明されています。 ロジスティック回帰で使われる分布は二項分布なので、いろいろなパラメータについて二項分布をプロットしてくれるコードを用意しました。 コードはRで書きました。

qs <- seq(0.1, 0.9, by = 0.2)
xl <- "y"
yl <- "Probability"
legends <- paste0("q = ", qs)
for (N in 1:20) {
  title <- paste0("Binomial distribution, p(y|N = ", N, ", q)")
  plot(0, 0, type = "n", xlim = c(0, N), ylim = c(0.0, 1.0), main = title, xlab = xl, ylab = yl)
  for (i in 1:5) {
    lines(0:N, dbinom(0:N, N, qs[i]), type = "l")
    points(0:N, dbinom(0:N, N, qs[i]), pch = i)
  }
  legend("topright", legend = legends, pch = 1:5)
}

実行するとグラフが20枚プロットされます。 例えば下図のようなグラフがプロットされます。

binomial_15.png

二項分布は確率 q で起きる事象が N 回中 y 回起きる確率の分布です。 この図は N = 15 の場合です。 q が小さいうちは事象は少ない回数起きる確率が高く、 q が大きくなるにしたがって多くの回数起きる確率が高くなるのが分かります。

このブログ記事は『データ解析のための統計モデリング入門』(久保拓弥 著、岩波書店)という、とても分かりやすい統計モデリングの入門書を、さらに分かりやすくするために読書メモをまとめたものです。 今回は第5章、GLMの尤度比検定と検定の非対称性についてのまとめの六回目です。

この章では尤度比検定などについて説明がされています。 尤度比検定の手法の一つとしてパラメトリックブートストラップ法が説明されています。 対立仮説が正しい場合でも、尤度比検定は必ずしも帰無仮説を棄却できるわけではありません。 そこで、パラメトリックブートストラップ法で帰無仮説を棄却できない場合がどれくらい現れるかを見るためのコードを用意しました。 コードはRで書きました。

createData <- function(n, b1, b2) {
  x <- runif(n, 1, 20)
  y <- rpois(n, lambda = exp(b1 + b2 * x))
  data.frame(x, y)
}

bootstrap <- function(b1) {
  simulatedData <- createData(100, b1, 0)
  fit1 <- glm(formula = y ~ 1, family = poisson, data = simulatedData)
  fit2 <- glm(formula = y ~ x, family = poisson, data = simulatedData)
  fit1$deviance - fit2$deviance
}

test <- function(n, b2) {
  expData <- createData(n, 0.1, b2)
  fit1 <- glm(formula = y ~ 1, family = poisson, data = expData)
  fit2 <- glm(formula = y ~ x, family = poisson, data = expData)
  dev <- fit1$deviance - fit2$deviance

  P <- mean(replicate(100, bootstrap(fit1$coefficient[1]) > dev))
  P > 0.05
}

xl <- "Data size"
yl <- "The ratio of failing to reject"
b2 <- c(0.05, 0.07, 0.1, 0.12)
legends <- paste0("b2 = ", b2)
plot(0, 0, type = "n", xlim = c(0, 20), ylim = c(0, 1), xlab = xl, ylab = yl)
n <- 1:20
for(i in 1:4) {
  r <- vector(length = 20)
  for (n in 1:20) {
    r[n] <- mean(replicate(100, test(n, b2[i])))
  }
  lines(1:20, r, type = "b", pch = i) 
}
legend("topright", legend = legends, pch = 1:4)

実行するとグラフを1枚プロットします。 最尤推定を 202 * 100 * 20 * 4 回やるので結構時間がかかります。

有意水準としては5%を設定しました。 なので、パラメトリックブートストラップ法を使ってP値を計算し、得られたP値が 0.05 より大きいと帰無仮説を棄却できません。 test は検定に使うデータサイズ n と、データをサンプリングする分布の中の説明変数 x の係数 b2 を受け取り、P値を計算して帰無仮説が棄却されなかった場合 true を返します。 あとは、各 nb2 について、それを100回繰り返して、true が帰ってくる確率を求め、プロットします。 プロットされるのは下図のような図です。

failing_to_reject.png

傾向として、データの数が多くなって x との相関が強いほど、帰無仮説は棄却されやすくなるのが分かります。

このブログ記事は『データ解析のための統計モデリング入門』(久保拓弥 著、岩波書店)という、とても分かりやすい統計モデリングの入門書を、さらに分かりやすくするために読書メモをまとめたものです。 今回は第5章、GLMの尤度比検定と検定の非対称性についてのまとめの五回目です。

この章では尤度比検定などについて説明がされています。 尤度比検定の手法の一つとしてパラメトリックブートストラップ法が説明されています。 なので、パラメトリックブートストラップ法でどのようにデータがシミュレートされているか見るためのコードを用意しました。 コードはRで書きました。

linkFunction <- function(b1, b2) {
  function(x) {
    exp(b1 + b2 * x)
  }
}

createData <- function(n, b1, b2) {
  x <- runif(n, 1, 20)
  y <- rpois(n, lambda = linkFunction(b1, b2)(x))
  data.frame(x, y)
}

plotGraph <- function(data, fit1, fit2, title) {
  plot(data$x, data$y, type = "p", main = title, xlab = "x", ylab = "y")
  lines(0:20, linkFunction(fit1$coefficient[1], 0)(0:20), lty = "solid")
  lines(0:20, linkFunction(fit2$coefficient[1], fit2$coefficient[2])(0:20), lty = "dashed")
  nl <- paste0("y = exp(", round(fit1$coefficient[1], 3), ")")
  xl <- paste0("y = exp(", round(fit2$coefficient[1], 3),
    " + ", round(fit2$coefficient[2], 3), " * x)")
  legend("topleft", legend = c(nl, xl), lty = c("solid", "dashed"))
}

bootstrap <- function(b1) {
  simulatedData <- createData(100, b1, 0)
  fit1 <- glm(formula = y ~ 1, family = poisson, data = simulatedData)
  fit2 <- glm(formula = y ~ x, family = poisson, data = simulatedData)
  fit1$deviance - fit2$deviance
}

for(i in 1:10) {
  expData <- createData(20, 0.1, 0.07)
  fit1 <- glm(formula = y ~ 1, family = poisson, data = expData)
  fit2 <- glm(formula = y ~ x, family = poisson, data = expData)
  dev <- fit1$deviance - fit2$deviance

  P <- mean(replicate(1000, bootstrap(fit1$coefficient[1])) > dev)
  title <- paste("P value", P)
  plotGraph(expData, fit1, fit2, title)
}

コードを実行するとグラフが10枚プロットされます。 それぞれにポアソン回帰で使われるデータと、推定されたモデルをプロットしてあります。 あと、パラメトリックブートストラップ法でP値を計算してタイトルに入れてあります。 典型的には下図のようなグラフがプロットされます。

P_value_0.png

ですが、時々は下図のようなグラフがプロットされます。

P_value_0191.png

有意水準として5%を設定するなら、P値が0.05以上だと帰無仮説を棄却できません。 棄却できる時は対立仮説は帰無仮説とはだいぶ形の違う曲線になりますが、棄却できない時は直線に近い形になるのが分かります。 データは x に相関のある分布から生成されていますが、たまたまデータが偏っていると最尤推定で x と相関のないようなパラメータを推定してしまう場合があるわけですね。

このブログ記事は『データ解析のための統計モデリング入門』(久保拓弥 著、岩波書店)という、とても分かりやすい統計モデリングの入門書を、さらに分かりやすくするために読書メモをまとめたものです。 今回は第5章、GLMの尤度比検定と検定の非対称性についてのまとめの四回目です。

この章では尤度比検定などについて説明がされています。 尤度比検定の手法の一つとしてパラメトリックブートストラップ法が説明されています。 前回はパラメトリックブートストラップ法でどのようにデータがシミュレートされているか見るためのコードを紹介しました。

コードを実行するとグラフが3枚プロットされます。 1枚目は以下のような図です。 createData を使ってサンプリングしたデータと、帰無仮説と対立仮説に対してポアソン回帰して得られた曲線がプロットされています。 このデータは何かの実験データのような、検定にかけたいデータだとします。 帰無仮説は formulay ~ 1 を使っている方で、対立仮説は y ~ x を使っている方です。 データは x と正の相関のあるデータなので、対立仮説は右肩上がりの曲線になっています。

ExperimentalData.png

2枚目は以下のような図です。 パラメトリックブートストラップ法で使う、ポアソン回帰で得られた帰無仮説のパラメータを使ってシミュレートしたデータと、そのデータに対して新たにポアソン回帰して得られた帰無仮説と対立仮説の曲線がプロットされています。 帰無仮説は x と相関のないモデルなので、これを真の分布だと仮定してシミュレートすると、得られるデータは x に相関のないデータになります。 すると、帰無仮説も対立仮説も、ほとんど同じような直線になってしまいます。

Bootstrap_null.png

3枚目は以下のような図です。 パラメトリックブートストラップ法ではこういうシミュレーションは行いませんが、参考までに帰無仮説ではなく、対立仮説の方のパラメータを使ってデータをシミュレートして、2枚目と同じようなプロットを作ってみました。 対立仮説は x と相関のあるモデルなので、そのパラメータを使うと x に相関のあるデータがシミュレートされます。 すると、ポアソン回帰で得られる曲線は、実験データの時と同じような曲線になります。

Bootstrap_x.png

このように、検定にかけたい元の実験データと、パラメトリックブートストラップで行うシミュレーションで生成されたデータは性質が違うデータです。 この違いを利用して尤度比検定は帰無仮説を棄却できるのですね。

このアーカイブについて

このページには、sato.yuichiroが最近書いたブログ記事が含まれています。

最近のコンテンツはインデックスページで見られます。過去に書かれたものはアーカイブのページで見られます。