TIM Labs

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

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

plotBootstrap <- function(b1, b2) {
  simulatedData <- createData(100, b1, b2)
  fit1 <- glm(formula = y ~ 1, family = poisson, data = simulatedData)
  fit2 <- glm(formula = y ~ x, family = poisson, data = simulatedData)
  title <- paste0("Bootstrap (b1 = ", b1, ", b2 = ", b2, ")")
  plotGraph(simulatedData, fit1, fit2, title)
}

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)
plotGraph(expData, fit1, fit2, "Experimental data")

plotBootstrap(fit1$coefficient[1], 0)
plotBootstrap(fit2$coefficient[1], fit2$coefficient[2])

コードを実行するとグラフが3枚プロットされます。 ぜひ実行してみてください。 実行結果の解説は次回にまわします。

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

この章では尤度比検定などについて説明がされています。 尤度比検定の手法の一つとしてパラメトリックブートストラップ法が説明されています。 前回はパラメトリックブートストラップ法がどのように動くのかを見るためのコードを紹介しました。 コードを実行すると結果は次のようになります。 何度か実行してみて、だいたいの傾向をつかんでみてください。

data size 5 , lambda = exp(0.01 + 0.03x) 
P = 0.956 
Failling to reject the null hypothesis.

data size 5 , lambda = exp(0.01 + 0.1x) 
P = 0.098 
Failling to reject the null hypothesis.

data size 20 , lambda = exp(0.01 + 0.03x) 
P = 0.128 
Failling to reject the null hypothesis.

data size 20 , lambda = exp(0.01 + 0.1x) 
P = 0 
The null hypothesis is rejected.

だいたいの傾向としてサンプル数が多くて x の相関が強いほど帰無仮説は棄却されやすくなります。 データを生成する真の分布は x と相関があるので、帰無仮説は棄却されて欲しいですが、そうならない場合があるわけです。 その様子をグラフにもしてあるので、紹介します。 まずは帰無仮説が棄却されて、無事正しいモデルが支持された場合です。

bootstrap_sample2.png

検定に使われたデータが丸、実線が帰無仮説、点線が対立仮設です。 対立仮設の x の係数は 0.098 と推定されています。 おかげで対立仮設は帰無仮説と比べて x と相関のあるデータによく当てはまり、二つのモデルの逸脱度の差が大きくなってP値が低くなり尤度比検定で正しいモデルが選択されたわけです。

一方、帰無仮説が棄却できなかった場合はこのような図になります。

bootstrap_sample1.png

一見して帰無仮説も対立仮設も違いが見えません。 対立仮設の x の係数は 0.008 と推定されていて 0 と違いはほとんどありません。 とすると、帰無仮説でも対立仮設でもデータへの当てはまりの良さに違いはほとんどなく、十分なP値が得られなかったということです。

だいたいの感覚として、帰無仮説が棄却されたりされなかったりする時はこういうことが起きています。

このブログ記事は『データ解析のための統計モデリング入門』(久保拓弥 著、岩波書店)という、とても分かりやすい統計モデリングの入門書を、さらに分かりやすくするために読書メモをまとめたものです。 今回は第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 の相関が強いほど帰無仮説は棄却されやすくなるのが分かると思います。 実行結果の解説は次回やります。

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

この章ではバイアスやAICについて説明がされています。 データを説明するモデルとして様々なモデルが考えられる時に、データへのあてはまりの良さだけを基準にモデルを選ぶことはできません。 前回はAICを使ってモデル選択の様子を見るためのコードの内容を解説しました。 今回は実行結果の内容を解説します。

コードを実行すると結果は次のようになります。

indipendent data
                    null        x           c           xc
log likelihood       -208.4236   -207.8885   -207.5525   -207.5126 
mean log likelihood  -209.1079   -209.7084   -209.7     -210.2147 
bias                 0.6843544   1.819908    2.14749     2.702091 
aic                  418.8471    419.7771    419.1049    421.0252 

correlated X data
                    null        x           c           xc
log likelihood       -244.4405   -191.2695   -243.681    -190.9884 
mean log likelihood  -246.9158   -193.4544   -248.0273   -193.9776 
bias                 2.475276    2.184949    4.346326    2.989236 
aic                 490.881      386.5389    491.3619    387.9768 

correlated C data
                    null        x           c           xc
log likelihood       -245.6007   -244.4792   -191.0656   -190.7746 
mean log likelihood  -246.9547   -247.9338   -193.4138   -193.9289 
bias                 1.354043    3.454623    2.348202    3.154229 
aic                 493.2014     492.9584    386.1313    387.5493 

correlated XC data
                    null        x           c           xc
log likelihood       -584.6693   -410.6131   -410.5213   -244.1359 
mean log likelihood  -591.7481   -421.7135   -421.5719   -246.7352 
bias                 7.078751    11.10039    11.05061    2.599335 
aic                  1171.339    825.2262    825.0427    494.2717 

結果の見方ですが、上から順に xc の両方に関係のないデータ、x と関係したデータ、c と関係したデータ、xc の両方と関係したデータについて計算した場合です。 それぞれの場合において、4つのモデルに対する最大対数尤度、平均対数尤度、バイアス、AICが計算されています。 モデルは左から順にリンク関数に exp(b1)exp(b1 + b2 * x)exp(b1 + b3 * c)exp(b1 + b2 * x + b3 * c) を使った場合です。

どの場合においても、対数尤度が最大になるのは、xc の両方を説明変数に含むモデルです。 ここには複雑なモデルほどデータによく当てはまるという一般的な性質が表れています。 データが xc の両方に関係する場合には望ましい結果ですが、実際にはデータが xc に関係のない場合にも、関係あるとしてモデルを作った方が特定のデータへの当てはまりは良くなるわけです。 しかし、平均対数尤度は大きくなっていません。 これは、無用に複雑なモデルは推定のために与えられたデータだけに過度に当てはまって、未知のデータへの当てはまりが悪くなっていることを意味します。 尤度だけを使ってモデルの良さを測ることはできないわけです。

そこで、AICの出番です。 AICはモデルのパラメータの個数を当てはまりの良さから割り引いて計算するので、無用なパラメータを増やしても良い数値になりません。 実行結果でも、AICが最小になるモデルを選べば正しくモデルを選択できることが分かります。

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

この章ではバイアスやAICについて説明がされています。 データを説明するモデルとして様々なモデルが考えられる時に、データへのあてはまりの良さだけを基準にモデルを選ぶことはできません。 前回はAICを使ってモデル選択の様子を見るためのコードを紹介しました。 今回はそのコードの内容を解説します。

このコードを実行すると、4種類ある分布からサンプリングしたデータに対して、4種類のモデルを使ってポアソン回帰をし、AICやバイアスを計算して表示します。 計算されたAICを使って複数あるモデルから最適な物を選択する様子を見ることができます。

回帰に使うデータを作っているのは以下のコードです。 xc0 から 10 までの一様分布をします。 yxccreateY を通して何らかの関係を持つ数です。 xc が説明変数で y が応答変数です。

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

関係としては、以下の4つがあります。 xc と関係のある無しで4種類です。

function(x, c) rpois(dataSize, lambda = 4)
function(x, c) rpois(dataSize, lambda = exp(0.1 + 0.2 * x))
function(x, c) rpois(dataSize, lambda = exp(0.1 + 0.2 * c))
function(x, c) rpois(dataSize, lambda = exp(0.1 + 0.2 * x + 0.2 * c))

ポアソン回帰をやっているのは以下のコードです。 f でリンク関数を指定して、 genData によって生成されたデータ d に対してポアソン回帰をやっています。 推定されたパラメータを使いやすくするために、n という名前を付けています。

calculateAic <- function(d, genData, f, n) {
  fit <- glm(formula = f, family = poisson, data = d)
  co <- append(fit$coefficients, replicate(2, 0))
  names(co) <- n

このようにして得たパラメータを使って、最大対数尤度、平均対数尤度、バイアス、AICを以下のように計算しています。

  meanLogLik <- meanLogLikelihood(genData, co["b1"], co["b2"], co["b3"])
  bias <- logLik(fit) - meanLogLik
  aic <- -2 * (logLik(fit) - length(fit$coefficients))
  result <- c(logLik(fit), meanLogLik, bias, aic)
  names(result) <- c("logLik", "meanLogLik", "bias", "aic")
  result
}

これで1つのモデルに対して、1回分の実験が終わります。 以下のコードで実験を何度も繰り返して、結果を平均しています。

calculateMeanAic <- function(genData, f, n) {
  genData <- match.fun(genData)
  m <- replicate(expSize, calculateAic(genData(), genData, f, n))
  apply(m, 1, mean)
}

以下のコードで4種類のモデルに対して上記の計算をしています。 リンク関数は exp(b1)exp(b1 + b2 * x)exp(b1 + b3 * c)exp(b1 + b2 * x + b3 * c) の4種類です。

showAic <- function(genData) {
  r1 <- calculateMeanAic(genData, y ~ 1, c("b1","b2","b3"))
  r2 <- calculateMeanAic(genData, y ~ x, c("b1","b2","b3"))
  r3 <- calculateMeanAic(genData, y ~ c, c("b1","b3","b2"))
  r4 <- calculateMeanAic(genData, y ~ x + c, c("b1","b2","b3"))

後はサンプリングのための関数をこの関数に渡して、結果を表示しています。 得られた結果の解説は次回に回します。

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

この章ではバイアスやAICについて説明がされています。 データを説明するモデルとして様々なモデルが考えられる時に、データへのあてはまりの良さだけを基準にモデルを選ぶことはできません。 そこでAICを使ったモデル選択の様子を見るためのコードを用意しました。 コードはRで書きました。

dataSize <- 100
expSize <- 1000

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

createIndipendentData <- function() {
  createData(function(x, c) rpois(dataSize, lambda = 4))
}

createCorrelatedXData <- function() {
  createData(function(x, c) rpois(dataSize, lambda = exp(0.1 + 0.2 * x)))
}

createCorrelatedCData <- function() {
  createData(function(x, c) rpois(dataSize, lambda = exp(0.1 + 0.2 * c)))
}

createCorrelatedXCData <- function() {
  createData(function(x, c) rpois(dataSize, lambda = exp(0.1 + 0.2 * x + 0.2 * c)))
}

logLikelihood <- function(genData, b1, b2, b3) {
  genData <- match.fun(genData)
  data <- genData()
  prob <- vector(length = dataSize)
  for (j in 1:dataSize) {
    prob[j] <- log(dpois(data$y[j], lambda = exp(b1 + b2 * data$x[j] + b3 * data$c[j])))
  }
  sum(prob)
}

meanLogLikelihood <- function(genData, b1, b2, b3) {
  mean(replicate(dataSize, logLikelihood(genData, b1, b2, b3)))
}

calculateAic <- function(d, genData, f, n) {
  fit <- glm(formula = f, family = poisson, data = d)
  co <- append(fit$coefficients, replicate(2, 0))
  names(co) <- n

  meanLogLik <- meanLogLikelihood(genData, co["b1"], co["b2"], co["b3"])
  bias <- logLik(fit) - meanLogLik
  aic <- -2 * (logLik(fit) - length(fit$coefficients))
  result <- c(logLik(fit), meanLogLik, bias, aic)
  names(result) <- c("logLik", "meanLogLik", "bias", "aic")
  result
}

calculateMeanAic <- function(genData, f, n) {
  genData <- match.fun(genData)
  m <- replicate(expSize, calculateAic(genData(), genData, f, n))
  apply(m, 1, mean)
}

showAic <- function(genData) {
  r1 <- calculateMeanAic(genData, y ~ 1, c("b1","b2","b3"))
  r2 <- calculateMeanAic(genData, y ~ x, c("b1","b2","b3"))
  r3 <- calculateMeanAic(genData, y ~ c, c("b1","b3","b2"))
  r4 <- calculateMeanAic(genData, y ~ x + c, c("b1","b2","b3"))

  cat("\t\t\tnull\t\tx\t\tc\t\txc\n")
  cat("log likelihood\t\t", r1["logLik"], "\t", r2["logLik"], "\t",
    r3["logLik"], "\t", r4["logLik"], "\n")
  cat("mean log likelihood\t", r1["meanLogLik"], "\t",
     r2["meanLogLik"], "\t", r3["meanLogLik"], "\t", r4["meanLogLik"], "\n")
  cat("bias\t\t\t", r1["bias"], "\t", r2["bias"], "\t", r3["bias"], "\t", r4["bias"], "\n")
  cat("aic\t\t\t", r1["aic"], "\t", r2["aic"], "\t", r3["aic"], "\t", r4["aic"], "\n\n")
}

cat("indipendent data\n")
showAic(createIndipendentData)

cat("correlated X data\n")
showAic(createCorrelatedXData)

cat("correlated C data\n")
showAic(createCorrelatedCData)

cat("correlated XC data\n")
showAic(createCorrelatedXCData)

少し長くなったので、解説は次回に回します。 ぜひ実行してみてください。

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

この章ではバイアスやAICについて説明がされています。 データを説明するモデルとして様々なモデルが考えられる時に、データへのあてはまりの良さだけを基準にモデルを選ぶことはできません。 それを理解するために、ポアソン回帰のバイアスを見るためのコードを用意しました。 コードはRで書きました。

createData <- function() {
  y <- rpois(10, lambda = 4)
  data.frame(y)
}

logLikelihood <- function(data, b1) {
  sum(log(dpois(data, lambda = exp(b1))))
}

printLogLikelihood <- function(data, i) {
  cat(paste0("Input data is data", i, "\n"))
  print(data[i]$y)

  fit <- glm(formula = y ~ 1, family = poisson, data = data[i])
  b1 <- fit$coefficients[1]

  cat("Estimated b1 =", b1, "\n")
  cat("Log likelihood for\n")
  cat("data1\t\tdata2\t\tdata3\t\n")
  cat(logLikelihood(data[1]$y, b1), "\t",
    logLikelihood(data[2]$y, b1), "\t",
    logLikelihood(data[3]$y, b1), "\n\n")
}

printData <- function(data, i) {
  cat(paste0("data", i))
  print(data[i]$y)
}

data <- c(createData(), createData(), createData())

for(i in 1:3) {
  printData(data, i)
}
cat("\n")
for(i in 1:3) {
  printLogLikelihood(data, i)
}

コードを実行すると、まずポアソン分布からサンプリングした 10 個のデータからなるデータセットを3個作ります。 それぞれのデータに対してポアソン回帰をして、推定されたパラメータを使って他の2つのデータセットに対して対数尤度を計算します。 対数尤度を計算しているコードは以下のコードです。

logLikelihood <- function(data, b1) {
  sum(log(dpois(data, lambda = exp(b1))))
}

b1 が最尤推定されたパラメータです。 data が対数尤度を計算する対象のデータです。 平均 exp(b1) を持つポアソン分布からデータが得られる確率の対数を計算しています。 実行結果は以下のようになります。

data1 [1] 5 1 2 4 9 3 3 2 3 4
data2 [1] 7 4 3 3 3 8 4 6 3 1
data3 [1] 9 1 5 6 1 3 3 8 3 2

Input data is data1
 [1] 5 1 2 4 9 3 3 2 3 4
Estimated b1 = 1.280934 
Log likelihood for
data1       data2       data3   
-20.59338    -21.43294   -24.32331 

Input data is data2
 [1] 7 4 3 3 3 8 4 6 3 1
Estimated b1 = 1.435085 
Log likelihood for
data1       data2       data3   
-21.04396    -20.95861   -24.00313 

Input data is data3
 [1] 9 1 5 6 1 3 3 8 3 2
Estimated b1 = 1.410987 
Log likelihood for
data1       data2       data3   
-20.91147    -20.97071   -23.99113 

傾向として、ポアソン回帰に使ったデータセットへの尤度が一番大きくて、他のデータセットへの尤度は小さくなります。 時々は、他のデータセットへの尤度が大きくなります。 最尤推定は与えられたデータに最もよく当てはまるパラメータを返すので、まだ見ぬデータにもよく当てはまるかどうかは分からないということです。

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

この章ではバイアスやAICについて説明がされています。 データを説明するモデルとして様々なモデルが考えられる時に、データへのあてはまりの良さだけを基準にモデルを選ぶことはできません。 それを理解するために、ポアソン回帰で推定されるパラメータがどれくらい真の値からずれるのかをヒストグラムにしてくれるコードを用意しました。 コードはRで書きました。

createHistogram <- function(n) {
  lambdas <- vector(length = 1000)
  for (i in 1:1000) {
    y <- rpois(n, lambda = 4)
    data <- data.frame(y)

    fit <- glm(formula = y ~ 1, family = poisson, data = data)
    lambdas[i] <- exp(fit$coefficients[1])
  }
  title <- paste("Number of data =", n)
  hist(lambdas[lambdas <= 10], breaks = seq(0, 10, by = 0.01), main = title, xlab = "Estimeted lambda")
}

for (i in c(1, 5, 10, 100, 200)) {
  createHistogram(i)
}

実行すると平均 4 のポアソン分布から n 個のサンプルを取って、ポアソン回帰をするということを1000回繰り返して、推定結果のヒストグラムを描きます。 リンク関数は定数です。 n10 の時は下図のようになります。

estimation_number_of_data10.png

n100 の時は下図のようになります。

estimation_number_of_data100.png

推定結果が真の値である 4 の周囲に分布しています。 データの数が多くなれば分布のピークは 4 に集中していきますが、データが 100 個あっても、やはりずれる時はずれるのが分かります。

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

この章ではバイアスやAICについて説明がされています。 データを説明するモデルとして様々なモデルが考えられる時に、データへのあてはまりの良さだけを基準にモデルを選ぶことはできません。 それを理解するために、ポアソン回帰で推定されるパラメータがどれくらい真の値からずれるのかを見るためのコードを用意しました。 コードはRで書きました。

legends <- c("ideal", "regression")
ltys <- c("dashed", "solid")
for (i in 1:10) {
  y <- rpois(5, lambda = 4)
  data <- data.frame(y)

  fit <- glm(formula = y ~ 1, family = poisson, data = data)
  l <- exp(fit$coefficients[1])

  title <- paste("Poisson regression, estimated lambda = ", l)
  plot(1:10, dpois(1:10, lambda = 4), ylim = c(0.0, 0.25), type = "l", lty = "dashed", main = title, xlab = "", ylab = "")
  lines(1:10, dpois(1:10, lambda = l), type = "l")
  legend("topright", legend = legends, lty = ltys)
}

平均 4 のポアソン分布から5つのサンプルをサンプリングして、そのデータに対してポアソン回帰をしています。 リンク関数は定数です。 正確に推定されるなら、推定されるパラメータは 4 に近い値になるはずです。 実際に実行してみると、5個しかサンプルが無いわりには結構正確に推定されるのが分かります。 しかし、だいたい10回も推定すると、1回くらいは目に見えておかしな推定をするのが分かります。 例えば下図のような推定をします。

poisson_regresion_estimated4.png

真の値は 4 なのに、推定結果は 5.4 です。 このように、与えられたデータに対して最尤推定するしかない以上、推定は推定でしかなくどうしても真の値からずれてしまいます。

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

この章ではポアソン回帰について説明がされています。 ポアソン回帰によって得られるのは、説明変数と応答変数の分布の関係です。 変数が変化すると別の変数の分布が変化するという、結構複雑な関係が得られるわけです。 分布そのものにも興味はありますが、やはり一番興味があるのは、分布の中でどこが一番確率が高くなるかです。 そこで、本の中で使われているサンプルデータを使ってポアソン回帰をやってみて、得られた分布の中の一番確率が高い場所をプロットしてみました。 コードはRで書きました。

d <- read.csv("data3a.csv")

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

b1 <- fit$coefficients[1]
b2 <- fit$coefficients[2]

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

fitted <- function(x, y) {
  dpois(y, lambda = linkFunction(x))
}

cat("p(y | exp(a + bx)) where p is Poisson distribution\n")
cat("a =", b1, "b =", b2, "\n\n")

y <- 1:20
max <- numeric(length = 20)
for (x in 1:20) {
  at <- fitted(x, y)
  cat("when x = ", x, ", y = ", which.max(at), " is most probable, probability is ", max(at), "\n")
  max[x] <- which.max(at)
}
title <- paste0("p(y|exp(", b1, " + ", b2, "x))")
xl <- "x"
yl <- "the most probable y"
plot(1:20, max, type = "b", main = title, xlab = xl, ylab = yl)

実行すると以下のようなグラフが得られます。

most_probable.png

サンプルデータの説明変数 x と応答変数 y には正の相関があるので、 x の増加に従って y の分布は大きな値を取るようになります。 グラフを見ると一番高い確率で得られる y の値が x の増加に従って増えているのが分かります。 正の相関の一面が見られますね。

リンク関数は指数関数なので、もっときれいに指数関数的増加が見られないかとプロットする x の範囲を広げてみました。

most_probable2.png

細かい増加は分かりにくくなりましたが、全体的な指数関数的増加がよく分かるようになりました。

このアーカイブについて

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

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