TIM Labs

データ解析のための統計モデリング入門 GLMの尤度比検定と検定の非対称性 読書メモ1

このブログ記事は『データ解析のための統計モデリング入門』(久保拓弥 著、岩波書店)という、とても分かりやすい統計モデリングの入門書を、さらに分かりやすくするために読書メモをまとめたものです。 今回は第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) {
  plot(data$x, data$y, type = "p", main = "Bootstrap", 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
}

showBootstrap <- function(dataSize, b1, b2) {
  d <- createData(dataSize, b1, b2)
  fit1 <- glm(formula = y ~ 1, family = poisson, data = d)
  fit2 <- glm(formula = y ~ x, family = poisson, data = d)
  dev <- fit1$deviance - fit2$deviance

  devs <- replicate(1000, bootstrap(fit1$coefficients[1]))
  cat("data size", dataSize, paste0(", lambda = exp(", b1, " + ", b2, "x)"), "\n")
  summary(d$y)
  P <- sum(devs >= dev) / 1000
  cat("P =", P, "\n")
  if (P < 0.05) {
    cat("The null hypothesis is rejected.\n\n")
  } else {
    cat("Failling to reject the null hypothesis.\n\n")
  }
  plotGraph(d, fit1, fit2)
}

showBootstrap(5, 0.01, 0.03)
showBootstrap(5, 0.01, 0.1)
showBootstrap(20, 0.01, 0.03)
showBootstrap(20, 0.01, 0.1)

コードを実行すると、ポアソン分布からサンプルデータをサンプリングして、二つのモデルでポアソン回帰をし、尤度比検定を行って帰無仮説が棄却できるかを見ます。 サンプルデータを用意して、ポアソン回帰をしているのは以下のコードです。 createData は平均 exp(b1 + b2 * x) を持つポアソン分布からデータをサンプリングします。 帰無仮説は formulay ~ 1 を使っている方で、対立仮説は y ~ x を使っている方です。 最後に逸脱度の差を計算しています。

  d <- createData(dataSize, b1, b2)
  fit1 <- glm(formula = y ~ 1, family = poisson, data = d)
  fit2 <- glm(formula = y ~ x, family = poisson, data = d)
  dev <- fit1$deviance - fit2$deviance

パラメトリックブートストラップ法では、推定した帰無仮説のパラメータを使ってデータをシミュレートしますが、それをやっているのは以下のコードです。 b1 には推定したパラメータが渡されてきます。 そして平均 exp(b1) を持つポアソン分布からデータをサンプリングして、そのデータについてポアソン回帰をやり直し、逸脱度の差を計算します。

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
}

あとはこのシミュレーションを繰り返し、最初の逸脱度の差がどれだけ珍しいか確率を計算して帰無仮説を棄却できるかどうか結論します。 コードを何回か実行してみると、だいたいの傾向としてサンプル数が多くてx の相関が強いほど帰無仮説は棄却されやすくなるのが分かると思います。 実行結果の解説は次回やります。

トラックバック(0)

トラックバックURL: http://labs.timedia.co.jp/mt/mt-tb.cgi/586

コメントする

このブログ記事について

このページは、sato.yuichiroが2017年3月17日 12:00に書いたブログ記事です。

ひとつ前のブログ記事は「Chainer:MNISTの手書き数字は使えそう」です。

次のブログ記事は「GPU:C言語開発環境CUDAのインストール」です。

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