Use of Kernel Densities to Evaluate Food Supply and Deserts in Dallas

Michael Tiefelsdorf & Laila Al Aswad

2020-02-23

Introduction

The key proposition of this analysis is that the spatial variability of food supply in an urban environment can be approximated by kernel densities around individual stores selling food. It is assumed that grocery stores offer nutritious food whereas convenience stores (including food aisles of gas stations etc.) predominately sell processed food comprising mainly of ‘empty’ calories.

The potential spatial extend of the store’s market areas is proportional to the bandwidth of the kernel function. The estimated sales volume of food products is used as weight of the kernel density estimator. Usually grocery stores have a larger market area and a higher food sales volume than convenience stores.

General access to food is modeled by the combination of both store types, whereas access to quality food is modeled by grocery stores and access to inferior food is modeled by convenience stores.

The concept of food deserts can either be modeled by the general access to food or more specifically by the access to quality food. This study aims to evaluate both concepts. Log relative risk ratios of convenience stores against grocery stores identify spatial imbalances between both food sources.

One needs to bear in mind that the location of stores also has to follow the potential demand. Locations of insufficient demand cannot support food stores.

Several factors are ignored in this study:

Evaluating the Concept of Weighted Kernel Densities

Kernel weights effectively increase or decrease the number of points or a fraction thereof at a location relative to the weights at other locations.

At any given distance from the point location just the intensity is modified by weights but not the spatial reach of the density estimator. Thus, weights have no impact the bandwidth.

The graph below compares the density at a location with a weight of 1 against that at a location with a weight of 3.

###################################################
### Impact of weights on Gaussian Kernel Density
###################################################
x <- c(-5,4.9,5,5.1)
nx <-length(x)
bw <-1
fx <- density(x , bw=bw, kernel="gaussian")
plot(fx,
     main="Impact on Weights on Gaussian Kernel Densities\nweight=1 versus weight=3 (equivalent to 3 points)")
abline(h=0)
points(x, rep(0, nx), col="red", cex=1.5, pch=20)
abline(v=c(-5-bw,-5+bw, 5-bw,5+bw), lty=2)

Derivation of the Log-Relative Risk Ratios

The following figure calculates the densities for two store types and the evaluates the log relative risk ratio of both densities. Since the convenience stores are in the numerator of the log risk ratio, the log relative risk ratio is positive if the density of the convenience stores is higher than that of the grocery stores.

Note that this example does not make use of differential weights and varying bandwidths.

######################################################
### Grovery versus Convenience Stores Kernel Densities
######################################################
bw <- 1
groc <- c(-1,-1.5,-2)
conv <- c(1,1.5,2)
both <- c(groc,conv)

f.both <- density(both , bw=bw, kernel="gaussian")
f.both$y <- f.both$y+0.006
f.groc <- density(groc , bw=bw, kernel="gaussian",
                  from=min(f.both$x), to=max(f.both$x))
f.groc$y <- f.groc$y/2
f.conv <- density(conv , bw=bw, kernel="gaussian",
                  from=min(f.both$x), to=max(f.both$x))
f.conv$y <- f.conv$y/2
f.risk <- log(f.conv$y/f.groc$y)/100

plot(f.both, ylim=c(-0.2,0.2), type="n",
     xlab="Store Locations & Cumulative Market Areas (bw=1)", ylab="Density & log(RR)",
     main="Market Areas by Store Types &\nlog(Relative Risks)")
lines(f.conv, lwd=4, col="red")
lines(f.groc, lwd=4, col="green")
lines(f.both, lwd=4, col="blue")
lines(f.both$x, f.risk, col="purple", lwd=4, lty=2)
abline(h=0, lwd=4)

points(conv, rep(0, length(conv)), col="red", cex=4, pch=20)
points(groc, rep(0, length(groc)), col="green", cex=4, pch=20)

legend("bottomright",
       legend = c("Grocery Stores",
                  "Convenience Stores",
                  "Joint Both Stores",
                  "log(Relative Risk)"),
       col=c("green","red","blue","purple"), lwd=rep(4,4),
       bty="o", title="Store Types & log(RR)")

Function densityKernelFct

Since kernel densities are calculate on a fixed raster grid, the pixel values need to be aggregated into pixelated census tracts. Both the grid dimensions and the pixel dimensions of the census tracts need to match. This function is a slow performer due to the inefficient use of a for loop.

densityKernelFct <- function(kernel, tess){
  ##
  ## Split a kernel density image into individual tessellations (also type image)
  ## and calculated for each tessellation descriptive statistics of the kernel
  ## density pixel values
  ##
  tract <- split(kernel, tess) # see Baddeley p 122
  #sapply(tract, integral.im)  # kernel sum in tract
  n <- length(tract)           # number of tracts
  ## Initialized data-frame of results
  df <- data.frame(SeqId=1:n, nOfCells=NA, Sum=NA, Density=NA, varDens=NA)
  ## Cycle over tracts. Data values are in $v and may include NAs.
  ## $v is a matrix.
  for (i in 1:n){
    Sum <- sum(tract[[i]]$v, na.rm = T)
    Cells <- length(na.omit(as.vector(tract[[i]]$v)))
    Density <- mean(as.vector(tract[[i]]$v), na.rm = T)
    varDens <- var(as.vector(tract[[i]]$v), na.rm = T)
    df[i,2:5] <- c(Cells, Sum, Density, varDens)
  }
  return(df)
} ## end::densityKernelFct

Projecting the Geographic Locations into the UTM System

The original maps are in the long-lat format. In order to perform isotropic distance calculations the maps need to be re-projected into the UTM coordinates system. Note that Dallas is in the UTM zone 14.

Lastly the point coordinates of the stores need to be converted into ppp-objects of the spatstat package.

## Check tracts unique sequence Id
plot(tractShp, axes=T, main="Check SeqId")
text(coordinates(tractShp), as.character(tractShp$SeqId), cex=0.5)


## Change the map coordinate system to UTM zone 14 (best fit for Texas)
proj4string(bndShp)                               # Current projection system
R> [1] "+proj=longlat +ellps=GRS80 +no_defs"
projUTM <- CRS("+proj=utm +zone=14  +units=m")    # New coordinate sytem definition
bndUTM <- spTransform(bndShp, projUTM)            # Re-project boundary
ptsUTM <- spTransform(foodStoresShp, projUTM)     # Re-project store locations
tractUTM <- spTransform(tractShp, projUTM)        # Re-project tracts

## Converting UTM shape-files to ppp object and assign marks
ptsDf <- as.data.frame(ptsUTM)
pts <- as.ppp(ptsUTM)
pts$marks <- NULL
pts$marks <- ptsDf$STORETYPE
#pts$marks <- factor(ptsDf$STORETYPE, labels=c("healhy","junk"))

Definition of Pixel Image Dimensions for spatstat

This code chunk sets the resolution of the pixels and converts the census tracts into a tessellation object. These tract tessellations are subsequently transformed into a raster image of matching dimensions.

## Define window with units and pixel resolution
win <- as.mask(as.owin(bndUTM), eps=200)   # Width and height of pixels = 200 m
unitname(win) <- list("meter","meters")
pts <- pts[win]        # assign owin to pts
summary(pts)
R> Marked planar point pattern:  1623 points
R> Average intensity 6.891603e-07 points per square meter
R> 
R> *Pattern contains duplicated points*
R> 
R> Coordinates are given to 1 decimal place
R> i.e. rounded to the nearest multiple of 0.1 meters
R> 
R> Multitype:
R>              frequency proportion    intensity
R> Healthy Food       616  0.3795441 2.615667e-07
R> Junk Food         1007  0.6204559 4.275936e-07
R> 
R> binary image mask
R> 246 x 244 pixel array (ny, nx)
R> pixel size: 200 by 200 meters
R> enclosing rectangle: [683852.5, 732590.5] x [3602947, 3652068] meters
R>                      (48740 x 49120 meters)
R> Window area = 2355040000 square meters
R> Unit of length: 1 meter
R> Fraction of frame area: 0.984

## Convert tracts to tesselation object
tractDf <- as.data.frame(tractUTM)          # Save the attribute informaton
tracts <- slot(tractUTM, "polygons")
tracts <- lapply(tracts, function(x) { SpatialPolygons(list(x)) })
tracts <- lapply(tracts, as.owin)
tractsTess <- tess(tiles=tracts)

## Converting tesselation to grid image with 200x200 meters cell size
tractsTessIm <- tess(tiles=tracts, window=win)

Evaluating Weighted Store Kernel Densities

The food supply densities in Dallas County are evaluated for grocery stores, convenience stores and the overall food supply by merging both store types. Each store is weighted by it food related sales volume.

It is assumed that each grocery store’s market area is larger than that of convenience stores. Thus the band width of grocery stores is 4000 meters and that of convenience stores only 2000 meters. The compromise band width of 3000 meters is selected when evaluating both store types together. Other bandwidths for the store’s market areas can be considered.

The Gaussian kernel function implies that approximately 68 % of the customers travel less than one band width unit to the closest store.

##############################################################################
## Sample code evaluating the food store market areas by food sales weighted
## kernel density estimates
##############################################################################

## subset cases and controls
groc <- subset(pts, marks=="Healthy Food")
conv <- subset(pts, marks=="Junk Food")
grocFoodSales <- ptsDf[ptsDf$STORETYPE=="Healthy Food","FOODSALES"]
convFoodSales <- ptsDf[ptsDf$STORETYPE=="Junk Food","FOODSALES"]

Grocery Stores

Convenience Stores

All Stores (Convenience and Grocery Stores)

Calculating the Weighted Log-Relative Risk Surface

Here the relative risk ratios per pixel and census tract are calculated and the significance of over- or under-supply of convenience stores relative to grocery stores are evaluated by a randomization algorithm.

Note that in the sparsely populated South-East sector of Dallas just a few convenience stores are located along I45 and no grocery store. This leads to an unrealistic large log-relative risk ratio. In contrast the Park Cities in the center of Dallas County lack convenience stores thus making the log-relative risk ratio significantly negative.

Log-Relative Risk Surface

Evaluation of Significance

Evaluating the General Store Densities Relative to Census Tract Attributes

The last code chunk models the general food supply within census tracts using a
set of socio-economic and demographic variables.

First, however, those census tracts without population need to be excluded from the analysis. These are the airports of Love Field and DFW International.

## Identify two tracts without night population and map them
tractShp@data[tractShp$NIGHTPOP <= 0, "SeqId" ]
R> [1] 154 248
tractShp$airport <- 0
tractShp$airport[c(154,248)] <- 1
tractShp$airport <- factor(tractShp$airport, labels=c("no","airport"))

mapColorQual(tractShp$airport, tractShp, legend.cex=0.7, legend.title = "Airport\nTracts",
             map.title = "Airport Census Tracts\nto be Excluded from the Analysis")
plot(lakesShp, col="blue", border="blue",add=T)
plot(hwyShp, col="red", lwd=3, add=T)


## Remove both records from tractShp
tractShp <- tractShp[-c(154,248),]

The following tasks are performed:

  1. A scatter-plot matrix explores the distribution of the variables.

  2. A preliminary regression model aims at explaining the distribution of the food supply density.

  3. Non-linearities are evaluated with the residualPlot function. It suggests that log(PCTDAYPOP) should be specified as a quadratic function.

  4. Effect plots visualize the non-linear relationships.

  5. The function influenceIndexPlot checks for outliers and influential observations.

##
## Save augmented (with Density Estimates) Census Tract shape file
##
# rgdal::writeOGR(tractShp, dsn=getwd(), layer="DallasFoodTracts",
#                 overwrite_layer=T, driver="ESRI Shapefile")
#
# tractShp <- rgdal::readOGR(dsn=".", layer="DallasFoodTracts",
#                            integer64="warn.loss")

## Preliminary regression model
car::scatterplotMatrix(~all3000D+log(POPDEN)+log(PCTDAYPOP)+log(PCTBLACK+1)+PCTUNIVDEG+
                         PCTBADENG +log(PCTHUVAC+1)+log(PCTNOVEH+1)+PCTNOHINS, data=tractShp,
                       main="Potential Variables Impacting the Store-Density",
                       pch=20, smooth=list(span = 0.35,lty.smooth=1, col.smooth="red", col.var="red"),
                       regLine=list(col="green"))

all3000.lm1 <- lm(all3000D~log(POPDEN)+log(PCTDAYPOP)+log(PCTBLACK+1)+PCTUNIVDEG+
                           PCTBADENG +log(PCTHUVAC+1)+log(PCTNOVEH+1)+PCTNOHINS, 
                           data=tractShp)
summary(all3000.lm1)
R> 
R> Call:
R> lm(formula = all3000D ~ log(POPDEN) + log(PCTDAYPOP) + log(PCTBLACK + 
R>     1) + PCTUNIVDEG + PCTBADENG + log(PCTHUVAC + 1) + log(PCTNOVEH + 
R>     1) + PCTNOHINS, data = tractShp)
R> 
R> Residuals:
R>     Min      1Q  Median      3Q     Max 
R> -3.7297 -0.7491 -0.0566  0.7129  4.1508 
R> 
R> Coefficients:
R>                    Estimate Std. Error t value Pr(>|t|)    
R> (Intercept)       -4.365164   0.772808  -5.648 2.67e-08 ***
R> log(POPDEN)        0.872422   0.087116  10.014  < 2e-16 ***
R> log(PCTDAYPOP)    -0.418317   0.165501  -2.528 0.011781 *  
R> log(PCTBLACK + 1) -0.501326   0.071020  -7.059 5.42e-12 ***
R> PCTUNIVDEG         0.039334   0.005340   7.366 6.98e-13 ***
R> PCTBADENG         -0.030244   0.008264  -3.660 0.000278 ***
R> log(PCTHUVAC + 1)  0.452220   0.099835   4.530 7.34e-06 ***
R> log(PCTNOVEH + 1)  0.534643   0.088928   6.012 3.46e-09 ***
R> PCTNOHINS          0.040413   0.010837   3.729 0.000213 ***
R> ---
R> Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
R> 
R> Residual standard error: 1.344 on 518 degrees of freedom
R> Multiple R-squared:  0.5714, Adjusted R-squared:  0.5647 
R> F-statistic: 86.31 on 8 and 518 DF,  p-value: < 2.2e-16
car::vif(all3000.lm1)
R>       log(POPDEN)    log(PCTDAYPOP) log(PCTBLACK + 1)        PCTUNIVDEG 
R>          1.492422          1.232416          1.851127          5.001884 
R>         PCTBADENG log(PCTHUVAC + 1) log(PCTNOVEH + 1)         PCTNOHINS 
R>          3.536014          1.237699          1.629145          4.634050
car::residualPlots(all3000.lm1)

R>                   Test stat Pr(>|Test stat|)    
R> log(POPDEN)          0.1605          0.87252    
R> log(PCTDAYPOP)      -4.1933        3.235e-05 ***
R> log(PCTBLACK + 1)    0.3175          0.75099    
R> PCTUNIVDEG           0.0815          0.93507    
R> PCTBADENG            0.0350          0.97210    
R> log(PCTHUVAC + 1)    0.5037          0.61470    
R> log(PCTNOVEH + 1)    1.8692          0.06216 .  
R> PCTNOHINS           -1.4602          0.14484    
R> Tukey test           2.3087          0.02096 *  
R> ---
R> Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
all3000.lm2 <- lm(all3000D~log(POPDEN)+log(PCTDAYPOP)+I(log(PCTDAYPOP)^2)+log(PCTBLACK+1)+
                           PCTUNIVDEG+PCTBADENG +log(PCTHUVAC+1)+log(PCTNOVEH+1)+PCTNOHINS, 
                           data=tractShp)
summary(all3000.lm2)
R> 
R> Call:
R> lm(formula = all3000D ~ log(POPDEN) + log(PCTDAYPOP) + I(log(PCTDAYPOP)^2) + 
R>     log(PCTBLACK + 1) + PCTUNIVDEG + PCTBADENG + log(PCTHUVAC + 
R>     1) + log(PCTNOVEH + 1) + PCTNOHINS, data = tractShp)
R> 
R> Residuals:
R>     Min      1Q  Median      3Q     Max 
R> -4.0738 -0.7509 -0.0558  0.6909  4.1917 
R> 
R> Coefficients:
R>                       Estimate Std. Error t value Pr(>|t|)    
R> (Intercept)         -22.927526   4.491596  -5.105 4.67e-07 ***
R> log(POPDEN)           0.953115   0.087887  10.845  < 2e-16 ***
R> log(PCTDAYPOP)        9.489680   2.368444   4.007 7.06e-05 ***
R> I(log(PCTDAYPOP)^2)  -1.350698   0.322111  -4.193 3.24e-05 ***
R> log(PCTBLACK + 1)    -0.456872   0.070709  -6.461 2.40e-10 ***
R> PCTUNIVDEG            0.038387   0.005261   7.296 1.12e-12 ***
R> PCTBADENG            -0.029891   0.008135  -3.674 0.000263 ***
R> log(PCTHUVAC + 1)     0.501473   0.098974   5.067 5.64e-07 ***
R> log(PCTNOVEH + 1)     0.472945   0.088766   5.328 1.48e-07 ***
R> PCTNOHINS             0.034231   0.010769   3.179 0.001568 ** 
R> ---
R> Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
R> 
R> Residual standard error: 1.323 on 517 degrees of freedom
R> Multiple R-squared:  0.5855, Adjusted R-squared:  0.5782 
R> F-statistic: 81.13 on 9 and 517 DF,  p-value: < 2.2e-16
plot(allEffects(all3000.lm2))


car::influenceIndexPlot(all3000.lm2, id=list(n=3))