第 67 章 贝叶斯logistic-binomial模型
library(tidyverse)
library(tidybayes)
library(rstan)
rstan_options(auto_write = TRUE)
options(mc.cores = parallel::detectCores())
theme_set(bayesplot::theme_default())67.1 企鹅案例
筛选出物种为”Gentoo”的企鹅，并构建gender变量，male 对应1，female对应0
library(palmerpenguins)
gentoo <- penguins %>%
filter(species == "Gentoo", !is.na(sex)) %>%
mutate(gender = if_else(sex == "male", 1, 0))
gentoo## # A tibble: 119 × 9
## species island bill_length_mm bill_depth_mm flipper_length_mm body_mass_g
## <fct> <fct> <dbl> <dbl> <int> <int>
## 1 Gentoo Biscoe 46.1 13.2 211 4500
## 2 Gentoo Biscoe 50 16.3 230 5700
## 3 Gentoo Biscoe 48.7 14.1 210 4450
## 4 Gentoo Biscoe 50 15.2 218 5700
## 5 Gentoo Biscoe 47.6 14.5 215 5400
## 6 Gentoo Biscoe 46.5 13.5 210 4550
## 7 Gentoo Biscoe 45.4 14.6 211 4800
## 8 Gentoo Biscoe 46.7 15.3 219 5200
## 9 Gentoo Biscoe 43.3 13.4 209 4400
## 10 Gentoo Biscoe 46.8 15.4 215 5150
## # … with 109 more rows, and 3 more variables: sex <fct>, year <int>,
## # gender <dbl>67.1.1 dotplots
借鉴ggdist的Logit dotplots 的画法，画出dotplot
gentoo %>%
ggplot(aes(x = body_mass_g, y = sex, side = ifelse(sex == "male", "bottom", "top"))) +
geom_dots(scale = 0.5) +
ggtitle(
"geom_dots(scale = 0.5)",
'aes(side = ifelse(sex == "male", "bottom", "top"))'
)
\[ \begin{align*} y_i & = \text{bernoulli}( p_i) \\ p_i & =\text{logit}^{-1}(X_i \beta) \end{align*} \]
67.1.2 bayesian logit模型
stan_program <- "
data {
int<lower=0> N;
vector[N] x;
int<lower=0,upper=1> y[N];
int<lower=0> M;
vector[M] new_x;
}
parameters {
real alpha;
real beta;
}
model {
// more efficient and arithmetically stable
y ~ bernoulli_logit(alpha + beta * x);
}
generated quantities {
vector[M] y_epred;
vector[M] mu = alpha + beta * new_x;
for(i in 1:M) {
y_epred[i] = inv_logit(mu[i]);
}
}
"
newdata <- data.frame(
body_mass_g = seq(min(gentoo$body_mass_g), max(gentoo$body_mass_g), length.out = 100)
)
stan_data <- list(
N = nrow(gentoo),
y = gentoo$gender,
x = gentoo$body_mass_g,
M = nrow(newdata),
new_x = newdata$body_mass_g
)
m <- stan(model_code = stan_program, data = stan_data)
fit <- m %>%
tidybayes::gather_draws(y_epred[i]) %>%
ggdist::mean_qi(.value)## Warning: `gather_()` was deprecated in tidyr 1.2.0.
## Please use `gather()` instead.
## This warning is displayed once every 8 hours.
## Call `lifecycle::last_lifecycle_warnings()` to see where this warning was generated.
fit## # A tibble: 100 × 8
## i .variable .value .lower .upper .width .point .interval
## <int> <chr> <dbl> <dbl> <dbl> <dbl> <chr> <chr>
## 1 1 y_epred 0.0000218 0.0000000114 0.000157 0.95 mean qi
## 2 2 y_epred 0.0000267 0.0000000169 0.000188 0.95 mean qi
## 3 3 y_epred 0.0000327 0.0000000252 0.000227 0.95 mean qi
## 4 4 y_epred 0.0000401 0.0000000369 0.000273 0.95 mean qi
## 5 5 y_epred 0.0000493 0.0000000546 0.000332 0.95 mean qi
## 6 6 y_epred 0.0000605 0.0000000814 0.000404 0.95 mean qi
## 7 7 y_epred 0.0000744 0.000000121 0.000492 0.95 mean qi
## 8 8 y_epred 0.0000916 0.000000178 0.000599 0.95 mean qi
## 9 9 y_epred 0.000113 0.000000264 0.000725 0.95 mean qi
## 10 10 y_epred 0.000139 0.000000389 0.000884 0.95 mean qi
## # … with 90 more rows两个图画在一起
fit %>%
bind_cols(newdata) %>%
ggplot(aes(x = body_mass_g)) +
geom_dots(
data = gentoo,
aes(y = gender, side = ifelse(sex == "male", "bottom", "top")),
scale = 0.4
) +
geom_lineribbon(
aes(y = .value, ymin = .lower, ymax = .upper),
alpha = 1/4,
fill = "#08306b"
) +
labs(
title = "logit dotplot: stat_dots() with stat_lineribbon()",
subtitle = 'aes(side = ifelse(sex == "male", "bottom", "top"))',
x = "Body mass (g) of Gentoo penguins",
y = "Pr(sex = male)"
)
67.2 篮球案例
我们模拟100个选手每人投篮20次，假定命中概率是身高的线性函数，案例来源chap15.3 of [Regression and Other Stories] (page270).
n <- 100
data <-
tibble(size = 20,
height = rnorm(n, mean = 72, sd = 3)) %>%
mutate(y = rbinom(n, size = size, p = 0.4 + 0.1 * (height - 72) / 3))
head(data)## # A tibble: 6 × 3
## size height y
## <dbl> <dbl> <int>
## 1 20 76.5 4
## 2 20 70.8 8
## 3 20 70.7 9
## 4 20 74.1 9
## 5 20 68.3 4
## 6 20 74.5 1267.2.1 常规做法
##
## Call: glm(formula = cbind(y, 20 - y) ~ height, family = binomial(link = "logit"),
## data = data)
##
## Coefficients:
## (Intercept) height
## -10.2702 0.1371
##
## Degrees of Freedom: 99 Total (i.e. Null); 98 Residual
## Null Deviance: 188.5
## Residual Deviance: 106.9 AIC: 443.267.2.2 stan 代码
\[ \begin{align*} y_i & = \text{Binomial}(n_i, p_i) \\ p_i & =\text{logit}^{-1}(X_i \beta) \end{align*} \]
stan_program <- "
data {
int<lower=0> N;
int<lower=0> K;
matrix[N, K] X;
int<lower=0> y[N];
int trials[N];
}
parameters {
vector[K] beta;
}
model {
for(i in 1:N) {
target += binomial_logit_lpmf(y[i] | trials[i], X[i] * beta);
}
}
"
stan_data <- data %>%
tidybayes::compose_data(
N = n,
K = 2,
y = y,
trials = size,
X = model.matrix(~ 1 + height)
)
fit <- stan(model_code = stan_program, data = stan_data)
fit## Inference for Stan model: anon_model.
## 4 chains, each with iter=2000; warmup=1000; thin=1;
## post-warmup draws per chain=1000, total post-warmup draws=4000.
##
## mean se_mean sd 2.5% 25% 50% 75% 97.5% n_eff Rhat
## beta[1] -10.24 0.05 1.18 -12.46 -11.02 -10.27 -9.49 -7.90 502 1
## beta[2] 0.14 0.00 0.02 0.10 0.13 0.14 0.15 0.17 502 1
## lp__ -220.62 0.04 1.03 -223.33 -221.06 -220.30 -219.88 -219.61 537 1
##
## Samples were drawn using NUTS(diag_e) at Sat Jul 16 11:17:06 2022.
## For each parameter, n_eff is a crude measure of effective sample size,
## and Rhat is the potential scale reduction factor on split chains (at
## convergence, Rhat=1).## Warning in rm(gentoo, fit, m, stan_data, stan_program, data, fit_glm, fit, :
## object 'fit' not found