TIM Labs

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

このブログ記事は『データ解析のための統計モデリング入門』(久保拓弥 著、岩波書店)という、とても分かりやすい統計モデリングの入門書を、さらに分かりやすくするために読書メモをまとめたものです。 今回は第8章、マルコフ連鎖モンテカルロ(MCMC)法とベイズ統計モデルについてのまとめの一回目です。

この章ではMCMCやメトロポリス法について説明がされています。 本の中ではMCMCの前に、まずはランダムにパラメータを更新していく手法が紹介されています。 なので、本の中で紹介されているような推定法を実際に動かしてみるコードを用意しました。 コードはRで書きました。

dataSize <- 20
N <- 8

data <- rbinom(dataSize, N, 0.45)
summary(data)
cat("Analytically, q =", sum(data) / (dataSize * N), "\n\n")

likelihood <- function (q) {
  sum(log(dbinom(data, N, q)))
}

develop <- function(q) {
  if (runif(1, 0, 1) < 0.5) {
    q + 0.01
  } else {
    q - 0.01
  }
}

q <- round(0.3 + runif(1, 0, 0.3), 2)
ql <- likelihood(q)
i <- 0
cat("start q =", q, "\n")
while(i < 30) {
  p <- develop(q)   
  pl <- likelihood(p)
  if (ql < pl) {
    q <- p
    ql <- pl
  }
  cat("q =", q, "likelihood", ql, "\n")
  i <- i + 1
}
cat("last q =", q, "\n")

サンプルデータとしては、本の中でやっているのと同じように、事象の発生確率が 0.45 で、事象の試行回数は 8 の二項分布から、20 個のサンプルを取って用意しました。 このデータに対して、事象の発生確率 q が分かっていないものとして推定していきます。

推定方法は q をランダムに変化させて、対数尤度を計算し、尤度が大きくなっていたら変化させた先の q に更新する、という方法です。 develop 関数の中で q0.01 足すか引くかして q を変化させています。 コードを実行すると次のような出力が得られます。

   Min. 1st Qu.  Median    Mean 3rd Qu.    Max.
      1       3       4       4       5       6
Analyticaly, q = 0.5

start q = 0.35
q = 0.36 likelihood -41.69397
q = 0.37 likelihood -40.76192
q = 0.38 likelihood -39.90848
q = 0.39 likelihood -39.13129
q = 0.39 likelihood -39.13129
q = 0.39 likelihood -39.13129
q = 0.39 likelihood -39.13129
q = 0.4 likelihood -38.42821
q = 0.4 likelihood -38.42821
q = 0.41 likelihood -37.79737
q = 0.42 likelihood -37.23712
q = 0.42 likelihood -37.23712
q = 0.43 likelihood -36.74602
q = 0.44 likelihood -36.32282
q = 0.45 likelihood -35.96647
q = 0.45 likelihood -35.96647
q = 0.45 likelihood -35.96647
q = 0.46 likelihood -35.67609
q = 0.47 likelihood -35.45097
q = 0.48 likelihood -35.29055
q = 0.49 likelihood -35.19445
q = 0.5 likelihood -35.16245
q = 0.5 likelihood -35.16245
q = 0.5 likelihood -35.16245
q = 0.5 likelihood -35.16245
q = 0.5 likelihood -35.16245
q = 0.5 likelihood -35.16245
q = 0.5 likelihood -35.16245
q = 0.5 likelihood -35.16245
q = 0.5 likelihood -35.16245
last q = 0.5

この例の場合は、最尤推定値を解析的に求めることができるので、まず最初にその値を出力しています。 q をランダムに変化させていくと、最終的にはその値に近づいていくのが分かります。

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

この章では複数の分布を混ぜて使う、一般化線形混合モデルついて説明がされています。 なので、サンプルデータを用意して、実際に一般化線形混合モデルを最尤推定してみるコードを用意しました。 コードはRで書きました。

library('glmmML')

dataSize <- 100
N <- 20

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

createData <- function() {
  x <- runif(dataSize, 0, 5)
  y <- vector(length = dataSize)
  for (i in 1:dataSize) {
    r <- rnorm(1, mean = 0, sd = 2)
    l <- linkFunction(-2, 1, x[i])
    y[i] <- rbinom(1, N, l(r))
  }
  id <- 1:dataSize
  data.frame(x, y, id)
}

data <- createData()
glmmML(cbind(y, N - y) ~ x, data = data, family = binomial, cluster = id)

コードを実行するには、glmmMLパッケージをインストールする必要があります。

リンク関数は logit(q) = b1 + b2 * x + r です。 r は平均 0 、標準偏差 s の正規分布から生成されます。

サンプルデータは、0 から 5 までの一様分布から x を生成し、標準偏差 s = 2 の正規分布から r を生成し、b1 = -2b2 = 1 のリンク関数を使って事象の発生確率 q を求め、試行回数 N = 20 の二項分布から値を生成しています。 なので、最尤推定すると、b1 = -2b2 = 1s = 2 に近い値を推定してくれるはずです。

コードを実行すると以下のような出力が出力されます。

Call:  glmmML(formula = cbind(y, N - y) ~ x, family = binomial, data = data,      cluster = id)


              coef se(coef)      z Pr(>|z|)
(Intercept) -2.163   0.4570 -4.732 2.22e-06
x            1.210   0.1604  7.539 4.73e-14

Scale parameter in mixing distribution:  2.144 gaussian
Std. Error:                              0.1654

        LR p-value for H_0: sigma = 0:  5.114e-131

Residual deviance: 311.6 on 97 degrees of freedom       AIC: 317.6

期待通りにパラメータが推定されているのが分かります。

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

この章では複数の分布を混ぜて使う、一般化線形混合モデルついて説明がされています。 その中で二項分布とと正規分布を混ぜて使う一般化線形混合モデルについて説明されています。 なので、サンプルデータを用意して、様々なパラメータにおける対数尤度の値をプロットしてくれるコードを用意しました。 コードはRで書きました。

dataSize <- 50
N <- 20

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

dmixture <- function(b1, b2, s) {
  function(y, x) {
    l <- linkFunction(b1, b2, x)
    sum <- 0
    for (r in seq(-100, 100, by = 0.1)) {
      sum <- sum + dbinom(y, N, l(r)) * dnorm(r, mean = 0, sd = s) * 0.1
    }
    sum
  }
}

likelihood <- function(b1, b2, s) {
  d <- dmixture(b1, b2, s)
  function (y, x) {
    d(y, x)
  }
}

logLikelihood <- function(data, b1, b2) {
  function (s) {
    sl <- 0
    ll <- likelihood(b1, b2, s)
    for (i in 1:dataSize) {
      sl <- sl + log(ll(data$y[i], data$x[i]))
    }
    sl
  }
}

createData <- function() {
  x <- runif(dataSize, 0, 5)
  y <- vector(length = dataSize)
  for (i in 1:dataSize) {
    r <- rnorm(1, mean = 0, sd = 1)
    l <- linkFunction(-2, 1, x[i])
    y[i] <- rbinom(1, N, l(r))
  }
  data.frame(x, y)
}

data <- createData()
xl <- "s"
yl <- "Lilelihood"
ss <- seq(0.2, 2, by = 0.2) 
b2s <- seq(-2, 2, by = 1)
legends <- paste0("b2 = ", b2s)
for (b1 in seq(-2, 2, by = 1)) {
  title <- paste0("logit(q) = ", b1, " + b2 * x + r, r ~ N(0, s)")
  plot(0, 0, type = "n", xlim = c(0, 2), ylim = c(-1000, 0), main = title, xlab = xl, ylab = yl)
  for (i in 1:5) {
    ls <- logLikelihood(data, b1, b2s[i])(ss)
    lines(ss, ls, type = "l")
    points(ss, ls, pch = i)
  }
  legend("topleft", legend = legends, pch = 1:5)
}

リンク関数は logit(q) = b1 + b2 * x + r です。 r は平均 0 、標準偏差 s の正規分布から生成されます。

サンプルデータを生成するには、まず 0 から 5 までの一様分布から x を、標準偏差 s = 1 の正規分布から r を生成します。 その後、b1 = -2b2 = 1 のリンク関数を使って事象の発生確率 q を求め、試行回数 N = 20 の二項分布から値を生成しています。 なので、様々なパラメータについてモデルを試すと、b1 = -2b2 = 1s = 1 の時に最もデータへの当てはまりがよくなるはずです。 つまり尤度が最大になるはずです。

コードを実行すると5枚のグラフがプロットされます。 例えば次のようなグラフが得られます。

mixture_likelihood.png

図を見ると、期待通りに b1 = -2b2 = 1s = 1 の場合に尤度が最大になっているのが分かります。 実際の問題では、これらのパラメータを最尤推定すればいいわけです。

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

この章では複数の分布を混ぜて使う、一般化線形混合モデルついて説明がされています。 その中でポアソン分布と正規分布を混ぜた分布について説明されています。 なので、ポアソン分布と正規分布を混ぜ合わせた混合分布をグラフにしてくれるコードを用意しました。 コードはRで書きました。

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

distribution <- function(s, linkFunction) {
  l <- match.fun(linkFunction)
  function(y) {
    sum <- 0
    for (r in seq(-100, 100, by = 0.1)) {
      sum <- sum + dpois(y, l(r)) * dnorm(r, mean = 0, sd = s) * 0.1
    }
    sum
  }
}

xl <- "y"
yl <- "Probability"
xs <- 1:5
ys <- 0:10
b1 <- 0.1
legends = paste0("x = ", xs)
for (b2 in seq(0.1, 0.5, by = 0.1)) {
  for (s in c(0.5, 1.0, 1.5, 2.0)) {
    title <- paste0("lambda = exp(", b1, " + ", b2, " * x + r, r ~ N(0, ", s, ")")
    plot(0, 0, type = "n", xlim = c(0, 10), ylim = c(0.0, 0.5), main = title, xlab = xl, ylab = yl)
    for (i in 1:5) {
      l <- linkFunction(b1, b2, xs[i])
      d <- distribution(s, l)
      lines(ys, d(ys), type = "l")
      points(ys, d(ys), pch = i)
    }
    legend("topright", legend = legends, pch = 1:5)
  }
}

リンク関数は λ = exp(b1 + b2 * x + r) です。 r は平均 0 、標準偏差 s の正規分布から生成されます。 r を生成する正規分布の標準偏差が0に近い数の時は、これはポアソン回帰のリンク関数とほとんど同じです。 混合分布はポアソン分布と正規分布をかけて r で積分したものです。

コードを実行すると20枚のグラフがプロットされます。 例えば次のようなグラフが得られます。

mixture_poisson.png

この図は b1 = 0.1b2 = 0.5s = 1 の場合です。 ポアソン分布よりも広く分布した分布が得られているのが分かります。

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

この章では複数の分布を混ぜて使う、一般化線形混合モデルついて説明がされています。 その中でポアソン分布と正規分布を混ぜた分布について説明されています。 なので、ポアソン分布と正規分布を掛けた値をグラフにしてくれるコードを用意しました。 コードはRで書きました。

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

distribution <- function(y, s, linkFunction) {
  l <- match.fun(linkFunction)
  function(r) {
    dpois(y, l(r)) * dnorm(r, mean = 0, sd = s)
  }
}

xl <- "r"
yl <- "Probability"
ys <- seq(0, 8, by = 2)
rs <- seq(-5, 5, by = 0.5)
b1 <- 0.1
legends = paste0("y = ", ys)
for (b2 in seq(0.1, 0.5, by = 0.1)) {
  for (s in c(0.5, 1.0, 1.5, 2.0)) {
    l <- linkFunction(b1, b2, 2)
    title <- paste0("lambda  = exp(", b1, " + ", b2, " * 2 +  r), r ~ N(0, ", s, ")")
    plot(0, 0, type = "n", xlim = c(-5, 5), ylim = c(0.0, 0.2), main = title, xlab = xl, ylab = yl)
    for (i in 1:5) {
      d <- distribution(ys[i], s, l)
      lines(rs, d(rs), type = "l")
      points(rs, d(rs), pch = i)
    }
    legend("topright", legend = legends, pch = 1:5)
  }
}

リンク関数は λ = exp(b1 + b2 * x + r) です。 r は平均 0、標準偏差 s の正規分布から生成される変数です。 あとは、正規分布から r が生成される確率と、平均 λ のポアソン分布を掛け合わせれば、欲しい値になります。

コードを実行すると20枚のグラフがプロットされます。 例えば次のようなグラフが得られます。

mixed_poission.png

y によってグラフの山の大きさが変わってくるのが分かります。 混合分布はこれを r で積分したものです。 なので、グラフ上で面積の大きな y が確率の高い分布になります。

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

この章では複数の分布を混ぜて使う、一般化線形混合モデルついて説明がされています。 なので、二項分布と正規分布を混ぜ合わせた混合分布をグラフにしてくれるコードを用意しました。 コードはRで書きました。

N <- 20

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

distribution <- function(s, linkFunction) {
  l <- match.fun(linkFunction)
  function(y) {
    sum <- 0
    for (r in seq(-100, 100, by = 0.1)) {
      sum <- sum + dbinom(y, N, l(r)) * dnorm(r, mean = 0, sd = s) * 0.1
    }
    sum
  }
}

xl <- "y"
yl <- "Probability"
xs <- c(-5,-2,0,2,5) 
ys <- 0:20
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)) {
    for (s in c(0.1, 0.5, 1, 1.5, 2)) {
      title <- paste0("logit(q) = ", b1, " + ", b2, " * x + r, r ~ N(0, ", s, ")")
      plot(0, 0, type = "n", xlim = c(0, 20), ylim = c(0.0, 0.5), main = title, xlab = xl, ylab = yl)
      for (i in 1:5) {
        l <- linkFunction(b1, b2, xs[i])
        d <- distribution(s, l)
        lines(ys, d(ys), type = "l")
        points(ys, d(ys), pch = i)
      }
      legend("topright", legend = legends, pch = 1:5)
    }
  }
}

リンク関数は logit(q) = b1 + b2 * x + r です。 r は平均 0 、標準偏差 s の正規分布から生成されます。 r を生成する正規分布の標準偏差が0に近い数の時は、これはロジスティック回帰のリンク関数とほとんど同じです。

混合分布は二項分布と正規分布をかけて r で積分したものです。 今は単純に r-100 から 100 まで 0.1 刻みで変化させて、面積の合計を求めておきました。 コードで言うと以下の部分です。

    for (r in seq(-100, 100, by = 0.1)) {
      sum <- sum + dbinom(y, N, l(r)) * dnorm(r, mean = 0, sd = s) * 0.1
    }

これで様々な b1b2s について混合分布のグラフをプロットできます。 コードを実行すると180枚のグラフがプロットされます。 例えば次のようなグラフが得られます。

mixture_distribution.png

この図は b1 = 0.5b2 = 0.5s = 1 の場合です。 二項分布とはたいぶ違った分布が得られているのが分かります。

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

この章では複数の分布を混ぜて使う、一般化線形混合モデルついて説明がされています。 なので、二項分布と正規分布を掛けた値をグラフにしてくれるコードを用意しました。 コードはRで書きました。

N <- 20

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

distribution <- function(y, linkFunction) {
  l <- match.fun(linkFunction)
  function(r) {
    dbinom(y, N, l(r)) * dnorm(r, mean = 0, sd = 2)
  }
}

xl <- "r"
yl <- "Probability"
ys <- seq(0, 20, by = 5)
rs <- seq(-5, 5, by = 0.5)
legends = paste0("y = ", ys)
for (b1 in seq(-0.5, 0.5, by = 0.2)) {
  for (b2 in seq(-0.5, 0.5, by = 0.2)) {
    for (x in 0:3) {
      l <- linkFunction(b1, b2, x)
      title <- paste0("logit(q) = ", b1, " + ", b2, " * ", x, " + r, r ~ N(0, 2)")
      plot(0, 0, type = "n", xlim = c(-5, 5), ylim = c(0.0, 0.1), main = title, xlab = xl, ylab = yl)
      for (i in 1:5) {
        d <- distribution(ys[i], l)
        lines(rs, d(rs), type = "l")
        points(rs, d(rs), pch = ys[i])
      }
      legend("topright", legend = legends, pch = ys)
    }
  }
}

リンク関数は logit(q) = b1 + b2 * x + r です。 r は平均 0、標準偏差 2 の正規分布から生成される変数です。 事象の試行回数は 20 にしました。 リンク関数を使って二項分布する事象の発生確率 q が求まれば、事象が y 回起こる確率を計算できます。 あとは、正規分布から r が生成される確率と、二項分布を掛け合わせれば、欲しい値になります。

様々な b1b2x を与えて、r に対して確率をプロットしてみました。 コードを実行すると144枚のグラフがプロットされます。 例えば次のようなグラフが得られます。

mixed_2.png

y によってグラフの山の大きさが変わってくるのが分かります。 混合分布はこれを r で積分したものです

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

この章では複数の分布を混ぜて使う、一般化線形混合モデルついて説明がされています。 まずは二項分布と正規分布を混ぜた混合分布について解説されています。 なので、混合分布をグラフにしてみたいところですが、やっぱり変数がたくさんあって複雑なので、まずは感覚をつかむために数値を羅列してくれるコードを用意しました。 コードはRで書きました。

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

l <- linkFunction(-0.5, 0.3, 2)

cat("r\ty=0\t\ty=1\t\ty=2\t\ty=3\n")
sum <- vector(length = 4)
for (r in seq(-5, 5, by = 0.5)) {
  q <- l(r)
  cat(r, "\t")
  for (y in 0:3) {
    b <- dbinom(y, 3, q)
    n <- dnorm(r, mean = 0, sd = 2)
    mixture <- b * n
    sum[y + 1] <- sum[y + 1] + mixture * 0.5
    cat(mixture, "\t")
  }
  cat("\n")
}
cat("sum\t", sum[1], "\t", sum[2], "\t", sum[3], "\t", sum[4], "\n")

リンク関数は xr を変数に含んでいます。 x2 に固定した場合について計算しました。 なので、あとは r を与えれば、二項分布している事象の発生確率 q が求まります。 r は正規分布する変数です。 今回はr は平均 0、標準偏差 2 の正規分布から得られるとしました。 事象の試行回数は 3 にして、0 回から 3 回事象が発生する場合の確率を計算して出力しています。 r-5 から 5 まで、0.5 刻みで計算しました。 コードを実行すると、次の出力が得られます。

r       y=0             y=1             y=2             y=3
-5      0.008571241     0.0001914794    1.425867e-06    3.539279e-09    
-4.5    0.01529937      0.0005635068    6.918364e-06    2.831304e-08    
-4      0.02542036      0.00154367      3.124683e-05    2.108318e-07    
-3.5    0.03909264      0.003913947     0.0001306212    1.453086e-06    
-3      0.05514584      0.009102904     0.0005008711    9.186514e-06    
-2.5    0.07038014      0.01915423      0.001737632     5.254481e-05    
-2      0.07963943      0.03573468      0.005344786     0.0002664708    
-1.5    0.07772425      0.05749969      0.01417925      0.00116552  
-1      0.06325714      0.0771553       0.031369        0.004251228     
-0.5    0.04148674      0.08342818      0.05592358      0.01249557  
0       0.02138051      0.07088734      0.07834263      0.02886066  
0.5     0.008601664     0.04701976      0.08567559      0.05203704  
1       0.002741934     0.02471168      0.07423798      0.07434107  
1.5     0.0007137071    0.01060504      0.05252712      0.08672284  
2       0.0001570974    0.003848653     0.03142876      0.08555086  
2.5     3.018185e-05    0.001219082     0.0164134       0.07366188  
3       5.187402e-06    0.0003454491    0.007668262     0.0567399   
3.5     8.116416e-07    8.911395e-05    0.003261413     0.03978732  
4       1.169644e-07    2.1173e-05      0.001277585     0.02569661  
4.5     1.564146e-08    4.66824e-06     0.0004644167    0.01540073  
5       1.950226e-09    9.596395e-07    0.0001574019    0.008605787     
sum      0.2548242       0.2235203       0.2303399      0.2828235 

r の増加に従って q も増加するので、r が小さいうちは事象の発生回数は少ない方が確率が高く、r が大きくなると多い方が確率が高くなるのが分かります。 混合分布では、r についての積分がとられるので、最後の行にざっくりとした積分の近似値を出力しておきました。 二項分布と正規分布の混合分布は、だいたいこういう数字になります。

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

この章では複数の分布を混ぜて使う、一般化線形混合モデルついて説明がされています。 なので、複数の分布を混ぜた分布をグラフにしてみたいところです。 ですが、変数の数が多くて複雑なので、まずは感覚をつかむために数値を羅列してくれるコードを用意しました。 コードはRで書きました。

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

l <- linkFunction(-0.5, 0.3)

cat("r\tq\t\tdbinom\t\tdnorm\t\tmixed\n")
for (r in seq(-1, 1, by = 0.1)) {
  q <- l(0.8, r)
  b <- dbinom(7, 20, q)
  n <- dnorm(r, mean = 0, sd = 1)
  mixed <- b * n
  cat(r, "\t", q, "\t", b, "\t", n, "\t", mixed, "\n")
}

このコードは二項分布と正規分布を掛けて得られる数値を出力します。 二項分布のパラメータである、事象の発生確率 q はリンク関数から得られます。 リンク関数には xr という変数があって、r の方が正規分布から生成されます。 計算したのは、dbinom(7, 20, q)dnorm(r, 0, 1) の積です。 コードを実行すると次の出力が得られます。

r       q           dbinom          dnorm       mixed
-1       0.2209739   0.07762054      0.2419707   0.0187819 
-0.9     0.2386673   0.0987262       0.2660852   0.02626959 
-0.8     0.2573095   0.1210895       0.2896916   0.03507862 
-0.7     0.2768782   0.1429499       0.3122539   0.04463668 
-0.6     0.2973393   0.1621339       0.3332246   0.05402699 
-0.5     0.3186463   0.1763681       0.3520653   0.06209308 
-0.4     0.3407396   0.1837016       0.3682701   0.06765181 
-0.3     0.3635475   0.1829363       0.3813878   0.06976967 
-0.2     0.3869858   0.1739386       0.3910427   0.06801741 
-0.1     0.4109596   0.1577226       0.3969525   0.06260837 
0        0.4353637   0.1362616       0.3989423   0.05436053 
0.1      0.4600851   0.1120753       0.3969525   0.04448858 
0.2      0.4850045   0.08771599      0.3910427   0.0343007 
0.3      0.5099987   0.06530726      0.3813878   0.02490739 
0.4      0.5349429   0.04625377      0.3682701   0.01703388 
0.5      0.5597136   0.03116972      0.3520653   0.01097378 
0.6      0.5841905   0.01999504      0.3332246   0.00666284 
0.7      0.608259    0.01221858      0.3122539   0.003815299 
0.8      0.6318124   0.007119175     0.2896916   0.002062365 
0.9      0.6547535   0.003959494     0.2660852   0.001053563 
1        0.6769959   0.00210485      0.2419707   0.000509312 

r に対して q は単調に増加しています。 二項分布の部分で計算しているのは、事象が 20 回中 7 回起きる確率なので、q が小さい時は小さい値を取ります。 q が増加すると、だんだん確率は上がりますが、増加し過ぎると 7 回より多い回数、事象が発生しやすくなるので、あるところで最大をとって、後は下がります。 正規分布は 0 を中心とした山なりの分布なので、r0 に近い値なら確率は高く、そこから離れるほど小さくなります。

この二つを掛け合わせた値は、だいたいどっちも大きな値を取るときに大きいですが、どちらの最大とも微妙にずれたところに最大がきているのが分かります。

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

この章では複数の分布を混ぜて使う、一般化線形混合モデルついて説明がされています。 まずはロジスティック回帰にパラメータを追加したモデルについて説明されています。 このモデルでは、リンク関数に新しいパラメータが追加されて、そのパラメータは正規分布から生成されます。 なので、いろいろなパラメータについて、その値をプロットしてくれるコードを用意しました。 コードはRで書きました。

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

xl <- "x"
yl <- "q"
xs <- runif(100, 0, 10)
ys <- function(r) {
  linkFunction(-0.5, 0.3)(xs, r)
}
for (sd in seq(0.1, 2.0, by = 0.1)) {
  rs <- rnorm(100, mean = 0, sd = sd)
  ps <- data.frame(x = xs, y = ys(rs))
  title <- paste0("logit(q) = -0.5 + 0.3 * x + r, r ~ N(0, ", sd, ")")
  plot(0, 0, type = "n", xlim = c(0, 10), ylim = c(0.0, 1.0), main = title, xlab = xl, ylab = yl)
  points(ps, pch = 1)
}

コードを実行すると図が20枚プロットされます。 横軸はリンク関数に使われている説明変数 x で、縦軸は事象の発生確率 q です。 q には x の他にも、平均 0 の正規分布から生成された乱数の影響があります。 この乱数が説明変数として考慮されていない影響を表します。 例えば以下のような図がプロットされます。

logit_2.png

この図は分散が 0.2 の正規分布から生成された乱数がリンク関数を通して q に影響している場合です。 x の増加に対して q が緩やかに上昇する傾向が見て取れます。 これは説明変数以外の要因の影響が少ない場合です。 一方、以下のような図もプロットされます。

logit_10.png

この図は分散が 1.0 の場合です。 xq の関係が分かりにくくなりました。 これは説明変数以外の要因の影響が多い場合です。 このように、パラメータを生成する正規分布の分散によって、説明変数として考慮されていない要因の影響が大きい場合や小さい場合を表現できるのが分かります。

このアーカイブについて

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

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