canopy_structure_analysis.Rmd

---
title: "Sorghum Canopy Structure Analysis"
author: "Kellyn Montgomery"
date: "June 17, 2019"
output: html_document
---

##Load CSMs and field data
```{r data-load}
library(raster)
library(rgdal)
library(data.table)
library(plyr)
library(lattice)
library(moments)
library(gridExtra)
library(spdep)
library(velox)
library(randomForest)
library(rgeos)
library(gstat)
library(usdm)
library(viridis)
library(car)
library(lidR)
library(diptest)
library(ggplot2)
library(rcompanion)
library(ggpubr)
library(olsrr)
library(spatialreg)
library(MuMIn)
library(MASS)

#Load raster data for all dates
#setwd("/media/Kellyn/F20E17B40E177139/kpmontgo@ncsu.edu/LkWheeler_Sorghum/LkWheeler_Fusarium_Sorghum")
setwd("G:/My Drive/LkWheeler_Sorghum/LkWheeler_Fusarium_Sorghum")

sorghum_area <- readOGR("sorghum_area", "sorghum_area", stringsAsFactors = F)

csm_1cm <- crop(raster('CSM/csm_adj_9_1_noDepth_1cm.tif'), sorghum_area)
csm_5cm <- crop(raster('CSM/csm_adj_9_1_noDepth_5cm.tif'), sorghum_area)
csm_10cm <- crop(raster('CSM/csm_adj_9_1_noDepth_10cm.tif'), sorghum_area)
csm_20cm <- crop(raster('CSM/csm_adj_9_1_noDepth_20cm.tif'), sorghum_area)

csm_1cm_velox <- velox(csm_1cm)
csm_5cm_velox <- velox(csm_5cm)
csm_10cm_velox <- velox(csm_10cm)
csm_20cm_velox <- velox(csm_20cm)

#ortho_9_1 <- crop(raster("orthos/ortho_9_1_18.tif"))

field_data <- fread("Sorghum 2018 JB Yield AUDPC.csv")
field_data$Plot <- as.integer(field_data$Plot)

# If normalizing:
# field_data_norm <- field_data[,3]
# field_data_norm$BuAc <- field_data$BuAc
# field_data_norm$AUDPC <- field_data$AUDPC
# field_data_norm <- as.data.frame(scale(field_data_norm))
# field_data_norm$Trt <- field_data$Trt
# field_data_norm$Plot <- field_data$Plot

plots <- readOGR("sorghum_plots", "sorghum_plots_small", stringsAsFactors = F)
plots <- spTransform(plots, CRS("+proj=lcc +lat_1=36.16666666666666 +lat_2=34.33333333333334 +lat_0=33.75 +lon_0=-79 +x_0=609601.22 +y_0=0 +datum=NAD83 +units=m +no_defs +ellps=GRS80 +towgs84=0,0,0"))
plots$id <- as.integer(plots$id)

plots_1cm <- plots
plots_5cm <- plots
plots_10cm <- plots
plots_20cm <- plots

plots_1cm <- merge(plots_1cm, field_data, by.x="id", by.y="Plot")
plots_5cm <- merge(plots_5cm, field_data, by.x="id", by.y="Plot")
plots_10cm <- merge(plots_10cm, field_data, by.x="id", by.y="Plot")
plots_20cm <- merge(plots_20cm, field_data, by.x="id", by.y="Plot")

```

##SIMWE water depth model
```{r simwe}

water <- crop(raster("simwe/water_depth.tif"), sorghum_area)

water_velox <- velox(water)
water_extract <- water_velox$extract(sp=plots)
names(water_extract) <- plots$id

shapiro.test(water_extract$`402`)

watermedian <- sapply(water_extract, median, na.rm=TRUE)
watermedian <- sapply(water_extract, mean, na.rm=TRUE)
wateriqr <- sapply(water_extract, IQR, na.rm=TRUE)
watersum <- sapply(water_extract, sum, na.rm=TRUE)

waterDF <- data.frame(plots=plots$id, watersum = watersum, wateriqr=wateriqr, watermedian=watermedian)
# norm_water <- as.data.frame(scale(waterDF[,2:4]))
# norm_water$plots <- waterDF$plots

plots_1cm <- merge(plots_1cm, waterDF, by.x="id", by.y="plots")
plots_5cm <- merge(plots_5cm, waterDF, by.x="id", by.y="plots")
plots_10cm <- merge(plots_10cm, waterDF, by.x="id", by.y="plots")
plots_20cm <- merge(plots_20cm, waterDF, by.x="id", by.y="plots")

```
## Water is not normally distributed. Use median and iqr.


##Canopy relief ratio
```{r crr}
crr <- function(x){
  (mean(x, na.rm=T)-min(x, na.rm=T))/(max(x, na.rm=T)-min(x, na.rm=T))
}

# calculate CRR for approx 0.5m moving window
crrRast_1cm <- focal(csm_1cm, w=matrix(1/2401, nc=49, nr=49), crr)
crrRast_5cm <- focal(csm_5cm, w=matrix(1/81, nc=9, nr=9), crr)
crrRast_10cm <- focal(csm_10cm, w=matrix(1/25, nc=5, nr=5), crr)
crrRast_20cm <- focal(csm_20cm, w=matrix(1/9, nc=3, nr=3), crr)
crrRast_list <- list(crrRast_1cm,crrRast_5cm,crrRast_10cm,crrRast_20cm)

# Check if crr within plot is normally distributed
crr_1cm_df <- as.data.frame(crrRast_1cm)
shapiro.test(sample(crr_1cm_df$layer, 5000))
histogram(crr_1cm_df$layer)
velox_crr_1cm <- velox(crrRast_1cm)
crr_extract_1cm <- velox_crr_1cm$extract(sp=plots)
names(crr_extract_1cm) <- plots$id

histogram(crr_extract_1cm$`302`)
shapiro.test(sample(crr_extract_1cm$`407`, 5000))
ggqqplot(sample(crr_extract_1cm$`302`, 5000))

# Test for multimodality in crr
library(diptest)
dip.test(crr_1cm_df$layer)
dip.test(crr_extract_1cm$`302`)

crr_metrics <- function(x){
  veloxRast <- velox(x)
  extract <- veloxRast$extract(sp=plots)
  names(extract) <- plots$id
  crrmedian <- sapply(extract, median, na.rm=T)
  crriqr <- sapply(extract, IQR, na.rm=T)
  crrDF <- data.frame(plots=plots$id, crrmedian=crrmedian, crriqr = crriqr)
  crrDF
}

crr_df_list <- lapply(crrRast_list, crr_metrics)

plots_1cm <- merge(plots_1cm, crr_df_list[[1]], by.x="id", by.y="plots")
plots_5cm <- merge(plots_5cm, crr_df_list[[2]], by.x="id", by.y="plots")
plots_10cm <- merge(plots_10cm, crr_df_list[[3]], by.x="id", by.y="plots")
plots_20cm <- merge(plots_20cm, crr_df_list[[4]], by.x="id", by.y="plots")

```
# CRR is not normally distributed. Use median and iqr to describe crr within each plot.
CRR is also multimodal.


##Crop height metrics
```{r height}

plot_extract_1cm <- csm_1cm_velox$extract(sp=plots)
plot_extract_5cm <- csm_5cm_velox$extract(sp=plots)
plot_extract_10cm <- csm_10cm_velox$extract(sp=plots)
plot_extract_20cm <- csm_20cm_velox$extract(sp=plots)
names(plot_extract_1cm) <- plots$id
names(plot_extract_5cm) <- plots$id
names(plot_extract_10cm) <- plots$id
names(plot_extract_20cm) <- plots$id

extract_list <- list(plot_extract_1cm, plot_extract_5cm, plot_extract_10cm, plot_extract_20cm)

# Check for multimodality
dip.test(plot_extract_20cm$`102`)
histogram(plot_extract_1cm$`302`)

multimodality_test <- function(x){
  plotNums <- plots$id
  dip_all <- numeric(length(plotNums))
  names(dip_all) <- "dip"
  count <- 1
  for (plot in x){
    test <- dip.test(plot)
    p <- test$p.value
    dip_all[count] <- p
    count <- count+1
  }
  dip_df <- data.frame(plots=plots$id, dip= dip_all)
  dip_df
}

dip_df_list <- lapply(extract_list, multimodality_test)

#merge into spatial polygons
plots_1cm <- merge(plots_1cm, dip_df_list[[1]], by.x="id", by.y="plots")
plots_5cm <- merge(plots_5cm, dip_df_list[[2]], by.x="id", by.y="plots")
plots_10cm <- merge(plots_10cm, dip_df_list[[3]], by.x="id", by.y="plots")
plots_20cm <- merge(plots_20cm, dip_df_list[[4]], by.x="id", by.y="plots")

# Check for normality
shapiro.test(plot_extract_20cm$`402`)



ch_metrics <- function(x){
  CHmedian <- sapply(x, median, na.rm=TRUE)
  CHmean <- sapply(x, mean, na.rm=TRUE)
  CHsd <- sapply(x, sd, na.rm=TRUE)
  CHiqr <- sapply(x, IQR, na.rm=TRUE)
  CHvar <- sapply(x, var, na.rm=TRUE)
  CHsum <- sapply(x, sum, na.rm=TRUE)
  CHskew <- sapply(x, skewness, na.rm=TRUE)
  CHkurt <- sapply(x, kurtosis, na.rm=TRUE)
  CHcrr <- sapply(x, crr)
  CHrange <- sapply(x, range, na.rm=TRUE)
  CHrange <- CHrange[2,]-CHrange[1,]
  cropHeightDF <- data.frame(plots=plots$id, median= CHmedian, mean=CHmean, sd = CHsd, iqr = CHiqr, var = CHvar, sum = CHsum, skew = CHskew, kurt= CHkurt, PlotCRR = CHcrr, chRange = CHrange)
  cropHeightDF
}

ch_df_list <- lapply(extract_list, ch_metrics)

#merge all metrics into spatial polygons
plots_1cm <- merge(plots_1cm, ch_df_list[[1]], by.x="id", by.y="plots")
plots_5cm <- merge(plots_5cm, ch_df_list[[2]], by.x="id", by.y="plots")
plots_10cm <- merge(plots_10cm, ch_df_list[[3]], by.x="id", by.y="plots")
plots_20cm <- merge(plots_20cm, ch_df_list[[4]], by.x="id", by.y="plots")

# 10th and 90th percentile yield
pct_10_90_df <- function(x){
  top_plots <- c('410', '401','407', '406', '408')
  top_yield <- x[top_plots]
  top_yield_all <- do.call(c, top_yield)
  top_yield_all <- as.data.frame(top_yield_all)
  names(top_yield_all) <- "height"
  top_yield_all$group <- "90th percentile yield"
  
  bottom_plots <- c('203', '202','205', '204', '212')
  bottom_yield <- x[bottom_plots]
  bottom_yield_all <- do.call(c, bottom_yield)
  bottom_yield_all <- as.data.frame(bottom_yield_all)
  names(bottom_yield_all) <- "height"
  bottom_yield_all$group <- "10th percentile yield"
  
  pct_10_90_yield <- rbind(top_yield_all, bottom_yield_all)
}

pct_10_90_df_list <- lapply(extract_list, pct_10_90_df)

# Histogram and denisty of top vs bottom yield
library(ggplot2)

hist_10_90 <- function(x){
  height_mean <- ddply(x, "group", summarise, height.mean=mean(height))
  height_median <- ddply(x, "group", summarise, height.median=median(height))
  ggplot(x, aes(x=height, fill=group, color=group)) +
    geom_histogram(binwidth=.03, position="identity", alpha=0.7) +
    geom_vline(data=height_median, aes(xintercept=height.median,  colour=group), linetype="solid", size=1) +
    scale_color_grey()+scale_fill_grey() +
    theme_classic() + labs(x="Crop Height (m)", y = "Count") + theme(legend.position= c(0.2,0.8), legend.spacing.x = unit(0.3, 'cm'), 
                                                                     legend.spacing.y = unit(0.2, 'cm'), legend.text = element_text(size=18),
                                                                     legend.title = element_blank(),
                                                                     axis.title.x = element_text(size=18, face="bold"),
                                                                     axis.text=element_text(size=14),
                                                                     axis.title.y = element_text(size=18, face="bold"))
}

lapply(pct_10_90_df_list, hist_10_90)
```
## Crop heights are not normally distributed. 
About half of the plots are bimodal (especially in row 1 plots with weeds).
Looking at the scatterplots, bimodality may be captured by some of the roughness metrics - iqr, rumple, kurtosis, crrsd.

```{r rumple}

library(lidR)

rumple <- function(x){
  plotNums <- plots$id
  rumple_all <- numeric(length(plotNums))
  names(rumple_all) <- "rumple"
  count <- 1
  for (num in plotNums){
    tmp <- plots@data
    position <- match(num, tmp$id)
    extract <- mask(x, plots[position,])
    extract_df <- rasterToPoints(extract)
    rumple <- rumple_index(extract_df[,1], extract_df[,2], extract_df[,3])
    rumple_all[count] <- rumple
    count <- count+1
  }
  rumple_df <- data.frame(plots=plots$id, rumple= rumple_all)
  rumple_df
}


csm_list <- list(csm_1cm, csm_5cm, csm_10cm, csm_20cm)
rumple_df_list <- lapply(csm_list, rumple)

#merge all metrics into spatial polygons
plots_1cm <- merge(plots_1cm, rumple_df_list[[1]], by.x="id", by.y="plots")
plots_5cm <- merge(plots_5cm, rumple_df_list[[2]], by.x="id", by.y="plots")
plots_10cm <- merge(plots_10cm, rumple_df_list[[3]], by.x="id", by.y="plots")
plots_20cm <- merge(plots_20cm, rumple_df_list[[4]], by.x="id", by.y="plots")

```

```{r sp-autocor}

library(ape)

library(ncf)

autocor_metrics <- function(x, w){
  plotNums <- plots$id
  morans_all <- numeric(length(plotNums))
  geary_all <- numeric(length(plotNums))
  count <- 1
  
  for (num in plotNums){
    tmp <- plots@data
    position <- match(num, tmp$id)
    extract <- crop(x, plots[position,])
    moran <- Moran(extract)
    morans_all[count] <- moran
    geary <- Geary(extract)
    geary_all[count] <- geary
    count <- count+1
  }
  
  names(morans_all) <- plots$id
  names(geary_all) <- plots$id
  autocor_df <- data.frame(plots=plots$id, moran= morans_all, geary=geary_all)
  autocor_df
}

#Using approx 30 cm x 30 cm weights matrix
autocor_1cm <- autocor_metrics(csm_1cm, matrix(c(1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,0,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1), nc=25, nr=25))
autocor_5cm <- autocor_metrics(csm_5cm, matrix(c(1,1,1,1,1,1,1,1,1,1,1,1,0,1,1,1,1,1,1,1,1,1,1,1,1), nc=5, nr=5))
autocor_10cm <- autocor_metrics(csm_10cm, matrix(c(1,1,1,1,0,1,1,1,1), nc=3, nr=3))
autocor_20cm <- autocor_metrics(csm_20cm, matrix(c(1,1,1,1,0,1,1,1,1), nc=3, nr=3))

plots_1cm <- merge(plots_1cm, autocor_1cm, by.x="id", by.y="plots")
plots_5cm <- merge(plots_5cm, autocor_5cm, by.x="id", by.y="plots")
plots_10cm <- merge(plots_10cm, autocor_10cm, by.x="id", by.y="plots")
plots_20cm <- merge(plots_20cm, autocor_20cm, by.x="id", by.y="plots")



# Try Moran.I from ape package
# ncf_cor_5cm_302 <-correlog(samp@coords[,1], samp@coords[,2], samp$csm_adj_9_1_noDepth_5cm, increment = 0.3)
# 
# w <- pointDistance(samp)
# diag(w) <- 0
# 
# tmp <- plots@data
# position <- match('302', tmp$id)
# extract <- crop(csm_5cm, plots[position,])
# samp <- sampleRandom(extract, 500, sp=T)
# samp_df <- as.data.frame(samp@data)
# moran <- Moran.I(samp[[1]], w)
# 
# autocor_metrics_ape <- function(x){
#   plotNums <- plots$id
#   morans_all <- numeric(length(plotNums))
#   #geary_all <- numeric(length(plotNums))
#   count <- 1
#   
#   for (num in plotNums){
#     tmp <- plots@data
#     position <- match(num, tmp$id)
#     extract <- crop(x, plots[position,])
#     samp <- sampleRandom(extract, 500, sp=T)
#     w <- pointDistance(samp)
#     diag(w) <- 0
#     samp_df <- as.data.frame(samp@data)
#     moran <- Moran.I(samp[[1]], w)
#     morans_all[count] <- moran[[1]]
#     #geary <- Geary(extract)
#     #geary_all[count] <- geary
#     count <- count+1
#   }
#   
#   names(morans_all) <- plots$id
#   #names(geary_all) <- plots$id
#   autocor_df <- data.frame(plots=plots$id, moran_ape= morans_all)
#   autocor_df
# }
# 
# autocor_1cm_ape <- autocor_metrics_ape(csm_1cm)
# plots_1cm <- merge(plots_1cm, autocor_1cm_ape, by.x="id", by.y="plots")

```

```{r var-norm}

library(rcompanion)

shapiro.test(plots_1cm$BuAc)
shapiro.test(plots_1cm$AUDPC)
shapiro.test(plots_1cm$median)
shapiro.test(plots_1cm$iqr)
shapiro.test(plots_1cm$skew)
shapiro.test(plots_1cm$kurt)
shapiro.test(plots_1cm$crriqr)
shapiro.test(plots_1cm$crrmedian)
shapiro.test(plots_1cm$geary)
shapiro.test(plots_1cm$moran)
shapiro.test(plots_1cm$watersum)
shapiro.test(plots_1cm$watermedian)
shapiro.test(plots_1cm$wateriqr)
shapiro.test(plots_1cm$PlotCRR)
shapiro.test(plots_1cm$rumple)
shapiro.test(plots_1cm$dip)
shapiro.test(plots_1cm$moran_ape)
shapiro.test(plots_1cm$dip)

shapiro.test(crrRast_1cm@data)

# NOTE the transformations below are not all natural log. Tukey transformation choses the best transformation method.
plots_1cm$BuAc_ln <- transformTukey(plots_1cm$BuAc)
plots_1cm$iqr_ln <- transformTukey(plots_1cm$iqr)
plots_1cm$AUDPC_ln <- transformTukey(plots_1cm$AUDPC)
plots_1cm$crrmedian_ln <- transformTukey(plots_1cm$crrmedian)
plots_1cm$watersum_ln <- transformTukey(plots_1cm$watersum)
plots_1cm$watermedian_ln <- transformTukey(plots_1cm$watermedian)
plots_1cm$wateriqr_ln <- transformTukey(plots_1cm$wateriqr)
plots_1cm$PlotCRR_ln <- transformTukey(plots_1cm$PlotCRR)
plots_1cm$rumple_ln <- transformTukey(plots_1cm$rumple)
#plots_1cm$moran_ape_ln <- transformTukey(plots_1cm$moran_ape)
#plots_1cm$dip_ln <- transformTukey(plots_1cm$dip)

plots_5cm$BuAc_ln <- transformTukey(plots_5cm$BuAc)
plots_5cm$iqr_ln <- transformTukey(plots_5cm$iqr)
plots_5cm$AUDPC_ln <- transformTukey(plots_5cm$AUDPC)
plots_5cm$crrmedian_ln <- transformTukey(plots_5cm$crrmedian)
plots_5cm$watersum_ln <- transformTukey(plots_5cm$watersum)
plots_5cm$watermedian_ln <- transformTukey(plots_5cm$watermedian)
plots_5cm$wateriqr_ln <- transformTukey(plots_5cm$wateriqr)
plots_5cm$PlotCRR_ln <- transformTukey(plots_5cm$PlotCRR)
plots_5cm$rumple_ln <- transformTukey(plots_5cm$rumple)

plots_10cm$BuAc_ln <- transformTukey(plots_10cm$BuAc)
plots_10cm$iqr_ln <- transformTukey(plots_10cm$iqr)
plots_10cm$AUDPC_ln <- transformTukey(plots_10cm$AUDPC)
plots_10cm$crrmedian_ln <- transformTukey(plots_10cm$crrmedian)
plots_10cm$watersum_ln <- transformTukey(plots_10cm$watersum)
plots_10cm$watermedian_ln <- transformTukey(plots_10cm$watermedian)
plots_10cm$wateriqr_ln <- transformTukey(plots_10cm$wateriqr)
plots_10cm$PlotCRR_ln <- transformTukey(plots_10cm$PlotCRR)
plots_10cm$rumple_ln <- transformTukey(plots_10cm$rumple)

plots_20cm$BuAc_ln <- transformTukey(plots_20cm$BuAc)
plots_20cm$iqr_ln <- transformTukey(plots_20cm$iqr)
plots_20cm$AUDPC_ln <- transformTukey(plots_20cm$AUDPC)
plots_20cm$crrmedian_ln <- transformTukey(plots_20cm$crrmedian)
plots_20cm$watersum_ln <- transformTukey(plots_20cm$watersum)
plots_20cm$watermedian_ln <- transformTukey(plots_20cm$watermedian)
plots_20cm$wateriqr_ln <- transformTukey(plots_20cm$wateriqr)
plots_20cm$PlotCRR_ln <- transformTukey(plots_20cm$PlotCRR)
plots_20cm$rumple_ln <- transformTukey(plots_20cm$rumple)

shapiro.test(plots_1cm$BuAc_ln)
shapiro.test(plots_1cm$AUDPC_ln)
shapiro.test(plots_1cm$iqr_ln)
shapiro.test(plots_1cm$crrmedian_ln)
shapiro.test(plots_1cm$watersum_ln)
shapiro.test(plots_1cm$watermedian_ln)
shapiro.test(plots_1cm$wateriqr_ln)
shapiro.test(plots_1cm$PlotCRR_ln)
shapiro.test(plots_1cm$rumple_ln)


```

## Non-normally distributed variables: BuAc, AUDPC, iqr, crrmedian, watersum, watermedian, wateriqr, PlotCRR, rumple
```{r var-cor}

library(ggpubr)

imp_vars <- c("BuAc", "median", "skew", "iqr", "kurt", "crrsd", "crrmean", "geary", "Trt", "AUDPC", "watersum", "watermedian", "wateriqr", "PlotCRR", "rumple") 

pairs(BuAc_ln~median+skew+PlotCRR_ln+crriqr+crrmedian_ln+kurt+Trt+AUDPC_ln+rumple_ln+geary+iqr_ln+dip,data=plots_5cm, lower.panel = NULL, 
      main="Simple Scatterplot Matrix")

par(mfrow=c(2,4), mar=c(5, 6, 2, 3))

cor.test(plots_5cm$AUDPC_ln, plots_5cm$BuAc_ln, method = "spearman")
ggscatter(plots_1cm@data, x = "AUDPC_ln", y = "BuAc_ln", 
          add = "reg.line", conf.int = TRUE, 
          cor.coef = TRUE, cor.method = "spearman",
          xlab = "log(AUDPC)", ylab = "log(Yield) (bu/ac)")

cor.test(plots_5cm$median, plots_5cm$BuAc_ln, method = "spearman")
ggscatter(plots_1cm@data, x = "median", y = "BuAc_ln", 
          add = "reg.line", conf.int = TRUE, 
          cor.coef = TRUE, cor.method = "spearman",
          xlab = "Median Height (m)", ylab = "log(Yield) (bu/ac)")

cor.test(plots_5cm$iqr_ln, plots_5cm$BuAc_ln, method = "spearman")
ggscatter(plots_1cm@data, x = "iqr_ln", y = "BuAc_ln", 
          add = "reg.line", conf.int = TRUE, 
          cor.coef = TRUE, cor.method = "spearman",
          xlab = "IQR (m)", ylab = "log(Yield) (bu/ac)")

cor.test(plots_5cm$rumple_ln, plots_5cm$BuAc_ln, method = "spearman")
ggscatter(plots_1cm@data, x = "rumple_ln", y = "BuAc_ln", 
          add = "reg.line", conf.int = TRUE, 
          cor.coef = TRUE, cor.method = "spearman",
          xlab = "Rumple index", ylab = "log(Yield) (bu/ac)")

cor.test(plots_5cm$skew, plots_5cm$BuAc_ln, method = "spearman")
ggscatter(plots_1cm@data, x = "skew", y = "BuAc_ln", 
          add = "reg.line", conf.int = TRUE, 
          cor.coef = TRUE, cor.method = "spearman",
          xlab = "Skewness", ylab = "log(Yield) (bu/ac)")

cor.test(plots_5cm$kurt, plots_5cm$BuAc_ln, method = "spearman")
ggscatter(plots_1cm@data, x = "kurt", y = "BuAc_ln", 
          add = "reg.line", conf.int = TRUE, 
          cor.coef = TRUE, cor.method = "spearman",
          xlab = "Kurtosis", ylab = "log(Yield) (bu/ac)")

cor.test(plots_5cm$crrmean, plots_5cm$BuAc_ln, method = "spearman")
ggscatter(plots_1cm@data, x = "crrmean", y = "BuAc_ln", 
          add = "reg.line", conf.int = TRUE, 
          cor.coef = TRUE, cor.method = "spearman",
          xlab = "CRR mean", ylab = "log(Yield) (bu/ac)")

cor.test(plots_5cm$crrsd_ln, plots_5cm$BuAc_ln, method = "spearman")
ggscatter(plots_1cm@data, x = "crrsd_ln", y = "BuAc_ln", 
          add = "reg.line", conf.int = TRUE, 
          cor.coef = TRUE, cor.method = "spearman",
          xlab = "log(CRR standard deviation)", ylab = "log(Yield) (bu/ac)")

cor.test(plots_5cm$PlotCRR_ln, plots_5cm$BuAc_ln, method = "spearman")
ggscatter(plots_1cm@data, x = "PlotCRR_ln", y = "BuAc_ln", 
          add = "reg.line", conf.int = TRUE, 
          cor.coef = TRUE, cor.method = "spearman",
          xlab = "log(Plot CRR)", ylab = "log(Yield) (bu/ac)")

cor.test(plots_1cm$geary, plots_5cm$BuAc_ln, method = "spearman")
ggscatter(plots_1cm@data, x = "geary", y = "BuAc_ln", 
          add = "reg.line", conf.int = TRUE, 
          cor.coef = TRUE, cor.method = "spearman",
          xlab = "Geary's C", ylab = "log(Yield) (bu/ac)")

cor.test(plots_1cm$moran, plots_5cm$BuAc_ln, method = "spearman")
ggscatter(plots_1cm@data, x = "moran", y = "BuAc_ln", 
          add = "reg.line", conf.int = TRUE, 
          cor.coef = TRUE, cor.method = "spearman",
          xlab = "Moran's I", ylab = "log(Yield) (bu/ac)")

cor.test(plots_5cm$watersum_ln, plots_5cm$BuAc_ln, method = "spearman")
ggscatter(plots_1cm@data, x = "watersum_ln", y = "BuAc_ln", 
          add = "reg.line", conf.int = TRUE, 
          cor.coef = TRUE, cor.method = "spearman",
          xlab = "log(Water Accumulation sum) (m)", ylab = "log(Yield) (bu/ac)")


par(mfrow=c(1,1))

# Plot histogram distribution of each variable
histogram(plots_5cm$BuAc_ln, nint=52)
histogram(plots_5cm$AUDPC_ln, nint=52)
histogram(plots_5cm$watersum_ln, nint=52)
histogram(plots_5cm$crrsd_ln, nint=52)
histogram(plots_5cm$crrmean, nint=20)
histogram(plots_5cm$median, nint=52)
histogram(plots_5cm$iqr_ln, nint=52)
histogram(plots_5cm$skew, nint=52)
histogram(plots_5cm$kurt, nint=30)
histogram(plots_5cm$PlotCRR_ln, nint=30)
histogram(plots_5cm$rumple_ln, nint=30)
histogram(plots_5cm$geary, nint=30)


```
## Visual inspection of scatterplots. Colinearity between:
skew & kurt (neg)
skew & iqr
plotcrr & crriqr & rumple & iqr
plotcrr & kurt
crriqr & crrmedian
crrmedian & rumple
all the water metrics


```{r stats}

library(olsrr)
library(spatialreg)
library(MuMIn)

# standardize data
norm_1cm <- as.data.frame(scale(plots_1cm@data[,c(21, 23, 26, 29, 30,34, 35, 36, 37, 38, 39, 40 ,41, 42, 43, 44)]), scale = F)
norm_1cm$Trt <- plots_1cm$Trt
norm_1cm$id <- plots_1cm$id

norm_5cm <- as.data.frame(scale(plots_5cm@data[,c(10,20,21, 23, 26, 29, 30,34, 35, 36, 37, 38, 39, 40 ,41, 42, 43, 44)]), center = T, scale = FALSE)
norm_5cm$Trt <- plots_5cm$Trt
norm_5cm$id <- plots_5cm$id

norm_10cm <- as.data.frame(scale(plots_10cm@data[,c(21, 23, 26, 29, 30,34, 35, 36, 37, 38, 39, 40 ,41, 42, 43, 44)]), scale = F)
norm_10cm$Trt <- plots_10cm$Trt
norm_10cm$id <- plots_10cm$id

norm_20cm <- as.data.frame(scale(plots_20cm@data[,c(21, 23, 26, 29, 30,34, 35, 36, 37, 38, 39, 40 ,41, 42, 43, 44)]), scale = F)
norm_20cm$Trt <- plots_20cm$Trt
norm_20cm$id <- plots_20cm$id

# Weights matrix, same for each resolution
w <- knn2nb(knearneigh(coordinates(plots), k=8))

# Function for calculating r^2 in SAR model, null model will be same for each resolution

null <- spautolm(BuAc_ln ~ 1, data = norm_1cm, listw=nb2listw(w), zero.policy = T, na.action = na.fail, family="SAR")
#null <- lm(BuAc ~ 1, data=norm_1cm, na.action=na.fail)
r2 <- function(model){
  r.squaredLR(model, null = null)[1]
}

adj_r2 <- function(model){
  r2 <- r.squaredLR(model, null = null)[1]
  n <- model$fit[[5]]
  p <- model$parameters
  1-(1-r2)*((n - 1)/(n - p - 1))
}

```

### All : median+skew+PlotCRR_ln+crrmedian_ln+crriqr+kurt+Trt+AUDPC_ln+watersum_ln+watermedian_ln+wateriqr_ln+rumple_ln+geary+iqr_ln

Post colinearity 
1cm: median+skew+PlotCRR_ln+crrmedian_ln+crriqr+Trt+AUDPC_ln+watermedian_ln+wateriqr_ln+geary+iqr_ln
5cm: median+skew+PlotCRR_ln+crrmedian_ln+crriqr+Trt+AUDPC_ln+watermedian_ln+wateriqr_ln+geary+iqr_ln
10cm: median+skew+PlotCRR_ln+crrmedian_ln+crriqr+Trt+AUDPC_ln+watermedian_ln+wateriqr_ln+geary+iqr_ln
20cm: median+PlotCRR_ln+crrmedian_ln+crriqr+Trt+AUDPC_ln+watermedian_ln+wateriqr_ln+rumple_ln+geary
use: median+PlotCRR_ln+crrmean+

Dip test pvalue was used in 2nd best model for 1cm, but was not significant. Omitting it dropped r2 from 0.4858 to 0.46.

Moran's I using ape package was not used in top models. Using 30cmx30cm neighborhood imporved relationship between moran/geary and yield slightly. Still not certain I'm using the neighborhood correctly though. Decided to continue using raster package autocor with 30x30 matrix.


```{r stats-1cm}

library(car)
ols_yield_1cm <- lm(BuAc_ln ~ median+skew+PlotCRR_ln+crrmedian_ln+crriqr+Trt+AUDPC_ln+watermedian_ln+wateriqr_ln+geary+iqr_ln, data=norm_1cm)
#dredge(ols_yield_1cm, beta = "none", trace = T, extra = c('r2', 'adj_r2'))
#top_ols_1cm <- lm(BuAc ~ crrmean+iqr+median+watermedian+wateriqr, data=norm_1cm)
#summary(top_ols_1cm)
vif(ols_yield_1cm)
norm_1cm$ols_res <- ols_yield_1cm$residuals
#summary(ols_yield_1cm)
moran.test(norm_1cm$ols_res, nb2listw(w))

sar_1cm <- spautolm(BuAc_ln ~ median+skew+PlotCRR_ln+crrmedian_ln+crriqr+Trt+AUDPC_ln+watermedian_ln+wateriqr_ln+geary+iqr_ln, data = norm_1cm, listw=nb2listw(w), zero.policy = T, na.action = na.fail, family="SAR")
summary(sar_1cm)
sar_resi <- sar_1cm$fit[[9]]
norm_1cm$sar_resi <- sar_resi
moran.test(norm_1cm$sar_resi, nb2listw(w))

modsel_1cm <- dredge(sar_1cm, beta = "none", trace = T, extra = c('r2', 'adj_r2'))
topmod_1cm <- spautolm(BuAc_ln ~ median + AUDPC_ln + crrmedian_ln, data = norm_1cm, listw=nb2listw(w), zero.policy = T, na.action = na.fail, family="SAR")
summary(topmod_1cm)

r2_1cm <- r2(topmod_1cm)
adj_r2_1cm <- adj_r2(topmod_1cm)

vars <- names(topmod_1cm$fit[[1]])[-1]
sd_all <- numeric(length(vars))
mean_all <- numeric(length(vars))
count <- 1
for (x in vars){
  values <- as.data.frame(plots_1cm)[x]
  sd_all[count] <- sd(values[,1])
  mean_all[count] <- mean(values[,1])
  count <- count+1
}

# backwards stepwise - median+AUDPC+geary
# sp_err_yi_1cm <- errorsarlm(BuAc ~ median+PlotCRR+crrsd+crrmean+kurt+Trt+AUDPC+watersum+geary+iqr, data = norm_1cm, listw = nb2listw(w), zero.policy = T)
# names(sp_err_yi_1cm)
#summary(sp_err_yi_1cm)
#norm_1cm$sp_err_resi <- sp_err_yi_1cm$residuals
#moran.test(norm_1cm$sp_err_resi, nb2listw(w))


```

```{r stats-5cm}

ols_yield_5cm <- lm(BuAc_ln ~ median+skew+PlotCRR_ln+crrmedian_ln+crriqr+Trt+AUDPC_ln+watermedian_ln+wateriqr_ln+geary+iqr_ln, data=norm_5cm)
car::vif(ols_yield_5cm)
#norm_5cm$ols_yi_res <- ols_yield_5cm$residuals
#summary(ols_yield_5cm)
#moran.test(norm_5cm$ols_yi_res, nb2listw(w))

sar_5cm <- spatialreg::spautolm(BuAc_ln ~ median+skew+PlotCRR_ln+crrmedian_ln+crriqr+Trt+AUDPC_ln+watermedian_ln+wateriqr_ln+geary+iqr_ln, data = norm_5cm,  listw=nb2listw(w), zero.policy = T, na.action = na.fail, family="SAR")
summary(sar_5cm)
sar_resi <- sar_5cm$fit[[9]]
norm_5cm$sar_resi <- sar_resi
moran.test(norm_5cm$sar_resi, nb2listw(w))

modsel_5cm <- dredge(sar_5cm, beta = "none", trace = T, extra = c('r2', 'adj_r2'))
modsel_5cm
topmod_5cm <- spatialreg::spautolm(BuAc_ln ~ median + crrmedian_ln + skew +Trt, data = norm_5cm, listw=nb2listw(w), zero.policy = T, na.action = na.fail, family="SAR")
summary(topmod_5cm)
r2_5cm <- r2(topmod_5cm)
adj_r2_5cm <- adj_r2(topmod_5cm)

vars <- names(topmod_5cm$fit[[1]])[-1]
sd_all <- numeric(length(vars))
mean_all <- numeric(length(vars))
count <- 1
for (x in vars){
  values <- as.data.frame(plots_5cm)[x]
  sd_all[count] <- sd(values[,1])
  mean_all[count] <- mean(values[,1])
  count <- count+1
}

# model selection results when using same dependent variables as 1cm models, ignoring collinearity analysis results


# backwards stepwise - median+Trt+geary
# sp_err_yi_5cm <- errorsarlm(BuAc ~median+crrsd+crrmedian+skew+Trt+AUDPC+watersum+geary+iqr, data = norm_5cm, listw = nb2listw(w), zero.policy = T)
# summary(sp_err_yi_5cm)
#norm_5cm$sp_err_resi <- sp_err_yi_5cm$residuals
#moran.test(norm_5cm$sp_err_resi, nb2listw(w))

```

```{r stats-10cm}

ols_yield_10cm <- lm(BuAc_ln ~ median+skew+PlotCRR_ln+crrmedian_ln+crriqr+Trt+AUDPC_ln+watermedian_ln+wateriqr_ln+geary+iqr_ln, data=norm_10cm)
car::vif(ols_yield_10cm)
#norm_10cm$ols_yi_res <- ols_yield_10cm$residuals
#summary(ols_yield_10cm)
#moran.test(norm_10cm$ols_yi_res, nb2listw(w))

sar_10cm <- spatialreg::spautolm(BuAc_ln ~ median+skew+PlotCRR_ln+crrmedian_ln+crriqr+Trt+AUDPC_ln+watermedian_ln+wateriqr_ln+geary+iqr_ln, data = norm_10cm, listw=nb2listw(w), zero.policy = T, na.action = na.fail, family="SAR")
summary(sar_10cm)
sar_resi <- sar_10cm$fit[[9]]
norm_10cm$sar_resi <- sar_resi
moran.test(norm_10cm$sar_resi, nb2listw(w))

modsel_10cm <- dredge(sar_10cm, beta = "none", trace = T, extra = c('r2', 'adj_r2'))
modsel_10cm
topmod_10cm <- spautolm(BuAc_ln ~ crrmedian_ln + median + skew + Trt, data = norm_10cm, listw=nb2listw(w), zero.policy = T, na.action = na.fail, family="SAR")
summary(topmod_10cm)
r2_10cm <- r2(topmod_10cm)
adj_r2_10cm <- adj_r2(topmod_10cm)

# model selection results when using same dependent variables as 1cm models, ignoring collinearity analysis results
tmp_modsel_10cm <- modsel_10cm


# backwards stepwise - median+crrsd+Trt+iqr
# sp_err_yi_10cm <- errorsarlm(BuAc ~ median+PlotCRR+crrsd+crrmean+kurt+Trt+AUDPC+watermedian+wateriqr+geary+iqr, data = norm_10cm, listw = nb2listw(w), zero.policy = T)
# summary(sp_err_yi_10cm)
#norm_10cm$sp_err_resi <- sp_err_yi_10cm$residuals
#moran.test(norm_10cm$sp_err_resi, nb2listw(w))

```

```{r stats-20cm}
ols_yield_20cm <- lm(BuAc_ln ~ median+PlotCRR_ln+crrmedian_ln+crriqr+Trt+AUDPC_ln+watermedian_ln+wateriqr_ln+rumple_ln+geary, data=norm_20cm)
car::vif(ols_yield_20cm)
# norm_20cm$ols_yi_res <- ols_yield_20cm$residuals
# summary(ols_yield_20cm)
# moran.test(norm_20cm$ols_yi_res, nb2listw(w))

sar_20cm <- spatialreg::spautolm(BuAc_ln ~ median+PlotCRR_ln+crrmedian_ln+crriqr+Trt+AUDPC_ln+watermedian_ln+wateriqr_ln+rumple_ln+geary, data = norm_20cm, listw=nb2listw(w), zero.policy = T, na.action = na.fail, family="SAR")
summary(sar_20cm)
sar_resi <- sar_20cm$fit[[9]]
norm_20cm$sar_resi <- sar_resi
moran.test(norm_20cm$sar_resi, nb2listw(w))

modsel_20cm <- dredge(sar_20cm, beta = "none", trace = T, extra = c('r2', 'adj_r2'))
modsel_20cm
topmod_20cm <- spatialreg::spautolm(BuAc_ln ~ median + crrmedian_ln + rumple_ln + Trt, data = norm_20cm, listw=nb2listw(w), zero.policy = T, na.action = na.fail, family="SAR")
summary(topmod_20cm)
r2_20cm <- r2(topmod_20cm)
adj_r2_20cm <- adj_r2(topmod_20cm)

# backwards stepwise - median+Trt
# sp_err_yi_20cm <- errorsarlm(BuAc ~ median+crrsd+crrmedian+skew+Trt+AUDPC+watersum+geary+iqr, data = norm_20cm, listw = nb2listw(w), zero.policy = T)
# summary(sp_err_yi_20cm)
# norm_20cm$sp_err_resi <- sp_err_yi_20cm$residuals
# moran.test(norm_20cm$sp_err_resi, nb2listw(w))

```

```{r tables}

extractResults <- function(x){
  coef <- coef(x)[-length(coef(x))]
  se <- summary(x)$rest.se
  p <- summary(x)$Coef[,4]
  results_var <- cbind(coef,se,p)
  mod_p <- summary(x)$LR1$p.value
  lambda <- summary(x)$lambda
  aic <- AIC(x)
  ll <- summary(x)$LL
  r2 <- r2(x)
  adjr2 <- adj_r2(x)
  results_mod <- cbind(ll, mod_p, lambda, aic, r2, adjr2)
  results <- list(results_var, results_mod)
}

sar_results_1cm <- extractResults(topmod_1cm)
sar_results_5cm <- extractResults(topmod_5cm)
sar_results_10cm <- extractResults(topmod_10cm)
sar_results_20cm <- extractResults(topmod_20cm)

write.csv(sar_results_1cm, "Q:/My Drive/Research/Canopy_Morphology/results/sar_results_1cm.csv")
write.csv(sar_results_5cm, "Q:/My Drive/Research/Canopy_Morphology/results/sar_results_5cm.csv")
write.csv(sar_results_10cm, "Q:/My Drive/Research/Canopy_Morphology/results/sar_results_10cm.csv")
write.csv(sar_results_20cm, "Q:/My Drive/Research/Canopy_Morphology/results/sar_results_20cm.csv")



```

```{r randomforest}

# ------------------------------------- Random Forest Classification ------------------------------------------------------
library(randomForest)
set.seed(42)
train_data_indices <- rep(FALSE, nrow(norm_1cm))
train_data_indices[sample(1:nrow(norm_1cm), round(0.7 * nrow(norm_1cm)))] <- TRUE # randomly select 80% of the data for training
rf_regression_1cm<- randomForest(BuAc_ln ~ median+skew+PlotCRR_ln+crrmedian_ln+crriqr+Trt+AUDPC_ln+geary+iqr_ln, data=norm_1cm[train_data_indices, ], importance=T)
rf_regression_1cm
varImpPlot(rf_regression_1cm)
pred_trt <- predict(rf_regression_1cm, norm_1cm[!train_data_indices,]) # predict the rings
table(pred_trt, norm_1cm$BuAc_ln[!train_data_indices])
plot(norm_1cm$BuAc_ln[!train_data_indices], pred_trt, xlab="Observed", ylab="Predicted")
#abline(a=0, b=1, lty=2, col=2)

set.seed(42)
train_data_indices <- rep(FALSE, nrow(norm_5cm))
train_data_indices[sample(1:nrow(norm_5cm), round(0.7 * nrow(norm_5cm)))] <- TRUE # randomly select 80% of the data for training
rf_regression_5cm<- randomForest(BuAc_ln ~ median+skew+PlotCRR_ln+crrmedian_ln+crriqr+Trt+AUDPC_ln+geary+iqr_ln, data=norm_5cm[train_data_indices, ], importance=T)
rf_regression_5cm
varImpPlot(rf_regression_5cm)
pred_trt <- predict(rf_regression_5cm, norm_5cm[!train_data_indices,]) # predict the rings
table(pred_trt, norm_5cm$BuAc_ln[!train_data_indices])
plot(norm_5cm$BuAc_ln[!train_data_indices], pred_trt, xlab="Observed", ylab="Predicted")
#abline(a=0, b=1, lty=2, col=2)

set.seed(42)
train_data_indices <- rep(FALSE, nrow(norm_10cm))
train_data_indices[sample(1:nrow(norm_10cm), round(0.7 * nrow(norm_10cm)))] <- TRUE # randomly select 80% of the data for training
rf_regression_10cm<- randomForest(BuAc_ln ~ median+skew+PlotCRR_ln+crrmedian_ln+crriqr+Trt+AUDPC_ln+geary+iqr_ln, data=norm_10cm[train_data_indices, ], importance=T)
rf_regression_10cm
varImpPlot(rf_regression_10cm)
pred_trt <- predict(rf_regression_10cm, norm_10cm[!train_data_indices,]) # predict the rings
table(pred_trt, norm_10cm$BuAc_ln[!train_data_indices])
plot(norm_10cm$BuAc_ln[!train_data_indices], pred_trt, xlab="Observed", ylab="Predicted")
#abline(a=0, b=1, lty=2, col=2)

set.seed(42)
train_data_indices <- rep(FALSE, nrow(norm_20cm))
train_data_indices[sample(1:nrow(norm_20cm), round(0.7 * nrow(norm_20cm)))] <- TRUE # randomly select 80% of the data for training
rf_regression_20cm<- randomForest(BuAc_ln ~ median+skew+PlotCRR_ln+crrmedian_ln+crriqr+Trt+AUDPC_ln+geary+iqr_ln, data=norm_20cm[train_data_indices, ], importance=T)
rf_regression_20cm
varImpPlot(rf_regression_20cm)
pred_trt <- predict(rf_regression_20cm, norm_20cm[!train_data_indices,]) # predict the rings
table(pred_trt, norm_20cm$BuAc_ln[!train_data_indices])
plot(norm_20cm$BuAc_ln[!train_data_indices], pred_trt, xlab="Observed", ylab="Predicted")
#abline(a=0, b=1, lty=2, col=2)

```

```{r figures}

spplot(plots_10cm, zcol="BuAc")
spplot(plots_10cm, zcol="AUDPC")

```