我还没有尝试挖掘我根据 Gelman 和 Hill (2006) 编写的预测代码，我似乎记得他们使用了模拟。我仍然打算这样做。在我有限的经验中，你的问题的一个方面似乎很独特，那就是我习惯于预测单个观察结果（在本例中是单个学生参加单个测试）。然而，您似乎想要预测两组预测之间的差异。换句话说，您想要预测如果进行一组 5 项简单考试而不是一组 5 项困难考试，将会有多少学生通过。
我不确定 Gelman 和 Hill (2006) 是否涵盖了这一点。您似乎也想用频率主义方法来做到这一点。
我认为，如果您可以预测单个观察结果，以便每个观察结果都有一个置信区间，那么也许您可以估计每个组内通过的加权平均概率并减去两个加权平均值。 Delta 方法可用于估计加权平均值及其差异的置信区间。
为了实施该方法，预测观测值之间的协方差可能必须假设为 0。
如果假设协方差为 0 不能令人满意，那么贝叶斯方法也许会更好。同样，我只熟悉预测单个观察结果。使用贝叶斯方法，我通过包含自变量（但不包含因变量）来预测单个观测值，以便预测观测值。我想你可以预测同一贝叶斯运行中的每个观察结果（预测每个学生的高和低）。每组通过测试的加权平均值和加权平均值的差异是派生参数，我怀疑可以直接包含在贝叶斯逻辑回归的代码中。然后，您将获得通过每组测试的概率以及通过每组测试的概率差异的点估计和方差估计。如果您想要通过每组测试的学生人数差异，也许也可以将其作为派生参数包含在贝叶斯代码中。
我意识到到目前为止，这个答案比预期的更具对话性。我只是制定了要尝试的策略，而还没有时间尝试实施这些策略。提供所有 R 和 WinBUGS 代码来实现这两个提议的策略可能需要我几天的时间。 （可以从 R 内部调用 WinBUGS 或 OpenBUGS。）我将在继续过程中将代码附加到这个答案中。如果有人认为我提出的策略和/或即将发布的代码不正确，我希望他们能够随意指出我的错误并提供更正。
EDIT
下面的代码生成虚假数据并使用频率论和贝叶斯方法分析该数据。我还没有添加代码来实现上述预测想法。我将尝试在接下来的 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) 有一个不使用计算机代码的有效示例，但该有效示例包含六个变量，其中包括 2 个交互作用，并且我在让我的贝叶斯估计与她的频率估计相匹配时遇到了一些麻烦。
由于有大量派生参数，程序需要很长时间才能运行。
我相信，即使没有预测代码，OpenBUGS 显然也经常崩溃。我不确定这是因为我做错了什么，还是因为最新版本的 R 的更改或最新版本的 R 软件包的更改，或者可能是因为我尝试使用 64 位 R 运行代码或其他原因别的。
我会尝试尽快发布预测代码，但是上述所有问题都减慢了我的速度。
