Bayesian logistic regression with mixed (two level) data

I have two datasets, the first on schools, and the second lists students in each school who have failed in standardized test (emphasis intentional). The first dataset has the structure below (I'm sorry I could not figure out how to insert a proper data snippet):

School_ID   Total_White   Total_Black   Total_Asian   School_Revenue


The second looks like this:

Student_ID   School_ID   Race


I am trying to estimate the probability of failure, given student race and school revenue. If I run a multinomial discrete choice model on the second dataset, I shall clearly be estimating P(Race | Fail=1). I obviously have to estimate the inverse, P(Fail=1 | Race), by merging the two datasets by School_ID. Since all the pieces of information are available in the two datasets (P(Fail), P(Race), Revenue), I see no reason why this can't be done. But I am stumped as to actually how to implement in R. Any pointer would be much appreciated. Thanks.

-
What exactly is stumping you, merging the datasets or how to set up the analysis after the merge, or...? –  jbowman Jan 16 '12 at 21:56
why multinomial? the dependent variable seems to be binary. –  Manoel Galdino Jan 16 '12 at 22:30

There are many good tutorials on logistic regression in R on the web. A couple tutorials that helped me are here and here. You want to merge the two data sets so that you end up with the number of passes and the number of failures for each school by race. The glm function will perform the logistic regression. The race variable will be interpreted as factors. A summary report will give you the coefficients and the significance of each variable. Halfway through the second tutorial there is a good demonstration on generating probability predictions.

Here is a stab at it with some random data. Additional school parameters can be included by updating the model formula.

# -- Generate some fake data

#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 failed.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) ) # -- end fake data #roll-up the number of failed students by school #count the number of students using the length function failed = aggregate(list(failed=students$student_id),by=list(school_id=students$school_id,race=students$race),FUN=length)

#merge the failure data with the school data
schools.test = merge(schools.data,failed,by.x="school_id",by.y="school_id")

#compute the proportion of students by race
schools.test$pct_white = schools.test$tot_white/schools.test$tot_students schools.test$pct_black = schools.test$tot_black/schools.test$tot_students
schools.test$pct_asian = schools.test$tot_asian/schools.test$tot_students #Get the number of passed students schools.test$passed=0
race.totals = c('tot_white','tot_black','tot_asian')
race.factors = c('white','black','asian')
for (i in 1:3){
mask = schools.test$race==race.factors[i] schools.test$passed[mask] = schools.test[,race.totals[i]][mask] - schools.test$failed[mask] } #modify the default base race factor schools.test$race = relevel(schools.test\$race,ref='white')

#compute the logistic regression
fit.formula = cbind(failed,passed)~race+school_rev+pct_white+pct_black+pct_asian

print(summary(fit))
print(anova(fit))

par(mfrow=c(2,2))
plot(fit)

-
Tharen... Many thanks for this simple and elegant solution. This approach is the exact opposite of my line of thought, which was to do a multinomial discrete choice at the individual level, weighted by the distribution of race at the school level (essentially a Bayesian inverse probability). But one problem I see with your solution is that I can't add other student-characteristics (say, poverty status). Or am I missing something? –  user3671 Jan 17 '12 at 6:01
@user3671 I modified the answer to allow including additional school variables. I don't see how you could add additional student variables unless you include students who passed as well. If you had that data as well, the school data could be bound to the individual students and their respective pass/fail would be the dependent variable. Also, I agree with the comment above that this looks like a binary problem as a student can only pass or fail. Finally, in hind-sight this is not a Bayesian solution as you had specified. –  tharen Jan 17 '12 at 8:06
I edited the answer to use percent of students in each race as the terms of the model. This seems more appropriate. –  tharen Jan 17 '12 at 8:44
Thanks, Tharen. Somehow I was thinking that there could be a way to find P(Fail|Race) as a product of P(Race|Fail), estimated at the individual level and the likelihood of race estimated at the school level. But I must be mistaken. Thanks a bunch. –  user3671 Jan 17 '12 at 8:57