rules.Rmd

---
title: "rules"
author: "Ilya"
date: "5/13/2018"
output: github_document
---

#####install packages
```{r}
#install package dplyr if not already installed
list.of.packages <- c("dplyr", "rgdal", "raster", "sp", "data.table", "reshape", "auk", "lubridate", "apriori", "arulesViz")
new.packages <- list.of.packages[!(list.of.packages %in% installed.packages()[,"Package"])]
if(length(new.packages)) install.packages(new.packages, repos = "http://cran.us.r-project.org")
print(new.packages)
library(dplyr)
library(auk)
library(lubridate)
library(reshape)
library(arules)
library(data.table)
library(arulesViz)
```

#####read in shapefile of protected areas and reproject and save as PAD_wgs84.Rdata -- commenting out this part
```{r}

#   library("rgdal")
#  library("raster") 
# library("sp")
# 
# #source: https://gapanalysis.usgs.gov/padus/data/download/
# PAD <- shapefile("PADUS1_4Combined.shp")
# #change datum
# PAD_wgs84 <- spTransform(PAD, CRS("+proj=longlat +datum=WGS84"))
# 
# #assign cooordinates (package: sp)
# PAD_coordinates = coordinates(PAD_wgs84)
# print(proc.time() - ptm)
# 
# PAD_wgs84$X = PAD_coordinates[,1]
# PAD_wgs84$Y = PAD_coordinates[,2]
# 
# save(PAD_wgs84, file = "PAD_wgs84.Rdata")

```

#####read in iNaturalist data and assign it to parks -- commenting out this part
```{r}

# rm(list = ls())
# library("raster")
# library("data.table")
# load("PAD_wgs84.Rdata")
# 
# #occurrence download available here: 
# #https://www.gbif.org/occurrence/download/0016851-180131172636756
# inat = fread("occurrence.txt",blank.lines.skip=TRUE)
# inat = as.data.frame(inat)
# inat$longitude <- inat$decimalLongitude
# inat$latitude <- inat$decimalLatitude
# inat.ok= filter(inat, !is.na(decimalLatitude) & !is.na(decimalLongitude))
# 
# #turn it into a SpatialPointsDataFrame
# coordinates(inat.ok) <- c("decimalLongitude","decimalLatitude")#package sp
# 
# save(inat.ok,file = "inat.ok.Rdata")
# 
#  # tell R that inat coordinates are in the same lat/lon reference system
# # as the parks data 
# projection(inat.ok) <- projection(PAD_wgs84)
# 
# # use 'over' with parks as a SpatialPolygonsDataFrame
# # object, to determine which park (if any) contains each sighting, and
# # store the park name as an attribute of the data
# inat.ok$park <- over(inat.ok, PAD_wgs84)$Unit_Nm#name of park
# inat.ok$X.park <- over(inat.ok, PAD_wgs84)$X
# inat.ok$Y.park <- over(inat.ok, PAD_wgs84)$Y
# inat.ok$park.acres <- over(inat.ok, PAD_wgs84)$GIS_Acres
# 
# save(inat.ok, file = "inat.ok.Rdata")

```

#####make small version of data -- comment out this part
```{r}
# load("inat.ok.Rdata")
# inat.ok$date =  as.Date(paste(inat.ok$year, inat.ok$month, inat.ok$day, sep='-'))
# inat.ok.sm = inat.ok[, c("park", "recordedBy", "species", "date", "X.park", "Y.park")]
# save(inat.ok.sm, file = "inat.ok.sm.Rdata")

```

#####subset data for NYC 
```{r}
load("inat.ok.sm.Rdata")
inat.nyc = subset(inat.ok.sm,X.park >=-74.2589 & X.park< -73.7004 & Y.park <=40.9176 & Y.park>=40.4774)
write.csv(inat.nyc, file = "inat.nyc.csv")
```

#####find association rules for rats in NYC; plot strongly supported associations
```{r}
inat.nyc = read.csv("inat.nyc.csv")
inat.ok = inat.nyc
inat.df = as.data.frame(inat.ok) 

#get count by park, date, and species for all data
inat.ct <- inat.df %>%
    group_by(park, date, species) %>%
  summarize(count = n())
inat.ct$count = 1#make all counts 1
 inat.ct = data.frame(inat.ct)
 inat.ct$park_date_recordedBy = factor(paste(inat.ct$park, inat.ct$date, inat.ct$recordedBy, sep = "_"))
  inat.ct$park_date = factor(paste(inat.ct$park, inat.ct$date, sep = "_"))
 i.df = inat.ct[, c("park_date_recordedBy", "park_date", "species", "count")]
 
 RaNo = subset(i.df, species == "Rattus norvegicus")
 dim(RaNo)#only seen 5 times
#write and then read back in again as transaction
  write.csv(i.df, file = "i.df.csv")
 
itrans = read.transactions("i.df.csv", format = "single", sep = ",", cols = c("park_date_recordedBy", "species"))

#find all rules involving Norway rats
irulesRaNo <- apriori(itrans, parameter = list(support = 0.0001, confidence = 0.2, maxlen=3), appearance = list(rhs="Rattus norvegicus"))# 
subrulesRaNo <- head(irulesRaNo, n = 5, by = "lift")
plot(subrulesRaNo, method = "graph", main = "Norway rat, NYC, N=5 rules")
inspect(subrulesRaNo)

#find all rules for black-capped chickadees:
irulesBlCa <- apriori(itrans, parameter = list(support = 0.0001, confidence = 0.2, maxlen=3), appearance = list(rhs="Poecile atricapillus"))#
subrulesBlCa <- head(irulesBlCa, n = 5, by = "lift")
plot(subrulesBlCa, method = "graph", main = "Black-capped Chickadee, NYC, N=5 rules")
inspect(subrulesBlCa)
    ```

#####find association rules across all iNaturalist data
```{r}

load("inat.ok.sm.Rdata")
inat.df = inat.ok.sm
inat.df = data.frame(inat.df)
#get count by park, date, and species for all data
inat.ct <- inat.df %>%
    group_by(park, date, species) %>%
  summarize(count = n())
inat.ct$count = 1#make all counts 1
 inat.ct = data.frame(inat.ct)
 inat.ct$park_date_recordedBy = factor(paste(inat.ct$park, inat.ct$date, inat.ct$recordedBy, sep = "_"))
  inat.ct$park_date = factor(paste(inat.ct$park, inat.ct$date, sep = "_"))
 i.df = inat.ct[, c("park_date_recordedBy", "park_date", "species", "count")]
#write and then read back in again as transaction
  write.csv(i.df, file = "i.df.csv")

itrans = read.transactions("i.df.csv", format = "single", sep = ",", cols = c("park_date_recordedBy", "species"))

#find all rat observations
rat=subset(i.df, species == "Rattus norvegicus")
#find all rules involving Norway rats
  irulesRaNo <- apriori(itrans, parameter = list(support = 0.0001, confidence = 0.05, maxlen=3), appearance = list(rhs="Rattus norvegicus"))#
subrulesRaNo <- head(irulesRaNo, n = 5, by = "lift")
plot(subrulesRaNo, method = "graph", main = "rats, all data, maxlen=3")
print("Norway rat rules")
inspect(subrulesRaNo)
#find all rules for black-capped chickadees:
irulesBlCa <- apriori(itrans, parameter = list(support = 0.0001, confidence = 0.1, maxlen=3), appearance = list(rhs="Poecile atricapillus"))#
subrulesBlCa <- head(irulesBlCa, n = 5, by = "lift")
plot(subrulesBlCa, method = "graph", main = "black-capped chickadees, all data, maxlen=3")
print("Black-capped Chickadee rules")
inspect(subrulesBlCa)

#find all rules, include maxlen of 3
irules <- apriori(itrans, parameter = list(support = 0.001, confidence = 0.1, maxlen=3))#
subrules <- head(irules, n = 10, by = "lift")
plot(subrules, method = "graph", main = "all data, maxlen=3")

```

#####read in basic (not reference) ebird data
```{r}

rm(list = ls())
input_file = "/Users/fischhoff/app_wild/bird_basic_NY_NY/ebd_US-NY-061_relFeb-2018.txt"

output_file = "/Users/fischhoff/app_wild/bird_basic_NY_NY/ebd_US-NY-061_relFeb-2018clean.txt"
auk_clean(input_file, output_file, overwrite=TRUE)

ebd_NYNY = read_ebd(output_file)

```


#####find association rules in eBird NYC data 
```{r}

write.csv(ebd_NYNY, file = "ebd_NYNY.csv")
etrans = read.transactions("ebd_NYNY.csv", format = "single", sep = ",", cols = c("checklist_id", "common_name"))

erulesBlCa <- apriori(etrans, parameter = list(support = 0.001, confidence = 0.6, maxlen=3), appearance = list(rhs="Black-capped Chickadee"))#could re-run with
subrulesBlCa <- head(erulesBlCa, n = 5, by = "lift")
plot(subrulesBlCa, method = "graph", main = "black-capped chickadee, eBird NYC, N=10 rules, maxlen=3")
inspect(subrulesBlCa)

erules <- apriori(etrans, parameter = list(support = 0.01, confidence = 0.6, maxlen=3))#could re-run with
subrules <- head(erules, n = 10, by = "lift")
plot(subrules, method = "graph", main = "all species, eBird NYC, N=10 rules, maxlen=3")
plot = plot(subrules, method = "grouped")
inspect(subrules)

```