## Zip codes!


Zip codes are another data playground. So many data sets are indexed (or include) zip code. For example, the [American Community Survey](https://en.wikipedia.org/wiki/American_Community_Survey) is easily accessible in R with the library [acs](https://cran.r-project.org/web/packages/acs/index.html). Locations of zip codes are available in the R library zipcode. 

```{r cache=TRUE}
set.seed(1)
source("http://pages.stat.wisc.edu/~karlrohe/netsci/code/loadData.R")
wi = DT[State == "WI"]
```



```{r}
library(data.table)  # so fast!
library(igraph)  # all the basic graph operations.
library(zipcode)
zip = wi$"Zip Code"
zip = substr(zip, start = 1, stop = 5)

data(zipcode)  # this contains the locations of zip codes
zipcode = as.data.table(zipcode); setkey(zipcode, zip)  # thanks data.table for making things so fast!  
loc =  zipcode[zip, c("latitude", "longitude"), with = F]
loc = loc[complete.cases(loc)]
loc = as.matrix(loc)
plot(loc)
plot(loc[,2], loc[,1])
```

We can make it prettier by finding an R library to help.

```{r}
library(maps); library(ggplot2)
plot(loc[,2], loc[,1])
map('state', region = c('wisconsin'), add = T)  # adds an outline

# or, you might like this version:
# install.packages("ggmap")
library(ggmap)
map <- get_map(location = 'Madison, Wisconsin', zoom=10)
ggmap(map)  +   geom_point(aes(x = loc[,2], y = loc[,1], alpha = .5))

```


Let's compute the geographic distance of some of the edges.  There are at least three different ways to do this.  The third way is with a "pipe" operator %>%.   While I will not use %>% here, this is the basic idea:  for each edge, get the NPI  %>% "Zip Code" (using DT) %>%   longitude, latitude (in dataset zipcode) %>% dist in distGeo.  The pipe operator is nice because the code becomes much more readable.  Below, you will see that everything must be "backwards".  First, take a simple random sample of rows of DT.  We will look at the edges coming out of these nodes.


```{r cache = T}
samp =   DT$NPI[sample(dim(DT)[1], 10000)]  # take a random sample of NPI's. 
DTsamp = DT[samp,mult ="first"]
dim(DTsamp)
DTsamp = DTsamp[complete.cases(DTsamp$"Zip Code")]
dim(DTsamp)
setkey(DTsamp, NPI)
tmp = Et[DTsamp$NPI]
Esamp = tmp[complete.cases(tmp)]  #lots of NA's.  Have not inspected why.
Esamp=as.matrix(Esamp)[,1:2] #igraph needs the edgelist to be in matrix format
```

Ok, first, the slow way using apply.

```{r cache=T}
library(geosphere)
edgeDistance  = function(x){
  from = x[1]; to = x[2]
  p1 = zipcode[substr(DTsamp[from]$"Zip Code", start = 1, stop =5), c("longitude", "latitude"), with = F] 
  p2 = zipcode[substr(DT[to, mult = "first"]$"Zip Code", start = 1, stop =5), c("longitude","latitude"), with = F] 
  if(mean(complete.cases(p2))==1 ) return(distGeo(p1,p2)/1000)
  return(NA)
}
# becuase this way is so slow, let's just do 100 edges.
edtmp = apply(Esamp[1:100,], 1, edgeDistance)
```

```{r}
# so ugly and so fast!
ed = distGeo(
  zipcode[
    substr(
      DT[Esamp[,1], mult = "first"]$"Zip Code" ,start = 1, stop = 5
    )
    , c("longitude", "latitude"), with = F] 
  , zipcode[substr(DT[Esamp[,2], mult = "first"]$"Zip Code" ,start = 1, stop = 5), c("longitude", "latitude"), with = F] 
)/1000
mean(ed ==0, na.rm = T)
hist(log(ed+1,10))

```

How do the distribution of referral distances vary between providers?  Let's study it using some characteristics of the provider's zip code.  Do providers in rural areas make referrals over longer distances?  You could use acs data for this.  For right now, I found some data on [population density and unemployment rate](http://blog.splitwise.com/2014/01/06/free-us-population-density-and-unemployment-rate-by-zip-code/) with a google search.  If this was for a publication, I would definitely want to find this data at a more citable place (e.g. ACS and BLS).  

```{r}
edgeZip = substr(DT[Esamp[,1], mult = "first"]$"Zip Code",1,5)
names(which.max(table(Esamp[,1])))
popd = fread("~/dataFiles/physicianReferral/popDensity.csv",sep = ",")
unemp = fread("~/dataFiles/physicianReferral/Unemployment.csv",sep = ",")
setkey(popd, Zip/ZCTA)
setkey(unemp, Zip)
str(unemp)  # careful!  zeros removed.

#unemployment has a %. remove it.
unemp$`Unemp. Rate`[nchar(unemp$`Unemp. Rate`)==2] = substr(unemp$`Unemp. Rate`[nchar(unemp$`Unemp. Rate`)==2],1,1)
unemp$`Unemp. Rate`[nchar(unemp$`Unemp. Rate`)==3] = substr(unemp$`Unemp. Rate`[nchar(unemp$`Unemp. Rate`)==3],1,2)
unemp$`Unemp. Rate`[nchar(unemp$`Unemp. Rate`)==4] = substr(unemp$`Unemp. Rate`[nchar(unemp$`Unemp. Rate`)==4],1,3)
unemp$`Unemp. Rate` = as.numeric(unemp$`Unemp. Rate`)
edgeUnemp = unemp$`Unemp. Rate`[as.numeric(edgeZip)]

edgeUnemp[edgeUnemp ==100] = NA
plot( edgeUnemp, log(ed+1,10))
abline( lm(log(ed+1,10)~edgeUnemp))
summary( lm(log(ed+1,10)~edgeUnemp))


edgeDen = popd[as.numeric(edgeZip)]$`Density Per Sq Mile`
plot( log(edgeDen+1,10), log(ed+1,10))
abline( lm(log(ed+1,10)~log(edgeDen+1,10)))
summary( lm(log(ed+1,10)~log(edgeDen+1,10)))

summary( lm(log(ed+1,10)~log(edgeDen+1,10) + log(edgeUnemp+1,10)))
```

There are so many problems with this analysis!  Let's count the ways.


