R：階層的データのベイジアンロジスティック回帰

Question

これはstats.stackexchangeからの再投稿です。ここで私は満足のいく回答を得られませんでした。私は2つのデータセットを持っています。最初は学校であり、2番目は各学校の学生です。標準化テストで失敗した人（重点は意図的です）。偽のデータセットは、（Tharenのおかげで）によって生成することができます：私は（失敗= 1 |学生の人種、学校収益）Pを推定しようとしています

#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（Race | Fail = 1）を見積もっているはずです。私は明らかにこれの逆数を推定する必要があります。すべての情報が2つのデータセット（P（Fail）、P（Race）、Revenue）で使用可能であるため、これを行うことができない理由はありません。しかし、私は実際にどのようにRで実装するのか困っています。どんなポインタも大いに感謝します。ありがとう。

Answer 1

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) 

（EDIT）あなたが渡されたもののために失敗した学生、 だけで集約されたデータに関する詳細な情報を持っている場合は、次のように あなたは、完全なデータセットを再作成することができます。

# 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) 

次に、頻繁にアプローチするためにlme4パッケージにglmer（）を使用できます。

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

ベイジアン見積もりが必要な場合は、MCMCglmmパッケージをご覧ください。