TIM Labs

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

このブログ記事は『データ解析のための統計モデリング入門』(久保拓弥 著、岩波書店)という、とても分かりやすい統計モデリングの入門書を、さらに分かりやすくするために読書メモをまとめたものです。 今回は第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 の関係が分かりにくくなりました。 これは説明変数以外の要因の影響が多い場合です。 このように、パラメータを生成する正規分布の分散によって、説明変数として考慮されていない要因の影響が大きい場合や小さい場合を表現できるのが分かります。

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

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

N <- 20

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

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

rs <- -2:2
ys <- 0:N
pchs <- 1:5
xl <- "y"
yl <- "Probability"
legends <- paste0("r = ", rs)
for (b1 in c(-0.5, 0.5, by = 0.2)) {
  for (b2 in c(-0.5, 0.5, by = 0.2)) {
    l <- linkFunction(b1, b2)
    title <- paste0("logit(q) = ", b1, " + ",b2," * 2"," + r")
    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) {
      d <- distribution(2, rs[i], l)
      lines(ys, d(ys), type = "l")
      points(ys, d(ys), pch = pchs[i])
    }
    legend("topright", legend = legends, pch = pchs)
  }
}

このモデルでは、二項分布のパラメータである事象の発生確率 q が以下のコードで与えられます。

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

GLMのモデルと比べると、新しい変数として r が追加されています。 この r が説明変数 x とは違った、何とは特定されていない、余分な事象の発生確率への影響を表します。

コードを実行すると様々な b1b2r の値に対してモデルをプロットします。 x2 で固定にしました。 事象の試行回数は N = 20 で、事象の発生回数は y です。 例えば下図のような図がプロットされます。

mixed.png

図を見ると、確かに r の増減に従って y の分布が影響を受けるのが分かります。

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

この章では様々なタイプのGLMについて説明がされています。 その中でガンマ回帰について説明がされています。 なので、サンプルデータを用意して、ガンマ回帰をやってみるコードを用意しました。 コードはRで書きました。

dataSize <- 1000
createData <- function(createY) {
  createY <- match.fun(createY)
  x <- runif(dataSize, 0.1, 0.9)
  y <- createY(x)
  data.frame(x, y)
}

dataGenerator <- function()
  createData(function(x) rgamma(dataSize, shape = 4, rate = 4 / exp(0.5 + log(x))))

data <- dataGenerator()
fit <- glm(y ~ log(x), data = data, family = Gamma(link = "log"))
b1 <- fit$coefficients[1]
b2 <- fit$coefficients[2]
phi <- summary(fit)$dispersion

cat("b1", b1, "\n")
cat("b2", b2, "\n")
cat("shape", 1 / phi, "\n")
cat(paste0("rate = ", round(1 / phi, 3), " / exp(", round(b1, 3), " + ", round(b2, 3), " * log(x))\n"))

rate <- function(x) {
  1 / (phi * exp(b1 + b2 * log(x)))
}

xs <- c(0.1, 1.0, 2.0, 3.0, 4.0)
ys <- 1:10
title <- "Gamma regression"
xl <- "y"
yl <- "Probability"
legends <- paste0("x = ", xs)
plot(0, 0, type = "n", xlim = c(1, 10), ylim = c(0.0, 1.0), main = title, xlab = xl, ylab = yl)
for (x in xs) {
  y <- dgamma(ys, shape = 1 / phi, rate = rate(x))
  lines(ys, y, type="l")
  points(ys, y, pch = x)
}
legend("topleft", legend = legends, pch = xs)

ガンマ回帰のリンク関数は平均を与えるだけで、ガンマ分布のパラメータそのものは得られません。 なので、推定された分布を計算するには、もうひと工夫必要になります。 どうすればいいかはサポートページに書いてありました。 以下のコードがガンマ回帰でリンク関数のパラメータと、ガンマ分布のパラメータを決定するために必要なパラメータを計算している部分です。

fit <- glm(y ~ log(x), data = data, family = Gamma(link = "log"))
b1 <- fit$coefficients[1]
b2 <- fit$coefficients[2]
phi <- summary(fit)$dispersion

phi はshapeパラメータと shape = 1 / phi という関係を持つパラメータです。 shapeパラメータが得られれば、rateパラメータは以下のように得られます。

rate <- function(x) {
  1 / (phi * exp(b1 + b2 * log(x)))
}

これで推定したガンマ分布をプロットすることができるようになりました。 コードを実行すると以下のような図がプロットされます。

gamma_regression.png

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

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

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

distribution <- function(x, s, linkFunction) {
  l <- match.fun(linkFunction)
  function(y) {
    dgamma(y, shape = s, rate = s / l(x))
  }
}

xs <- c(0.1, 1.0, 2.0, 3.0, 4.0)
ys <- seq(0.5, 10.0, by = 0.5)
pchs <- 1:5
xl <- "y"
yl <- "Probability"
legends <- paste0("x = ", xs)
for (b1 in seq(0.5, 2.5, by = 0.5)) {
  for (b2 in seq(0.5, 2.5, by = 0.5)) {
    for (s in 1:5) {
      l <- linkFunction(b1, b2)
      title <- paste0("shape = ", s, ", rate = ", s, " / exp(", b1, " + ", b2, "log(x))")
      plot(0, 0, type = "n", xlim = c(0, 10), ylim = c(0.0, 1.0), main = title, xlab = xl, ylab = yl)
      for (i in 1:5) {
        d <- distribution(xs[i], s, l)
        lines(ys, d(ys), type = "l")
        points(ys, d(ys), pch = pchs[i])
      }
      legend("topright", legend = legends, pch = pchs)
    }
  }
}

リンク関数は x の単調増加関数です。 コードでは以下の部分です。

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

これがガンマ分布の平均を与えます。 ガンマ分布にはshapeパラメータとrateパラメータがありますが、平均は shape / rate なのでリンク関数のパラメータだけを与えても、どんなガンマ分布になるか分かりません。 なので、リンク関数のパラメータ b1b2 の他にshapeパラメータも与えました。

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

gamma_model_4.png

これは b1 = 1b2 = 0.5shape = 4 の場合です。 ガンマ分布の平均は shape / rate = exp(1 + 0.5 * log(x)) なので、rate = 4 / exp(1 + 0.5 * log(x)) になります。 平均が x の単調増加関数なので、x が大きな値を取るようになるほど y の分布も大きな値を取るように変化するのが分かります。

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

このアーカイブについて

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

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