異なる条件のデータに対するフィッティングカーブ重ね書きです。
仮想データは二項分布、その確率が平均と分散の異なるガウス分布に従います（リンク関数がプロビット）。 まあ心理物理でよくあるデータです。今回はガンマとラムダは無視。
set.seed(42)
library(plyr)
d <- expand.grid(C = 1:3, x = sample(10, 1000, TRUE))
m <- c(3, 5, 7)
s <- 1:3
d <- ddply(d, .(C, x),
function(a) data.frame(r = rbinom(nrow(a), 1, prob=pnorm(a$x, m[a$C], s[a$C]))))
ベースグラフィックスで真のカーブを描いてみます。
for (i in 1:3) {
curve(pnorm(x, m[i], s[i]), 0, 10, col =i, add=if (i == 1) FALSE else TRUE)
}
ggplot2でデータの平均をプロットしてみます。
library(ggplot2)
p1 <- ggplot(d, aes(x, r, fill = factor(C), colour = factor(C))) +
stat_summary(fun.y = mean, geom = "point")
p1 + stat_summary(fun.y = mean, geom = "line")
フィッティングしてパラメータを取り出します。
fit <- dlply(d, .(C), function(a) {
f <- glm(r~x, data=a, family=binomial(link="probit"))
c0<-coef(f)[1]
c1<-coef(f)[2]
f$mu<- unname(-c0/c1)
f$sigma<- unname(1/abs(c1))
f
})
params <- ldply(fit, function(x) data.frame(m = x$mu, s = x$sigma))
paramsは条件ごとのパラメータ推定値です。
params
## C m s
## 1 1 2.875 0.9484
## 2 2 4.856 1.9371
## 3 3 6.651 2.9699
このパラメータを使ってフィッティングカーブを書きます。
stat_functionを使います。
p1 +
dlply(params, .(C), function(a)
stat_function(fun = pnorm, args = list(mean = a$m, sd = a$s),
data = data.frame(C = a$C),
mapping = aes(x = NULL, y = NULL)))
