Bayesian logistic regression with mixed (two level) data

Question

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.

Answer 1

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
fit = glm(fit.formula, data=schools.test, family=binomial(link='logit'))

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

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