在glm()中的逻辑回归中预测(0,1)。

Question

我正在尝试使用二元logit模型对“假设”情况进行建模。我正在估计在性别控制下，通过测试的概率，给定测试难度的水平（1=最容易，5=最困难）。 （数据在这里）。学生们参加的测试通常很难（数据中为“HIGH”）。从这个模型中，我们可以估计测试难度对通过考试的可能性的影响：

model = glm(PASS ~ as.factor(SEX) + as.factor(HIGH), family=binomial(link="logit"), data=df)
summary(model)

我们还可以使用以下代码获取通过考试的预测概率：

predict.high = predict(model, type="response")

问题是如果进行“LOW”测试怎么办？为了得到新的概率，我们可以执行以下操作：

newdata = rename.vars(subset(df, select=c(-HIGH)), 'LOW','HIGH')
predict.low = predict(model, newdata=newdata, type="response")

但是我怎么知道在这种情况下有多少额外的学生会通过考试呢？glm() 中是否有一个明显的开关我没有注意到？

- user702432

---

马克，我基本上正在进行模拟，即用“LOW”替换“HIGH”的实际值。 - user702432

好的，看起来HIGH代表5个测试。在估计通过这5个测试的概率后，您正在使用相同的模型来预测通过LOW代表的5个新测试的概率。我在安装了“gregmisc”包后使用了您的数据运行了您的代码。然而，我没有我的Gelman和Hill（2006）副本。我有他们的代码注释版本，如果您的问题没有被其他人首先回答，我可以在周末挖出来并尝试找到答案。 - Mark Miller

提前感谢，马克。既然您提到了G＆H，我还想向您指出“arm”包中的bayesglm（）函数。它非常灵活。 - user702432

Ben...你可以将“HIGH”和“LOW”视为同一“难度”变量的两个备选数字抽取。我使用“HIGH”值来估计系数值。如果我将“HIGH”值插入到估计方程中，我将得到一组预测概率。如果我在相同的估计方程中插入“LOW”值，我将得到另一组预测值。我可以通过叠加两个预测密度图来了解差异。但我的问题是是否可能以数字术语获得差异。我想说使用“HIGH”值 - user702432

“which went into estimating the coefficients”，我预测有800名学生通过。现在，如果我插入“LOW”值，那么就会有1000名学生通过。这在R中是否可行？ - user702432

Answer 1

我还没有尝试挖掘我基于Gelman和Hill（2006）编写的预测代码，他们似乎使用了模拟。我仍然打算这样做。你问题中似乎独特的一个方面是，我习惯于为单个观察值（在这种情况下是单个学生参加单个考试）进行预测。然而，你似乎想要预测两组预测之间的差异。换句话说，你想预测如果给定5个简单考试的一组而不是5个困难考试的一组，将有多少个学生通过。

我不确定Gelman和Hill（2006）是否涵盖了这一点。你似乎也想以频率学派的方法来做这件事。

我认为，如果你可以为单个观察值进行预测，使得每个观察值都有置信区间，那么也许你可以估计每个组内通过的加权平均概率，并减去两个加权平均值。Delta方法可用于估计加权平均值及其差异的置信区间。

为了实施这种方法，可能必须假定预测观察之间的协方差为0。

如果假设协方差为0不令人满意，那么也许贝叶斯方法会更好。同样，我只熟悉为单个观察值进行预测。在贝叶斯方法中，我通过包括独立变量而不是因变量来预测单个观察值。我想你可以在同一贝叶斯运行中为每个观察值进行预测（在HIGH和LOW中预测每个学生）。每组测试通过的加权平均值和加权平均值之间的差异是派生参数，我怀疑可以直接将它们包含在贝叶斯逻辑回归的代码中。然后，你将得到每组测试通过概率及其差异的点估计和方差估计。如果你想要每组测试通过的学生人数的差异，也许这也可以作为一个派生参数包含在贝叶斯代码中。

我意识到，到目前为止，这个答案比较口语化，可能不如期望的那么正式。我只是在没有时间尝试实施这些策略的情况下制定了这些策略。提供实现两种策略的所有R和WinBUGS代码可能需要我几天时间。（WinBUGS或OpenBUGS可以从R中调用。）我会在我继续编写代码的同时附加代码到这个答案中。如果有人认为我的建议策略和/或即将发布的代码不正确，我希望他们可以指出我的错误并提供更正。

编辑

下面是生成虚假数据并使用频率学派和贝叶斯方法分析该数据的代码。我还没有添加实现上述预测想法的代码。我将尝试在接下来的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日

我没有忘记这个。相反，我遇到了几个问题：

1. 用逻辑回归可以预测：a）给定组中学生通过考试的概率p和b）给定学生参加考试的结果（0或1）。然后对所有0和1求平均值。我不确定要使用哪个。预测的p的点估计和SD与已知测试结果的估计p相同。预测的0和1的平均值的点估计略有不同，而平均0和1的SD要大得多。我相信我想要b，预测的0和1的平均值。但是，我正在尝试检查各种网站和书籍以确保。Collett（1991）有一个未使用计算机代码的示例，但该示例包括半打变量，包括2个交互作用，并且我在让我的贝叶斯估计匹配她的频率估计时遇到了一些麻烦。

2. 由于有很多派生参数，程序运行时间很长。

3. 显然，OpenBUGS经常崩溃，我认为即使没有预测代码也是如此。我不确定这是因为我做错了什么，还是因为最近版本的R或R软件包中的更改，或者因为我正在尝试使用64位R运行代码，还是其他原因。

我将尽快发布预测代码，但以上所有问题都使我放缓了速度。

- Mark Miller

---

非常感谢你的时间和努力，马克。我觉得现在我的头脑开始清晰了。你对问题的解释非常准确。我的问题是，如果学生们随机获得一组平均难度比他们实际获得的测试要容易的测试，那该怎么办？步骤1：通过回归难度水平和实际通过率的变化来估计通过概率。从中，我们可以得到一个预测密度曲线。步骤2：插入反事实测试集的值，并获得另一个预测概率密度曲线。差异 - user702432

在这两条曲线中，给出了预测的概率差异。顺便说一下...我实际上是在bayesglm(arm)中使用贝叶斯结构和弱先验。我没有在问题中提到的原因是glm()和bayesglm()非常相似，我不想让人们感到困惑。 - user702432

Answer 2

您可以轻松使用此方法来找到截止点：

cutoff <- runif(length(predicted_probabilities)) 

这是基于Metropolis-Hastings的确定性决策。

- user3430235