Interpolating and aggregating Median Income from binned Census data

Question

I'm working with the most recent ACS Census data for New York City. I'm trying to calculate median income for Neighborhood Tabulation Areas (NTAs), which consist of several Census tracts.

The Census provides Median Family Income, Mean Family Income and binned income data, at the tract level - I'm trying to find a way to aggregate this to the NTA level.

I know which Tracts are in each NTA, so my approach was to simply combine the binned data, as such:

         NTA                     Income Bin           n    min    max  cumsum
1  Neighborhood 1     Less.than..10.000_Households   540      0   9999    540
2  Neighborhood 1   X.10.000.to..14.999_Households   431  10000  14999    971
3  Neighborhood 1   X.15.000.to..24.999_Households   617  15000  24999   1588
4  Neighborhood 1   X.25.000.to..34.999_Households   755  25000  34999   2343
5  Neighborhood 1   X.35.000.to..49.999_Households   961  35000  49999   3304
6  Neighborhood 1   X.50.000.to..74.999_Households  1439  50000  74999   4743
7  Neighborhood 1   X.75.000.to..99.999_Households  1652  75000  99999   6395
8  Neighborhood 1 X.100.000.to..149.999_Households  2517 100000 149999   8912
9  Neighborhood 1 X.150.000.to..199.999_Households  1534 150000 199999  10446
10 Neighborhood 1     X.200.000.or.more_Households  3816 200000         14262

I then used the approach mentioned in this post:

Some books give formulas by which to calculate an estimate of the median which amount to doing pretty much what I just did, such as a formula like this: median = L+w(n2−c)fL+w(n2−c)f (where LL is the lower limit of the bin containing the median, ww is the width of that bin, nn is the total population, cc is the cumulative count (cumulative frequency) up to LL (the end of the previous bin), and ff is the count (frequency) in the median bin does pretty much the same thing (aside from the (n+1)/2 vs n/2, it is the same).

However, when I try it on 2010 data, which has already been interpolated to the NTA level by the city, my results do not match at all.

This is the (bad) code in R I used to produce my results:

last_real_cum <- nrow(branch_1)
branch_cum <- branch_1$cumsum[last_real_cum]
med_value <- (branch_cum + 1)/2
for (i in 1:nrow(branch_1)) {
if (med_value>branch_1$cumsum [i] & med_value<branch_1$cumsum[i+1])
{
  branch_1$group[i+1] <- "Y"
  var_L <- branch_1$min[i+1]
  var_F <- branch_1$value[i+1]
  var_W <- branch_1$max[i+1] - branch_1$min[i+1]
  var_C <- branch_1$cumsum[i]
}
else {
  branch_1$group[i] <- "N"
}

    var_N <- branch_cum
    est_med <- var_L + (var_W * ((var_N/2)-var_C)) / var_F
    return(est_med)

Is there a better approach to doing this which would produce more accurate results?

Answer 1

A bit late but try this:

#The quantile.bins function takes 3 inputs:
# binedges is a numerical vector representing the TOP value of each bin. So for incomes 1000-14999, binedge should be 14999. 
#         binedges must be sorted in an ascending fashion. The top value will be ingored by the function,
#         so you can make the top value one dollar more than the top bin (e.g. 150,000) or just leave it NA
# bincounts is a numerical vector of people in that income group
# quantile is the value you're looking for - usually .5 for medians. 

library(binsmooth)


##see bottom for a little tutorial
median.approxfun<-function(x,probability){
  jane <- 0
  tim <- 0
  while(jane < probability){
    jane <- x(tim)
    tim <- tim + 1000
  }
  while(jane > probability){
    jane <- x(tim)
    tim <- tim - 100
  }
  while(jane < probability){
    jane <- x(tim)
    tim <- tim + 10
  }
  while(jane > probability){
    jane <- x(tim)
    tim <- tim - 1
  }
  while(jane < probability){
    jane <- x(tim)
    tim <- tim + .1
  }
  return(tim)}

quantile.bins<-function(binedges,bincounts,quantile=.5){
  splb <- splinebins(binedges, bincounts)
  output<-median.approxfun(splb$splineCDF,quantile)
  output}

Just turn the cum/sum into raw count data, then you can use the tidyverse code below:

output<- data %>% group_by(NTA) %>% summarise(median=quantile.bins(max,count))

Answer 2

For any income $x$, you can estimate the population in a tract with income $\le x$ by linearly interpolating the given data, eg:

i = sum(df$min <= x)
frac = (x - df$min[i]) / (df$max[i] - df$min[i])
pop = (1 - frac)*df$n[i] + frac*df$n[i+1]

Then find an $x$ so that the total population with income $\le x$ is half the total population over all the tracts. For instance, you can start with $x$ equal to the mean income, and then examine lower incomes as needed.