Р: Байесовская логистическая регрессия для иерархических данных

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 (Ошибка = 1 | Раса учеников, Доход школы).Если я использую полиномиальную модель дискретного выбора в наборе студенческих данных, я точно буду оценивать P (Race | Fail = 1).Я, очевидно, должен оценить обратное.Поскольку все фрагменты информации доступны в двух наборах данных (P (Fail), P (Race), Доход), я не вижу причин, почему это невозможно сделать.Но я нахожусь в тупик, как на самом деле, как реализовать в 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, если вам нужны байесовские оценки.