我还没有尝试挖掘出我根据 Gelman 和 Hill (2006) 编写的预测代码，我似乎记得使用过模拟。我仍然打算这样做。在我有限的经验中，您的问题的一个方面似乎是独一无二的，那就是我习惯于预测一个单一的观察结果（在这种情况下，一个学生参加一个测试）。但是，您似乎想要预测两组预测之间的差异。换句话说，如果给定一组 5 门简单的考试而不是一组 5 门难的考试，您想预测有多少学生将通过。
我不确定 Gelman 和 Hill (2006) 是否涵盖了这一点。你似乎也想用一种常客的方法来做到这一点。
我在想，如果您可以预测单个观察结果，以便您对每个观察结果都有一个置信区间，那么也许您可以估计每个组内通过的加权平均概率并减去两个加权平均值。delta 方法可用于估计加权平均值及其差异的置信区间。
预测观测值之间的协方差可能必须假定为 0 才能实施该方法。
如果假设协方差 0 不令人满意，那么贝叶斯方法可能会更好。同样，我只熟悉预测单个观察值。使用贝叶斯方法，我通过包含自变量而不是因变量来预测单个观察结果，以便预测观察结果。我想您可以预测同一贝叶斯运行中的每个观察结果（预测每个学生的高和低）。每组通过测试的加权平均值和加权平均值的差异是派生参数，我怀疑可以直接包含在贝叶斯逻辑回归的代码中。然后，您将对通过每组测试的概率以及通过每组测试的概率差异进行点估计和方差估计。
我意识到，到目前为止，这个答案比预期的更具对话性。我只是在制定策略来尝试，而没有时间尝试实施这些策略。提供所有 R 和 WinBUGS 代码来实现这两种建议的策略可能需要我几天时间。（可以从 R 中调用 WinBUGS 或 OpenBUGS。）我将在此过程中将代码附加到此答案中。如果有人认为我提出的策略和/或即将发布的代码不正确，我希望他们随时指出我的错误并提供更正。
编辑
下面是生成假数据并使用常客和贝叶斯方法分析该数据的代码。我还没有添加代码来实现上述预测的想法。我会在接下来的 1-2 天内尝试添加贝叶斯预测代码。我只使用了三个测试而不是五个。通过下面的代码编写方式，您可以将学生人数 n 更改为任何可分为 6 个相等整数的非零数。
# Bayesian_logistic_regression_June2012.r
# June 24, 2012
library(R2WinBUGS)
library(arm)
library(BRugs)
set.seed(3234)
# create fake data for n students and three tests
n <- 1200
# create factors for n/6 students in each of 6 categories
gender <- c(rep(0, (n/2)), rep(1, (n/2)))
test2 <- c(rep(0, (n/6)), rep(1, (n/6)), rep(0, (n/6)),
rep(0, (n/6)), rep(1, (n/6)), rep(0, (n/6)))
test3 <- c(rep(0, (n/6)), rep(0, (n/6)), rep(1, (n/6)),
rep(0, (n/6)), rep(0, (n/6)), rep(1, (n/6)))
# assign slopes to factors
B0 <- 0.4
Bgender <- -0.2
Btest2 <- 0.6
Btest3 <- 1.2
# estimate probability of passing test
p.pass <- ( exp(B0 + Bgender * gender +
Btest2 * test2 +
Btest3 * test3) /
(1 + exp(B0 + Bgender * gender +
Btest2 * test2 +
Btest3 * test3)))
# identify which students passed their test, 0 = fail, 1 = pass
passed <- rep(0, n)
r.passed <- runif(n,0,1)
passed[r.passed <= p.pass] = 1
# use frequentist approach in R to estimate probability
# of passing test
m.freq <- glm(passed ~ as.factor(gender) +
as.factor(test2) +
as.factor(test3) ,
family = binomial)
summary(m.freq)
# predict(m.freq, type = "response")
# use OpenBUGS to analyze same data set
# Define model
sink("Bayesian.logistic.regression.txt")
cat("
model {
# Priors
alpha ~ dnorm(0,0.01)
bgender ~ dnorm(0,0.01)
btest2 ~ dnorm(0,0.01)
btest3 ~ dnorm(0,0.01)
# Likelihood
for (i in 1:n) {
passed[i] ~ dbin(p[i], 1)
logit(p[i]) <- (alpha + bgender * gender[i] +
btest2 * test2[i] +
btest3 * test3[i])
}
# Derived parameters
p.g.t1 <- exp(alpha) / (1 + exp(alpha))
p.b.t1 <- exp(alpha + bgender) / (1 + exp(alpha + bgender))
p.g.t2 <- ( exp(alpha + btest2) /
(1 + exp(alpha + btest2)))
p.b.t2 <- ( exp(alpha + bgender + btest2) /
(1 + exp(alpha + bgender + btest2)))
p.g.t3 <- ( exp(alpha + btest3) /
(1 + exp(alpha + btest3)))
p.b.t3 <- ( exp(alpha + bgender + btest3) /
(1 + exp(alpha + bgender + btest3)))
}
", fill = TRUE)
sink()
my.data <- list(passed = passed,
gender = gender,
test2 = test2,
test3 = test3,
n = length(passed))
# Inits function
inits <- function(){ list(alpha = rlnorm(1),
bgender = rlnorm(1),
btest2 = rlnorm(1),
btest3 = rlnorm(1)) }
# Parameters to estimate
params <- c("alpha", "bgender", "btest2", "btest3",
"p.g.t1", "p.b.t1", "p.g.t2", "p.b.t2",
"p.g.t3", "p.b.t3")
# MCMC settings
nc <- 3
ni <- 2000
nb <- 500
nt <- 2
# Start Gibbs sampling
out <- bugs(data = my.data, inits = inits,
parameters.to.save = params,
"c:/users/Mark W Miller/documents/Bayesian.logistic.regression.txt",
program = 'OpenBUGS',
n.thin = nt, n.chains = nc,
n.burnin = nb, n.iter = ni, debug = TRUE)
print(out, dig = 5)
在我尝试实施加权平均预测方法之前，我想说服自己它可能有效。因此，我编写了以下代码，这似乎表明它可能：
# specify number of girls taking each test and
# number of boys taking each test
g.t1 <- rep(0,400)
b.t1 <- rep(0,120)
g.t2 <- rep(0,1200)
b.t2 <- rep(0,50)
g.t3 <- rep(0,1000)
b.t3 <- rep(0,2000)
# specify probability of individuals in each of the
# 6 groups passing their test
p.g1.t1 <- 0.40
p.b1.t1 <- 0.30
p.g1.t2 <- 0.60
p.b1.t2 <- 0.50
p.g1.t3 <- 0.80
p.b1.t3 <- 0.70
# identify which individuals in each group passed their test
g.t1[1:(p.g1.t1 * length(g.t1))] = 1
sum(g.t1)
b.t1[1:(p.b1.t1 * length(b.t1))] = 1
sum(b.t1)
g.t2[1:(p.g1.t2 * length(g.t2))] = 1
sum(g.t2)
b.t2[1:(p.b1.t2 * length(b.t2))] = 1
sum(b.t2)
g.t3[1:(p.g1.t3 * length(g.t3))] = 1
sum(g.t3)
b.t3[1:(p.b1.t3 * length(b.t3))] = 1
sum(b.t3)
# determine the weighted average probability of passing
# on test day for all individuals as a class
wt.ave.p <- ((p.g1.t1 * length(g.t1) + p.b1.t1 * length(b.t1) +
p.g1.t2 * length(g.t2) + p.b1.t2 * length(b.t2) +
p.g1.t3 * length(g.t3) + p.b1.t3 * length(b.t3) ) /
(length(g.t1) + length(b.t1) + length(g.t2) +
length(b.t2) + length(g.t3) + length(b.t3)))
wt.ave.p
# determine the expected number of individuals passing
# their test in the class as a whole
exp.num.pass <- wt.ave.p * (length(g.t1) + length(b.t1) +
length(g.t2) + length(b.t2) +
length(g.t3) + length(b.t3))
exp.num.pass
# determine the number of individuals passing
num.passing <- (sum(g.t1) + sum(b.t1) +
sum(g.t2) + sum(b.t2) +
sum(g.t3) + sum(b.t3) )
num.passing
# the expected number of students passing, exp.num.pass,
# should equal the observed number of students passing,
# num.passing regardless of the number of students in each
# group and regardless of the probability of passing a
# given test, within rounding error
identical(round(exp.num.pass), round(num.passing))
希望在接下来的几天里，我可以尝试将预测代码添加到上述贝叶斯代码中。
编辑 - 2012 年 6 月 27 日
我没有忘记这件事。相反，我遇到了几个问题：
使用逻辑回归可以预测：a）给定组中的学生通过考试的概率 p 和 b）给定学生参加考试的结果（0 或 1）。然后对所有的 0 和 1 进行平均。我不确定要使用其中的哪一个。预测 p 的点估计和 SD 与已知测试结果的估计 p 相同。预测的 0 和 1 的平均值的点估计略有不同，平均 0 和 1 的 SD 要大得多。我相信我想要 b，预测的 0 和 1 的平均值。但是，我正在尝试检查各种网站和书籍以确定。Collett (1991) 有一个不使用计算机代码的工作示例，
由于有很多派生参数，该程序需要很长时间才能运行。
显然 OpenBUGS 经常崩溃，我相信，即使没有预测代码。我不确定这是因为我做错了什么，还是因为最近版本的 R 的变化或最近版本的 R 包的变化，或者可能是因为我试图用 64 位 R 或其他东西运行代码别的。
我会尽快发布预测代码，但上述所有问题都让我放慢了速度。