R：分层数据的贝叶斯逻辑回归

Question

这是来自stats.stackexchange 的转贴，我没有得到满意的回复。我有两个数据集，第一个关于学校，第二个列出每所学校在标准化考试中未通过的学生（强调是故意的）。可以通过以下方式生成假数据集（感谢Tharen）：

#random school data for 30 schools
schools.num = 30
schools.data = data.frame(school_id=seq(1,schools.num)
                         ,tot_white=sample(100:300,schools.num,TRUE)
                         ,tot_black=sample(100:300,schools.num,TRUE)
                         ,tot_asian=sample(100:300,schools.num,TRUE)
                         ,school_rev=sample(4e6:6e6,schools.num,TRUE)
                         )

#total students in each school
schools.data$tot_students = schools.data$tot_white + schools.data$tot_black + schools.data$tot_asian
#sum of all students all schools
tot_students = sum(schools.data$tot_white, schools.data$tot_black, schools.data$tot_asian)
#generate some random failing students
fail.num = as.integer(tot_students * 0.05)

students = data.frame(student_id=sample(seq(1:tot_students), fail.num, FALSE)
                      ,school_id=sample(1:schools.num, fail.num, TRUE)
                      ,race=sample(c('white', 'black', 'asian'), fail.num, TRUE)
                      )

我正在尝试估算 P(Fail=1 | Student Race, School Revenue)。如果我在学生数据集上运行多项式离散选择模型，我显然会估计 P(Race | Fail=1)。我显然必须估计这个的倒数。由于两个数据集（P（失败）、P（种族）、收入）中的所有信息都可用，我认为没有理由不能这样做。但我对如何在 R 中实现感到困惑。任何指针都将不胜感激。谢谢。

Answer 1

如果你有一个单独的 data.frame 会更容易。

library(reshape2)
library(plyr)
d1 <- ddply(
  students, 
  c("school_id", "race"), 
  summarize,
  fail=length(student_id)
) 
d2 <- with( schools.data, data.frame( 
  school_id = school_id, 
  white = tot_white, 
  black = tot_black, 
  asian = tot_asian, 
  school_rev = school_rev 
) )
d2 <- melt(d2, 
  id.vars=c("school_id", "school_rev"), 
  variable.name="race", 
  value.name="total"
)
d <- merge( d1, d2, by=c("school_id", "race") )
d$pass <- d$total - d$fail

然后你可以查看数据

library(lattice)
xyplot( d$fail / d$total ~ school_rev | race, data=d )

或者计算任何你想要的东西。

r <- glm(
  cbind(fail,pass) ~ race + school_rev, 
  data=d, 
  family=binomial() # Logistic regression (not bayesian)
)
summary(r)

（编辑）如果您有更多关于失败学生的信息， 但只有通过的汇总数据， 您可以按如下方式重新创建完整的数据集。

# Unique student_id for the passed students
d3 <- ddply( d, 
  c("school_id", "race"), 
  summarize, student_id=1:pass 
)
d3$student_id <- - seq_len(nrow(d3))
# All students
d3$result <- "pass"
students$result <- "fail"
d3 <- merge( # rather than rbind, in case there are more columns
  d3, students, 
  by=c("student_id", "school_id", "race", "result"), 
  all=TRUE 
)
# Students and schools in a single data.frame
d3 <- merge( d3, schools.data, by="school_id", all=TRUE )
# Check that the results did not change
r <- glm(
  (result=="fail") ~ race + school_rev, 
  data=d3, 
  family=binomial()
)
summary(r)

Answer 2

您需要一个包含所有学生信息的数据集。都失败了，都通过了。

schools.num = 30
schools.data = data.frame(school_id=seq(1,schools.num)
                          ,tot_white=sample(100:300,schools.num,TRUE)
                          ,tot_black=sample(100:300,schools.num,TRUE)
                          ,tot_asian=sample(100:300,schools.num,TRUE)
                          ,school_rev=sample(4e6:6e6,schools.num,TRUE)
                          )

library(plyr)
fail_ratio <- 0.05
dataset <- ddply(schools.data, .(school_id, school_rev), function(x){
  data.frame(Fail = rbinom(sum(x$tot_white, x$tot_asian, x$tot_black), size = 1, prob = fail_ratio), Race = c(rep("white", x$tot_white), rep("asian", x$tot_asian), rep("black", x$tot_black)))
})
dataset$Race <- factor(dataset$Race)

然后，您可以将 glmer() 用于 lme4 包，以获得频繁的方法。

library(lme4)
glmer(Fail ~ school_rev + Race + (1|school_id), data = dataset, family = binomial)

如果您需要贝叶斯估计，请查看 MCMCglmm 包。