---
title: "Limiting nutrients drive mountain stream ecosystem processes along an elevation gradients"
author: "Beck, Rugenski, and Poff 2021"
date: "March 11, 2021"
output: html_document
---

```{r setup, include=FALSE}
knitr::opts_chunk$set(echo = TRUE)
```

##Preparing Data

```{r}
#load packages
library(dplyr)
library(tidyr)
library(ggplot2)
library(devtools)
library(ggfortify)
library('vtreat')
library('WVPlots') 
library(ggpubr)

#read data
Data<-read.csv("Data Summary.csv",header=T)

#use first column to assign rownames 
row.names(Data)<-Data$PCA.Name

#remove first column
Data<-Data[,-1]

#Take out correlated variables 
#Take out the gradients themselves (Elevation & Season) and watershed factors (watershed area, forest cover) to focus on mechanistic environmental factors
#Take out response variables
Data_nocorr <- Data%>%select(-NO3,-T_Min, -T_Max, -Q, -pH, -Site, -Elevation, -Season, -WA, -Forest, -N_LRR, -NP_LRR, -P_LRR, -Accrual)

```

##PCA and Initial Plots

```{r}

#run PCA
Data_pca<-prcomp(Data_nocorr,scale.=T, center=T) #scaling variables since they have different units

#view output
summary(Data_pca) #return proportions of variance explained

Data_pca$rotation #print eigenvectors, or variable loadings

#Plot using ggplot
Second=autoplot(Data_pca, loadings=TRUE, loadings.label=TRUE, 
                loadings.label.hjust=1.2, 
                loadings.label.vjust=-1.1, 
loadings.label.size=3, loadings.label.colour="black", loadings.colour="black")+
theme(axis.text.x=element_text(colour="black"), panel.grid.minor = element_blank(), 
panel.grid.major = element_blank(), panel.background=element_blank(), 
axis.text.y=element_text(colour="black"), axis.line.x = element_line(color="black"), axis.line.y = element_line(color="black"))+
  xlim(-0.5, 0.5)+
  ylim(-0.5, 0.5)+
  ggtitle(bquote(bold(A)))

Third=autoplot(Data_pca, x=1, y=3, loadings=TRUE, loadings.label=TRUE, loadings.label.hjust=0.65, loadings.label.vjust=-1.5, loadings.label.size=3, loadings.label.colour="black", loadings.colour="black")+
  theme(axis.text.x=element_text(colour="black"), panel.grid.minor = element_blank(), panel.grid.major = element_blank(), panel.background=element_blank(), axis.text.y=element_text(colour="black"), axis.line.x = element_line(color="black"), axis.line.y = element_line(color="black"))+
  xlim(-0.5, 0.5)+
  ylim(-0.5, 0.5)+
  ggtitle(bquote(bold(B)))

figure <- ggarrange(Second, Third,
                    ncol = 1, nrow = 2)

#Save file to desktop folder
ggsave("PCA_plot.png", width=13, height=20, units="cm")





```

##Examine PCA variance explained to determine how many PCA components to include in regression models 
##Visualize variable loadings

```{r}

# extract the principal components
extractProjection <- function(ndim,princ) {
  proj <- princ$rotation[,1:ndim]}

rsq <- function(x,y) {
  1 - sum((y-x)^2)/sum((y-mean(y))^2)
}
PCs <- extractProjection(4 , Data_pca) #This is specifying that we are extracting the first 4 PCs to use as predictors for the multiple regression.

# prepare the data to plot the variable loadings
PCdf = as.data.frame(PCs)
PCdf$varName = rownames(PCdf)
PCdflong = gather(PCdf, "PC", "loading",
                       starts_with("PC"))
PCdflong$vartype = ifelse(grepl("noise", 
                                     PCdflong$varName),
                               "noise", "signal")

# plot the singular values
dotplot_identity = function(frame, xvar, yvar, colorvar=NULL) {
  if(is.null(colorvar)) {
    gplot = ggplot(frame, aes_string(x=xvar, y=yvar, ymax=yvar))
  } else {
    gplot = ggplot(frame, aes_string(x=xvar, y=yvar, ymax=yvar, color=colorvar))
  }
  
  gplot + geom_point() + geom_linerange(aes(ymin=0))
}

dotplot_identity(frame = data.frame(pc=1:length(Data_pca$sdev), 
                                    magnitude=Data_pca$sdev), 
                 xvar="pc",yvar="magnitude") +
  ggtitle("Magnitudes of singular values")

#Plot variable loadings
dotplot_identity(PCdflong, "varName", "loading", "vartype") + 
  facet_wrap(~PC,nrow=1) + coord_flip() + 
  ggtitle("x scaled variable loadings")

ggsave("Variable Loadings.png", width=24, height=16, units="cm")


```

##N-LRR Model Selection

```{r}

#Add Data for response variable
Resp=read.csv("Data Summary.csv")

#Scaling the response variable and adding it to the dataframe with the PCs
Resp_scaled <- 
  as.data.frame(scale(Resp[,colnames(Resp)=='N_LRR'],
                      center=TRUE,scale=TRUE),stringsAsFactors = FALSE)
projectedTrain <- as.data.frame(predict(Data_pca,Data_nocorr),
                                 stringsAsFactors = FALSE)
projectedTrain$y <- Resp_scaled$V1

#AICc for all possible models
model1=lm(y~PC1, data=projectedTrain)
model2=lm(y~PC2, data=projectedTrain)
model3=lm(y~PC3, data=projectedTrain)
model4=lm(y~PC4, data=projectedTrain)
model5=lm(y~PC1+PC2, data=projectedTrain)
model6=lm(y~PC1+PC3, data=projectedTrain)
model7=lm(y~PC1+PC4, data=projectedTrain)
model8=lm(y~PC2+PC3, data=projectedTrain)
model9=lm(y~PC3+PC4, data=projectedTrain)
model10=lm(y~PC2+PC4, data=projectedTrain)
model11=lm(y~PC1+PC2+PC3, data=projectedTrain)
model12=lm(y~PC1+PC2+PC4, data=projectedTrain)
model13=lm(y~PC1+PC3+PC4, data=projectedTrain)
model14=lm(y~PC2+PC3+PC4, data=projectedTrain)
model15=lm(y~PC1+PC2+PC3+PC4, data=projectedTrain)

AICc(model1)
AICc(model2)
AICc(model3)
AICc(model4)
AICc(model5)
AICc(model6)
AICc(model7)
AICc(model8)
AICc(model9)
AICc(model10)
AICc(model11)
AICc(model12)
AICc(model13)
AICc(model14)
AICc(model15)

#Examining summary of the model with the lowest AICc
summary(model6)
```

##P-LRR Model Selection

```{r}
#Scaling the response variable and adding it to the dataframe with the PCs
Resp_scaled <- 
  as.data.frame(scale(Resp[,colnames(Resp)=='P_LRR'],
                      center=TRUE,scale=TRUE),stringsAsFactors = FALSE)
projectedTrain <- as.data.frame(predict(Data_pca,Data_nocorr),
                                 stringsAsFactors = FALSE)
projectedTrain$y <- Resp_scaled$V1

#AICc for all possible models

model1=lm(y~PC1, data=projectedTrain)
model2=lm(y~PC2, data=projectedTrain)
model3=lm(y~PC3, data=projectedTrain)
model4=lm(y~PC4, data=projectedTrain)
model5=lm(y~PC1+PC2, data=projectedTrain)
model6=lm(y~PC1+PC3, data=projectedTrain)
model7=lm(y~PC1+PC4, data=projectedTrain)
model8=lm(y~PC2+PC3, data=projectedTrain)
model9=lm(y~PC3+PC4, data=projectedTrain)
model10=lm(y~PC2+PC4, data=projectedTrain)
model11=lm(y~PC1+PC2+PC3, data=projectedTrain)
model12=lm(y~PC1+PC2+PC4, data=projectedTrain)
model13=lm(y~PC1+PC3+PC4, data=projectedTrain)
model14=lm(y~PC2+PC3+PC4, data=projectedTrain)
model15=lm(y~PC1+PC2+PC3+PC4, data=projectedTrain)

AICc(model1)
AICc(model2)
AICc(model3)
AICc(model4)
AICc(model5)
AICc(model6)
AICc(model7)
AICc(model8)
AICc(model9)
AICc(model10)
AICc(model11)
AICc(model12)
AICc(model13)
AICc(model14)
AICc(model15)

#Examining summary of the model with the lowest AICc
summary(model2)
```

##NP-LRR Model Selection

```{r}

#Scaling the response variable and adding it to the dataframe with the PCs
Resp_scaled <- 
  as.data.frame(scale(Resp[,colnames(Resp)=='NP_LRR'],
                      center=TRUE,scale=TRUE),stringsAsFactors = FALSE)
projectedTrain <- as.data.frame(predict(Data_pca,Data_nocorr),
                                 stringsAsFactors = FALSE)
projectedTrain$y <- Resp_scaled$V1

#AICc for all possible models

model1=lm(y~PC1, data=projectedTrain)
model2=lm(y~PC2, data=projectedTrain)
model3=lm(y~PC3, data=projectedTrain)
model4=lm(y~PC4, data=projectedTrain)
model5=lm(y~PC1+PC2, data=projectedTrain)
model6=lm(y~PC1+PC3, data=projectedTrain)
model7=lm(y~PC1+PC4, data=projectedTrain)
model8=lm(y~PC2+PC3, data=projectedTrain)
model9=lm(y~PC3+PC4, data=projectedTrain)
model10=lm(y~PC2+PC4, data=projectedTrain)
model11=lm(y~PC1+PC2+PC3, data=projectedTrain)
model12=lm(y~PC1+PC2+PC4, data=projectedTrain)
model13=lm(y~PC1+PC3+PC4, data=projectedTrain)
model14=lm(y~PC2+PC3+PC4, data=projectedTrain)
model15=lm(y~PC1+PC2+PC3+PC4, data=projectedTrain)

AICc(model1)
AICc(model2)
AICc(model3)
AICc(model4)
AICc(model5)
AICc(model6)
AICc(model7)
AICc(model8)
AICc(model9)
AICc(model10)
AICc(model11)
AICc(model12)
AICc(model13)
AICc(model14)
AICc(model15)

#Examining summary for the model with the lowest AICc
summary(model6)

```

##Accrual Model Selection

```{r}

#Scaling the response variable and adding it to the dataframe with the PCs
Resp_scaled <- 
  as.data.frame(scale(Resp[,colnames(Resp)=='Accrual'],
                      center=TRUE,scale=TRUE),stringsAsFactors = FALSE)
projectedTrain <- as.data.frame(predict(Data_pca,Data_nocorr),
                                 stringsAsFactors = FALSE)
projectedTrain$y <- Resp_scaled$V1

#AICc for all possible models

model1=lm(y~PC1, data=projectedTrain)
model2=lm(y~PC2, data=projectedTrain)
model3=lm(y~PC3, data=projectedTrain)
model4=lm(y~PC4, data=projectedTrain)
model5=lm(y~PC1+PC2, data=projectedTrain)
model6=lm(y~PC1+PC3, data=projectedTrain)
model7=lm(y~PC1+PC4, data=projectedTrain)
model8=lm(y~PC2+PC3, data=projectedTrain)
model9=lm(y~PC3+PC4, data=projectedTrain)
model10=lm(y~PC2+PC4, data=projectedTrain)
model11=lm(y~PC1+PC2+PC3, data=projectedTrain)
model12=lm(y~PC1+PC2+PC4, data=projectedTrain)
model13=lm(y~PC1+PC3+PC4, data=projectedTrain)
model14=lm(y~PC2+PC3+PC4, data=projectedTrain)
model15=lm(y~PC1+PC2+PC3+PC4, data=projectedTrain)

AICc(model1)
AICc(model2)
AICc(model3)
AICc(model4)
AICc(model5)
AICc(model6)
AICc(model7)
AICc(model8)
AICc(model9)
AICc(model10)
AICc(model11)
AICc(model12)
AICc(model13)
AICc(model14)
AICc(model15)

#Examining summary for model with lowest AICc

summary(model6)

```

