Food Access Grand Rapids

Data Modifications

Loading packages

#install.packages("readxl")
library(reshape2)
library(readxl)

## Warning: package 'readxl' was built under R version 4.3.1

library(corrplot)

## corrplot 0.92 loaded

library(plotly)

## Loading required package: ggplot2

## 
## Attaching package: 'plotly'

## The following object is masked from 'package:ggplot2':
## 
##     last_plot

## The following object is masked from 'package:stats':
## 
##     filter

## The following object is masked from 'package:graphics':
## 
##     layout

library(ggthemes)
library(ggdist)
library(ggplot2)
library(psych)

## 
## Attaching package: 'psych'

## The following objects are masked from 'package:ggplot2':
## 
##     %+%, alpha

library(epiDisplay)

## Loading required package: foreign

## Loading required package: survival

## Loading required package: MASS

## 
## Attaching package: 'MASS'

## The following object is masked from 'package:plotly':
## 
##     select

## Loading required package: nnet

## 
## Attaching package: 'epiDisplay'

## The following objects are masked from 'package:psych':
## 
##     alpha, cs, lookup

## The following object is masked from 'package:ggplot2':
## 
##     alpha

library(dplyr)

## 
## Attaching package: 'dplyr'

## The following object is masked from 'package:MASS':
## 
##     select

## The following objects are masked from 'package:stats':
## 
##     filter, lag

## The following objects are masked from 'package:base':
## 
##     intersect, setdiff, setequal, union

setwd("C:/Users/ianmt/Documents/R")
ACSOutletJoin <- read.csv("Census_Outlets_Join.csv")
#View(ACSOutletJoin)
ACSOutletRecode <- read.csv("Census_RecodeType_Join.csv")
#View(ACSOutletRecode)
reACSOutletJoin <- dcast(data = ACSOutletJoin, formula = GEOID ~ FoodOutletType, value.var = "Join_Count")

## Aggregation function missing: defaulting to length

#View(reACSOutletJoin)
#write.csv(reACSOutletJoin, "C:/Users/ianmt/Documents/R\\ACSOutletCast.csv", row.names = FALSE)
reACSOutletRecode <- dcast(data = ACSOutletRecode, formula = GEOID ~ Type_Recode, value.var = "Join_Count")

## Aggregation function missing: defaulting to length

#View(reACSOutletRecode)
#write.csv(reACSOutletRecode, "C:/Users/ianmt/Documents/R\\ACSRecodeCast.csv", row.names = FALSE)

Adding percentages & categories

ACSOutletZip <- read.csv("ACSOutlet_Zip.csv")
#View(ACSOutletZip)
PWhite = (ACSOutletZip$White/ACSOutletZip$TotPop)*100
PBlack = (ACSOutletZip$Black/ACSOutletZip$TotPop)*100
PNativeAmerican = (ACSOutletZip$NativeAmerican/ACSOutletZip$TotPop)*100
PAsian = (ACSOutletZip$Asian/ACSOutletZip$TotPop)*100
PHawaiianPI = (ACSOutletZip$HawaiianPI/ACSOutletZip$TotPop)*100
POtherRace = (ACSOutletZip$OtherRace/ACSOutletZip$TotPop)*100
PHispanic = (ACSOutletZip$Hispanic/ACSOutletZip$TotPop)*100
PPOC = ((ACSOutletZip$TotPop-ACSOutletZip$White)/ACSOutletZip$TotPop)*100
Incomeadj = (ACSOutletZip$MedianHHIncome)^(1/2)
PSNAP = (ACSOutletZip$HHSNAP/ACSOutletZip$TotalHouseholds)*100
PRegHS = (ACSOutletZip$RegHS/ACSOutletZip$X25pTotal)*100
PGED = (ACSOutletZip$GED/ACSOutletZip$X25pTotal)*100
PBachelor = (ACSOutletZip$Bachelor/ACSOutletZip$X25pTotal)*100
PProfessional = (ACSOutletZip$Professional/ACSOutletZip$X25pTotal)*100
PopDensity = ACSOutletZip$TotPop/(ACSOutletZip$Shape_Area/1000)
PClosure = (ACSOutletZip$Closures/(ACSOutletZip$Total+1))*100
GRData=data.frame(ACSOutletZip,PWhite,PBlack,PNativeAmerican,PAsian,PHawaiianPI,POtherRace,PHispanic,PPOC,Incomeadj,PSNAP,PRegHS,PGED,PBachelor,PProfessional,PopDensity,PClosure)
GRData$BlackCat[GRData$PBlack<=13.11] <- "under"
GRData$BlackCat[GRData$PBlack>13.11] <- "over"
GRData$HispanicCat[GRData$PHispanic<=13.21] <- "under"
GRData$HispanicCat[GRData$PHispanic>13.21] <- "over"
##GRData$EduCat[GRData$PBachelor<=30] <- "low"
##GRData$EduCat[GRData$PBachelor>30] <- "high"
GRData$POCCat[GRData$PPOC<=15] <- "VL"
GRData$POCCat[((GRData$PPOC>15)&(GRData$PPOC<=30))] <- "L"
GRData$POCCat[((GRData$PPOC>30)&(GRData$PPOC<=45))] <- "H"
GRData$POCCat[GRData$PPOC>45] <- "VH"
#View(GRData)
ACSreTypeZip <- read.csv("ACSRecodeType_Zip.csv")
PWhite = (ACSreTypeZip$White/ACSreTypeZip$TotPop)*100
PBlack = (ACSreTypeZip$Black/ACSreTypeZip$TotPop)*100
PNativeAmerican = (ACSreTypeZip$NativeAmerican/ACSreTypeZip$TotPop)*100
PAsian = (ACSreTypeZip$Asian/ACSreTypeZip$TotPop)*100
PHawaiianPI = (ACSreTypeZip$HawaiianPI/ACSreTypeZip$TotPop)*100
POtherRace = (ACSreTypeZip$OtherRace/ACSreTypeZip$TotPop)*100
PHispanic = (ACSreTypeZip$Hispanic/ACSreTypeZip$TotPop)*100
PPOC = ((ACSreTypeZip$TotPop-ACSreTypeZip$White)/ACSreTypeZip$TotPop)*100
Incomeadj = (ACSreTypeZip$MedianHHIncome)^(1/2)
PSNAP = (ACSreTypeZip$HHSNAP/ACSreTypeZip$TotalHouseholds)*100
PRegHS = (ACSreTypeZip$RegHS/ACSreTypeZip$X25pTotal)*100
PGED = (ACSreTypeZip$GED/ACSreTypeZip$X25pTotal)*100
PBachelor = (ACSreTypeZip$Bachelor/ACSreTypeZip$X25pTotal)*100
PProfessional = (ACSreTypeZip$Professional/ACSreTypeZip$X25pTotal)*100
PopDensity = ACSreTypeZip$TotPop/(ACSreTypeZip$Shape_Area/1000)
PClosure = (ACSreTypeZip$Closures/(ACSreTypeZip$Total+1))*100
GRDataRe=data.frame(ACSreTypeZip,PWhite,PBlack,PNativeAmerican,PAsian,PHawaiianPI,POtherRace,PHispanic,PPOC,Incomeadj,PSNAP,PRegHS,PGED,PBachelor,PProfessional,PopDensity,PClosure)
GRDataRe$BlackCat[GRData$PBlack<=13.11] <- "under"
GRDataRe$BlackCat[GRData$PBlack>13.11] <- "over"
GRDataRe$HispanicCat[GRData$PHispanic<=13.21] <- "under"
GRDataRe$HispanicCat[GRData$PHispanic>13.21] <- "over"
##GRDataRe$EduCat[GRData$PBachelor<=30] <- "low"
##GRDataRe$EduCat[GRData$PBachelor>30] <- "high"
GRDataRe$POCCat[GRData$PPOC<=15] <- "VL"
GRDataRe$POCCat[((GRData$PPOC>15)&(GRData$PPOC<=30))] <- "L"
GRDataRe$POCCat[((GRData$PPOC>30)&(GRData$PPOC<=45))] <- "H"
GRDataRe$POCCat[GRData$PPOC>45] <- "VH"
#View(GRDataRe)

Census aggregations

tp <- GRData$TotPop
hi <- GRData$MedianHHIncome
x25 <- GRData$X25pTotal
bt <- GRData$Black
ht <- GRData$Hispanic
wt<- GRData$White
bat<- GRData$Bachelor
sh<- GRData$HHSNAP
th <- GRData$TotalHouseholds
pocc <- GRData$POCCat
bc <- GRData$BlackCat
hc <- GRData$HispanicCat

agg1a <- aggregate(tp ~ pocc , FUN=sum)
agg1a

##   pocc     tp
## 1    H 101796
## 2    L  94150
## 3   VH  81859
## 4   VL 147188

agg1b <- aggregate(tp ~ bc , FUN=sum)
agg1b

##      bc     tp
## 1  over 146540
## 2 under 278453

agg1c <- aggregate(bt ~ hc , FUN=sum)
agg1c

##      hc    bt
## 1  over 20312
## 2 under 34458

agg2a <- aggregate(bt ~ pocc , FUN=sum)
agg2a

##   pocc    bt
## 1    H 17793
## 2    L  8507
## 3   VH 25352
## 4   VL  3118

agg2b <- aggregate(bt ~ bc , FUN=sum)
agg2b

##      bc    bt
## 1  over 41069
## 2 under 13701

agg2c <- aggregate(ht ~ hc , FUN=sum)
agg2c

##      hc    ht
## 1  over 41720
## 2 under 13564

agg2d <- aggregate(ht ~ pocc , FUN=sum)
agg2d

##   pocc    ht
## 1    H 15603
## 2    L 10582
## 3   VH 24061
## 4   VL  5038

agg2e <- aggregate(wt ~ pocc , FUN =sum)
agg2e

##   pocc     wt
## 1    H  65875
## 2    L  73223
## 3   VH  34724
## 4   VL 135540

agg3a <- aggregate(hi ~ pocc , FUN =mean)
agg3a

##   pocc       hi
## 1    H 55666.25
## 2    L 55467.62
## 3   VH 43746.90
## 4   VL 83694.55

agg3b <- aggregate(hi ~ bc , FUN =mean)
agg3b

##      bc       hi
## 1  over 51640.89
## 2 under 68384.60

agg3c <- aggregate(hi ~ hc , FUN =mean)
agg3c

##      hc       hi
## 1  over 46228.91
## 2 under 70272.31

agg4a <- aggregate(x25 ~ pocc , FUN =sum)
agg4a

##   pocc    x25
## 1    H  67232
## 2    L  62563
## 3   VH  47425
## 4   VL 101161

agg4b <- aggregate(x25 ~ bc , FUN =sum)
agg4b

##      bc    x25
## 1  over  93492
## 2 under 184889

agg4c <- aggregate(x25 ~ hc , FUN =sum)
agg4c

##      hc    x25
## 1  over  79225
## 2 under 199156

agg5a <- aggregate(bat ~ pocc , FUN =sum)
agg5a

##   pocc   bat
## 1    H 14619
## 2    L 16724
## 3   VH  7516
## 4   VL 30530

agg5b <- aggregate(bat ~ bc , FUN =sum)
agg5b

##      bc   bat
## 1  over 22105
## 2 under 47284

agg5c <- aggregate(bat ~ hc , FUN =sum)
agg5c

##      hc   bat
## 1  over 11310
## 2 under 58079

agg6a <- aggregate(th ~ pocc , FUN =sum)
agg6a

##   pocc    th
## 1    H 39266
## 2    L 39419
## 3   VH 26825
## 4   VL 55935

agg6b <- aggregate(th ~ bc , FUN =sum)
agg6b

##      bc     th
## 1  over  55655
## 2 under 105790

agg6c <- aggregate(th ~ hc , FUN =sum)
agg6c

##      hc     th
## 1  over  45961
## 2 under 115484

agg7a <- aggregate(sh ~ pocc , FUN =sum)
agg7a

##   pocc   sh
## 1    H 4711
## 2    L 4752
## 3   VH 6677
## 4   VL 2604

agg7b <- aggregate(sh ~ bc , FUN =sum)
agg7b

##      bc   sh
## 1  over 9716
## 2 under 9028

agg7c <- aggregate(sh ~ hc , FUN =sum)
agg7c

##      hc   sh
## 1  over 8818
## 2 under 9926

Food Outlets

Types & count

OutletType <- read.csv("Outlets_Type.csv")
#View(OutletType)
OTgg1 <- ggplot(OutletType,
       aes(x = OutletCode,
           y = FREQUENCY,
           fill = FoodOutletType))+
  geom_col(stat = "identity")+
  scale_fill_viridis_d(name = "Group", option = "turbo")+
  labs(title = "Outlet", y = "Count", x = "Group")

## Warning in geom_col(stat = "identity"): Ignoring unknown parameters: `stat`

OTgg1 + theme_dark()

OutletType2 <- read.csv("RecodeType_Count.csv")
OTgg2 <- ggplot(OutletType2,
       aes(x = Type_Recode,
           y = FREQUENCY,
           fill = Label))+
  geom_col(stat = "identity")+
  scale_fill_viridis_d(name = "Group", option = "turbo")+
  labs(title = "Outlet", y = "Count", x = "Group")

## Warning in geom_col(stat = "identity"): Ignoring unknown parameters: `stat`

OTgg2 + theme_dark()

Outlet histograms

summary(GRDataRe$CLG)

##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
##   0.000   0.250   2.000   2.363   3.750  13.000

sd(GRDataRe$CLG)

## [1] 2.311355

sum(GRDataRe$CLG)

## [1] 241

GRDRh1<- ggplot(GRDataRe,
       aes(x=CLG))+
  geom_histogram(binwidth=1,color="white",fill="#13bfab")+
  labs(title = "Convenience, Corner Stores, Small Groceries, Liquor Stores,
and Gas Stations by Census Tract in Grand Rapids, MI", y = "Count", x = "Outlets")
GRDRh1 + theme_dark()

summary(GRDataRe$Restaurant)

##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
##   0.000   0.000   2.000   4.804   6.000  42.000

sd(GRDataRe$Restaurant)

## [1] 6.82459

sum(GRDataRe$Restaurant)

## [1] 490

GRDRh2<- ggplot(GRDataRe,
       aes(x=Restaurant))+
  geom_histogram(binwidth=3,color="white",fill="#13bfab")+
  labs(title = "Fast Food, Full Service, and Take Out Restaurants
by Census Tract in Grand Rapids, MI", y = "Count", x = "Outlets")
GRDRh2 + theme_dark()

summary(GRDataRe$Specialty)

##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
##    0.00    0.00    0.00    1.48    2.00    9.00

sd(GRDataRe$Specialty)

## [1] 2.289582

sum(GRDataRe$Specialty)

## [1] 151

GRDRh3<- ggplot(GRDataRe,
       aes(x=Specialty))+
  geom_histogram(binwidth=1,color="white",fill="#13bfab")+
  labs(title = "Specialty Food Outlets by Census Tract in Grand Rapids, MI", y = "Count", x = "Outlets")
GRDRh3 + theme_dark()

summary(GRData$Convenience)

##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
##  0.0000  0.0000  0.0000  0.9608  1.0000  9.0000

sd(GRData$Convenience)

## [1] 1.427603

sum(GRData$Convenience)

## [1] 98

GRDh4<- ggplot(GRData,
       aes(x=Convenience))+
  geom_histogram(binwidth=1,color="white",fill="#13bfab")+
  labs(title = "Convenience, Corner Stores, and Small Groceries
by Census Tract in Grand Rapids, MI", y = "Count", x = "Outlets")
GRDh4 + theme_dark()

summary(GRData$Pharmacy)

##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
##  0.0000  0.0000  0.0000  0.5196  1.0000  5.0000

sd(GRData$Pharmacy)

## [1] 0.9516051

sum(GRData$Pharmacy)

## [1] 53

GRDh5<- ggplot(GRData,
       aes(x=Pharmacy))+
  geom_histogram(binwidth=1,color="white",fill="#13bfab")+
  labs(title = "Pharmacies and Drug Stores by Census Tract in Grand Rapids, MI", y = "Count", x = "Outlets")
GRDh5 + theme_dark()

summary(GRData$Coffee)

##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
##  0.0000  0.0000  0.0000  0.7353  1.0000  7.0000

sd(GRData$Coffee)

## [1] 1.413699

sum(GRData$Coffee)

## [1] 75

GRDh6<- ggplot(GRData,
       aes(x=Coffee))+
  geom_histogram(binwidth=1,color="white",fill="#13bfab")+
  labs(title = "Coffee, Tea, and Juice Shops by Census Tract
in Grand Rapids, MI", y = "Count", x = "Outlets")
GRDh6 + theme_dark()

summary(GRData$Restaurant)

##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
##    0.00    0.00    1.00    2.99    4.00   35.00

sd(GRData$Restaurant)

## [1] 4.944234

sum(GRData$Restaurant)

## [1] 305

GRDh7<- ggplot(GRData,
       aes(x=Restaurant))+
  geom_histogram(binwidth=3,color="white",fill="#13bfab")+
  labs(title = "Full Service Restaurants in Grand Rapids, MI", y = "Count", x = "Outlets")
GRDh7 + theme_dark()

summary(GRData$FastFood)

##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
##   0.000   0.000   1.000   1.814   2.750  10.000

sd(GRData$FastFood)

## [1] 2.555252

sum(GRData$FastFood)

## [1] 185

GRDh8<- ggplot(GRData,
       aes(x=FastFood))+
  geom_histogram(binwidth=1,color="white",fill="#13bfab")+
  labs(title = "Fast Food Restaurants in Grand Rapids, MI", y = "Count", x = "Outlets")
GRDh8 + theme_dark()

summary(GRData$Bars)

##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
##  0.0000  0.0000  0.0000  0.6765  0.0000 15.0000

sd(GRData$Bars)

## [1] 1.904459

sum(GRData$Bars)

## [1] 69

GRDh9<- ggplot(GRData,
       aes(x=Bars))+
  geom_histogram(binwidth=1,color="white",fill="#13bfab")+
  labs(title = "Bars in Grand Rapids, MI", y = "Count", x = "Outlets")
GRDh9 + theme_dark()

summary(GRData$Gas)

##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
##  0.0000  0.0000  0.0000  0.7941  1.0000  4.0000

sd(GRData$Gas)

## [1] 1.065556

sum(GRData$Gas)

## [1] 81

GRDh10<- ggplot(GRData,
       aes(x=Gas))+
  geom_histogram(binwidth=1,color="white",fill="#13bfab")+
  labs(title = "Gas Stations in Grand Rapids, MI", y = "Count", x = "Outlets")
GRDh10 + theme_dark()

summary(GRData$Liquor)

##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
##   0.000   0.000   0.000   0.598   1.000   3.000

sd(GRData$Liquor)

## [1] 0.7992547

sum(GRData$Liquor)

## [1] 61

GRDh10<- ggplot(GRData,
       aes(x=Liquor))+
  geom_histogram(binwidth=1,color="white",fill="#13bfab")+
  labs(title = "Liquor Stores in Grand Rapids, MI", y = "Count", x = "Outlets")
GRDh10 + theme_dark()

Census Tract Demographics

Demographic histograms

summary(GRData$PPOC)

##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
##    0.00   11.70   25.62   27.88   38.06   73.97

sd(GRData$PPOC)

## [1] 18.90794

GRDha<- ggplot(GRData,
       aes(x=PPOC))+
  geom_histogram(binwidth=10,color="white",fill="#8410de")+
  labs(title = "Pct People of Color by Census Tract in Grand Rapids, MI", y = "Count", x = "Percent")
GRDha + theme_dark()

summary(GRData$PBlack)

##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
##   0.000   2.804   8.437  13.110  18.669  57.987

sd(GRData$PBlack)

## [1] 13.37966

GRDh1<- ggplot(GRData,
       aes(x=PBlack))+
  geom_histogram(binwidth=10,color="white",fill="#8410de")+
  labs(title = "Pct Black & African American Residents
by Census Tract in Grand Rapids, MI", y = "Count", x = "Percent")
GRDh1 + theme_dark()

summary(GRData$PHispanic)

##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
##  0.2112  3.3832  5.8707 13.2116 17.6121 82.8249

sd(GRData$PHispanic)

## [1] 16.83476

GRDh2<- ggplot(GRData,
       aes(x=PHispanic))+
  geom_histogram(binwidth=10,color="white",fill="#8410de")+
  labs(title = "Pct Latinx & Hispanic Residents
by Census Tract in Grand Rapids, MI", y = "Count", x = "Percent")
GRDh2 + theme_dark()

summary(GRData$PSNAP)

##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
##   0.000   4.095   8.452  12.535  18.757  43.348

sd(GRData$PSNAP)

## [1] 10.9205

GRDhb<- ggplot(GRData,
       aes(x=PSNAP))+
  geom_histogram(binwidth=10,color="white",fill="#8410de")+
  labs(title = "Pct Household SNAP Recipients by Census Tract
in Grand Rapids, MI", y = "Count", x = "Percent")
GRDhb + theme_dark()

summary(GRData$MedianHHIncome)

##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max.    NA's 
##   19572   46210   59444   62417   68486  154000       1

sd(GRData$MedianHHIncome, na.rm = TRUE)

## [1] 25292.82

GRDh3<- ggplot(GRData,
       aes(x=MedianHHIncome))+
  geom_histogram(binwidth=5000,color="white",fill="#8410de")+
  labs(title = "Median Household Income by Census Tract in Grand Rapids, MI", y = "Count", x = "Household Income ($)")
GRDh3 + theme_dark()

## Warning: Removed 1 rows containing non-finite values (`stat_bin()`).

summary(GRData$Incomeadj)

##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max.    NA's 
##   139.9   215.0   243.8   245.5   261.7   392.4       1

sd(GRData$Incomeadj)

## [1] NA

summary(GRData$PBachelor)

##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
##  0.8003 15.5341 23.1643 23.7676 32.7417 50.0000

sd(GRData$PBachelor)

## [1] 11.01472

GRDh4<- ggplot(GRData,
       aes(x=PBachelor))+
  geom_histogram(binwidth=7,color="white",fill="#8410de")+
  labs(title = "Pct Bachelor Attainment by Census Tract in Grand Rapids, MI", y = "Count", x = "Percent")
GRDh4 + theme_dark()

summary(GRData$TotalHouseholds)

##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
##     366    1208    1602    1583    1918    2976

sd(GRData$TotalHouseholds)

## [1] 538.7381

GRDh5<- ggplot(GRData,
       aes(x=TotalHouseholds))+
  geom_histogram(binwidth=250,color="white",fill="#8410de")+
  labs(title = "Total Households by Census Tract in Grand Rapids, MI", y = "Count", x = "Households")
GRDh5 + theme_dark()

summary(GRData$TotPop)

##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
##    1510    3215    4166    4167    5124    8397

sd(GRData$TotPop)

## [1] 1457.22

GRDh6<- ggplot(GRData,
       aes(x=TotPop))+
  geom_histogram(binwidth=750,color="white",fill="#8410de")+
  labs(title = "Total Population by Census Tract in Grand Rapids, MI", y = "Count", x = "Population")
GRDh6 + theme_dark()

Trivariate plots

p1 <- GRData %>% 
  plot_ly(x = ~Shape_Area, 
          z = ~Incomeadj, 
          y = ~TotPop,
          color = ~PPOC,
          colors = c("aquamarine3", "goldenrod1")) %>% 
  layout(title = 'Income by Area and Total Population',
         scene = list(xaxis = list(title = 'Area'),
                     yaxis = list(title = 'Total Population'),
                     zaxis = list(title = 'Income Adj ^1/2')))
p1

## No trace type specified:
##   Based on info supplied, a 'scatter3d' trace seems appropriate.
##   Read more about this trace type -> https://plotly.com/r/reference/#scatter3d

## No scatter3d mode specifed:
##   Setting the mode to markers
##   Read more about this attribute -> https://plotly.com/r/reference/#scatter-mode

## Warning: Ignoring 1 observations

p2 <- GRData %>% 
  plot_ly(x = ~PBlack, 
          z = ~Incomeadj, 
          y = ~PBachelor,
          color = ~TotPop,
          colors = c("green3", "goldenrod1")) %>% 
  layout(title = 'Income by Pct Black and African American Residents
and Pct Bachelor Degrees',
         scene = list(xaxis = list(title = 'Pct Black and African American'),
                     yaxis = list(title = 'Pct Bachor Degrees'),
                     zaxis = list(title = 'Income Adjusted ($^1/2)')))
p2

## No trace type specified:
##   Based on info supplied, a 'scatter3d' trace seems appropriate.
##   Read more about this trace type -> https://plotly.com/r/reference/#scatter3d
## No scatter3d mode specifed:
##   Setting the mode to markers
##   Read more about this attribute -> https://plotly.com/r/reference/#scatter-mode

## Warning: Ignoring 1 observations

p3 <- GRData %>% 
  plot_ly(x = ~PHispanic, 
          z = ~Incomeadj, 
          y = ~PBachelor,
          color = ~TotPop,
          colors = c("mediumorchid3", "goldenrod1")) %>% 
  layout(title = 'Income by Pct Hispanic  Residents and Pct Bachelor Degrees',
         scene = list(xaxis = list(title = 'Pct Hispanic'),
                     yaxis = list(title = 'Pct Bachor Degrees'),
                     zaxis = list(title = 'Income Adj ^1/2')))
p3

## No trace type specified:
##   Based on info supplied, a 'scatter3d' trace seems appropriate.
##   Read more about this trace type -> https://plotly.com/r/reference/#scatter3d
## No scatter3d mode specifed:
##   Setting the mode to markers
##   Read more about this attribute -> https://plotly.com/r/reference/#scatter-mode

## Warning: Ignoring 1 observations

POC categorizations

fwHHinPOC <- ggplot(GRData,
       aes(x = PPOC,
           y = MedianHHIncome,
           color = POCCat))+
  geom_jitter()+
  scale_color_manual(values = c("yellow", "cyan", "orangered", "darkblue"))+
  facet_wrap(~POCCat)+
  labs(title = "Median Household Income ~ People of Color (%)",
       x = "Percent",
       y = "Income ($)")
fwHHinPOC + theme_dark()

## Warning: Removed 1 rows containing missing values (`geom_point()`).

anova1 <- aov(Incomeadj~POCCat, data = GRData)
summary(anova1)

##             Df Sum Sq Mean Sq F value   Pr(>F)    
## POCCat       3  89113   29704   22.68 3.31e-11 ***
## Residuals   97 127066    1310                     
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 1 observation deleted due to missingness

fwSNAPPOC <- ggplot(GRData,
       aes(x = PPOC,
           y = PSNAP,
           color = POCCat))+
  geom_jitter()+
  scale_color_manual(values = c("yellow", "cyan", "orangered", "darkblue"))+
  facet_wrap(~POCCat)+
  labs(title = "Household SNAP Recipients (%) ~ People of Color (%)",
       x = "Percent POC",
       y = "Percent SNAP")
fwSNAPPOC + theme_dark()

anova2 <- aov(PSNAP~POCCat, data = GRData)
summary(anova2)

##             Df Sum Sq Mean Sq F value  Pr(>F)    
## POCCat       3   5215  1738.2   24.94 4.5e-12 ***
## Residuals   98   6830    69.7                    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

fwSNAPPOC <- ggplot(GRData,
       aes(x = PPOC,
           y = PBachelor,
           color = POCCat))+
  geom_jitter()+
  scale_color_manual(values = c("yellow", "cyan", "orangered", "darkblue"))+
  facet_wrap(~POCCat)+
  labs(title = "Bachelor Degree Attainment (%) ~ People of Color (%)",
       x = "Percent POC",
       y = "Percent Bachelor Degree Attainment")
fwSNAPPOC + theme_dark()

anova3 <- aov(PBachelor~POCCat, data = GRData)
summary(anova3)

##             Df Sum Sq Mean Sq F value   Pr(>F)    
## POCCat       3   2887   962.4   10.07 7.63e-06 ***
## Residuals   98   9367    95.6                     
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Correlations

Correlation matrices

GRData0 <- GRData[complete.cases(GRData),]
#View(GRData0)
GRData1 <- GRData0[c(5,8,24,25,28,31,36,40,50,58,66,67,69,70,75:78,81,83)]
CorGRData1 <- cor(GRData1)
corrplot(CorGRData1, is.corr = FALSE, method = "square")

GRDataRe0 <- GRDataRe[complete.cases(GRDataRe),]
#View(GRDataRe0)
GRDataRe1 <- GRDataRe0[c(3:13,21,29,30,32,33,38:41,44,46)]
CorGRDataRe1 <- cor(GRDataRe1)
corrplot(CorGRDataRe1, is.corr = FALSE, method = "square")

Correlation testing

#cor.test(GRData$PHispanic, GRData$FastFood, use = "complete obs", method = "pearson")
#GRDcor1 <- ggplot(GRData,
#       aes(x = PHispanic,
#           y = FastFood,
#           color = TotPop))+
#  geom_point(alpha = 0.6, color="white",)+
#  geom_smooth(color="#8410de")+
# labs(title = "Fast Food
#~ Pct Latinx & Hispanic Residents within Census Tracts",
#       x = "Latinx & Hispanic Residents (%)",
#       y = "Outlets")
#GRDcor1 + theme_dark()

#cor.test(GRData$PHispanic, GRData$Bars, use = "complete obs", method = "pearson")
#GRDcor2 <- ggplot(GRData,
#       aes(x = PHispanic,
#           y = Bars,
#           color = TotPop))+
#  geom_point(alpha = 0.6, color="white",)+
#  geom_smooth(color="#8410de")+
# labs(title = "Bars
#~ Pct Latinx & Hispanic Residents within Census Tracts",
#       x = "Latinx & Hispanic Residents (%)",
#       y = "Outlets")
#GRDcor2 + theme_dark()

#cor.test(GRData$PHispanic, GRData$Coffee, use = "complete obs", method = "pearson")
#GRDcor3 <- ggplot(GRData,
#       aes(x = PHispanic,
#           y = Coffee,
#           color = TotPop))+
#  geom_point(alpha = 0.6, color="white",)+
#  geom_smooth(color="#8410de")+
# labs(title = "Coffee, Tea, and Juice Shops
#~ Pct Latinx & Hispanic Residents within Census Tracts",
#       x = "Latinx & Hispanic Residents (%)",
#       y = "Outlets")
#GRDcor3 + theme_dark()

#cor.test(GRData$PHispanic, GRData$Restaurant, use = "complete obs", method = "pearson")
#GRDcor4 <- ggplot(GRData,
#       aes(x = PHispanic,
#           y = Restaurant,
#           color = TotPop))+
#  geom_point(alpha = 0.6, color="white",)+
#  geom_smooth(color="#8410de")+
# labs(title = "Full Service Restaurants
#~ Pct Latinx & Hispanic Residents within Census Tracts",
#       x = "Latinx & Hispanic Residents (%)",
#       y = "Outlets")
#GRDcor4 + theme_dark()

#cor.test(GRData$PHispanic, GRData$Gas, use = "complete obs", method = "pearson")
#GRDcor5 <- ggplot(GRData,
#       aes(x = PHispanic,
#           y = Gas,
#           color = TotPop))+
#  geom_point(alpha = 0.6, color="white",)+
#  geom_smooth(color="#8410de")+
# labs(title = "Gas Stations
#~ Pct Latinx & Hispanic Residents within Census Tracts",
#       x = "Latinx & Hispanic Residents (%)",
#       y = "Outlets")
#GRDcor5 + theme_dark()

#cor.test(GRData$PHispanic, GRData$Hotel, use = "complete obs", method = "pearson")
#GRDcor6 <- ggplot(GRData,
#       aes(x = PHispanic,
#           y = Hotel,
#           color = TotPop))+
#  geom_point(alpha = 0.6, color="white",)+
#  geom_smooth(color="#8410de")+
# labs(title = "Hotels
#~ Pct Latinx & Hispanic Residents within Census Tracts",
#       x = "Latinx & Hispanic Residents (%)",
#       y = "Outlets")
#GRDcor6 + theme_dark()

#cor.test(GRData$PHispanic, GRData$Liquor, use = "complete obs", method = "pearson")
#GRDcor7 <- ggplot(GRData,
#       aes(x = PHispanic,
#           y = Liquor,
#           color = TotPop))+
#  geom_point(alpha = 0.6, color="white",)+
#  geom_smooth(color="#8410de")+
# labs(title = "Liquor and Party Stores
#~ Pct Latinx & Hispanic Residents within Census Tracts",
#       x = "Latinx & Hispanic Residents (%)",
#       y = "Outlets")
#GRDcor7 + theme_dark()

cor.test(GRData$PHispanic, GRData$Convenience, use = "complete obs", method = "pearson")

## 
##  Pearson's product-moment correlation
## 
## data:  GRData$PHispanic and GRData$Convenience
## t = 3.6721, df = 100, p-value = 0.0003885
## alternative hypothesis: true correlation is not equal to 0
## 95 percent confidence interval:
##  0.1610225 0.5053044
## sample estimates:
##       cor 
## 0.3447031

GRDcor8 <- ggplot(GRData,
       aes(x = PHispanic,
           y = Convenience,
           color = TotPop))+
  geom_point(alpha = 0.6, color="white",)+
  geom_smooth(color="#8410de")+
 labs(title = "Convenience, Corner Stores, and Small Groceries
~ Pct Latinx & Hispanic Residents within Census Tracts",
       x = "Latinx & Hispanic Residents (%)",
       y = "Outlets")
GRDcor8 + theme_dark()

## `geom_smooth()` using method = 'loess' and formula = 'y ~ x'

cor.test(GRData$PHispanic, GRData$Pharmacy, use = "complete obs", method = "pearson")

## 
##  Pearson's product-moment correlation
## 
## data:  GRData$PHispanic and GRData$Pharmacy
## t = -2.18, df = 100, p-value = 0.0316
## alternative hypothesis: true correlation is not equal to 0
## 95 percent confidence interval:
##  -0.3912646 -0.0193224
## sample estimates:
##        cor 
## -0.2129969

GRDcor9 <-ggplot(GRData,
       aes(x = PHispanic,
           y = Pharmacy,
           color = TotPop))+
  geom_point(alpha = 0.6, color="white")+
  geom_smooth(color="#8410de")+
 labs(title = "Pharmacies and Drug Stores ~ Pct Latinx & Hispanic
Residents within Census Tracts",
       x = "Latinx & Hispanic Residents (%)",
       y = "Outlets")
GRDcor9 + theme_dark()

## `geom_smooth()` using method = 'loess' and formula = 'y ~ x'

#cor.test(GRData$PBlack, GRData$FastFood, use = "complete obs", method = "pearson")
#GRDcor10 <- ggplot(GRData,
#       aes(x = PBlack,
#           y = FastFood,
#           color = TotPop))+
#  geom_point(alpha = 0.6, color="white",)+
#  geom_smooth(color="#8410de")+
# labs(title = "Fast Food
#~ Pct Black & African American Residents within Census Tracts",
#       x = "Black & African American Residents (%)",
#       y = "Outlets")
#GRDcor10 + theme_dark()

#cor.test(GRData$PBlack, GRData$Bars, use = "complete obs", method = "pearson")
#GRDcor11 <- ggplot(GRData,
#       aes(x = PBlack,
#           y = Bars,
#           color = TotPop))+
#  geom_point(alpha = 0.6, color="white",)+
#  geom_smooth(color="#8410de")+
# labs(title = "Bars
#~ Pct Black & African American Residents within Census Tracts",
#       x = "Black & African American Residents (%)",
#       y = "Outlets")
#GRDcor11 + theme_dark()

#cor.test(GRData$PBlack, GRData$Coffee, use = "complete obs", method = "pearson")
#GRDcor12 <- ggplot(GRData,
#       aes(x = PBlack,
#           y = Coffee,
#           color = TotPop))+
#  geom_point(alpha = 0.6, color="white",)+
#  geom_smooth(color="#8410de")+
# labs(title = "Coffee, Juice, and Tea Shops
#~ Pct Black & African American Residents within Census Tracts",
#       x = "Black & African American Residents (%)",
#       y = "Outlets")
#GRDcor12 + theme_dark()

#cor.test(GRData$PBlack, GRData$Restaurant, use = "complete obs", method = "pearson")
#GRDcor13 <- ggplot(GRData,
#       aes(x = PBlack,
#           y = Restaurant,
#           color = TotPop))+
#  geom_point(alpha = 0.6, color="white",)+
#  geom_smooth(color="#8410de")+
# labs(title = "Full Service Restaurants
#~ Pct Black & African American Residents within Census Tracts",
#       x = "Black & African American Residents (%)",
#       y = "Outlets")
#GRDcor13 + theme_dark()

#cor.test(GRData$PBlack, GRData$Gas, use = "complete obs", method = "pearson")
#GRDcor14 <- ggplot(GRData,
#       aes(x = PBlack,
#           y = Gas,
#           color = TotPop))+
#  geom_point(alpha = 0.6, color="white",)+
#  geom_smooth(color="#8410de")+
# labs(title = "Gas Stations
#~ Pct Black & African American Residents within Census Tracts",
#       x = "Black & African American Residents (%)",
#       y = "Outlets")
#GRDcor14 + theme_dark()

#cor.test(GRData$PBlack, GRData$Hotel, use = "complete obs", method = "pearson")
#GRDcor15 <- ggplot(GRData,
#       aes(x = PBlack,
#           y = Hotel,
#           color = TotPop))+
#  geom_point(alpha = 0.6, color="white",)+
#  geom_smooth(color="#8410de")+
# labs(title = "Hotels
#~ Pct Black & African American Residents within Census Tracts",
#       x = "Black & African American Residents (%)",
#       y = "Outlets")
#GRDcor15 + theme_dark()

#cor.test(GRData$PBlack, GRData$Liquor, use = "complete obs", method = "pearson")
#GRDcor16 <- ggplot(GRData,
#       aes(x = PBlack,
#           y = Liquor,
#           color = TotPop))+
#  geom_point(alpha = 0.6, color="white",)+
#  geom_smooth(color="#8410de")+
# labs(title = "Liquor and Party Stores
#~ Pct Black & African American Residents within Census Tracts",
#       x = "Black & African American Residents (%)",
#       y = "Outlets")
#GRDcor16 + theme_dark()

#cor.test(GRData$PBlack, GRData$Convenience, use = "complete obs", method = "pearson")
#GRDcor17 <- ggplot(GRData,
#       aes(x = PBlack,
#           y = Convenience,
#           color = TotPop))+
#  geom_point(alpha = 0.6, color="white",)+
#  geom_smooth(color="#8410de")+
# labs(title = "Convenience, Corner Stores, Small Groceries
#~ Pct Black & African American Residents within Census Tracts",
#       x = "Black & African American Residents (%)",
#       y = "Outlets")
#GRDcor17 + theme_dark()

#cor.test(GRData$PBlack, GRData$Pharmacy, use = "complete obs", method = "pearson")
#GRDcor18 <- ggplot(GRData,
#       aes(x = PBlack,
#           y = Pharmacy,
#           color = TotPop))+
#  geom_point(alpha = 0.6, color="white",)+
#  geom_smooth(color="#8410de")+
# labs(title = "Pharmacies and Drug Stores
#~ Pct Black & African American Residents within Census Tracts",
#       x = "Black & African American Residents (%)",
#       y = "Outlets")
#GRDcor18 + theme_dark()

#cor.test(GRData$PPOC, GRData$FastFood, use = "complete obs", method = "pearson")
#GRDcor20 <- ggplot(GRData,
#       aes(x = PPOC,
#           y = FastFood,
#           color = TotPop))+
#  geom_point(alpha = 0.6, color="white",)+
#  geom_smooth(color="#8410de")+
# labs(title = "Fast Food
#~ Pct People of Color within Census Tracts",
#       x = "People of Color (%)",
#       y = "Outlets")
#GRDcor20 + theme_dark()

#cor.test(GRData$PPOC, GRData$Bars, use = "complete obs", method = "pearson")
#GRDcor21 <- ggplot(GRData,
#       aes(x = PPOC,
#           y = Bars,
#           color = TotPop))+
#  geom_point(alpha = 0.6, color="white",)+
#  geom_smooth(color="#8410de")+
# labs(title = "Bars
#~ Pct People of Color within Census Tracts",
#       x = "People of Color (%)",
#       y = "Outlets")
#GRDcor21 + theme_dark()

cor.test(GRData$PPOC, GRData$Coffee, use = "complete obs", method = "pearson")

## 
##  Pearson's product-moment correlation
## 
## data:  GRData$PPOC and GRData$Coffee
## t = -2.0599, df = 100, p-value = 0.04201
## alternative hypothesis: true correlation is not equal to 0
## 95 percent confidence interval:
##  -0.381267433 -0.007574095
## sample estimates:
##        cor 
## -0.2017518

GRDcor21 <- ggplot(GRData,
       aes(x = PPOC,
           y = Coffee,
           color = TotPop))+
  geom_point(alpha = 0.6, color="white",)+
  geom_smooth(color="#8410de")+
 labs(title = "Coffee, Juice, and Tea Shops
~ Pct People of Color within Census Tracts",
       x = "People of Color (%)",
       y = "Outlets")
GRDcor21 + theme_dark()

## `geom_smooth()` using method = 'loess' and formula = 'y ~ x'

#cor.test(GRData$PPOC, GRData$Restaurant, use = "complete obs", method = "pearson")
#GRDcor22 <- ggplot(GRData,
#       aes(x = PPOC,
#           y = Restaurant,
#           color = TotPop))+
#  geom_point(alpha = 0.6, color="white",)+
#  geom_smooth(color="#8410de")+
# labs(title = "Full Service Restaurants
#~ Pct People of Color within Census Tracts",
#       x = "People of Color (%)",
#       y = "Outlets")
#GRDcor22 + theme_dark()

#cor.test(GRData$PPOC, GRData$Gas, use = "complete obs", method = "pearson")
#GRDcor23 <- ggplot(GRData,
#       aes(x = PPOC,
#           y = Gas,
#           color = TotPop))+
#  geom_point(alpha = 0.6, color="white",)+
#  geom_smooth(color="#8410de")+
# labs(title = "Gas Stations
#~ Pct People of Color within Census Tracts",
#       x = "People of Color (%)",
#       y = "Outlets")
#GRDcor23 + theme_dark()

#cor.test(GRData$PPOC, GRData$Hotel, use = "complete obs", method = "pearson")
#GRDcor24 <- ggplot(GRData,
#       aes(x = PPOC,
#           y = Hotel,
#           color = TotPop))+
#  geom_point(alpha = 0.6, color="white",)+
#  geom_smooth(color="#8410de")+
# labs(title = "Hotels
#~ Pct People of Color within Census Tracts",
#       x = "People of Color (%)",
#       y = "Outlets")
#GRDcor24 + theme_dark()

#cor.test(GRData$PPOC, GRData$Liquor, use = "complete obs", method = "pearson")
#GRDcor25 <- ggplot(GRData,
#       aes(x = PPOC,
#           y = Liquor,
#           color = TotPop))+
#  geom_point(alpha = 0.6, color="white",)+
#  geom_smooth(color="#8410de")+
# labs(title = "Liquor Stores
#~ Pct People of Color within Census Tracts",
#       x = "People of Color (%)",
#       y = "Outlets")
#GRDcor25 + theme_dark()

cor.test(GRData$PPOC, GRData$Convenience, use = "complete obs", method = "pearson")

## 
##  Pearson's product-moment correlation
## 
## data:  GRData$PPOC and GRData$Convenience
## t = 2.6951, df = 100, p-value = 0.008255
## alternative hypothesis: true correlation is not equal to 0
## 95 percent confidence interval:
##  0.06925734 0.43279922
## sample estimates:
##       cor 
## 0.2602273

GRDcor26 <- ggplot(GRData,
       aes(x = PPOC,
           y = Convenience,
           color = TotPop))+
  geom_point(alpha = 0.6, color="white",)+
  geom_smooth(color="#8410de")+
 labs(title = "Convenience, Corner Stores, Small Groceries
~ Pct People of Color within Census Tracts",
       x = "People of Color (%)",
       y = "Outlets")
GRDcor26 + theme_dark()

## `geom_smooth()` using method = 'loess' and formula = 'y ~ x'

#cor.test(GRData$PPOC, GRData$Pharmacy, use = "complete obs", method = "pearson")
#GRDcor27 <- ggplot(GRData,
#       aes(x = PPOC,
#           y = Pharmacy,
#           color = TotPop))+
#  geom_point(alpha = 0.6, color="white",)+
#  geom_smooth(color="#8410de")+
# labs(title = "Pharmacies & Drug Stores
#~ Pct People of Color within Census Tracts",
#       x = "People of Color (%)",
#       y = "Outlets")
#GRDcor27 + theme_dark()

Codependent variables

fit1 <- lm(GRData$Coffee ~ GRData$PPOC + GRData$PBachelor + GRData$Incomeadj)
summary(fit1)

## 
## Call:
## lm(formula = GRData$Coffee ~ GRData$PPOC + GRData$PBachelor + 
##     GRData$Incomeadj)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -1.8757 -0.7816 -0.2667  0.4414  5.1822 
## 
## Coefficients:
##                   Estimate Std. Error t value Pr(>|t|)    
## (Intercept)       2.852394   1.053370   2.708 0.008004 ** 
## GRData$PPOC      -0.016084   0.008939  -1.799 0.075077 .  
## GRData$PBachelor  0.061877   0.015470   4.000 0.000124 ***
## GRData$Incomeadj -0.012714   0.003979  -3.195 0.001887 ** 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 1.293 on 97 degrees of freedom
##   (1 observation deleted due to missingness)
## Multiple R-squared:  0.1939, Adjusted R-squared:  0.169 
## F-statistic: 7.779 on 3 and 97 DF,  p-value: 0.0001044

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

x1 <- GRData[c(58,66,69,70,75:78,81,83)]
#install.packages("GGally")
library(GGally)

## Registered S3 method overwritten by 'GGally':
##   method from   
##   +.gg   ggplot2

ggpairs(x1)

## Warning in ggally_statistic(data = data, mapping = mapping, na.rm = na.rm, :
## Removing 1 row that contained a missing value

## Warning in ggally_statistic(data = data, mapping = mapping, na.rm = na.rm, :
## Removing 1 row that contained a missing value

## Warning in ggally_statistic(data = data, mapping = mapping, na.rm = na.rm, :
## Removing 1 row that contained a missing value

## Warning in ggally_statistic(data = data, mapping = mapping, na.rm = na.rm, :
## Removing 1 row that contained a missing value

## Warning in ggally_statistic(data = data, mapping = mapping, na.rm = na.rm, :
## Removing 1 row that contained a missing value

## Warning in ggally_statistic(data = data, mapping = mapping, na.rm = na.rm, :
## Removing 1 row that contained a missing value

## Warning: Removed 1 rows containing missing values (`geom_point()`).
## Removed 1 rows containing missing values (`geom_point()`).
## Removed 1 rows containing missing values (`geom_point()`).
## Removed 1 rows containing missing values (`geom_point()`).
## Removed 1 rows containing missing values (`geom_point()`).
## Removed 1 rows containing missing values (`geom_point()`).

## Warning: Removed 1 rows containing non-finite values (`stat_density()`).

## Warning in ggally_statistic(data = data, mapping = mapping, na.rm = na.rm, :
## Removing 1 row that contained a missing value

## Warning in ggally_statistic(data = data, mapping = mapping, na.rm = na.rm, :
## Removing 1 row that contained a missing value

## Warning in ggally_statistic(data = data, mapping = mapping, na.rm = na.rm, :
## Removing 1 row that contained a missing value

## Warning: Removed 1 rows containing missing values (`geom_point()`).
## Removed 1 rows containing missing values (`geom_point()`).
## Removed 1 rows containing missing values (`geom_point()`).

library(ggfortify)
x2 <- GRData[c(7,8,24,36,40,58,66,69,70,75:77,81,83)]
x3 <- x2[complete.cases(x2),]
pca_df <- prcomp(x3, scale. = TRUE)
autoplot(pca_df)

autoplot(pca_df, data = x3, label = TRUE)

biplot(pca_df)
#write.csv(GRData, "C:/Users/ianmt/Documents\\GRData.csv")

Stat Testing

T-tests

#ttest1a <- t.test(Pharmacy~BlackCat, data=GRData, paired=F)
#ttest1a
#GRDRt1a <- ggplot(GRData,
#       aes(x = BlackCat,
#           y = Pharmacy,
#           fill = BlackCat))+
#  geom_col(stat = "identity")+
#  scale_fill_viridis_d(name = "Group", option = "turbo")+
#  labs(title = "Pharmacies and Drug Stores in Census Tracts with
#Pct Black & African American Residents +/- average", y = "Count", x = "Group")
#GRDRt1a+ theme_dark()

#ttest2a <- t.test(Convenience~BlackCat, data=GRData, paired=F)
#ttest2a
#GRDRt2a <- ggplot(GRData,
#       aes(x = BlackCat,
#           y = Convenience,
#           fill = BlackCat))+
#  geom_col(stat = "identity")+
#  scale_fill_viridis_d(name = "Group", option = "turbo")+
#  labs(title = "Convenience, Corner Stores, and Small Groceries in
#Tracts with Pct Black & African American Residents +/- average", y = "Count", x = "Group")
#GRDRt2a+ theme_dark()

#ttest3a <- t.test(Coffee~BlackCat, data=GRData, paired=F)
#ttest3a
#GRDRt3a <- ggplot(GRData,
#       aes(x = BlackCat,
#           y = Coffee,
#           fill = BlackCat))+
#  geom_col(stat = "identity")+
#  scale_fill_viridis_d(name = "Group", option = "turbo")+
#  labs(title = "Coffee, Tea, and Juice Shops in Tracts with 
#Pct Black & African American Residents +/- average", y = "Count", x = "Group")
#GRDRt3a+ theme_dark()

#ttest4a <- t.test(Bars~BlackCat, data=GRData, paired=F)
#ttest4a
#GRDRt4a <- ggplot(GRData,
#       aes(x = BlackCat,
#           y = Bars,
#           fill = BlackCat))+
#  geom_col(stat = "identity")+
#  scale_fill_viridis_d(name = "Group", option = "turbo")+
#  labs(title = "Bars in Tracts with Pct Black & African American Residents +/- Average", y = "Count", x = "Group")
#GRDRt4a+ theme_dark()

#ttest5a <- t.test(FastFood~BlackCat, data=GRData, paired=F)
#ttest5a
#GRDRt5a <- ggplot(GRData,
#       aes(x = BlackCat,
#           y = FastFood,
#           fill = BlackCat))+
#  geom_col(stat = "identity")+
#  scale_fill_viridis_d(name = "Group", option = "turbo")+
#  labs(title = "Fast Food in Tracts with 
#Pct Black & African American Residents +/- average", y = "Count", x = "Group")
#GRDRt5a+ theme_dark()

#ttest6a <- t.test(Restaurant~BlackCat, data=GRData, paired=F)
#ttest6a
#GRDRt6a <- ggplot(GRData,
#       aes(x = BlackCat,
#           y = Restaurants,
#           fill = BlackCat))+
#  geom_col(stat = "identity")+
#  scale_fill_viridis_d(name = "Group", option = "turbo")+
#  labs(title = "Full Service Restaurants in Tracts with 
#Pct Black & African American Residents +/- average", y = "Count", x = "Group")
#GRDRt6a+ theme_dark()

#ttest7a <- t.test(Gas~BlackCat, data=GRData, paired=F)
#ttest7a
#GRDRt7a <- ggplot(GRData,
#       aes(x = BlackCat,
#           y = Gas,
#           fill = BlackCat))+
#  geom_col(stat = "identity")+
#  scale_fill_viridis_d(name = "Group", option = "turbo")+
#  labs(title = "Gas Stations in Tracts with 
#Pct Black & African American Residents +/- average", y = "Count", x = "Group")
#GRDRt7a+ theme_dark()

#ttest8a <- t.test(Hotel~BlackCat, data=GRData, paired=F)
#ttest8a
#GRDRt8a <- ggplot(GRData,
#       aes(x = BlackCat,
#           y = Hotel,
#           fill = BlackCat))+
#  geom_col(stat = "identity")+
#  scale_fill_viridis_d(name = "Group", option = "turbo")+
#  labs(title = "Hotels in Tracts with 
#Pct Black & African American Residents +/- average", y = "Count", x = "Group")
#GRDRt8a+ theme_dark()

#ttest9a <- t.test(Liquor~BlackCat, data=GRData, paired=F)
#ttest9a
#GRDRt9a <- ggplot(GRData,
#       aes(x = BlackCat,
#           y = Liquor,
#           fill = BlackCat))+
#  geom_col(stat = "identity")+
#  scale_fill_viridis_d(name = "Group", option = "turbo")+
#  labs(title = "Liquor Stores in Tracts with 
#Pct Black & African American Residents +/- average", y = "Count", x = "Group")
#GRDRt9a+ theme_dark()

ttest1b <- t.test(Pharmacy~HispanicCat, data=GRData, paired=F)
ttest1b

## 
##  Welch Two Sample t-test
## 
## data:  Pharmacy by HispanicCat
## t = -2.7683, df = 96.825, p-value = 0.006752
## alternative hypothesis: true difference in means between group over and group under is not equal to 0
## 95 percent confidence interval:
##  -0.7804357 -0.1286552
## sample estimates:
##  mean in group over mean in group under 
##           0.2121212           0.6666667

GRDRt1b <- ggplot(GRData,
       aes(x = HispanicCat,
           y = Pharmacy,
           fill = HispanicCat))+
  geom_col(stat = "identity")+
  scale_fill_viridis_d(name = "Group", option = "turbo")+
  labs(title = "Pharmacies and Drug Stores in Census Tracts with
Pct Latinx & Hispanic Residents +/- average", y = "Count", x = "Group")

## Warning in geom_col(stat = "identity"): Ignoring unknown parameters: `stat`

GRDRt1b+ theme_dark()

ttest2b <- t.test(Convenience~HispanicCat, data=GRData, paired=F)
ttest2b

## 
##  Welch Two Sample t-test
## 
## data:  Convenience by HispanicCat
## t = 3.1153, df = 54.504, p-value = 0.002929
## alternative hypothesis: true difference in means between group over and group under is not equal to 0
## 95 percent confidence interval:
##  0.3401264 1.5676470
## sample estimates:
##  mean in group over mean in group under 
##           1.6060606           0.6521739

GRDRt2b <- ggplot(GRData,
       aes(x = HispanicCat,
           y = Convenience,
           fill = HispanicCat))+
  geom_col(stat = "identity")+
  scale_fill_viridis_d(name = "Group", option = "turbo")+
  labs(title = "Convenience, Corner Stores, and Small Groceries in
Tracts with Pct Latinx & Hispanic Residents +/- average", y = "Count", x = "Group")

## Warning in geom_col(stat = "identity"): Ignoring unknown parameters: `stat`

GRDRt2b+ theme_dark()

ttest3b <- t.test(Coffee~HispanicCat, data=GRData, paired=F)
ttest3b

## 
##  Welch Two Sample t-test
## 
## data:  Coffee by HispanicCat
## t = -2.213, df = 97.965, p-value = 0.02922
## alternative hypothesis: true difference in means between group over and group under is not equal to 0
## 95 percent confidence interval:
##  -0.95711541 -0.05210725
## sample estimates:
##  mean in group over mean in group under 
##           0.3939394           0.8985507

GRDRt3b <- ggplot(GRData,
       aes(x = HispanicCat,
           y = Coffee,
           fill = HispanicCat))+
  geom_col(stat = "identity")+
  scale_fill_viridis_d(name = "Group", option = "turbo")+
  labs(title = "Coffee, Tea, and Juice Shops in Tracts with 
Pct Latinx & Hispanic Residents +/- average", y = "Count", x = "Group")

## Warning in geom_col(stat = "identity"): Ignoring unknown parameters: `stat`

GRDRt3b+ theme_dark()

#ttest4b <- t.test(Bars~HispanicCat, data=GRData, paired=F)
#ttest4b
#GRDRt4b <- ggplot(GRData,
#       aes(x = HispanicCat,
#           y = Bars,
#           fill = HispanicCat))+
#  geom_col(stat = "identity")+
#  scale_fill_viridis_d(name = "Group", option = "turbo")+
#  labs(title = "Bars in Tracts with Pct Latinx & Hispanic Residents +/- Average", y = "Count", x = "Group")
#GRDRt4b+ theme_dark()

#ttest5b <- t.test(FastFood~HispanicCat, data=GRData, paired=F)
#ttest5b
#GRDRt5b <- ggplot(GRData,
#       aes(x = HispanicCat,
#           y = FastFood,
#           fill = HispanicCat))+
#  geom_col(stat = "identity")+
#  scale_fill_viridis_d(name = "Group", option = "turbo")+
#  labs(title = "Fast Food in Tracts with 
#Pct Latinx & Hispanic Residents +/- average", y = "Count", x = "Group")
#GRDRt5b+ theme_dark()

#ttest6b <- t.test(Restaurant~HispanicCat, data=GRData, paired=F)
#ttest6b
#GRDRt6b <- ggplot(GRData,
#       aes(x = HispanicCat,
#           y = Restaurants,
#           fill = HispanicCat))+
#  geom_col(stat = "identity")+
#  scale_fill_viridis_d(name = "Group", option = "turbo")+
#  labs(title = "Full Service Restaurants in Tracts with 
#Pct Latinx & Hispanic Residents +/- average", y = "Count", x = "Group")
#GRDRt6b+ theme_dark()

#ttest7b <- t.test(Gas~HispanicCat, data=GRData, paired=F)
#ttest7b
#GRDRt7b <- ggplot(GRData,
#       aes(x = HispanicCat,
#           y = Gas,
#           fill = HispanicCat))+
#  geom_col(stat = "identity")+
#  scale_fill_viridis_d(name = "Group", option = "turbo")+
#  labs(title = "Gas Stations in Tracts with 
#Pct Latinx & Hispanic Residents +/- average", y = "Count", x = "Group")
#GRDRt7b+ theme_dark()

#ttest8b <- t.test(Hotel~HispanicCat, data=GRData, paired=F)
#ttest8b
#GRDRt8b <- ggplot(GRData,
#       aes(x = HispanicCat,
#           y = Hotel,
#           fill = HispanicCat))+
#  geom_col(stat = "identity")+
#  scale_fill_viridis_d(name = "Group", option = "turbo")+
#  labs(title = "Hotels in Tracts with 
#Pct Latinx & Hispanic Residents +/- average", y = "Count", x = "Group")
#GRDRt8b+ theme_dark()

#ttest9b <- t.test(Liquor~HispanicCat, data=GRData, paired=F)
#ttest9b
#GRDRt9b <- ggplot(GRData,
#       aes(x = HispanicCat,
#           y = Liquor,
#           fill = HispanicCat))+
#  geom_col(stat = "identity")+
#  scale_fill_viridis_d(name = "Group", option = "turbo")+
#  labs(title = "Liquor Stores in Tracts with 
#Pct Latinx & Hispanic Residents +/- average", y = "Count", x = "Group")
#GRDRt9b+ theme_dark()

Kruskal-Wallis tests & multiple pairwise-comparisons

#View(GRData)
GRData$POCCat <- as.factor(GRData$POCCat)
levels(GRData$POCCat)

## [1] "H"  "L"  "VH" "VL"

GRData$POCCat <- ordered(GRData$POCCat)

group_by(GRData, POCCat) %>%
  summarise(
    count = n(),
    mean = mean(Bars, na.rm = TRUE),
    sd = sd(Bars, na.rm = TRUE),
    median = median(Bars, na.rm=TRUE),
    IQR = IQR(Bars, na.rm=TRUE)
  )

## # A tibble: 4 × 6
##   POCCat count  mean    sd median   IQR
##   <ord>  <int> <dbl> <dbl>  <dbl> <dbl>
## 1 H         25 0.32  0.852    0     0  
## 2 L         24 2     3.46     0.5   2.5
## 3 VH        20 0.2   0.523    0     0  
## 4 VL        33 0.273 0.626    0     0

kruskal.test(Bars ~ POCCat, data = GRData)

## 
##  Kruskal-Wallis rank sum test
## 
## data:  Bars by POCCat
## Kruskal-Wallis chi-squared = 12.631, df = 3, p-value = 0.005506

pairwise.wilcox.test(GRData$Bars, GRData$POCCat,
                 p.adjust.method = "BH")

## Warning in wilcox.test.default(xi, xj, paired = paired, ...): cannot compute
## exact p-value with ties

## Warning in wilcox.test.default(xi, xj, paired = paired, ...): cannot compute
## exact p-value with ties

## Warning in wilcox.test.default(xi, xj, paired = paired, ...): cannot compute
## exact p-value with ties

## Warning in wilcox.test.default(xi, xj, paired = paired, ...): cannot compute
## exact p-value with ties

## Warning in wilcox.test.default(xi, xj, paired = paired, ...): cannot compute
## exact p-value with ties

## Warning in wilcox.test.default(xi, xj, paired = paired, ...): cannot compute
## exact p-value with ties

## 
##  Pairwise comparisons using Wilcoxon rank sum test with continuity correction 
## 
## data:  GRData$Bars and GRData$POCCat 
## 
##    H     L     VH   
## L  0.022 -     -    
## VH 0.896 0.022 -    
## VL 0.896 0.022 0.896
## 
## P value adjustment method: BH

group_by(GRData, POCCat) %>%
  summarise(
    count = n(),
    mean = mean(Coffee, na.rm = TRUE),
    sd = sd(Coffee, na.rm = TRUE),
    median = median(Coffee, na.rm=TRUE),
    IQR = IQR(Coffee, na.rm=TRUE)
  )

## # A tibble: 4 × 6
##   POCCat count  mean    sd median   IQR
##   <ord>  <int> <dbl> <dbl>  <dbl> <dbl>
## 1 H         25 0.2   0.5        0  0   
## 2 L         24 1.5   2.09       1  2   
## 3 VH        20 0.25  0.444      0  0.25
## 4 VL        33 0.879 1.43       0  1

kruskal.test(Coffee ~ POCCat, data = GRData)

## 
##  Kruskal-Wallis rank sum test
## 
## data:  Coffee by POCCat
## Kruskal-Wallis chi-squared = 12.444, df = 3, p-value = 0.006007

pairwise.wilcox.test(GRData$Coffee, GRData$POCCat,
                 p.adjust.method = "BH")

## Warning in wilcox.test.default(xi, xj, paired = paired, ...): cannot compute
## exact p-value with ties

## Warning in wilcox.test.default(xi, xj, paired = paired, ...): cannot compute
## exact p-value with ties

## Warning in wilcox.test.default(xi, xj, paired = paired, ...): cannot compute
## exact p-value with ties

## Warning in wilcox.test.default(xi, xj, paired = paired, ...): cannot compute
## exact p-value with ties

## Warning in wilcox.test.default(xi, xj, paired = paired, ...): cannot compute
## exact p-value with ties

## Warning in wilcox.test.default(xi, xj, paired = paired, ...): cannot compute
## exact p-value with ties

## 
##  Pairwise comparisons using Wilcoxon rank sum test with continuity correction 
## 
## data:  GRData$Coffee and GRData$POCCat 
## 
##    H      L      VH    
## L  0.0096 -      -     
## VH 0.5214 0.0312 -     
## VL 0.1105 0.2340 0.2463
## 
## P value adjustment method: BH

#ns
group_by(GRData, POCCat) %>%
  summarise(
    count = n(),
    mean = mean(FastFood, na.rm = TRUE),
    sd = sd(FastFood, na.rm = TRUE),
    median = median(FastFood, na.rm=TRUE),
    IQR = IQR(FastFood, na.rm=TRUE)
  )

## # A tibble: 4 × 6
##   POCCat count  mean    sd median   IQR
##   <ord>  <int> <dbl> <dbl>  <dbl> <dbl>
## 1 H         25  2.04  3.01    0    3   
## 2 L         24  2.17  2.39    2    3   
## 3 VH        20  1.05  1.47    0.5  1.25
## 4 VL        33  1.85  2.82    0    2

kruskal.test(FastFood ~ POCCat, data = GRData)

## 
##  Kruskal-Wallis rank sum test
## 
## data:  FastFood by POCCat
## Kruskal-Wallis chi-squared = 2.4512, df = 3, p-value = 0.4842

#ns
group_by(GRData, POCCat) %>%
  summarise(
    count = n(),
    mean = mean(Restaurant, na.rm = TRUE),
    sd = sd(Restaurant, na.rm = TRUE),
    median = median(Restaurant, na.rm=TRUE),
    IQR = IQR(Restaurant, na.rm=TRUE)
  )

## # A tibble: 4 × 6
##   POCCat count  mean    sd median   IQR
##   <ord>  <int> <dbl> <dbl>  <dbl> <dbl>
## 1 H         25  2.16  2.56      1  4   
## 2 L         24  5.12  7.58      3  5.25
## 3 VH        20  1.45  1.50      1  2.25
## 4 VL        33  3     4.91      1  4

kruskal.test(Restaurant ~ POCCat, data = GRData)

## 
##  Kruskal-Wallis rank sum test
## 
## data:  Restaurant by POCCat
## Kruskal-Wallis chi-squared = 4.8809, df = 3, p-value = 0.1807

#ns
group_by(GRData, POCCat) %>%
  summarise(
    count = n(),
    mean = mean(Gas, na.rm = TRUE),
    sd = sd(Gas, na.rm = TRUE),
    median = median(Gas, na.rm=TRUE),
    IQR = IQR(Gas, na.rm=TRUE)
  )

## # A tibble: 4 × 6
##   POCCat count  mean    sd median   IQR
##   <ord>  <int> <dbl> <dbl>  <dbl> <dbl>
## 1 H         25 0.72  0.936      0  1   
## 2 L         24 0.708 0.955      0  1.25
## 3 VH        20 0.85  1.18       0  1.25
## 4 VL        33 0.879 1.19       0  2

kruskal.test(Gas ~ POCCat, data = GRData)

## 
##  Kruskal-Wallis rank sum test
## 
## data:  Gas by POCCat
## Kruskal-Wallis chi-squared = 0.15175, df = 3, p-value = 0.985

group_by(GRData, POCCat) %>%
  summarise(
    count = n(),
    mean = mean(Hotel, na.rm = TRUE),
    sd = sd(Hotel, na.rm = TRUE),
    median = median(Hotel, na.rm=TRUE),
    IQR = IQR(Hotel, na.rm=TRUE)
  )

## # A tibble: 4 × 6
##   POCCat count  mean    sd median   IQR
##   <ord>  <int> <dbl> <dbl>  <dbl> <dbl>
## 1 H         25 0.88   1.81      0     1
## 2 L         24 0.875  1.36      0     1
## 3 VH        20 0      0         0     0
## 4 VL        33 0.667  2.01      0     1

kruskal.test(Hotel ~ POCCat, data = GRData)

## 
##  Kruskal-Wallis rank sum test
## 
## data:  Hotel by POCCat
## Kruskal-Wallis chi-squared = 11.536, df = 3, p-value = 0.009154

pairwise.wilcox.test(GRData$Hotel, GRData$POCCat,
                 p.adjust.method = "BH")

## Warning in wilcox.test.default(xi, xj, paired = paired, ...): cannot compute
## exact p-value with ties

## Warning in wilcox.test.default(xi, xj, paired = paired, ...): cannot compute
## exact p-value with ties

## Warning in wilcox.test.default(xi, xj, paired = paired, ...): cannot compute
## exact p-value with ties

## Warning in wilcox.test.default(xi, xj, paired = paired, ...): cannot compute
## exact p-value with ties

## Warning in wilcox.test.default(xi, xj, paired = paired, ...): cannot compute
## exact p-value with ties

## Warning in wilcox.test.default(xi, xj, paired = paired, ...): cannot compute
## exact p-value with ties

## 
##  Pairwise comparisons using Wilcoxon rank sum test with continuity correction 
## 
## data:  GRData$Hotel and GRData$POCCat 
## 
##    H     L     VH   
## L  0.547 -     -    
## VH 0.019 0.004 -    
## VL 0.589 0.206 0.024
## 
## P value adjustment method: BH

#ns
group_by(GRData, POCCat) %>%
  summarise(
    count = n(),
    mean = mean(Liquor, na.rm = TRUE),
    sd = sd(Liquor, na.rm = TRUE),
    median = median(Liquor, na.rm=TRUE),
    IQR = IQR(Liquor, na.rm=TRUE)
  )

## # A tibble: 4 × 6
##   POCCat count  mean    sd median   IQR
##   <ord>  <int> <dbl> <dbl>  <dbl> <dbl>
## 1 H         25 0.52  0.823      0     1
## 2 L         24 0.708 0.955      0     1
## 3 VH        20 0.65  0.671      1     1
## 4 VL        33 0.545 0.754      0     1

kruskal.test(Liquor ~ POCCat, data = GRData)

## 
##  Kruskal-Wallis rank sum test
## 
## data:  Liquor by POCCat
## Kruskal-Wallis chi-squared = 1.225, df = 3, p-value = 0.747

#ns
group_by(GRData, POCCat) %>%
  summarise(
    count = n(),
    mean = mean(Pharmacy, na.rm = TRUE),
    sd = sd(Pharmacy, na.rm = TRUE),
    median = median(Pharmacy, na.rm=TRUE),
    IQR = IQR(Pharmacy, na.rm=TRUE)
  )

## # A tibble: 4 × 6
##   POCCat count  mean    sd median   IQR
##   <ord>  <int> <dbl> <dbl>  <dbl> <dbl>
## 1 H         25 0.36  0.757      0     0
## 2 L         24 0.583 1.10       0     1
## 3 VH        20 0.25  0.639      0     0
## 4 VL        33 0.758 1.09       0     1

kruskal.test(Pharmacy ~ POCCat, data = GRData)

## 
##  Kruskal-Wallis rank sum test
## 
## data:  Pharmacy by POCCat
## Kruskal-Wallis chi-squared = 5.9272, df = 3, p-value = 0.1152

group_by(GRData, POCCat) %>%
  summarise(
    count = n(),
    mean = mean(Convenience, na.rm = TRUE),
    sd = sd(Convenience, na.rm = TRUE),
    median = median(Convenience, na.rm=TRUE),
    IQR = IQR(Convenience, na.rm=TRUE)
  )

## # A tibble: 4 × 6
##   POCCat count  mean    sd median   IQR
##   <ord>  <int> <dbl> <dbl>  <dbl> <dbl>
## 1 H         25 1.12  1.54       1     2
## 2 L         24 1.25  1.15       1     2
## 3 VH        20 1.4   2.14       1     2
## 4 VL        33 0.364 0.653      0     1

kruskal.test(Convenience ~ POCCat, data = GRData)

## 
##  Kruskal-Wallis rank sum test
## 
## data:  Convenience by POCCat
## Kruskal-Wallis chi-squared = 11.71, df = 3, p-value = 0.008445

pairwise.wilcox.test(GRData$Convenience, GRData$POCCat,
                 p.adjust.method = "BH")

## Warning in wilcox.test.default(xi, xj, paired = paired, ...): cannot compute
## exact p-value with ties

## Warning in wilcox.test.default(xi, xj, paired = paired, ...): cannot compute
## exact p-value with ties

## Warning in wilcox.test.default(xi, xj, paired = paired, ...): cannot compute
## exact p-value with ties

## Warning in wilcox.test.default(xi, xj, paired = paired, ...): cannot compute
## exact p-value with ties

## Warning in wilcox.test.default(xi, xj, paired = paired, ...): cannot compute
## exact p-value with ties

## Warning in wilcox.test.default(xi, xj, paired = paired, ...): cannot compute
## exact p-value with ties

## 
##  Pairwise comparisons using Wilcoxon rank sum test with continuity correction 
## 
## data:  GRData$Convenience and GRData$POCCat 
## 
##    H      L      VH    
## L  0.5176 -      -     
## VH 0.7519 0.7269 -     
## VL 0.0649 0.0046 0.0576
## 
## P value adjustment method: BH

#ns
group_by(GRData, POCCat) %>%
  summarise(
    count = n(),
    mean = mean(Bakery, na.rm = TRUE),
    sd = sd(Bakery, na.rm = TRUE),
    median = median(Bakery, na.rm=TRUE),
    IQR = IQR(Bakery, na.rm=TRUE)
  )

## # A tibble: 4 × 6
##   POCCat count  mean    sd median   IQR
##   <ord>  <int> <dbl> <dbl>  <dbl> <dbl>
## 1 H         25 0.24  0.597      0     0
## 2 L         24 0.292 0.859      0     0
## 3 VH        20 0.15  0.366      0     0
## 4 VL        33 0.455 0.794      0     1

kruskal.test(Bakery ~ POCCat, data = GRData)

## 
##  Kruskal-Wallis rank sum test
## 
## data:  Bakery by POCCat
## Kruskal-Wallis chi-squared = 3.0155, df = 3, p-value = 0.3892

pairwise.wilcox.test(GRData$Bakery, GRData$POCCat,
                 p.adjust.method = "BH")

## Warning in wilcox.test.default(xi, xj, paired = paired, ...): cannot compute
## exact p-value with ties

## Warning in wilcox.test.default(xi, xj, paired = paired, ...): cannot compute
## exact p-value with ties

## Warning in wilcox.test.default(xi, xj, paired = paired, ...): cannot compute
## exact p-value with ties

## Warning in wilcox.test.default(xi, xj, paired = paired, ...): cannot compute
## exact p-value with ties

## Warning in wilcox.test.default(xi, xj, paired = paired, ...): cannot compute
## exact p-value with ties

## Warning in wilcox.test.default(xi, xj, paired = paired, ...): cannot compute
## exact p-value with ties

## 
##  Pairwise comparisons using Wilcoxon rank sum test with continuity correction 
## 
## data:  GRData$Bakery and GRData$POCCat 
## 
##    H    L    VH  
## L  0.99 -    -   
## VH 0.99 0.99 -   
## VL 0.48 0.48 0.48
## 
## P value adjustment method: BH

group_by(GRData, POCCat) %>%
  summarise(
    count = n(),
    mean = mean(Variety, na.rm = TRUE),
    sd = sd(Variety, na.rm = TRUE),
    median = median(Variety, na.rm=TRUE),
    IQR = IQR(Variety, na.rm=TRUE)
  )

## # A tibble: 4 × 6
##   POCCat count  mean    sd median   IQR
##   <ord>  <int> <dbl> <dbl>  <dbl> <dbl>
## 1 H         25 0.4   0.645      0     1
## 2 L         24 0.458 0.588      0     1
## 3 VH        20 0.55  0.759      0     1
## 4 VL        33 0.121 0.331      0     0

kruskal.test(Variety ~ POCCat, data = GRData)

## 
##  Kruskal-Wallis rank sum test
## 
## data:  Variety by POCCat
## Kruskal-Wallis chi-squared = 7.9749, df = 3, p-value = 0.04653

pairwise.wilcox.test(GRData$Variety, GRData$POCCat,
                 p.adjust.method = "BH")

## Warning in wilcox.test.default(xi, xj, paired = paired, ...): cannot compute
## exact p-value with ties

## Warning in wilcox.test.default(xi, xj, paired = paired, ...): cannot compute
## exact p-value with ties

## Warning in wilcox.test.default(xi, xj, paired = paired, ...): cannot compute
## exact p-value with ties

## Warning in wilcox.test.default(xi, xj, paired = paired, ...): cannot compute
## exact p-value with ties

## Warning in wilcox.test.default(xi, xj, paired = paired, ...): cannot compute
## exact p-value with ties

## Warning in wilcox.test.default(xi, xj, paired = paired, ...): cannot compute
## exact p-value with ties

## 
##  Pairwise comparisons using Wilcoxon rank sum test with continuity correction 
## 
## data:  GRData$Variety and GRData$POCCat 
## 
##    H     L     VH   
## L  0.714 -     -    
## VH 0.714 0.861 -    
## VL 0.113 0.043 0.043
## 
## P value adjustment method: BH

#ns
group_by(GRData, POCCat) %>%
  summarise(
    count = n(),
    mean = mean(IceCream, na.rm = TRUE),
    sd = sd(IceCream, na.rm = TRUE),
    median = median(IceCream, na.rm=TRUE),
    IQR = IQR(IceCream, na.rm=TRUE)
  )

## # A tibble: 4 × 6
##   POCCat count  mean    sd median   IQR
##   <ord>  <int> <dbl> <dbl>  <dbl> <dbl>
## 1 H         25 0.04  0.2        0     0
## 2 L         24 0.208 0.509      0     0
## 3 VH        20 0.1   0.308      0     0
## 4 VL        33 0.152 0.442      0     0

kruskal.test(IceCream ~ POCCat, data = GRData)

## 
##  Kruskal-Wallis rank sum test
## 
## data:  IceCream by POCCat
## Kruskal-Wallis chi-squared = 2.187, df = 3, p-value = 0.5345

pairwise.wilcox.test(GRData$IceCream, GRData$POCCat,
                 p.adjust.method = "BH")

## Warning in wilcox.test.default(xi, xj, paired = paired, ...): cannot compute
## exact p-value with ties

## Warning in wilcox.test.default(xi, xj, paired = paired, ...): cannot compute
## exact p-value with ties

## Warning in wilcox.test.default(xi, xj, paired = paired, ...): cannot compute
## exact p-value with ties

## Warning in wilcox.test.default(xi, xj, paired = paired, ...): cannot compute
## exact p-value with ties

## Warning in wilcox.test.default(xi, xj, paired = paired, ...): cannot compute
## exact p-value with ties

## Warning in wilcox.test.default(xi, xj, paired = paired, ...): cannot compute
## exact p-value with ties

## 
##  Pairwise comparisons using Wilcoxon rank sum test with continuity correction 
## 
## data:  GRData$IceCream and GRData$POCCat 
## 
##    H    L    VH  
## L  0.77 -    -   
## VH 0.77 0.77 -   
## VL 0.77 0.77 0.80
## 
## P value adjustment method: BH

#ns
group_by(GRData, POCCat) %>%
  summarise(
    count = n(),
    mean = mean(Supermarket, na.rm = TRUE),
    sd = sd(Supermarket, na.rm = TRUE),
    median = median(Supermarket, na.rm=TRUE),
    IQR = IQR(Supermarket, na.rm=TRUE)
  )

## # A tibble: 4 × 6
##   POCCat count   mean    sd median   IQR
##   <ord>  <int>  <dbl> <dbl>  <dbl> <dbl>
## 1 H         25 0.04   0.2        0  0   
## 2 L         24 0.0833 0.282      0  0   
## 3 VH        20 0.3    0.571      0  0.25
## 4 VL        33 0.273  0.517      0  0

kruskal.test(Supermarket ~ POCCat, data = GRData)

## 
##  Kruskal-Wallis rank sum test
## 
## data:  Supermarket by POCCat
## Kruskal-Wallis chi-squared = 6.7751, df = 3, p-value = 0.07942

pairwise.wilcox.test(GRData$Supermarket, GRData$POCCat,
                 p.adjust.method = "BH")

## Warning in wilcox.test.default(xi, xj, paired = paired, ...): cannot compute
## exact p-value with ties

## Warning in wilcox.test.default(xi, xj, paired = paired, ...): cannot compute
## exact p-value with ties

## Warning in wilcox.test.default(xi, xj, paired = paired, ...): cannot compute
## exact p-value with ties

## Warning in wilcox.test.default(xi, xj, paired = paired, ...): cannot compute
## exact p-value with ties

## Warning in wilcox.test.default(xi, xj, paired = paired, ...): cannot compute
## exact p-value with ties

## Warning in wilcox.test.default(xi, xj, paired = paired, ...): cannot compute
## exact p-value with ties

## 
##  Pairwise comparisons using Wilcoxon rank sum test with continuity correction 
## 
## data:  GRData$Supermarket and GRData$POCCat 
## 
##    H    L    VH  
## L  0.66 -    -   
## VH 0.13 0.20 -   
## VL 0.13 0.20 0.93
## 
## P value adjustment method: BH

#ns
group_by(GRData, POCCat) %>%
  summarise(
    count = n(),
    mean = mean(FarmMarkets, na.rm = TRUE),
    sd = sd(FarmMarkets, na.rm = TRUE),
    median = median(FarmMarkets, na.rm=TRUE),
    IQR = IQR(FarmMarkets, na.rm=TRUE)
  )

## # A tibble: 4 × 6
##   POCCat count  mean    sd median   IQR
##   <ord>  <int> <dbl> <dbl>  <dbl> <dbl>
## 1 H         25 0.12  0.332      0     0
## 2 L         24 0.208 0.658      0     0
## 3 VH        20 0.1   0.308      0     0
## 4 VL        33 0.273 0.876      0     0

kruskal.test(FarmMarkets ~ POCCat, data = GRData)

## 
##  Kruskal-Wallis rank sum test
## 
## data:  FarmMarkets by POCCat
## Kruskal-Wallis chi-squared = 0.1135, df = 3, p-value = 0.9902

pairwise.wilcox.test(GRData$FarmMarkets, GRData$POCCat,
                 p.adjust.method = "BH")

## Warning in wilcox.test.default(xi, xj, paired = paired, ...): cannot compute
## exact p-value with ties

## Warning in wilcox.test.default(xi, xj, paired = paired, ...): cannot compute
## exact p-value with ties

## Warning in wilcox.test.default(xi, xj, paired = paired, ...): cannot compute
## exact p-value with ties

## Warning in wilcox.test.default(xi, xj, paired = paired, ...): cannot compute
## exact p-value with ties

## Warning in wilcox.test.default(xi, xj, paired = paired, ...): cannot compute
## exact p-value with ties

## Warning in wilcox.test.default(xi, xj, paired = paired, ...): cannot compute
## exact p-value with ties

## 
##  Pairwise comparisons using Wilcoxon rank sum test with continuity correction 
## 
## data:  GRData$FarmMarkets and GRData$POCCat 
## 
##    H L VH
## L  1 - - 
## VH 1 1 - 
## VL 1 1 1 
## 
## P value adjustment method: BH

#ns
group_by(GRData, POCCat) %>%
  summarise(
    count = n(),
    mean = mean(School, na.rm = TRUE),
    sd = sd(School, na.rm = TRUE),
    median = median(School, na.rm=TRUE),
    IQR = IQR(School, na.rm=TRUE)
  )

## # A tibble: 4 × 6
##   POCCat count  mean    sd median   IQR
##   <ord>  <int> <dbl> <dbl>  <dbl> <dbl>
## 1 H         25 0.04  0.2        0     0
## 2 L         24 0.292 0.690      0     0
## 3 VH        20 0.05  0.224      0     0
## 4 VL        33 0.182 0.465      0     0

kruskal.test(School ~ POCCat, data = GRData)

## 
##  Kruskal-Wallis rank sum test
## 
## data:  School by POCCat
## Kruskal-Wallis chi-squared = 4.6566, df = 3, p-value = 0.1987

pairwise.wilcox.test(GRData$School, GRData$POCCat,
                 p.adjust.method = "BH")

## Warning in wilcox.test.default(xi, xj, paired = paired, ...): cannot compute
## exact p-value with ties

## Warning in wilcox.test.default(xi, xj, paired = paired, ...): cannot compute
## exact p-value with ties

## Warning in wilcox.test.default(xi, xj, paired = paired, ...): cannot compute
## exact p-value with ties

## Warning in wilcox.test.default(xi, xj, paired = paired, ...): cannot compute
## exact p-value with ties

## Warning in wilcox.test.default(xi, xj, paired = paired, ...): cannot compute
## exact p-value with ties

## Warning in wilcox.test.default(xi, xj, paired = paired, ...): cannot compute
## exact p-value with ties

## 
##  Pairwise comparisons using Wilcoxon rank sum test with continuity correction 
## 
## data:  GRData$School and GRData$POCCat 
## 
##    H    L    VH  
## L  0.34 -    -   
## VH 0.90 0.34 -   
## VL 0.34 0.70 0.39
## 
## P value adjustment method: BH

#ns
group_by(GRData, POCCat) %>%
  summarise(
    count = n(),
    mean = mean(Retirement, na.rm = TRUE),
    sd = sd(Retirement, na.rm = TRUE),
    median = median(Retirement, na.rm=TRUE),
    IQR = IQR(Retirement, na.rm=TRUE)
  )

## # A tibble: 4 × 6
##   POCCat count   mean    sd median   IQR
##   <ord>  <int>  <dbl> <dbl>  <dbl> <dbl>
## 1 H         25 0.48   0.714      0  1   
## 2 L         24 0.0833 0.282      0  0   
## 3 VH        20 0.35   0.671      0  0.25
## 4 VL        33 0.515  0.906      0  1

kruskal.test(Retirement ~ POCCat, data = GRData)

## 
##  Kruskal-Wallis rank sum test
## 
## data:  Retirement by POCCat
## Kruskal-Wallis chi-squared = 6.2434, df = 3, p-value = 0.1004

pairwise.wilcox.test(GRData$Retirement, GRData$POCCat,
                 p.adjust.method = "BH")

## Warning in wilcox.test.default(xi, xj, paired = paired, ...): cannot compute
## exact p-value with ties

## Warning in wilcox.test.default(xi, xj, paired = paired, ...): cannot compute
## exact p-value with ties

## Warning in wilcox.test.default(xi, xj, paired = paired, ...): cannot compute
## exact p-value with ties

## Warning in wilcox.test.default(xi, xj, paired = paired, ...): cannot compute
## exact p-value with ties

## Warning in wilcox.test.default(xi, xj, paired = paired, ...): cannot compute
## exact p-value with ties

## Warning in wilcox.test.default(xi, xj, paired = paired, ...): cannot compute
## exact p-value with ties

## 
##  Pairwise comparisons using Wilcoxon rank sum test with continuity correction 
## 
## data:  GRData$Retirement and GRData$POCCat 
## 
##    H     L     VH   
## L  0.071 -     -    
## VH 0.648 0.248 -    
## VL 0.889 0.071 0.648
## 
## P value adjustment method: BH

group_by(GRData, POCCat) %>%
  summarise(
    count = n(),
    mean = mean(ChildCare, na.rm = TRUE),
    sd = sd(ChildCare, na.rm = TRUE),
    median = median(ChildCare, na.rm=TRUE),
    IQR = IQR(ChildCare, na.rm=TRUE)
  )

## # A tibble: 4 × 6
##   POCCat count  mean    sd median   IQR
##   <ord>  <int> <dbl> <dbl>  <dbl> <dbl>
## 1 H         25 0.16  0.374      0     0
## 2 L         24 0.542 0.721      0     1
## 3 VH        20 0.15  0.366      0     0
## 4 VL        33 0.545 0.905      0     1

kruskal.test(ChildCare ~ POCCat, data = GRData)

## 
##  Kruskal-Wallis rank sum test
## 
## data:  ChildCare by POCCat
## Kruskal-Wallis chi-squared = 7.8403, df = 3, p-value = 0.04943

pairwise.wilcox.test(GRData$ChildCare, GRData$POCCat,
                 p.adjust.method = "BH")

## Warning in wilcox.test.default(xi, xj, paired = paired, ...): cannot compute
## exact p-value with ties

## Warning in wilcox.test.default(xi, xj, paired = paired, ...): cannot compute
## exact p-value with ties

## Warning in wilcox.test.default(xi, xj, paired = paired, ...): cannot compute
## exact p-value with ties

## Warning in wilcox.test.default(xi, xj, paired = paired, ...): cannot compute
## exact p-value with ties

## Warning in wilcox.test.default(xi, xj, paired = paired, ...): cannot compute
## exact p-value with ties

## Warning in wilcox.test.default(xi, xj, paired = paired, ...): cannot compute
## exact p-value with ties

## 
##  Pairwise comparisons using Wilcoxon rank sum test with continuity correction 
## 
## data:  GRData$ChildCare and GRData$POCCat 
## 
##    H     L     VH   
## L  0.088 -     -    
## VH 0.942 0.088 -    
## VL 0.163 0.713 0.163
## 
## P value adjustment method: BH

#ns
group_by(GRData, POCCat) %>%
  summarise(
    count = n(),
    mean = mean(Wholesale, na.rm = TRUE),
    sd = sd(Wholesale, na.rm = TRUE),
    median = median(Wholesale, na.rm=TRUE),
    IQR = IQR(Wholesale, na.rm=TRUE)
  )

## # A tibble: 4 × 6
##   POCCat count  mean    sd median   IQR
##   <ord>  <int> <dbl> <dbl>  <dbl> <dbl>
## 1 H         25 0.12  0.6        0  0   
## 2 L         24 0.333 0.637      0  0.25
## 3 VH        20 0.3   0.801      0  0   
## 4 VL        33 0.333 0.990      0  0

kruskal.test(Wholesale ~ POCCat, data = GRData)

## 
##  Kruskal-Wallis rank sum test
## 
## data:  Wholesale by POCCat
## Kruskal-Wallis chi-squared = 3.7165, df = 3, p-value = 0.2938

pairwise.wilcox.test(GRData$Wholesale, GRData$POCCat,
                 p.adjust.method = "BH")

## Warning in wilcox.test.default(xi, xj, paired = paired, ...): cannot compute
## exact p-value with ties

## Warning in wilcox.test.default(xi, xj, paired = paired, ...): cannot compute
## exact p-value with ties

## Warning in wilcox.test.default(xi, xj, paired = paired, ...): cannot compute
## exact p-value with ties

## Warning in wilcox.test.default(xi, xj, paired = paired, ...): cannot compute
## exact p-value with ties

## Warning in wilcox.test.default(xi, xj, paired = paired, ...): cannot compute
## exact p-value with ties

## Warning in wilcox.test.default(xi, xj, paired = paired, ...): cannot compute
## exact p-value with ties

## 
##  Pairwise comparisons using Wilcoxon rank sum test with continuity correction 
## 
## data:  GRData$Wholesale and GRData$POCCat 
## 
##    H    L    VH  
## L  0.30 -    -   
## VH 0.46 0.62 -   
## VL 0.46 0.62 1.00
## 
## P value adjustment method: BH

ANOVA w/ categorized scatter plots

anova1b <- aov(Specialty~POCCat, data=GRDataRe)
summary(anova1b)

##             Df Sum Sq Mean Sq F value Pr(>F)  
## POCCat       3   52.6  17.537   3.604 0.0161 *
## Residuals   98  476.9   4.866                 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

a1b <- ggplot(GRDataRe,
       aes(x = PPOC,
           y = Specialty,
           color = POCCat))+
  geom_jitter()+
  scale_color_manual(values = c("yellow", "cyan", "orangered", "darkblue"))+
  facet_wrap(~POCCat)+
  labs(title = "Specialty Outlets ~ People of Color (%)",
       x = "Percent",
       y = "Outlets")
a1b + theme_dark()

#anova2b <- aov(Restaurant~POCCat, data=GRDataRe)
#summary(anova2b)
#a2b <- ggplot(GRDataRe,
#       aes(x = PPOC,
#           y = Restaurant,
#           color = POCCat))+
#  geom_jitter()+
#  scale_color_manual(values = c("yellow", "cyan", "orangered", "darkblue"))+
#  facet_wrap(~POCCat)+
#  labs(title = "Full Service, Fast Food, and Take Out Restaurants ~ People of Color (%)",
#       x = "Percent",
#       y = "Outlets")
#a2b + theme_dark()

#anova3b <- aov(CLG~POCCat, data=GRDataRe)
#summary(anova3b)
#a3b <- ggplot(GRDataRe,
#       aes(x = PPOC,
#           y = CLG,
#           color = POCCat))+
#  geom_jitter()+
#  scale_color_manual(values = c("yellow", "cyan", "orangered", "darkblue"))+
#  facet_wrap(~POCCat)+
#  labs(title = "Convenience, Corner Stores, Small Groceries, Liquor, and Gas Stations ~ People of Color (%)",
#       x = "Percent",
#       y = "Outlets")
#a3b + theme_dark()

anova4b <- aov(Convenience~POCCat, data=GRData)
summary(anova4b)

##             Df Sum Sq Mean Sq F value Pr(>F)  
## POCCat       3  18.27   6.089   3.181 0.0273 *
## Residuals   98 187.58   1.914                 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

a4b <- ggplot(GRData,
       aes(x = PPOC,
           y = Convenience,
           color = POCCat))+
  geom_jitter()+
  scale_color_manual(values = c("yellow", "cyan", "orangered", "darkblue"))+
  facet_wrap(~POCCat)+
  labs(title = "Convenience, Corner Stores, and Small groceries
~ People of Color (%)",
       x = "Percent",
       y = "Outlets")
a4b + theme_dark()

anova5b <- aov(Restaurant~POCCat, data=GRData)
summary(anova5b)

##             Df Sum Sq Mean Sq F value Pr(>F)  
## POCCat       3  174.1   58.02   2.478 0.0658 .
## Residuals   98 2294.9   23.42                 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

a5b <- ggplot(GRData,
       aes(x = PPOC,
           y = Restaurant,
           color = POCCat))+
  geom_jitter()+
  scale_color_manual(values = c("yellow", "cyan", "orangered", "darkblue"))+
  facet_wrap(~POCCat)+
  labs(title = "Restaurants ~ People of Color (%)",
       x = "Percent",
       y = "Outlets")
a5b + theme_dark()

anova6b <- aov(Coffee~POCCat, data=GRData)
summary(anova6b)

##             Df Sum Sq Mean Sq F value  Pr(>F)   
## POCCat       3  26.59   8.863   4.956 0.00303 **
## Residuals   98 175.27   1.788                   
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

a6b <- ggplot(GRData,
       aes(x = PPOC,
           y = Coffee,
           color = POCCat))+
  geom_jitter()+
  scale_color_manual(values = c("yellow", "cyan", "orangered", "darkblue"))+
  facet_wrap(~POCCat)+
  labs(title = "Coffee, Tea, and Juice Shops ~ People of Color (%)",
       x = "Percent",
       y = "Outlets")
a6b + theme_dark()

anova7b <- aov(Bars~POCCat, data=GRData)
summary(anova7b)

##             Df Sum Sq Mean Sq F value Pr(>F)   
## POCCat       3  55.14  18.379   5.788 0.0011 **
## Residuals   98 311.19   3.175                  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

a7b <- ggplot(GRData,
       aes(x = PPOC,
           y = Bars,
           color = POCCat))+
  geom_jitter()+
  scale_color_manual(values = c("yellow", "cyan", "orangered", "darkblue"))+
  facet_wrap(~POCCat)+
  labs(title = "Bars ~ People of Color (%)",
       x = "Percent",
       y = "Outlets")
a7b + theme_dark()

#anova8b <- aov(FastFood~POCCat, data=GRData)
#summary(anova8b)
#a8b <- ggplot(GRData,
#       aes(x = PPOC,
#           y = FastFood,
#           color = POCCat))+
#  geom_jitter()+
#  scale_color_manual(values = c("yellow", "cyan", "orangered", "darkblue"))+
#  facet_wrap(~POCCat)+
#  labs(title = "Fast Food ~ People of Color (%)",
#       x = "Percent",
#       y = "Outlets")
#a8b + theme_dark()

#anova9b <- aov(Gas~POCCat, data=GRData)
#summary(anova9b)
#a9b <- ggplot(GRData,
#       aes(x = PPOC,
#           y = Gas,
#           color = POCCat))+
#  geom_jitter()+
#  scale_color_manual(values = c("yellow", "cyan", "orangered", "darkblue"))+
#  facet_wrap(~POCCat)+
#  labs(title = "Gas Stations ~ People of Color (%)",
#       x = "Percent",
#       y = "Outlets")
#a9b + theme_dark()

anova10b <- aov(Hotel~POCCat, data=GRData)
summary(anova10b)

##             Df Sum Sq Mean Sq F value Pr(>F)
## POCCat       3  10.98   3.660   1.431  0.238
## Residuals   98 250.60   2.557

a10b <- ggplot(GRData,
       aes(x = PPOC,
           y = Hotel,
           color = POCCat))+
  geom_jitter()+
  scale_color_manual(values = c("yellow", "cyan", "orangered", "darkblue"))+
  facet_wrap(~POCCat)+
  labs(title = "Hotels ~ People of Color (%)",
       x = "Percent",
       y = "Outlets")
a10b + theme_dark()

#anova11b <- aov(Pharmacy~POCCat, data=GRData)
#summary(anova11b)
#a11b <- ggplot(GRData,
#       aes(x = PPOC,
#           y = Pharmacy,
#           color = POCCat))+
#  geom_jitter()+
#  scale_color_manual(values = c("yellow", "cyan", "orangered", "darkblue"))+
#  facet_wrap(~POCCat)+
#  labs(title = "Pharmacies & Drug Stores ~ People of Color (%)",
#       x = "Percent",
#       y = "Outlets")
#a11b + theme_dark()

anova12b <- aov(ChildCare~POCCat, data=GRData)
summary(anova12b)

##             Df Sum Sq Mean Sq F value Pr(>F)  
## POCCat       3   3.79  1.2643   2.813 0.0433 *
## Residuals   98  44.05  0.4495                 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

a12b <- ggplot(GRData,
       aes(x = PPOC,
           y = ChildCare,
           color = POCCat))+
  geom_jitter()+
  scale_color_manual(values = c("yellow", "cyan", "orangered", "darkblue"))+
  facet_wrap(~POCCat)+
  labs(title = "Child Care w Food ~ People of Color (%)",
       x = "Percent",
       y = "Outlets")
a12b + theme_dark()

anova13b <- aov(Variety~POCCat, data=GRData)
summary(anova13b)

##             Df Sum Sq Mean Sq F value Pr(>F)  
## POCCat       3   2.87  0.9569   2.892 0.0392 *
## Residuals   98  32.42  0.3309                 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

a13b <- ggplot(GRData,
       aes(x = PPOC,
           y = Variety,
           color = POCCat))+
  geom_jitter()+
  scale_color_manual(values = c("yellow", "cyan", "orangered", "darkblue"))+
  facet_wrap(~POCCat)+
  labs(title = "Dollar & Variety Stores ~ People of Color (%)",
       x = "Percent",
       y = "Outlets")
a13b + theme_dark()