R Notebook

This is an R Markdown Notebook. When you execute code within the notebook, the results appear beneath the code.

Try executing this chunk by clicking the Run button within the chunk or by placing your cursor inside it and pressing Cmd+Shift+Enter.

plot(cars)

Add a new chunk by clicking the Insert Chunk button on the toolbar or by pressing Cmd+Option+I.

When you save the notebook, an HTML file containing the code and output will be saved alongside it (click the Preview button or press Cmd+Shift+K to preview the HTML file).

The preview shows you a rendered HTML copy of the contents of the editor. Consequently, unlike Knit, Preview does not run any R code chunks. Instead, the output of the chunk when it was last run in the editor is displayed.

summary(cars)

     speed           dist       
 Min.   : 4.0   Min.   :  2.00  
 1st Qu.:12.0   1st Qu.: 26.00  
 Median :15.0   Median : 36.00  
 Mean   :15.4   Mean   : 42.98  
 3rd Qu.:19.0   3rd Qu.: 56.00  
 Max.   :25.0   Max.   :120.00  

library(terra)
library(landscapemetrics)
library(sf)
library(tidyverse)
library(ggplot2)
library(tigris)
library(rgdal)
library(raster)
library(kimisc)
library(AICcmodavg)

## Read in Data

states=states(cb=FALSE)
va=filter(states, NAME=="Virginia") ##Virginia boundary
##huc12=st_read("HUC12_shapefile") ##HUC12 tributary systems
watershed=st_read("watershed_shapefile") ##Chesapeaky Bay Watershed boundary

Reading layer `Chesapeake_Bay_Watershed_Boundary' from data source 
  `/Users/joshuasementelli/Downloads/Final Project/knitting_project/watershed_shapefile' 
  using driver `ESRI Shapefile'
Simple feature collection with 1 feature and 3 fields
Geometry type: POLYGON
Dimension:     XY
Bounding box:  xmin: -80.54099 ymin: 36.66942 xmax: -74.58096 ymax: 42.98733
Geodetic CRS:  WGS 84

##wq=read.csv("waterquality1.csv") ##WAter quality data

##Project data

proj=crs(va)
va_watershed=st_transform(watershed, proj)
crs(va_watershed)

CRS arguments: +proj=longlat +datum=NAD83 +no_defs 

##Intersecting Chesapeake Bay boundary with VA state boundary

va_watershed=st_intersection(va_watershed, va) ##Intersecting watershed & VA data

Warning: attribute variables are assumed to be spatially constant throughout all geometries

plot(st_geometry(va_watershed))

##Pull in raster data

varast=rast("nlcd1.tiff")
plot(varast)

crs(varast)

[1] "PROJCRS[\"Albers_Conical_Equal_Area\",\n    BASEGEOGCRS[\"WGS 84\",\n        DATUM[\"World Geodetic System 1984\",\n            ELLIPSOID[\"WGS 84\",6378137,298.257223563,\n                LENGTHUNIT[\"metre\",1]]],\n        PRIMEM[\"Greenwich\",0,\n            ANGLEUNIT[\"degree\",0.0174532925199433]],\n        ID[\"EPSG\",4326]],\n    CONVERSION[\"Albers Equal Area\",\n        METHOD[\"Albers Equal Area\",\n            ID[\"EPSG\",9822]],\n        PARAMETER[\"Latitude of false origin\",23,\n            ANGLEUNIT[\"degree\",0.0174532925199433],\n            ID[\"EPSG\",8821]],\n        PARAMETER[\"Longitude of false origin\",-96,\n            ANGLEUNIT[\"degree\",0.0174532925199433],\n            ID[\"EPSG\",8822]],\n        PARAMETER[\"Latitude of 1st standard parallel\",29.5,\n            ANGLEUNIT[\"degree\",0.0174532925199433],\n            ID[\"EPSG\",8823]],\n        PARAMETER[\"Latitude of 2nd standard parallel\",45.5,\n            ANGLEUNIT[\"degree\",0.0174532925199433],\n            ID[\"EPSG\",8824]],\n        PARAMETER[\"Easting at false origin\",0,\n            LENGTHUNIT[\"metre\",1],\n            ID[\"EPSG\",8826]],\n        PARAMETER[\"Northing at false origin\",0,\n            LENGTHUNIT[\"metre\",1],\n            ID[\"EPSG\",8827]]],\n    CS[Cartesian,2],\n        AXIS[\"easting\",east,\n            ORDER[1],\n            LENGTHUNIT[\"metre\",1,\n                ID[\"EPSG\",9001]]],\n        AXIS[\"northing\",north,\n            ORDER[2],\n            LENGTHUNIT[\"metre\",1,\n                ID[\"EPSG\",9001]]]]"

rastproj=crs(varast)

##Reclassify data
#NLCD classes:
#11 (open water)
#21 (Developed open space - <20% imperv)
#22 (low intensity development, 20-49% imperv)
#23 (high intensity development, 50-79% imperv)
#24 (very high intesnity development, >80%imperv)
#31 (barren)
#41 (deciduous forest)
#42 (evergreen forest)
#43 (mixed forest)
#52 (shrub)
#71 (grassland) 
#81 (pasture/hay)
#82 (cultivated crop)
#90 (woody wetland)
#95 (emergent wetland)

## Reclassify raster to only include Developed land, forest cover, and agriculture

reclass=c(rep(0,2), rep(21,4), rep(0,1), rep(41,3), rep(0,2), rep(81,2), rep(0,2))
reclass

 [1]  0  0 21 21 21 21  0 41 41 41  0  0 81 81  0  0

reclass.mat=cbind(unique(varast)[[1]], reclass)
reclass.mat

         reclass
 [1,]  0       0
 [2,] 11       0
 [3,] 21      21
 [4,] 22      21
 [5,] 23      21
 [6,] 24      21
 [7,] 31       0
 [8,] 41      41
 [9,] 42      41
[10,] 43      41
[11,] 52       0
[12,] 71       0
[13,] 81      81
[14,] 82      81
[15,] 90       0
[16,] 95       0

classtype19=classify(varast, reclass.mat)

=====|====|====|====|====|====|====|====|====|====
--------------------------------------------------
                                                  

plot(classtype19)
plot(st_geometry(va_watershed), add=TRUE)

##Water Quality Data

## Read in water quality data and match projections with raster

wq=st_read("waterqualitygood", crs=3968) ##Assigngning crs data

Reading layer `Huc12_subset' from data source 
  `/Users/joshuasementelli/Downloads/Final Project/knitting_project/waterqualitygood' 
  using driver `ESRI Shapefile'
Simple feature collection with 36 features and 12 fields
Geometry type: POLYGON
Dimension:     XY
Bounding box:  xmin: -40573.28 ymin: 128871.5 xmax: 260078.6 ymax: 365919.6
Projected CRS: NAD83 / Virginia Lambert

wq=st_transform(wq, rastproj)
plot(classtype19)
plot(st_geometry(wq), add=TRUE)

##Transforming va watershed crs to match raster projection

va_watershed1=st_transform(va_watershed, rastproj)
plot(st_geometry(va_watershed1))
plot(st_geometry(wq), add=TRUE)

### Class level landscape metrics: percent of the landcape (pland), number of patches (np), largest patch index (lpi), patch density (pd), edge density (ed)

cl_metrics_wq<-sample_lsm(landscape = classtype19, y = wq, what = c("lsm_c_pland", "lsm_c_np", "lsm_c_lpi", "lsm_c_pd", "lsm_c_ed"))
## Filter out 0's that existed within class types

cl_metrics_wq=filter(cl_metrics_wq, class != "0")

## Create and assign plot id column values to merge with landscape metrics
wq <- tibble::rowid_to_column(wq, "plot_id")

cl_metrics19<- merge(cl_metrics_wq, wq, by="plot_id")


cl_metrics19$class=as.factor(cl_metrics19$class)
summary(cl_metrics19$class)

 21  41  81 
180 180 180 

## Spread out the dataframe for organization   

cl_metrics19_wide=cl_metrics19 %>% group_by(plot_id)  %>% 
  pivot_wider(names_from = c(metric, class), values_from=value)

#Linear Models of Water Quality Pollutants and Forest level landscape metrics

# Dissolved Oxygen 
DO_F_np<- lm(DO~np_41, data = cl_metrics19_wide)
DO_F_pland<- lm(DO~pland_41, data = cl_metrics19_wide)
DO_F_lpi<- lm(DO~lpi_41, data = cl_metrics19_wide)
summary(DO_F_lpi)

Call:
lm(formula = DO ~ lpi_41, data = cl_metrics19_wide)

Residuals:
     Min       1Q   Median       3Q      Max 
-2.36751 -0.55388  0.08003  0.67178  2.47673 

Coefficients:
            Estimate Std. Error t value Pr(>|t|)    
(Intercept)  9.01773    0.30460  29.605  < 2e-16 ***
lpi_41       0.05323    0.01718   3.099  0.00389 ** 
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 1.206 on 34 degrees of freedom
Multiple R-squared:  0.2202,    Adjusted R-squared:  0.1973 
F-statistic: 9.602 on 1 and 34 DF,  p-value: 0.003887

DO_F_pd<- lm(DO~pd_41, data = cl_metrics19_wide)
DO_F_ed<- lm(DO~ed_41, data = cl_metrics19_wide)

## Total Nitrogen
TN_F_np<- lm(TN~np_41, data = cl_metrics19_wide)
TN_F_pland<- lm(TN~pland_41, data = cl_metrics19_wide)
TN_F_lpi<- lm(TN~lpi_41, data = cl_metrics19_wide)
TN_F_pd<- lm(TN~pd_41, data = cl_metrics19_wide)
TN_F_ed<- lm(TN~ed_41, data = cl_metrics19_wide)

## Total Phosphorus
TP_F_np<- lm(TP~np_41, data = cl_metrics19_wide)
TP_F_pland<- lm(TP~pland_41, data = cl_metrics19_wide)
summary(TP_F_pland)

Call:
lm(formula = TP ~ pland_41, data = cl_metrics19_wide)

Residuals:
     Min       1Q   Median       3Q      Max 
-0.05942 -0.02340 -0.01013  0.01214  0.08965 

Coefficients:
              Estimate Std. Error t value Pr(>|t|)    
(Intercept)  0.1217783  0.0193586   6.291  8.4e-07 ***
pland_41    -0.0010619  0.0003392  -3.131  0.00405 ** 
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 0.03753 on 28 degrees of freedom
  (6 observations deleted due to missingness)
Multiple R-squared:  0.2593,    Adjusted R-squared:  0.2328 
F-statistic: 9.802 on 1 and 28 DF,  p-value: 0.004054

TP_F_lpi<- lm(TP~lpi_41, data = cl_metrics19_wide)
TP_F_pd<- lm(TP~pd_41, data = cl_metrics19_wide)
TP_F_ed<- lm(TP~ed_41, data = cl_metrics19_wide)

## Total Suspended Solids
TSS_F_np<- lm(TSS~np_41, data = cl_metrics19_wide)
TSS_F_pland<- lm(TSS~pland_41, data = cl_metrics19_wide)
TSS_F_lpi<- lm(TSS~lpi_41, data = cl_metrics19_wide)
TSS_F_pd<- lm(TSS~pd_41, data = cl_metrics19_wide)
TSS_F_ed<- lm(TSS~ed_41, data = cl_metrics19_wide)

## Linear Models of Water Quality Pollutants and impervious surface level landscape metrics

# Dissolved Oxygen 
DO_imperv_np<- lm(DO~np_21, data = cl_metrics19_wide)
DO_imperv_pland<- lm(DO~pland_21, data = cl_metrics19_wide)
summary(DO_imperv_pland)

Call:
lm(formula = DO ~ pland_21, data = cl_metrics19_wide)

Residuals:
    Min      1Q  Median      3Q     Max 
-2.5052 -0.8683  0.1239  1.0665  1.7734 

Coefficients:
             Estimate Std. Error t value Pr(>|t|)    
(Intercept) 10.293061   0.281388  36.580  < 2e-16 ***
pland_21    -0.027159   0.009323  -2.913  0.00628 ** 
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 1.221 on 34 degrees of freedom
Multiple R-squared:  0.1997,    Adjusted R-squared:  0.1762 
F-statistic: 8.486 on 1 and 34 DF,  p-value: 0.006282

DO_imperv_lpi<- lm(DO~lpi_21, data = cl_metrics19_wide)
DO_imperv_pd<- lm(DO~pd_21, data = cl_metrics19_wide)
DO_imperv_ed<- lm(DO~ed_21, data = cl_metrics19_wide)

## Total Nitrogen
TN_imperv_np<- lm(TN~np_21, data = cl_metrics19_wide)
TN_imperv_pland<- lm(TN~pland_21, data = cl_metrics19_wide)
TN_imperv_lpi<- lm(TN~lpi_21, data = cl_metrics19_wide)
TN_imperv_pd<- lm(TN~pd_21, data = cl_metrics19_wide)
TN_imperv_ed<- lm(TN~ed_21, data = cl_metrics19_wide)

## Total Phosphorus
TP_imperv_np<- lm(TP~np_21, data = cl_metrics19_wide)
TP_imperv_pland<- lm(TP~pland_21, data = cl_metrics19_wide)
TP_imperv_lpi<- lm(TP~lpi_21, data = cl_metrics19_wide)
TP_imperv_pd<- lm(TP~pd_21, data = cl_metrics19_wide)
TP_imperv_ed<- lm(TP~ed_21, data = cl_metrics19_wide)

## Total Suspended Solids
TSS_imperv_np<- lm(TSS~np_21, data = cl_metrics19_wide)
TSS_imperv_pland<- lm(TSS~pland_21, data = cl_metrics19_wide)
summary(TSS_imperv_pland)

Call:
lm(formula = TSS ~ pland_21, data = cl_metrics19_wide)

Residuals:
    Min      1Q  Median      3Q     Max 
-33.600 -16.646 -11.934   8.033  72.346 

Coefficients:
            Estimate Std. Error t value Pr(>|t|)   
(Intercept)  16.4163     6.3640   2.580  0.01543 * 
pland_21      0.7280     0.2575   2.827  0.00858 **
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 26.36 on 28 degrees of freedom
  (6 observations deleted due to missingness)
Multiple R-squared:  0.222, Adjusted R-squared:  0.1942 
F-statistic: 7.991 on 1 and 28 DF,  p-value: 0.00858

TSS_imperv_lpi<- lm(TSS~lpi_21, data = cl_metrics19_wide)
summary(TSS_imperv_lpi)

Call:
lm(formula = TSS ~ lpi_21, data = cl_metrics19_wide)

Residuals:
    Min      1Q  Median      3Q     Max 
-33.161 -17.503 -12.282   8.781  72.392 

Coefficients:
            Estimate Std. Error t value Pr(>|t|)   
(Intercept)  18.6737     5.9220   3.153  0.00383 **
lpi_21        0.6981     0.2513   2.778  0.00966 **
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 26.46 on 28 degrees of freedom
  (6 observations deleted due to missingness)
Multiple R-squared:  0.216, Adjusted R-squared:  0.188 
F-statistic: 7.715 on 1 and 28 DF,  p-value: 0.009661

TSS_imperv_pd<- lm(TSS~pd_21, data = cl_metrics19_wide)
TSS_imperv_ed<- lm(TSS~ed_21, data = cl_metrics19_wide)

## Linear Models of Water Quality pollutants and Agriculture Level Landscape metrics 

# Dissolved Oxygen 
DO_agr_np<- lm(DO~np_81, data = cl_metrics19_wide)
DO_agr_pland<- lm(DO~pland_81, data = cl_metrics19_wide)
summary(DO_agr_pland)

Call:
lm(formula = DO ~ pland_81, data = cl_metrics19_wide)

Residuals:
     Min       1Q   Median       3Q      Max 
-2.67605 -0.65212  0.08758  0.76989  2.25047 

Coefficients:
            Estimate Std. Error t value Pr(>|t|)    
(Intercept)  9.15582    0.28778  31.815  < 2e-16 ***
pland_81     0.03421    0.01211   2.825  0.00785 ** 
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 1.229 on 34 degrees of freedom
Multiple R-squared:  0.1901,    Adjusted R-squared:  0.1663 
F-statistic: 7.981 on 1 and 34 DF,  p-value: 0.007853

DO_agr_lpi<- lm(DO~lpi_81, data = cl_metrics19_wide)
summary(DO_agr_lpi)

Call:
lm(formula = DO ~ lpi_81, data = cl_metrics19_wide)

Residuals:
     Min       1Q   Median       3Q      Max 
-2.84658 -0.71494  0.09639  0.94008  2.05013 

Coefficients:
            Estimate Std. Error t value Pr(>|t|)    
(Intercept)  9.31496    0.25646  36.321   <2e-16 ***
lpi_81       0.20896    0.07719   2.707   0.0105 *  
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 1.238 on 34 degrees of freedom
Multiple R-squared:  0.1773,    Adjusted R-squared:  0.1531 
F-statistic: 7.328 on 1 and 34 DF,  p-value: 0.01054

DO_agr_pd<- lm(DO~pd_81, data = cl_metrics19_wide)
DO_agr_ed<- lm(DO~ed_81, data = cl_metrics19_wide)

## Total Nitrogen
TN_agr_np<- lm(TN~np_81, data = cl_metrics19_wide)
TN_agr_pland<- lm(TN~pland_81, data = cl_metrics19_wide)
summary(TN_agr_pland)

Call:
lm(formula = TN ~ pland_81, data = cl_metrics19_wide)

Residuals:
    Min      1Q  Median      3Q     Max 
-1.6168 -0.4603 -0.1458  0.2108  4.6545 

Coefficients:
            Estimate Std. Error t value Pr(>|t|)    
(Intercept)  0.09182    0.30004   0.306 0.761842    
pland_81     0.04723    0.01171   4.033 0.000385 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 1.145 on 28 degrees of freedom
  (6 observations deleted due to missingness)
Multiple R-squared:  0.3674,    Adjusted R-squared:  0.3448 
F-statistic: 16.26 on 1 and 28 DF,  p-value: 0.000385

TN_agr_lpi<- lm(TN~lpi_81, data = cl_metrics19_wide)
summary(TN_agr_lpi)

Call:
lm(formula = TN ~ lpi_81, data = cl_metrics19_wide)

Residuals:
    Min      1Q  Median      3Q     Max 
-2.6477 -0.3836 -0.1134  0.1869  4.7828 

Coefficients:
            Estimate Std. Error t value Pr(>|t|)    
(Intercept)  0.33125    0.26713   1.240 0.225259    
lpi_81       0.28610    0.07391   3.871 0.000593 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 1.162 on 28 degrees of freedom
  (6 observations deleted due to missingness)
Multiple R-squared:  0.3486,    Adjusted R-squared:  0.3253 
F-statistic: 14.98 on 1 and 28 DF,  p-value: 0.0005932

TN_agr_pd<- lm(TN~pd_81, data = cl_metrics19_wide)
TN_agr_ed<- lm(TN~ed_81, data = cl_metrics19_wide)

## Total Phosphorus
TP_agr_np<- lm(TP~np_81, data = cl_metrics19_wide)
TP_agr_pland<- lm(TP~pland_81, data = cl_metrics19_wide)
TP_agr_lpi<- lm(TP~lpi_81, data = cl_metrics19_wide)
TP_agr_pd<- lm(TP~pd_81, data = cl_metrics19_wide)
TP_agr_ed<- lm(TP~ed_81, data = cl_metrics19_wide)

## Total Suspended Solids
TSS_agr_np<- lm(TSS~np_81, data = cl_metrics19_wide)
TSS_agr_pland<- lm(TSS~pland_81, data = cl_metrics19_wide)
TSS_agr_lpi<- lm(TSS~lpi_81, data = cl_metrics19_wide)
TSS_agr_pd<- lm(TSS~pd_81, data = cl_metrics19_wide)
TSS_agr_ed<- lm(TSS~ed_81, data = cl_metrics19_wide)

#create an AIC table to see which model is best

## Dissolved Oxygen tables
DO_F_list <- tibble::lst(DO_F_ed, DO_F_lpi, DO_F_np, DO_F_pd, DO_F_pland) 
DO_agr_list=tibble::lst(DO_agr_ed, DO_agr_lpi, DO_agr_np, DO_agr_pd, DO_agr_pland) 
DO_imperv_list=tibble::lst(DO_imperv_ed, DO_imperv_lpi, DO_imperv_np, DO_imperv_pd, DO_imperv_pland)

DO_all_list=tibble::lst(DO_F_ed, DO_F_lpi, DO_F_np, DO_F_pd, DO_F_pland,
                        DO_agr_ed, DO_agr_lpi, DO_agr_np, DO_agr_pd, DO_agr_pland,
                        DO_imperv_ed, DO_imperv_lpi, DO_imperv_np, DO_imperv_pd, DO_imperv_pland)

DO_all_aictable=aictab(DO_all_list)
DO_all_aictable

Model selection based on AICc:

                K   AICc Delta_AICc AICcWt Cum.Wt     LL
DO_F_lpi        3 120.32       0.00   0.27   0.27 -56.78
DO_imperv_pland 3 121.25       0.93   0.17   0.44 -57.25
DO_agr_pland    3 121.68       1.36   0.14   0.58 -57.47
DO_agr_lpi      3 122.25       1.93   0.10   0.69 -57.75
DO_imperv_lpi   3 122.87       2.55   0.08   0.76 -58.06
DO_F_ed         3 123.63       3.31   0.05   0.82 -58.44
DO_agr_ed       3 123.97       3.65   0.04   0.86 -58.61
DO_F_pd         3 124.46       4.14   0.03   0.89 -58.85
DO_F_np         3 125.01       4.69   0.03   0.92 -59.13
DO_F_pland      3 125.03       4.71   0.03   0.95 -59.14
DO_agr_np       3 126.06       5.74   0.02   0.96 -59.66
DO_agr_pd       3 126.06       5.75   0.02   0.98 -59.66
DO_imperv_pd    3 126.90       6.58   0.01   0.99 -60.07
DO_imperv_ed    3 127.44       7.13   0.01   0.99 -60.35
DO_imperv_np    3 128.25       7.93   0.01   1.00 -60.75

##Total Nitrogen tables
TN_F_list <- tibble::lst(TN_F_ed, TN_agr_lpi, TN_agr_np, TN_agr_pd, TN_agr_pland) 
TN_agr_list=tibble::lst(TN_agr_ed, TN_agr_lpi, TN_agr_np, TN_agr_pd, TN_agr_pland) 
TN_imperv_list=tibble::lst(TN_imperv_ed, TN_imperv_lpi, TN_imperv_np, TN_imperv_pd, TN_imperv_pland) 

##Combining all TN~LSM models metrics into one table
TN_all_list=tibble::lst(TN_F_ed, TN_F_lpi, TN_F_np, TN_F_pd, TN_F_pland,
                        TN_agr_ed, TN_agr_lpi, TN_agr_np, TN_agr_pd, TN_agr_pland,
                        TN_imperv_ed, TN_imperv_lpi, TN_imperv_np, TN_imperv_pd, TN_imperv_pland)


#Create AIC table from list
TN_all_aictable=aictab(TN_all_list)
TN_all_aictable

Model selection based on AICc:

                K   AICc Delta_AICc AICcWt Cum.Wt     LL
TN_agr_pland    3  98.11       0.00   0.57   0.57 -45.59
TN_agr_lpi      3  98.99       0.88   0.37   0.95 -46.03
TN_F_pland      3 103.57       5.46   0.04   0.98 -48.32
TN_agr_ed       3 107.10       8.99   0.01   0.99 -50.09
TN_F_lpi        3 109.06      10.95   0.00   0.99 -51.07
TN_imperv_ed    3 110.00      11.89   0.00   0.99 -51.54
TN_imperv_pd    3 110.07      11.96   0.00   0.99 -51.57
TN_F_ed         3 111.05      12.94   0.00   1.00 -52.06
TN_F_pd         3 111.12      13.01   0.00   1.00 -52.10
TN_agr_np       3 111.12      13.02   0.00   1.00 -52.10
TN_imperv_lpi   3 111.63      13.52   0.00   1.00 -52.35
TN_imperv_pland 3 111.67      13.56   0.00   1.00 -52.37
TN_F_np         3 111.82      13.72   0.00   1.00 -52.45
TN_agr_pd       3 111.82      13.72   0.00   1.00 -52.45
TN_imperv_np    3 111.84      13.74   0.00   1.00 -52.46

##Total Phosphorus tables

TP_F_list <- tibble::lst(TP_F_ed, TP_agr_lpi, TP_agr_np, TP_agr_pd, TP_agr_pland)
TP_agr_list=tibble::lst(TP_agr_ed, TP_agr_lpi, TP_agr_np, TP_agr_pd, TP_agr_pland) 
TP_imperv_list=tibble::lst(TP_imperv_ed, TP_imperv_lpi, TP_imperv_np, TP_imperv_pd, TP_imperv_pland) 

##Combining all TP~LSM models metrics into one table
TP_all_list=tibble::lst(TP_F_ed, TP_F_lpi, TP_F_np, TP_F_pd, TP_F_pland,
                        TP_agr_ed, TP_agr_lpi, TP_agr_np, TP_agr_pd, TP_agr_pland,
                        TP_imperv_ed, TP_imperv_lpi, TP_imperv_np, TP_imperv_pd, TP_imperv_pland)

#Create AIC table from list
TP_all_aictable=aictab(TP_all_list)
TP_all_aictable

Model selection based on AICc:

                K    AICc Delta_AICc AICcWt Cum.Wt    LL
TP_F_pland      3 -106.97       0.00   0.59   0.59 56.95
TP_F_lpi        3 -104.75       2.23   0.20   0.79 55.83
TP_F_pd         3 -101.37       5.61   0.04   0.82 54.15
TP_imperv_pland 3 -100.79       6.18   0.03   0.85 53.86
TP_imperv_ed    3 -100.65       6.32   0.03   0.88 53.79
TP_imperv_lpi   3 -100.44       6.53   0.02   0.90 53.68
TP_F_ed         3 -100.24       6.74   0.02   0.92 53.58
TP_F_np         3  -99.97       7.00   0.02   0.94 53.45
TP_agr_pd       3  -99.04       7.94   0.01   0.95 52.98
TP_agr_pland    3  -99.01       7.96   0.01   0.96 52.97
TP_imperv_pd    3  -98.48       8.49   0.01   0.97 52.70
TP_agr_lpi      3  -98.44       8.53   0.01   0.98 52.68
TP_agr_np       3  -98.38       8.59   0.01   0.99 52.65
TP_agr_ed       3  -98.33       8.65   0.01   0.99 52.63
TP_imperv_np    3  -98.01       8.97   0.01   1.00 52.46

##Total Suspended Solids tables

TSS_F_list <- tibble::lst(TSS_F_ed, TSS_agr_lpi, TSS_agr_np, TSS_agr_pd, TSS_agr_pland)
TSS_agr_list=tibble::lst(TSS_agr_ed, TSS_agr_lpi, TSS_agr_np, TSS_agr_pd, TSS_agr_pland) 
TSS_imperv_list=tibble::lst(TSS_imperv_ed, TSS_imperv_lpi, TSS_imperv_np, TSS_imperv_pd, TSS_imperv_pland)

##Combining all TSS~LSM models metrics into one table
TSS_all_list=tibble::lst(TSS_F_ed, TSS_F_lpi, TSS_F_np, TSS_F_pd, TSS_F_pland,
                        TSS_agr_ed, TSS_agr_lpi, TSS_agr_np, TSS_agr_pd, TSS_agr_pland,
                        TSS_imperv_ed, TSS_imperv_lpi, TSS_imperv_np, TSS_imperv_pd, TSS_imperv_pland)

#Create AIC table from list
TSS_all_aictable=aictab(TSS_all_list)
TSS_all_aictable

Model selection based on AICc:

                 K   AICc Delta_AICc AICcWt Cum.Wt      LL
TSS_imperv_pland 3 286.29       0.00   0.32   0.32 -139.68
TSS_imperv_lpi   3 286.52       0.23   0.28   0.60 -139.80
TSS_F_pd         3 288.30       2.01   0.12   0.71 -140.69
TSS_F_np         3 288.81       2.52   0.09   0.80 -140.94
TSS_imperv_ed    3 289.06       2.77   0.08   0.88 -141.07
TSS_F_pland      3 290.99       4.70   0.03   0.91 -142.03
TSS_F_ed         3 292.14       5.85   0.02   0.93 -142.61
TSS_F_lpi        3 292.68       6.39   0.01   0.94 -142.88
TSS_agr_ed       3 293.03       6.74   0.01   0.95 -143.05
TSS_imperv_pd    3 293.57       7.28   0.01   0.96 -143.33
TSS_imperv_np    3 293.61       7.32   0.01   0.97 -143.34
TSS_agr_pland    3 293.61       7.32   0.01   0.98 -143.34
TSS_agr_lpi      3 293.63       7.34   0.01   0.99 -143.36
TSS_agr_np       3 293.80       7.51   0.01   0.99 -143.44
TSS_agr_pd       3 293.82       7.53   0.01   1.00 -143.45

## Additive Models of delta_AIC values <2 


## Total Suspended Solids
TSS_addit_model<- lm(TSS~pland_21+lpi_21, data = cl_metrics19_wide)
summary(TSS_addit_model)

TSS_addit_model1<- lm(TSS~pland_21+lpi_21+pd_41, data = cl_metrics19_wide)
summary(TSS_addit_model1)
TSS_addit_aictable=aictab(TSS_all_list)

###Neither models produced significant p-values of <.05 when added together


## Total Phosphorus ran only one significant modeled relationship

## Total Nitrogen
TN_addit_model<- lm(TN~pland_81+lpi_81, data = cl_metrics19_wide)
summary(TN_addit_model)

##No significant p-values either

## Dissolved Oxygen

DO_addit_model<- lm(DO~pland_81+lpi_81, data = cl_metrics19_wide)
summary(DO_addit_model)

##No significant p-values here

DO_addit_model3<- lm(DO~lpi_41+pland_81, data = cl_metrics19_wide)
summary(DO_addit_model3)
## Adjusted R-squared:  0.4998

### Here, forested lpi and agricultural pland both display significant p-values, while the others do not

##ggplot of significant results

## Significant results from Dissolved Oxygen models: forest lpi, impervious surface pland, agriculture pland & lpi

ggplot(cl_metrics19_wide, aes(lpi_41, DO))+
  geom_point()+
  geom_smooth(method=lm)+
  ylab("Dissolved Oxygen (mg/L)")+
  xlab("Largest Forest Patch Index")

`geom_smooth()` using formula 'y ~ x'

### Adjusted R-squared:  0.1973

ggplot(cl_metrics19_wide, aes(pland_21, DO))+
  geom_point()+
  geom_smooth(method=lm)+
  ylab("Dissolved Oxygen (mg/L")+
  xlab("Percent Impervious Surface")

`geom_smooth()` using formula 'y ~ x'

##Adjusted R-squared:  0.1762

ggplot(cl_metrics19_wide, aes(pland_81, DO))+
  geom_point()+
  geom_smooth(method=lm)+
  ylab("Dissolved Oxygen (mg/L)")+
  xlab("Percent Agriculture")

`geom_smooth()` using formula 'y ~ x'

## Adjusted R-squared:  0.1663

ggplot(cl_metrics19_wide, aes(lpi_81, DO))+
  geom_point()+
  geom_smooth(method=lm)+
  ylab("Dissolved Oxygen (mg/L)")+
  xlab("Largest Agricultural Patch Index")

`geom_smooth()` using formula 'y ~ x'

## Adjusted R-squared:  0.1531 

## additive model

summary(DO_addit_model3)

Call:
lm(formula = DO ~ lpi_41 + pland_81, data = cl_metrics19_wide)

Residuals:
    Min      1Q  Median      3Q     Max 
-1.8351 -0.3927  0.1058  0.5363  1.7284 

Coefficients:
            Estimate Std. Error t value Pr(>|t|)    
(Intercept) 8.078937   0.314168  25.715  < 2e-16 ***
lpi_41      0.067679   0.013913   4.865 2.75e-05 ***
pland_81    0.044685   0.009624   4.643 5.26e-05 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 0.9517 on 33 degrees of freedom
Multiple R-squared:  0.5283,    Adjusted R-squared:  0.4998 
F-statistic: 18.48 on 2 and 33 DF,  p-value: 4.12e-06

ggplot(cl_metrics19_wide, aes(pland_81, DO, color=lpi_41))+
  geom_point()+
  geom_smooth(method=lm)+
  ylab("Dissolved Oxygen (mg/L)")+
  xlab("Percent Agriculture")

`geom_smooth()` using formula 'y ~ x'

## Where forested patch index is lighest (higher) percent agriculture is also low. Furthermore, higher dissolved oxygen values are when forested patch is at its largest
# Adjusted R-squared:  0.4998

## Significant results from Total Nitrogen models: agriculture pland & lpi

ggplot(cl_metrics19_wide, aes(lpi_81, TN))+
  geom_point()+
  geom_smooth(method=lm)+
  ylab("Total Nitrogen (mg/L)")+
  xlab("Largest Agriculture Patch Index")

`geom_smooth()` using formula 'y ~ x'
Warning: Removed 6 rows containing non-finite values (stat_smooth).
Warning: Removed 6 rows containing missing values (geom_point).

## Adjusted R-squared:  0.3253 

ggplot(cl_metrics19_wide, aes(pland_81, TN))+
  geom_point()+
  geom_smooth(method=lm)+
  ylab("Total Nitrogen (mg/L)")+
  xlab("Percent Agricultural Landscape")

`geom_smooth()` using formula 'y ~ x'
Warning: Removed 6 rows containing non-finite values (stat_smooth).
Warning: Removed 6 rows containing missing values (geom_point).

## Adjusted R-squared:  0.3448 

## Increased percent agriculture shows higher total nitrogen

## Significant results from Total Phosphorus models: forest pland

ggplot(cl_metrics19_wide, aes(pland_41, TP))+
  geom_point()+
  geom_smooth(method=lm)+
  ylab("Total Phosphorus (mg/L)")+
  xlab("Percent Forest")

`geom_smooth()` using formula 'y ~ x'
Warning: Removed 6 rows containing non-finite values (stat_smooth).
Warning: Removed 6 rows containing missing values (geom_point).

## Adjusted R-squared:  0.2328 

## Negative linear relationship, more percent forest lower phosphorus volume

## Significant results from Total Suspended Solids models: impervious surface pland & lpi

ggplot(cl_metrics19_wide, aes(lpi_21, TSS))+
  geom_point()+
  geom_smooth(method=lm)+
  ylab("Total Suspended Solids (mg/L)")+
  xlab("Largest Impervious Surface Patch Index")

`geom_smooth()` using formula 'y ~ x'
Warning: Removed 6 rows containing non-finite values (stat_smooth).
Warning: Removed 6 rows containing missing values (geom_point).

## Adjusted R-squared:  0.188

ggplot(cl_metrics19_wide, aes(pland_21, TSS))+
  geom_point()+
  geom_smooth(method=lm)+
  ylab("Total Suspended Solids (mg/L)")+
  xlab("Percent Impervious Surface")

`geom_smooth()` using formula 'y ~ x'
Warning: Removed 6 rows containing non-finite values (stat_smooth).
Warning: Removed 6 rows containing missing values (geom_point).

## Adjusted R-squared:  0.1942 

### As percent of impervious surface increase, suspended solid volume also increases