MIVI.Rmd

---
title: "MIVI"
author: "Will Bennett"
date: "12/14/2023"
output: pdf_document
---

# Effects of Latitude and Elevation on the Phenology of an Invasive Grass
## Japanese Stiltgrass (Microstegium vimineum)

Data source: iNaturalist ([link](https://www.inaturalist.org/taxa/116710-Microstegium-vimineum))  
Total observations (North America): 13,782  
Observations needing phenology annotation: 5,270 ([link for reviewing](https://www.inaturalist.org/observations/identify?reviewed=any&quality_grade=research&order_by=random&geoprivacy=open&place_id=97394&taxon_id=116710&verifiable=true&acc_below_or_unknown=28000&without_term_id=12))

### Setup

```{R, message=FALSE}
# set lwd
setwd("~/Documents/MIVI/")

library(dplyr)
library(ggplot2)
library(elevatr)
library(corrtable)
library(table1)
library(knitr)
```


# Loading data

## Export iNaturalist Data

1. Export all _Microstegium vimineum_ observations from iNat with columns (`id`, `observed_on`, `latitude`, `longitude`, `place_state_name`, `place_country_name`):
[Link](https://www.inaturalist.org/observations/export?quality_grade=research&identifications=any&geoprivacy=open&place_id=97394&taxon_id=116710&verifiable=true&acc_below_or_unknown=28000)

2. Export all _M. vimineum_ observations with phenology 'No Evidence of Flowering' with column `id`:
[Link](https://www.inaturalist.org/observations/export?quality_grade=research&identifications=any&geoprivacy=open&place_id=97394&taxon_id=116710&verifiable=true&acc_below_or_unknown=28000&term_id=12&term_value_id=21)

2. Export all _M. vimineum_ observations with phenology 'Flowering' with column `id`:
[Link](https://www.inaturalist.org/observations/export?quality_grade=research&identifications=any&geoprivacy=open&place_id=97394&taxon_id=116710&verifiable=true&acc_below_or_unknown=28000&term_id=12&term_value_id=13)

2. Export all _M. vimineum_ observations with phenology 'Fruiting' with column `id`:
[Link](https://www.inaturalist.org/observations/export?quality_grade=research&identifications=any&geoprivacy=open&place_id=97394&taxon_id=116710&verifiable=true&acc_below_or_unknown=28000&term_id=12&term_value_id=14)

<!--
Note: `term_value=15` is Flower Budding, unused
-->

```{R, eval=FALSE}
# 1
mivi_all <- read.csv("./MIVI-ALL.csv") %>%
    mutate(date=as.Date(observed_on, format="%Y-%m-%d")) %>% select(-observed_on)

# 2
mivi_young <- read.csv("MIVI-YOUNG.csv") %>%
    mutate(stage="Vegetation")

# 3
mivi_flowering <- read.csv("./MIVI-FLOWERING.csv") %>%
    mutate(stage="Flowering")
mivi_flowering <- mivi_flowering %>% left_join(mivi_all, by="id")

# 4
mivi_fruiting <- read.csv("./MIVI-FRUITING.csv") %>%
    mutate(stage="Fruiting")


# join each based on id
mivi_all <- mivi_all %>% left_join(mivi_young, by="id")
mivi_all <- mivi_all %>% left_join(mivi_fruiting, by="id") %>%
    mutate(stage = coalesce(stage.x, stage.y)) %>% select(-stage.x, -stage.y)
mivi_all <- rbind(mivi_all, mivi_flowering)

# memory cleanup
rm(mivi_young, mivi_flowering, mivi_fruiting)
```

### Retrieve Elevation Information

```{R, eval=FALSE}
mivi_annotated <- mivi_all %>% filter(!is.na(stage))

coords <- data.frame(x=mivi_annotated$longitude,
                     y=mivi_annotated$latitude, ele_id=mivi_annotated$id)
# pipe to 'slice(1:100)' to get a subset for reducing retrieval time

# retrieve elevation from USGS (takes a while)
elevations <- get_elev_point(coords, prj=4326, src="epqs")

mivi_all <- mivi_all %>% left_join(elevations, by=join_by("id" == "ele_id")) %>%
    select(-elev_units, -geometry)

rm(coords, elevations)

write.csv(mivi_all, file="./MIVI-PROCESSED.csv", na='')
```

## Or load from processed file

```{R}
mivi_all <- read.csv("./MIVI-PROCESSED.csv") %>% select(-X) %>%
    mutate(date=as.Date(date, format="%Y-%m-%d"))
```


# Data processing

```{R}
# remove known incorrect records
mivi_all <- subset(mivi_all, id != "130398055")

# Get Julian day
mivi_all <- mivi_all %>% mutate(julian = as.integer(strftime(date, format="%j")))

# make phenology a factor type (not necessary)
mivi_all$stage <- factor(mivi_all$stage, ordered=TRUE,
                         levels=c("Vegetation", "Flowering", "Fruiting"))

# Select only observations with phenology data
mivi_annotated <- mivi_all %>% filter(!is.na(stage))

# Group into quartiles by latitude
mivi_annotated$group <- ntile(mivi_annotated$latitude, 4)
```


# Descriptive Tables

```{R}
# All observations by country
table1(~ place_country_name, data=mivi_all)

# All observations by state (includes Ontario)
table1(~ place_state_name, data=mivi_all)

table1(~ latitude + longitude + elevation + julian + place_country_name | stage, data=mivi_all %>% mutate(stage=addNA(stage)))
```


# Basic Descriptive Plots

```{R}
# Density of all observations by lat/lon
ggplot(mivi_all, aes(x=longitude,y=latitude)) + geom_point(alpha=0.25) +
    labs(title="Density of observations by location")

# Histogram of all observations by lat
ggplot(mivi_all, aes(x=latitude)) + geom_histogram(binwidth=0.2) +
    labs(title="Histogram of observations by latitude")

# Histogram with city labels
ggplot(mivi_all, aes(x=latitude)) + geom_histogram(binwidth=0.2) +
    geom_vline(xintercept=33.75) + annotate("text", x=33.5, y=900, label="Atlanta", angle=90) +
    geom_vline(xintercept=35.84) + annotate("text", x=35.59, y=900, label="Raleigh", angle=90) +
    geom_vline(xintercept=37.54) + annotate("text", x=37.29, y=900, label="Richmond", angle=90) +
    geom_vline(xintercept=38.9) + annotate("text", x=38.65, y=900, label="Washington", angle=90) +
    geom_vline(xintercept=40.75) + annotate("text", x=41, y=900, label="New York", angle=90) +
    ylab("count") + labs(title="Histogram of observations by latitude")

# Boxplots of phenology by quartile (notches break hinges for Flowering)
# Slight earlier trend in Fruiting visible
ggplot(mivi_annotated, aes(julian, stage)) + geom_boxplot() + facet_grid(~group) +
    coord_flip() + xlab("Julian day") + labs(title="Phenology stage, quartiles by latitude")
```

## Time Series Plots
```{R}
# Latitude against Julian day
ggplot(mivi_annotated, aes(julian, latitude)) + geom_point(aes(color=stage), alpha=0.5) +
    scale_color_hue() + xlab("Julian day") +
    labs(title="Annotated Observations by Latitude and Julian Day")

# Latitude against Julian day
# ggplot(mivi_all, aes(julian, latitude)) + geom_point(aes(color=stage), alpha=0.5) +
#     scale_color_hue() + xlab("Julian day") +
#     labs(title="Observations by Latitude and Julian Day")

# Time series by latitude, color by stage
ggplot(mivi_all, aes(date, latitude)) + geom_point(aes(color=stage), alpha=0.5) +
    scale_x_date(date_breaks = "1 year", date_labels = "%Y") + scale_color_hue() +
    labs(title="Observations by latitude over time")

# Zoom in on recent data
timeclip <- c(as.Date("2019-02-01"), as.Date("2023-11-30"))
ggplot(mivi_all, aes(date, latitude)) + geom_point(aes(color=stage), alpha=0.5) +
    scale_x_date(limits=timeclip, date_breaks = "1 year", date_labels = "%Y") + scale_color_hue() +
    labs(title="Observations by latitude over time (2019-2023)")
```


# Analysis!

```{R}
# Note: returns Inf if there are none in the selection
first_fruit <- function(df) {
    df <- df %>% filter(stage == "Fruiting")
    if(nrow(df)==0) {return (Inf)}
    return (min(df$julian))
}
```

## Latitude

```{R}
get_lat_quants <- function(df, n) {
    # create groups
    df$group <- ntile(df$latitude, n)
    
    # calculate latitude variables for each group
    a <- df %>% group_by(group) %>% summarize(avglat=mean(latitude),
                                              minlat=min(latitude),
                                              maxlat=max(latitude))
    
    # first fruiting date in each group
    b <- df %>% group_by(group) %>% group_map(~first_fruit(.x))
    
    a <- bind_cols(a, do.call(rbind.data.frame, b)[,1], .name_repair = "unique_quiet") %>%
        mutate(firstfruit = ...5) %>% select(-...5)
    a <- remove_missing(a, finite=TRUE, na.rm=TRUE)  # remove any Inf's
    
    return (a)
}

# Loop to determine best quantile amount for correlation
quants = data.frame()

for(i in 10:300) {
    # print(i)
    a <- get_lat_quants(mivi_annotated, i)
    quants <- rbind(quants, data.frame(i, cor(a$firstfruit, a$avglat)))
}
rm(i, a)

colnames(quants) <- c("n quantiles", "correlation")
plot(quants)

y <- which.min(quants$corr)
n <- quants[y,1]  # select n with strongest correlation
paste0("n quantiles for best correlation is: ", n)

rm(quants, y)


# Latitude linear model
data <- get_lat_quants(mivi_annotated, n)
model_lat <- lm(firstfruit~avglat, data=data)

# Plot linear model
print(
    ggplot(data, aes(avglat, firstfruit)) + geom_point() + geom_smooth(method='lm') +
        ylab("Julian day") + xlab("mean latitude") + labs(title="First Fruiting Day by Latitude")
)

# Q-Q residual plot
res <- resid(model_lat)
qqnorm(res)
qqline(res)

paste0("Average absolute residual: ", format(mad(res), digits=6))
paste0("avglat coefficient is: ", format(coefficients(model_lat)["avglat"], digits=5))

# Pearson's correlation test
cor.test(data$firstfruit, data$avglat, alternative="less")

correlation_matrix(data, use="lower")

# kable(correlation_matrix(data, use="lower"), booktabs=TRUE, format="latex")
# kable(correlation_matrix(data, use="lower"), booktabs=TRUE)

rm(data, model_lat, res)
```

## Elevation

```{R}
get_ele_quants <- function(df, n) {
    # Create groups
    df$group <- ntile(mivi_annotated$elevation, n)
    
    # mean lat for each group
    a <- df %>% group_by(group) %>% summarize(avgele=mean(elevation))
    
    # first fruiting date in each group
    b <- df %>% group_by(group) %>% group_map(~first_fruit(.x))
    
    a <- bind_cols(a, do.call(rbind.data.frame, b)[,1], .name_repair = "unique_quiet") %>%
        mutate(firstfruit = ...3) %>% select(-...3)
    a <- remove_missing(a, finite=TRUE, na.rm=TRUE)  # remove any Inf's
    
    return (a)
}

# Loop to determine best quantile amount for correlation
quants = data.frame()

for(i in 10:300) {
    # print(i)
    a <- get_ele_quants(mivi_annotated, i)
    quants <- rbind(quants, data.frame(i, cor(a$firstfruit, a$avgele)))
}
rm(i, a)

colnames(quants) <- c("n quantiles", "correlation")
plot(quants)

y <- which.min(quants$corr)
n <- quants[y,1]  # select n with strongest correlation
paste0("n quantiles for best correlation is: ", n)

rm(quants, y)


# Elevation linear model
data = get_ele_quants(mivi_annotated, n)
model_ele <- lm(firstfruit~avgele, data=data)

# Plot linear model
ggplot(data, aes(avgele, firstfruit)) + geom_point() + geom_smooth(method='lm') +
    ylab("Julian day") + xlab("mean elevation") + labs(title="First Fruiting Day by Elevation")

# Q-Q residual plot
res <- resid(model_ele)
qqnorm(res)
qqline(res)

paste0("Average absolute residual: ", format(mad(res), digits=6))
paste0("avgele coefficient is: ", format(coefficients(model_ele)["avgele"], digits=5))

# Pearson's correlation test
cor.test(data$firstfruit, data$avgele, alternative="less")

rm(data, model_ele, res)
```

<!--
data <- data.frame("category" = c('Annotated', 'Not Annotated'), "amount" = c(2952, 12066))
ggplot(data, aes(x="", y=amount, fill=category)) +
    geom_bar(stat="identity", width=1) +
    coord_polar("y", start=0)
-->