読者です 読者をやめる 読者になる 読者になる

StatModeling Memorandum

StanとRとPythonでベイズ統計モデリングします. たまに書評.

分布から見た線形モデル・GLM・GLMM

久保さんのみどりぼん勉強会もせっかく催されていることだし、それにちなんだ記事を書きたいと思っていました。ここまでいい加減にGLMとGLMMをすっ飛ばして紹介して、さっさとBUGS/Stanのラビリンスパラダイスへいざないたいなぁという心境をスライドにしました。

僕が勉強し始めた頃、GLMやGLMMがとっつきにくく感じる時がありました。しかし今は、できあがる分布(と説明変数Xの値を変えた時どうなるか)をイメージすることがまずは大切と思っています。それに親しんでいれば自然と適用範囲も分かります。そしてモデルの数式(もしくはBUGS/Stanコード)を併せて見ながら「GLMMの場合、この項は個体差を考えていることに相当するんだな」などと理解できればよいと思っています。言いたいことはそれに尽きています。

スライドの最後から2番目のページの各解析へのリンクは以下になります。

最後にスライドを作成するにあたって使ったRコードとパラメータを載せます。

library(ggplot2)
library(gamlss.dist)

xaxis <- -50:50/10
amp <- 2
rot.m <- matrix(c(0, 1, -1, 0), 2, 2, byrow=T)

beta0 <- 1
sigma <- 1
beta1 <- 1
x.v <- c(0,2,4)
input <- data.frame(x=x.v, y=beta0+beta1*x.v)

df <- data.frame()
for(i in 1:nrow(input)){
   x0 <- xaxis + input[i,]$y
   y0 <- dNO(x0, mu=input[i,]$y, sigma=sigma)
   m <- as.matrix(data.frame(x=x0, y=amp*y0))
   z <- m %*% rot.m
   d <- data.frame(x=z[,1]+input[i,]$x, y=z[,2], group=i)
   df <- rbind(df, d)
}

df.point <- data.frame(x=x.v, y=input$y, group=0)

p <- ggplot(df, aes(x=x, y=y, group=group))
p <- p + geom_polygon(fill="#000000", alpha=3/5)
p <- p + geom_abline(intercept = beta0, slope = beta1, size=1.5)
p <- p + geom_point(data=df.point, size=5)
p <- p + theme(
   axis.title.x = element_text(size=16),
   axis.title.y = element_text(size=16, angle=0),
   axis.text.x = element_text(size=14),
   axis.text.y = element_text(size=14)
)
p <- p + scale_x_continuous(breaks=x.v)
ggsave(file = "output/lm.png", plot = p, dpi = 300, width = 4, height = 3)
  • 5行目: 分布の図の引き伸ばし具合です。見やすくなるように増幅しています。
  • 6行目: 分布の図を90度回転させる行列です。
  • 17行目: 分布自体はひいきにしている{gamlss.dist}dNOを使っています。
  • 27行目: geom_polygonで分布の図を描きます。
library(ggplot2)
library(gamlss.dist)

logistic <- function(z) 1 / (1 + exp(-z))
N <- 8
xaxis <- 0:N
amp <- 2
rot.m <- matrix(c(0, 1, -1, 0), 2, 2, byrow=T)

beta0 <- -4.0
beta1 <- 1.6
x.v <- c(0,2,4)
input <- data.frame(x=x.v, y=beta0+beta1*x.v)

df <- data.frame()
for(i in 1:nrow(input)){
   y0 <- dBI(xaxis, bd=N, mu=logistic(input[i,]$y))
   m <- as.matrix(data.frame(x=xaxis, y=amp*y0))
   z <- m %*% rot.m
   d <- data.frame(x=0, y=0, xmin=z[,1]+input[i,]$x, xmax=input[i,]$x, ymin=z[,2]-0.3, ymax=z[,2]+0.3, group=i)
   df <- rbind(df, d)
}

df.line <- data.frame(x=0:40/10, y=N*logistic(beta0+beta1*(0:40/10)), xmin=0, xmax=0, ymin=0, ymax=0, group=0)
df.point <- data.frame(x=x.v, y=N*logistic(input$y), xmin=0, xmax=0, ymin=0, ymax=0, group=0)

p <- ggplot(df, aes(x=x, y=y, xmin=xmin, xmax=xmax, ymin=ymin, ymax=ymax, group=group))
p <- p + geom_rect(fill="#F8766D", alpha=.9)
p <- p + geom_line(data=df.line, size=1.5)
p <- p + geom_point(data=df.point, size=5)
p <- p + theme(
   axis.title.x = element_text(size=16),
   axis.title.y = element_text(size=16, angle=0),
   axis.text.x = element_text(size=14),
   axis.text.y = element_text(size=14)
)
p <- p + scale_x_continuous(breaks=x.v)
ggsave(file = "output/glm_binomial.png", plot = p, dpi = 300, width = 4, height = 3)
  • 28行目: 連続型の分布ではなく離散型の分布の場合はgeom_rectで描きます。
library(ggplot2)
library(gamlss.dist)

logistic <- function(z) 1 / (1 + exp(-z))
N <- 8
xaxis <- 0:N
amp <- 2
rot.m <- matrix(c(0, 1, -1, 0), 2, 2, byrow=T)

f.gaussian.binom <- function(r, x, size, fixed, sd)
   dbinom(x=x, size=size, prob=logistic(fixed + r)) * dnorm(x=r, mean=0, sd=sd)

d.gaussian.binom <- function(xaxis, size, fixed, sd){
   sapply(xaxis, function(x)
      integrate(
         f = f.gaussian.binom, lower = -sd * 10, upper = sd * 10,
         # for f.gaussian.binom
         x = x, size = size, fixed = fixed, sd = sd
      )$value
   )
}

beta0 <- -4.0
beta1 <- 1.6
x.v <- c(0,2,4)
sd <- 3   # (0.1, 1, 3)
input <- data.frame(x=x.v, y=beta0+beta1*x.v)

df <- data.frame()
for(i in 1:nrow(input)){
   y0 <- d.gaussian.binom(xaxis, size=N, fixed=input[i,]$y, sd=sd)
   m <- as.matrix(data.frame(x=xaxis, y=amp*y0))
   z <- m %*% rot.m
   d <- data.frame(x=0, y=0, xmin=z[,1]+input[i,]$x, xmax=input[i,]$x, ymin=z[,2]-0.3, ymax=z[,2]+0.3, group=i)
   df <- rbind(df, d)
}

df.line <- data.frame(x=0:40/10, y=N*logistic(beta0+beta1*(0:40/10)), xmin=0, xmax=0, ymin=0, ymax=0, group=0)
df.point <- data.frame(x=x.v, y=N*logistic(input$y), xmin=0, xmax=0, ymin=0, ymax=0, group=0)

p <- ggplot(df, aes(x=x, y=y, xmin=xmin, xmax=xmax, ymin=ymin, ymax=ymax, group=group))
p <- p + geom_rect(fill="#F8766D", alpha=.9)
p <- p + geom_line(data=df.line, size=1.5)
p <- p + geom_point(data=df.point, size=5)
p <- p + theme(
   axis.title.x = element_text(size=16),
   axis.title.y = element_text(size=16, angle=0),
   axis.text.x = element_text(size=14),
   axis.text.y = element_text(size=14)
)
p <- p + scale_x_continuous(breaks=x.v)
ggsave(file = "output/glmm_binomial.png", plot = p, dpi = 300, width = 4, height = 3)
  • 11-21行目: みどりぼんのサポートサイトから拝借しました。GLMMの場合は数値積分が入ります。詳しくはみどりぼん参照。

poisson分布を使ったGLM・GLMMの時のパラメータは beta0=0.1, beta1=0.6 で {gamlss.dist} の dPO(xaxis, mu=exp(input[i,]$y))を使いました。gamma分布を使ったGLM・GLMMの時のパラメータは beta0=0.1, beta1=0.3, sigma=0.8 で {gamlss.dist}dGA(xaxis, mu=exp(input[i,]$y), sigma=sigma) を使いました。