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

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

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

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

スライドの最後から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) を使いました。

StatModeling Memorandum

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

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