Macrophyte Niches in MN
Mike Verhoeven, Wes Glisson, Dan Larkin
2021-01-29
# header ------------------------------------------------------------------Preface
DRUM submission version (uses cleaned up dataset and filepaths, adds footer)
This code will pull in dataset and analyze the niche occupancy of aquatic plants in three dimensions– depth, light availability, phenology–and will produce statistics tests, and visualizations used in the accompanying manuscript.
Load libraries:
# +warning=FALSE, message=FALSE
library(knitr)## Warning: package 'knitr' was built under R version 3.6.3 library(ggplot2) ## Warning: package 'ggplot2' was built under R version 3.6.3 library(ezknitr)
library(data.table)
library(TPD)## Loading required package: ks## Warning: package 'ks' was built under R version 3.6.3 library(gridExtra)
library(ggnewscale)## Warning: package 'ggnewscale' was built under R version 3.6.3 library(GGally)## Warning: package 'GGally' was built under R version 3.6.3## Registered S3 method overwritten by 'GGally':
## method from
## +.gg ggplot2 library(ngram)
library(sf)## Warning: package 'sf' was built under R version 3.6.3## Linking to GEOS 3.8.0, GDAL 3.0.4, PROJ 6.3.1 library(maps)
library(rgdal)## Warning: package 'rgdal' was built under R version 3.6.3## Loading required package: sp## Warning: package 'sp' was built under R version 3.6.3## rgdal: version: 1.5-18, (SVN revision 1082)
## Geospatial Data Abstraction Library extensions to R successfully loaded
## Loaded GDAL runtime: GDAL 3.0.4, released 2020/01/28
## Path to GDAL shared files: C:/Users/verh0064/Documents/R/R-3.6.1/library/rgdal/gdal
## GDAL binary built with GEOS: TRUE
## Loaded PROJ runtime: Rel. 6.3.1, February 10th, 2020, [PJ_VERSION: 631]
## Path to PROJ shared files: C:/Users/verh0064/Documents/R/R-3.6.1/library/rgdal/proj
## Linking to sp version:1.4-4
## To mute warnings of possible GDAL/OSR exportToProj4() degradation,
## use options("rgdal_show_exportToProj4_warnings"="none") before loading rgdal. library(ggsn)## Loading required package: grid library(moments)
# load in data ------------------------------------------------------------Refresh data folder:
# run this line to update your local copy of project data folder files
# source("scripts/data_folder_pull.R")Load in Data:
#plant surveys
plants <- fread(file = "plants_DRUM.csv", drop = 1)
# explore dataset ---------------------------------------------------------
str(plants)## Classes 'data.table' and 'data.frame': 724228 obs. of 41 variables:
## $ DOWLKNUM : chr "1000100" "1000100" "1000100" "1000100" ...
## $ SURVEY_DATE : chr "6/25/2008" "6/25/2008" "6/25/2008" "6/25/2008" ...
## $ LAKE_NAME : chr "Pine" "Pine" "Pine" "Pine" ...
## $ DATASOURCE : chr "DNR Fisheries" "DNR Fisheries" "DNR Fisheries" "DNR Fisheries" ...
## $ STA_NBR : chr "1" "1" "1" "1" ...
## $ DEPTH_FT : num 2.3 2.3 2.3 2.3 2.3 ...
## $ SURVEYOR : chr "DNR Fisheries" "DNR Fisheries" "DNR Fisheries" "DNR Fisheries" ...
## $ NO_VEG_FOUND : logi FALSE FALSE FALSE FALSE FALSE FALSE ...
## $ TAXON : chr "Equisetum" NA "Lemna trisulca" "Nuphar" ...
## $ POINT_LVL_SECCHI : num NA NA NA NA NA NA NA NA NA NA ...
## $ SURVEY_ID.x : int 282 282 282 282 282 282 282 282 282 282 ...
## $ POINT_ID : int 19107 19107 19107 19107 19107 19107 19107 19107 19107 19108 ...
## $ OBS_ID : int 40234 40235 40236 40237 40238 40239 40240 40241 40242 40243 ...
## $ YEAR.SURVEY : int 2008 2008 2008 2008 2008 2008 2008 2008 2008 2008 ...
## $ MONTH.SURVEY : int 6 6 6 6 6 6 6 6 6 6 ...
## $ SURVEY_ID.y : int 282 282 282 282 282 282 282 282 282 282 ...
## $ Secchi_m.mean : num 2.47 2.47 2.47 2.47 2.47 ...
## $ YEAR.SECCHI.mean : num 2008 2008 2008 2008 2008 ...
## $ MONTH.SECCHI.mean : num 7.91 7.91 7.91 7.91 7.91 ...
## $ Secchi_m.min : num 1.22 1.22 1.22 1.22 1.22 1.22 1.22 1.22 1.22 1.22 ...
## $ YEAR.SECCHI.min : int 2007 2007 2007 2007 2007 2007 2007 2007 2007 2007 ...
## $ MONTH.SECCHI.min : int 7 7 7 7 7 7 7 7 7 7 ...
## $ Secchi_m.max : num 3.81 3.81 3.81 3.81 3.81 3.81 3.81 3.81 3.81 3.81 ...
## $ YEAR.SECCHI.max : int 2009 2009 2009 2009 2009 2009 2009 2009 2009 2009 ...
## $ MONTH.SECCHI.max : int 9 9 9 9 9 9 9 9 9 9 ...
## $ Secchi_m.sd : num 0.727 0.727 0.727 0.727 0.727 ...
## $ Secchi_m.length : int 34 34 34 34 34 34 34 34 34 34 ...
## $ Secchi_m.se : num 0.853 0.853 0.853 0.853 0.853 ...
## $ chillR_code : chr "727450_14913" "727450_14913" "727450_14913" "727450_14913" ...
## $ STATION.NAME : chr "DULUTH INTERNATIONAL AIRPORT" "DULUTH INTERNATIONAL AIRPORT" "DULUTH INTERNATIONAL AIRPORT" "DULUTH INTERNATIONAL AIRPORT" ...
## $ gdd : num 493 493 493 493 493 ...
## $ center_utm : num 494097 494097 494097 494097 494097 ...
## $ center_u_1 : num 5114741 5114741 5114741 5114741 5114741 ...
## $ long : num -93.1 -93.1 -93.1 -93.1 -93.1 ...
## $ lat : num 46.2 46.2 46.2 46.2 46.2 ...
## $ Lat_WX : num 46.8 46.8 46.8 46.8 46.8 ...
## $ Long_WX : num -92.2 -92.2 -92.2 -92.2 -92.2 ...
## $ distance : num 99.7 99.7 99.7 99.7 99.7 ...
## $ Overlap_years : int 19 19 19 19 19 19 19 19 19 19 ...
## $ Perc_interval_covered: int 100 100 100 100 100 100 100 100 100 100 ...
## $ proplight : num 0.521 0.521 0.521 0.521 0.521 ...
## - attr(*, ".internal.selfref")=<externalptr>Data Exploration and Prep
Ensure complete on all vars we’ll be using.
#whats in the dataset?
plants[ , .N , POINT_ID ] #353,148 sampled locations## POINT_ID N
## 1: 19107 9
## 2: 19108 2
## 3: 19109 2
## 4: 19110 1
## 5: 19111 6
## ---
## 353144: 353248 5
## 353145: 353249 1
## 353146: 353250 3
## 353147: 353251 4
## 353148: 353252 1 plants[, .N ,SURVEY_ID.x ] # Surveys## SURVEY_ID.x N
## 1: 282 334
## 2: 1 63
## 3: 2 126
## 4: 3 160
## 5: 4 147
## ---
## 3400: 3400 31
## 3401: 3401 109
## 3402: 3402 175
## 3403: 3403 328
## 3404: 3404 8 plants[, .N ,SURVEY_ID.y ] # Surveys with covariates## SURVEY_ID.y N
## 1: 282 334
## 2: 1 63
## 3: 2 126
## 4: 3 160
## 5: NA 125325
## ---
## 2580: 3386 55
## 2581: 3388 111
## 2582: 3395 336
## 2583: 3398 20
## 2584: 3400 31 plants[ , .N , .(DOWLKNUM, YEAR.SURVEY) ] #Lake-years## DOWLKNUM YEAR.SURVEY N
## 1: 1000100 2008 334
## 2: 1001600 2011 63
## 3: 1003400 2011 126
## 4: 1003500 2014 160
## 5: 1005300 NA 147
## ---
## 2510: W0655001 2005 31
## 2511: W0844001 NA 109
## 2512: W0854002 NA 175
## 2513: W0889002 NA 328
## 2514: W0900701 NA 8 plants[ , .N, DOWLKNUM]# Lakes## DOWLKNUM N
## 1: 1000100 334
## 2: 1001600 63
## 3: 1003400 126
## 4: 1003500 160
## 5: 1005300 147
## ---
## 1522: W0655001 31
## 1523: W0844001 109
## 1524: W0854002 175
## 1525: W0889002 328
## 1526: W0900701 8 plants[ , .N, POINT_ID]# Survey points## POINT_ID N
## 1: 19107 9
## 2: 19108 2
## 3: 19109 2
## 4: 19110 1
## 5: 19111 6
## ---
## 353144: 353248 5
## 353145: 353249 1
## 353146: 353250 3
## 353147: 353251 4
## 353148: 353252 1 plants[!is.na(TAXON) , .N, OBS_ID]## OBS_ID N
## 1: 40234 1
## 2: 40236 1
## 3: 40237 1
## 4: 40238 1
## 5: 40239 1
## ---
## 564034: 723011 1
## 564035: 723012 1
## 564036: 723013 1
## 564037: 723014 1
## 564038: 723015 1 surveypoints <- plants[ , length(unique(POINT_ID)), SURVEY_ID.x][ , V1,] #survey points per taxon
hist(surveypoints, breaks = 250, xlim = c(0,500))
mean(surveypoints)## [1] 103.745 sd(surveypoints)## [1] 107.8648 skewness(surveypoints)## [1] 3.044416 kurtosis(surveypoints)## [1] 17.41512 observations <- plants[ , length(unique(OBS_ID)), POINT_ID][ , V1,] #survey points per taxon
hist(observations)
mean(observations)## [1] 2.047343 sd(observations)## [1] 1.615402 skewness(observations)## [1] 2.154917 kurtosis(observations)## [1] 9.590638 #survey point per survey
#retain only locs for which we have all vars to estimate on
plants[is.na(DEPTH_FT) == F &
is.na(gdd) == F &
is.na(proplight) == F, .N, POINT_ID] # of sampled locs, 292,824 have all data## POINT_ID N
## 1: 19107 9
## 2: 19108 2
## 3: 19109 2
## 4: 19110 1
## 5: 19111 6
## ---
## 292820: 353130 4
## 292821: 353131 4
## 292822: 353132 4
## 292823: 353133 5
## 292824: 353134 6 #percent coverage by sample location for complete cases
nrow(plants[is.na(DEPTH_FT) == F &
is.na(gdd) == F &
is.na(proplight) == F, .N, POINT_ID])/nrow(plants[ , .N , POINT_ID ])## [1] 0.8291821 #percent coverage by species observationssample location for
plants[is.na(DEPTH_FT) == F &
is.na(gdd) == F &
is.na(proplight) == F, .N,]/plants[ , .N ,]## [1] 0.8249087Include only the complete cases
#retain only complete cases
plants <- plants[is.na(DEPTH_FT) == F &
is.na(gdd) == F &
is.na(proplight) == F, ,]
plants[ , .N, POINT_ID]# Survey points## POINT_ID N
## 1: 19107 9
## 2: 19108 2
## 3: 19109 2
## 4: 19110 1
## 5: 19111 6
## ---
## 292820: 353130 4
## 292821: 353131 4
## 292822: 353132 4
## 292823: 353133 5
## 292824: 353134 6 plants[!is.na(TAXON) , .N, OBS_ID]# plant obs## OBS_ID N
## 1: 40234 1
## 2: 40236 1
## 3: 40237 1
## 4: 40238 1
## 5: 40239 1
## ---
## 462114: 722391 1
## 462115: 722392 1
## 462116: 722393 1
## 462117: 722394 1
## 462118: 722395 1WX distance from sample lakes
# grab each lake's distance (take mean w/in lake)
WXlk_dist <- plants[!is.na(distance) , mean(distance) , DOWLKNUM] # distance of lakes from WX
# calc distance stats across lakes
WXlk_dist[ ,c(meandist=mean(V1), maxdist = max(V1), mindist = min(V1), sddist = sd(V1)), ,]## meandist maxdist mindist sddist
## 64.00260 171.02000 1.49000 30.65851Survey Effort
Check out the data. We want to have a look at the overall distribution of survey effort, and all of the niche parameters we hope to estimate on
# locations sampled (done by including POINT_ID in "by")
# summary(plants)
locssamp <- plants[ , .(Secchi = mean(Secchi_m.mean) , DEPTH_FT = mean(DEPTH_FT), proplight = mean(proplight), gdd = mean(gdd), year = mean(YEAR.SURVEY)) , .(DOWLKNUM, DATASOURCE, SURVEY_DATE, POINT_ID, SURVEY_ID.x) ]
hist(locssamp$proplight)
hist(locssamp$DEPTH_FT, breaks = seq(0,1000,1), xlim = c(0,30))
hist(locssamp$gdd, breaks = seq(0,3500,100))
#sample effort by lakes, years, lake year, and surveys
locssamp[ , .N, .(DOWLKNUM, year) ]## DOWLKNUM year N
## 1: 1000100 2008 86
## 2: 1001600 2011 41
## 3: 1003400 2011 66
## 4: 1003500 2014 37
## 5: 1006000 2009 52
## ---
## 1932: 87010200 2010 28
## 1933: 87010200 2013 40
## 1934: 87011600 2013 86
## 1935: R1961 2005 35
## 1936: W0655001 2005 7 locssamp[ , .N, .(DOWLKNUM) ]## DOWLKNUM N
## 1: 1000100 86
## 2: 1001600 41
## 3: 1003400 66
## 4: 1003500 37
## 5: 1006000 52
## ---
## 1106: 87006000 93
## 1107: 87010200 97
## 1108: 87011600 86
## 1109: R1961 35
## 1110: W0655001 7 locssamp[ , .N, .(year) ]## year N
## 1: 2008 33398
## 2: 2011 32607
## 3: 2014 21797
## 4: 2009 34987
## 5: 2002 4040
## 6: 2012 25432
## 7: 2010 34827
## 8: 2013 25084
## 9: 2015 19586
## 10: 2003 6289
## 11: 2006 16701
## 12: 2016 7691
## 13: 2017 6692
## 14: 2018 5449
## 15: 2005 9878
## 16: 2004 7174
## 17: 2007 14929
## 18: 2001 273 locssamp[ , .N, SURVEY_ID.x]## SURVEY_ID.x N
## 1: 282 86
## 2: 1 41
## 3: 2 66
## 4: 3 37
## 5: 5 52
## ---
## 2574: 3385 28
## 2575: 3386 40
## 2576: 3388 86
## 2577: 3395 35
## 2578: 3400 7PDF of sampling:
What does the probability distribution fn of our sampling look like?
#kernel density plots
kde_depths <- density(locssamp$DEPTH_FT)
plot(kde_depths, main = "TPD depth")
plot(kde_depths, main = "TPD depth", xlim = c(0,30))
kde_light <- density(locssamp$proplight)
plot(kde_light, main = "TPD light")
kde_gdd <- density(locssamp$gdd)
plot(kde_gdd, main = "TPD gdd")
Niche Axes Orthogonality
What about covariation in our niche axes:
# plot covariance and correllation coefs (r) from the niche axes (uncomment for visulaization thats computationally expensive)
ggpairs(locssamp[ ,.(proplight, DEPTH_FT, gdd) ,], lower = list(continuous = "smooth"), ggplot2::aes(fill = NULL))
cor.test(locssamp[ ,c(proplight) ,], locssamp[ ,c(DEPTH_FT) ,])##
## Pearson's product-moment correlation
##
## data: locssamp[, c(proplight), ] and locssamp[, c(DEPTH_FT), ]
## t = -253.32, df = 306832, p-value < 2.2e-16
## alternative hypothesis: true correlation is not equal to 0
## 95 percent confidence interval:
## -0.4188133 -0.4129607
## sample estimates:
## cor
## -0.4158913 cor.test(locssamp[ ,c(gdd) ,], locssamp[ ,c(proplight) ,])##
## Pearson's product-moment correlation
##
## data: locssamp[, c(gdd), ] and locssamp[, c(proplight), ]
## t = -6.871, df = 306832, p-value = 6.386e-12
## alternative hypothesis: true correlation is not equal to 0
## 95 percent confidence interval:
## -0.015940931 -0.008865387
## sample estimates:
## cor
## -0.01240331 cor.test(locssamp[ ,c(gdd) ,], locssamp[ ,c(DEPTH_FT) ,])##
## Pearson's product-moment correlation
##
## data: locssamp[, c(gdd), ] and locssamp[, c(DEPTH_FT), ]
## t = -26.528, df = 306832, p-value < 2.2e-16
## alternative hypothesis: true correlation is not equal to 0
## 95 percent confidence interval:
## -0.05136642 -0.04430598
## sample estimates:
## cor
## -0.0478368Let’s take a better look at those light-depth relationships:
ggplot(locssamp, aes(DEPTH_FT,proplight))+
geom_point(alpha = .01)+
xlim(c(0,23))## Warning: Removed 8510 rows containing missing values (geom_point).
hist(locssamp$Secchi)
summary(locssamp$Secchi)## Min. 1st Qu. Median Mean 3rd Qu. Max.
## 0.03048 0.83278 1.37333 1.81496 2.60588 6.76667We can see from this that we do have some non-orthogonality in the light- depth relationship (Correlation Coef 0.416). It doesn’t seem too high to me, but by the pearson correllation we can reject the null that there’s no correlation. In additon, we can see that the light-depth relationship is well represented considering that light exponentially decays with depth.Thus, the variation we see w/ light at depth is due to variation in water clarity– which we consider good because it provides us with a range of light vals for any given depth. I don’t see guidance in the lit about bounds on reasonable covariance for this application, so we are going to run with it.
Plant Occurrences
We want to keep only the species observations for niche models. Clean up the Taxon == NA. These are cases where no veg was found.
plants[is.na(TAXON) == T , .N ,] # no veg found cases## [1] 135272 plants[is.na(TAXON) == F , .N ,] # plant observations## [1] 462150 plantocc <- plants[is.na(TAXON) == F , ,]Are the plants occurring at reasonable values of light avialability? Ususally consider submerged plants to be adapted to 2-10% of incoming light.
plantocc[proplight < .10 & proplight > 0.02, .N , ]/plantocc[,.N,]## [1] 0.2277075 plantocc[proplight > 0.02, .N , ]/plantocc[,.N,]## [1] 0.8520697 hist(plantocc$proplight)
Depth outliers
Need to trim some of the depth outliers. Let’s have a look at what those outliers look like.
#explore outliers
plantocc[, summary(DEPTH_FT) , ]## Min. 1st Qu. Median Mean 3rd Qu. Max.
## 0.000 3.000 4.265 5.235 6.800 230.000 hist(plantocc$DEPTH_FT, breaks = 200)
hist(plantocc$DEPTH_FT, xlim = c(0,25), breaks = 200)
boxplot(plantocc$DEPTH_FT, ylim = c(0,75))
sd(plantocc$DEPTH_FT)## [1] 3.542636Lets call everything > 5 Standard Deviations form the mean an outlier for depth (Muthukrishnan et al 2018)
outlier_cut <- mean(plantocc$DEPTH_FT)+5*sd(plantocc$DEPTH_FT) #conveniently this corresponds to 6.998m
#who gets cut?
plantocc[DEPTH_FT > outlier_cut |
DEPTH_FT==0 , .N,]## [1] 1844 plantocc[DEPTH_FT > outlier_cut |
DEPTH_FT==0, c(1:15) ]## DOWLKNUM SURVEY_DATE LAKE_NAME DATASOURCE STA_NBR DEPTH_FT
## 1: 1000100 6/25/2008 Pine DNR Fisheries 35 0.0000
## 2: 1000100 6/25/2008 Pine DNR Fisheries 35 0.0000
## 3: 1000100 6/25/2008 Pine DNR Fisheries 35 0.0000
## 4: 1000100 6/25/2008 Pine DNR Fisheries 35 0.0000
## 5: 1000100 6/25/2008 Pine DNR Fisheries 35 0.0000
## ---
## 1840: 86023400 8/8/2016 bass James Johnson 62 23.2951
## 1841: 86023400 8/8/2016 bass James Johnson 69 23.9513
## 1842: 86023400 7/12/2018 bass James Johnson 65 25.2637
## 1843: 86023400 7/12/2018 bass James Johnson 67 24.9356
## 1844: 86023400 7/12/2018 bass James Johnson 68 24.9356
## SURVEYOR NO_VEG_FOUND TAXON POINT_LVL_SECCHI
## 1: DNR Fisheries FALSE Elodea canadensis NA
## 2: DNR Fisheries FALSE Lemna minor NA
## 3: DNR Fisheries FALSE Nuphar NA
## 4: DNR Fisheries FALSE Potamogeton crispus NA
## 5: DNR Fisheries FALSE Ranunculus NA
## ---
## 1840: FALSE Nitella NA
## 1841: FALSE Nitella NA
## 1842: FALSE Nitella NA
## 1843: FALSE Nitella NA
## 1844: FALSE Nitella NA
## SURVEY_ID.x POINT_ID OBS_ID YEAR.SURVEY MONTH.SURVEY
## 1: 282 19162 40438 2008 6
## 2: 282 19162 40440 2008 6
## 3: 282 19162 40441 2008 6
## 4: 282 19162 40442 2008 6
## 5: 282 19162 40443 2008 6
## ---
## 1840: 3342 348238 712002 2016 8
## 1841: 3342 348245 712009 2016 8
## 1842: 3343 348369 712367 2018 7
## 1843: 3343 348371 712369 2018 7
## 1844: 3343 348372 712370 2018 7 plantocc[ , .N,]## [1] 462150 #proportion
plantocc[DEPTH_FT > outlier_cut |
DEPTH_FT==0 , .N,]/plantocc[ , .N,]## [1] 0.003990047 plantocc <- plantocc[!DEPTH_FT > outlier_cut , ,]
plantocc <- plantocc[!DEPTH_FT==0]
# occurrence data by species ------------------------------------------------------Descriptive Stats by Species
Calculate H Index Values for these species
Nsurveys, Noccurrance, H-index by species
hind <- plants[ , .N, .(SURVEY_DATE, DATASOURCE, DOWLKNUM, TAXON) ] #n obs per species per survey
hind$TAXON <- as.factor(hind$TAXON)
setorder(hind, TAXON, -N)
hind[, index := 1:.N, by=TAXON]
hind[N > index , hindex:= .N, by = TAXON]
hind[ , nsurv := .N , by = TAXON]
hind[is.na(hindex) == T , hindex:= 0,]
hind[ , nocc:= sum(N) , by = .(TAXON)]
hind2 <- hind[ , .(hindex=max(hindex)) , by = .(TAXON,nsurv,nocc)]
setorder(hind2, -nocc)
#label taxonomic level
hind2[ , TAXON := as.character(TAXON), ]
# strips brackets from non-species taxa
hind2[, taxacl := gsub("\\s*\\([^\\)]+\\)","",TAXON),]
hind2$taxawords <- sapply(hind2$taxacl, wordcount)
hind2[taxawords>1, species := T,]
hind2[is.na(species), species:= F,]
hind2[ , .N , species]## species N
## 1: FALSE 56
## 2: TRUE 163 hind2[ nocc>10000, .N , species]## species N
## 1: FALSE 2
## 2: TRUE 10 # check out results:
hind2[1:75]## TAXON nsurv nocc hindex
## 1: <NA> 2235 135272 161
## 2: Ceratophyllum demersum 2074 74849 114
## 3: Chara 1443 40029 89
## 4: Potamogeton crispus 1127 35579 92
## 5: Potamogeton zosteriformis 1225 26358 65
## 6: Stuckenia pectinata 1636 22551 57
## 7: Elodea canadensis 1149 22152 63
## 8: Myriophyllum sibiricum 1174 21497 66
## 9: Najas flexilis 957 15298 54
## 10: Myriophyllum spicatum 493 15030 65
## 11: Lemna trisulca 646 13979 63
## 12: Nymphaea odorata 1105 13822 40
## 13: Vallisneria americana 715 9685 44
## 14: Potamogeton robbinsii 318 8693 52
## 15: Utricularia vulgaris 743 8434 38
## 16: Zizania palustris 435 7720 41
## 17: Potamogeton praelongus 709 7707 38
## 18: Potamogeton richardsonii 774 7249 35
## 19: Nuphar variegata 856 6890 29
## 20: Potamogeton amplifolius 560 5551 32
## 21: Najas guadalupensis 188 5442 41
## 22: Potamogeton (narrow) 570 5335 32
## 23: Najas 400 5302 37
## 24: Heteranthera dubia 646 5155 28
## 25: Potamogeton illinoensis 552 5100 26
## 26: Potamogeton friesii 479 4367 31
## 27: Lemna minor 426 3823 30
## 28: Schoenoplectus acutus 416 3823 33
## 29: Potamogeton pusillus 310 3269 29
## 30: Potamogeton natans 556 3212 21
## 31: Schoenoplectus 429 3188 25
## 32: Potamogeton 265 2801 24
## 33: Potamogeton gramineus 371 2688 22
## 34: Myriophyllum 42 2601 16
## 35: Typha angustifolia 463 2111 17
## 36: Potamogeton foliosus 216 2102 25
## 37: Nymphaea 140 2083 27
## 38: Bidens beckii 335 2046 19
## 39: Elodea 132 2046 21
## 40: Brasenia schreberi 234 1996 22
## 41: Spirodela polyrhiza 296 1894 22
## 42: Nitella 193 1709 19
## 43: Utricularia 102 1692 22
## 44: Typha 348 1442 14
## 45: Characeae 39 1403 21
## 46: Drepanocladus 251 1403 16
## 47: Ranunculus 242 1232 17
## 48: Sagittaria 380 1189 13
## 49: Nuphar 140 1186 20
## 50: Ranunculus aquatilis 134 1127 15
## 51: Zannichellia palustris 167 1091 16
## 52: Wolffia 79 957 17
## 53: Wolffia columbiana 81 775 13
## 54: Scirpus 87 721 14
## 55: Lemna 115 692 14
## 56: Sparganium (floating) 76 636 15
## 57: Eleocharis 186 605 10
## 58: Eleocharis acicularis 186 577 9
## 59: Najas marina 62 540 14
## 60: Sparganium 122 489 11
## 61: Zizania 43 476 15
## 62: Nuphar advena 66 464 12
## 63: Utricularia minor 113 460 10
## 64: Carex 145 405 7
## 65: Schoenoplectus subterminalis 76 397 9
## 66: Schoenoplectus tabernaemontani 128 379 7
## 67: Utricularia intermedia 70 356 10
## 68: Pontederia cordata 42 346 8
## 69: Ranunculus longirostris 47 292 8
## 70: Potamogeton nodosus 70 236 7
## 71: Equisetum 78 230 6
## 72: Potamogeton strictifolius 19 209 8
## 73: Typha latifolia 76 197 6
## 74: Nelumbo lutea 33 192 9
## 75: Potamogeton epihydrus 46 189 7
## TAXON nsurv nocc hindex
## taxacl taxawords species
## 1: <NA> 1 FALSE
## 2: Ceratophyllum demersum 2 TRUE
## 3: Chara 1 FALSE
## 4: Potamogeton crispus 2 TRUE
## 5: Potamogeton zosteriformis 2 TRUE
## 6: Stuckenia pectinata 2 TRUE
## 7: Elodea canadensis 2 TRUE
## 8: Myriophyllum sibiricum 2 TRUE
## 9: Najas flexilis 2 TRUE
## 10: Myriophyllum spicatum 2 TRUE
## 11: Lemna trisulca 2 TRUE
## 12: Nymphaea odorata 2 TRUE
## 13: Vallisneria americana 2 TRUE
## 14: Potamogeton robbinsii 2 TRUE
## 15: Utricularia vulgaris 2 TRUE
## 16: Zizania palustris 2 TRUE
## 17: Potamogeton praelongus 2 TRUE
## 18: Potamogeton richardsonii 2 TRUE
## 19: Nuphar variegata 2 TRUE
## 20: Potamogeton amplifolius 2 TRUE
## 21: Najas guadalupensis 2 TRUE
## 22: Potamogeton 1 FALSE
## 23: Najas 1 FALSE
## 24: Heteranthera dubia 2 TRUE
## 25: Potamogeton illinoensis 2 TRUE
## 26: Potamogeton friesii 2 TRUE
## 27: Lemna minor 2 TRUE
## 28: Schoenoplectus acutus 2 TRUE
## 29: Potamogeton pusillus 2 TRUE
## 30: Potamogeton natans 2 TRUE
## 31: Schoenoplectus 1 FALSE
## 32: Potamogeton 1 FALSE
## 33: Potamogeton gramineus 2 TRUE
## 34: Myriophyllum 1 FALSE
## 35: Typha angustifolia 2 TRUE
## 36: Potamogeton foliosus 2 TRUE
## 37: Nymphaea 1 FALSE
## 38: Bidens beckii 2 TRUE
## 39: Elodea 1 FALSE
## 40: Brasenia schreberi 2 TRUE
## 41: Spirodela polyrhiza 2 TRUE
## 42: Nitella 1 FALSE
## 43: Utricularia 1 FALSE
## 44: Typha 1 FALSE
## 45: Characeae 1 FALSE
## 46: Drepanocladus 1 FALSE
## 47: Ranunculus 1 FALSE
## 48: Sagittaria 1 FALSE
## 49: Nuphar 1 FALSE
## 50: Ranunculus aquatilis 2 TRUE
## 51: Zannichellia palustris 2 TRUE
## 52: Wolffia 1 FALSE
## 53: Wolffia columbiana 2 TRUE
## 54: Scirpus 1 FALSE
## 55: Lemna 1 FALSE
## 56: Sparganium 1 FALSE
## 57: Eleocharis 1 FALSE
## 58: Eleocharis acicularis 2 TRUE
## 59: Najas marina 2 TRUE
## 60: Sparganium 1 FALSE
## 61: Zizania 1 FALSE
## 62: Nuphar advena 2 TRUE
## 63: Utricularia minor 2 TRUE
## 64: Carex 1 FALSE
## 65: Schoenoplectus subterminalis 2 TRUE
## 66: Schoenoplectus tabernaemontani 2 TRUE
## 67: Utricularia intermedia 2 TRUE
## 68: Pontederia cordata 2 TRUE
## 69: Ranunculus longirostris 2 TRUE
## 70: Potamogeton nodosus 2 TRUE
## 71: Equisetum 1 FALSE
## 72: Potamogeton strictifolius 2 TRUE
## 73: Typha latifolia 2 TRUE
## 74: Nelumbo lutea 2 TRUE
## 75: Potamogeton epihydrus 2 TRUE
## taxacl taxawords species ggplot(hind,aes(hindex))+
geom_density()
ggplot(hind2, aes(hindex))+
geom_histogram()## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
plot(hind2$nsurv, hind2$hindex)
plot(hind2$nocc, hind2$hindex)
plot(hind2$nocc, hind2$nsurv)
Threshold for inclusion
Then we simply choose a cutoff value and press forward. We will use 1000 occurrences, and only species level IDs
# Set cutoff values for inclusion -----------------------------------------
# drop at some threshold
ptrm <- plantocc[ TAXON %in% hind2[nocc > 1000 &
species==T, TAXON, ], , ]
ptrm[,.N , TAXON]# by taxon, n obs## TAXON N
## 1: Lemna trisulca 13938
## 2: Potamogeton zosteriformis 26272
## 3: Schoenoplectus acutus 3820
## 4: Spirodela polyrhiza 1893
## 5: Elodea canadensis 22073
## 6: Ceratophyllum demersum 74653
## 7: Nuphar variegata 6876
## 8: Potamogeton amplifolius 5532
## 9: Potamogeton crispus 35330
## 10: Potamogeton richardsonii 7196
## 11: Myriophyllum sibiricum 21461
## 12: Stuckenia pectinata 22532
## 13: Lemna minor 3822
## 14: Potamogeton natans 3207
## 15: Vallisneria americana 9667
## 16: Potamogeton robbinsii 8665
## 17: Nymphaea odorata 13788
## 18: Typha angustifolia 2108
## 19: Zizania palustris 7713
## 20: Potamogeton praelongus 7681
## 21: Brasenia schreberi 1994
## 22: Utricularia vulgaris 8426
## 23: Najas flexilis 15285
## 24: Potamogeton illinoensis 5098
## 25: Potamogeton gramineus 2685
## 26: Bidens beckii 2023
## 27: Potamogeton friesii 4362
## 28: Heteranthera dubia 5021
## 29: Potamogeton pusillus 3267
## 30: Najas guadalupensis 5426
## 31: Myriophyllum spicatum 14934
## 32: Ranunculus aquatilis 1117
## 33: Zannichellia palustris 1091
## 34: Potamogeton foliosus 2097
## TAXON N ptrm[,.N ,]# total n obs## [1] 371053# TPD ----------------------------------------------------------------Modeling Niches
Now that we have trimmed our dataset to 1) only locs where plants occurred, 2) excluding observations that are depth outliers, 3) only taxa with more than 1000 occurrences, we can move on to modeling the species’ niches.
Total niche model
# calculating probability densities for multiple species
TPDs_plants <- TPDs(species = ptrm$TAXON, traits = ptrm[,.(scale(DEPTH_FT),scale(proplight),scale(gdd))] )## -------Calculating densities for One population_Multiple species-----------## Be careful, the integral of the pdf of some cases differ from 1.
## They have been reescaled, but you should consider increasing
## 'trait_ranges'
## Ceratophyllum demersum / Nuphar variegata / Potamogeton crispus / Stuckenia pectinata / Lemna minor / Potamogeton robbinsii / Nymphaea odorata / Zizania palustris / Najas guadalupensisCheckout evaluation grid. Each axis was divided into 50 intervals, creating a 3-dimensional matrix (50^3 = 125k cells) with a crossed set of those 50 part axes. Each species’ TPD ought to sum to one. Lets have a look at those things.
head(TPDs_plants$data$evaluation_grid) #first chunk of matrix only varies "V1"## V1 V1.1 V1.2
## 1 -2.595960 -1.814714 -2.856738
## 2 -2.413409 -1.814714 -2.856738
## 3 -2.230857 -1.814714 -2.856738
## 4 -2.048305 -1.814714 -2.856738
## 5 -1.865754 -1.814714 -2.856738
## 6 -1.683202 -1.814714 -2.856738 nrow(TPDs_plants$data$evaluation_grid) #size looks right## [1] 125000 names(TPDs_plants$TPDs) #Estimated for all plants## [1] "Lemna trisulca" "Potamogeton zosteriformis"
## [3] "Schoenoplectus acutus" "Spirodela polyrhiza"
## [5] "Elodea canadensis" "Ceratophyllum demersum"
## [7] "Nuphar variegata" "Potamogeton amplifolius"
## [9] "Potamogeton crispus" "Potamogeton richardsonii"
## [11] "Myriophyllum sibiricum" "Stuckenia pectinata"
## [13] "Lemna minor" "Potamogeton natans"
## [15] "Vallisneria americana" "Potamogeton robbinsii"
## [17] "Nymphaea odorata" "Typha angustifolia"
## [19] "Zizania palustris" "Potamogeton praelongus"
## [21] "Brasenia schreberi" "Utricularia vulgaris"
## [23] "Najas flexilis" "Potamogeton illinoensis"
## [25] "Potamogeton gramineus" "Bidens beckii"
## [27] "Potamogeton friesii" "Heteranthera dubia"
## [29] "Potamogeton pusillus" "Najas guadalupensis"
## [31] "Myriophyllum spicatum" "Ranunculus aquatilis"
## [33] "Zannichellia palustris" "Potamogeton foliosus" sum(TPDs_plants$TPDs$'Zanichellia palustris') ## [1] 0 # Do any not sum to one?
sapply(TPDs_plants$TPDs, sum)[sapply(TPDs_plants$TPDs, sum) != 1 ] #NOPE## Zannichellia palustris
## 1# full TPD dissimilarity table --------------------------------------------Niche Dissimilarity
#dissimilarity across in 3 dimensions
dissim_TPDs_plants <- dissim(TPDs_plants)## Calculating dissimilarities between 34 populations. It might take some time inv_dissim <- setDT(as.data.frame(dissim_TPDs_plants$populations$dissimilarity[ , c('Myriophyllum spicatum', 'Potamogeton crispus')]), keep.rownames = "Taxon")
inv_dissim[ , clpdis := `Potamogeton crispus`-`Myriophyllum spicatum` ,]
# relative contribution of abundance in shared space vs. spaces only occupied by that species
inv_pshared <- setDT(as.data.frame(dissim_TPDs_plants$populations$P_shared[ , c('Myriophyllum spicatum', 'Potamogeton crispus')]), keep.rownames = "Taxon")
inv_pnonshared <- setDT(as.data.frame(dissim_TPDs_plants$populations$P_non_shared[ , c('Myriophyllum spicatum', 'Potamogeton crispus')]), keep.rownames = "Taxon")
invsumm <- cbind(cbind(inv_dissim,inv_pnonshared[,2:3]),inv_pshared[, 2:3])
setnames(invsumm, c("Taxon", "MspicatumDissim", "PcrispusDissim", "PClessMS", "MspicatumNonShared", "PcrispusNonShared", "MspicatumShared", "PcrispusShared"))Test for proportional contrib of non-shared in CLP v EWM:
# nonshared prop of niche differentiation in CLP exceeds EWM (# of native species comparisons)
summary(0>invsumm[ , MspicatumNonShared-PcrispusNonShared , ])## Mode FALSE TRUE NA's
## logical 7 25 2 # proportional contribution of shared exceeds nonshared in each?
summary(0>invsumm[ , PcrispusShared-PcrispusNonShared , ]) # CLP## Mode FALSE TRUE NA's
## logical 24 9 1 summary(0>invsumm[ , MspicatumShared-MspicatumNonShared , ]) # EWM## Mode FALSE NA's
## logical 33 1 # Changing shared and non-shared to absolute dissims (i.e., sum to total dissim rather than relativized to sum to 1)
invsumm[, MspicatumNonShared := MspicatumNonShared*MspicatumDissim]
invsumm[, MspicatumShared := MspicatumShared*MspicatumDissim]
invsumm[, PcrispusNonShared := PcrispusNonShared*PcrispusDissim]
invsumm[, PcrispusShared := PcrispusShared*PcrispusDissim]
# Changing column orders and dropping "PClessMS"
setcolorder(invsumm, c("Taxon", "PcrispusDissim", "PcrispusNonShared", "PcrispusShared",
"MspicatumDissim", "MspicatumNonShared", "MspicatumShared"))
invsumm[,PClessMS:=NULL]
# Adding, and ordering by, mean dissimilarity of each species to all other species
sppMeanDissim <- rowSums(dissim_TPDs_plants$populations$dissimilarity)/(dim(dissim_TPDs_plants$populations$dissimilarity)[2]-1)
invsumm[, sppMeanDissim := sppMeanDissim]
setcolorder(invsumm, c("Taxon", "sppMeanDissim"))
setorder(invsumm, sppMeanDissim)
invsumm[is.na(invsumm)] <- 0
# Export Table 1:
write.csv(invsumm, file = "invsumm.csv", row.names = FALSE)
# Creating heatmap-style figure from table
invsumm.m <- melt(invsumm, id.vars = c("Taxon"))
invsumm.m$Taxon <- factor(invsumm.m$Taxon, levels = invsumm$Taxon)
levels(invsumm.m$variable) <- c(
"All species\nTotal",
"P. crispus\nTotal",
"P. crispus\nNon-shared",
"P. crispus\nShared",
"M. spicatum\nTotal",
"M. spicatum\nNon-shared",
"M. spicatum\nShared")Dissimilarity Table
dissim_table <- ggplot(data = invsumm.m, aes(x=variable, y=Taxon, fill=value)) +
geom_tile() +
scale_fill_gradient(low = "white", high = "steelblue", name = "Dissimilarity") +
geom_text(aes(label = sprintf("%0.2f", invsumm.m$value)), size=3) +
labs(x="Comparison") +
scale_x_discrete(expand = c(0, 0), position = "top") +
theme(axis.ticks = element_blank())+
theme(axis.text.y = element_text(face= "italic"))+
theme(legend.position = "none")
# Write to file:
tiff(file = "Fig_DissimTable.tif", width=8.5, height=10, units="in", res = 600, compression = "lzw")
dissim_table
dev.off()## png
## 2# TPD niche axes one at a time --------------------------------------------Decompose by niche axis
We will want to 1) TPD each axis individually 2) Calculate dissimilarity for each axis 3) Decompose each axis’ dissimilarity into shared and non-shared niche space
traits_plants_d1 <- ptrm[, c("gdd")]
sp_plants <- ptrm$TAXON
TPDs_plants_d1_a1 <- TPDs(species = sp_plants, traits = traits_plants_d1, alpha= 0.95, n_divisions = 100)## -------Calculating densities for One population_Multiple species----------- # plotTPD(TPD = TPDs_plants_d1_a1, nRowCol = c(10,5))
traits_plants_d2 <- ptrm[, c(proplight)]
sp_plants <- ptrm$TAXON
TPDs_plants_d2_a2 <- TPDs(species = sp_plants, traits = traits_plants_d2, alpha=0.95, n_divisions = 100)## -------Calculating densities for One population_Multiple species-----------## Be careful, the integral of the pdf of some cases differ from 1.
## They have been reescaled, but you should consider increasing
## 'trait_ranges'
## Potamogeton crispus / Stuckenia pectinata # plotTPD(TPD = TPDs_plants_d2_a2, nRowCol = c(10,5) )
ptrm[TAXON=="Potamogeton crispus", summary(proplight)]## Min. 1st Qu. Median Mean 3rd Qu. Max.
## 0.000000 0.003089 0.024849 0.079315 0.098837 0.911636 traits_plants_d3 <- ptrm[, c("DEPTH_FT")]
sp_plants <- ptrm$TAXON
TPDs_plants_d3_a3 <- TPDs(species = sp_plants, traits = traits_plants_d3, alpha=0.95, n_divisions = 27)## -------Calculating densities for One population_Multiple species-----------## Be careful, the integral of the pdf of some cases differ from 1.
## They have been reescaled, but you should consider increasing
## 'trait_ranges'
## Lemna trisulca / Schoenoplectus acutus / Nuphar variegata / Potamogeton amplifolius / Potamogeton natans / Nymphaea odorata / Typha angustifolia / Zizania palustris / Brasenia schreberi / Utricularia vulgaris / Najas flexilis # plotTPD(TPD = TPDs_plants_d3_a3, nRowCol = c(10,5))
# calc dissimilarity and decompose by each axis ------------------------------Single niche dissimilarity
#dissimilarity across for each dimension
dissim_TPDs_gdd <- dissim(TPDs_plants_d1_a1)## Calculating dissimilarities between 34 populations. It might take some time dissim_TPDs_light <- dissim(TPDs_plants_d2_a2)## Calculating dissimilarities between 34 populations. It might take some time dissim_TPDs_depth <- dissim(TPDs_plants_d3_a3)## Calculating dissimilarities between 34 populations. It might take some time# dimensional contribution to dissim --------------------------------------------------We can make a plot that shows the relative contrib from each niche dim:
# Grab mean dissimilarity of each species to all other species for total and each niche axis
dimdissim <- data.table(Total = rowSums(dissim_TPDs_plants$populations$dissimilarity)/(dim(dissim_TPDs_plants$populations$dissimilarity)[2]-1))
dimdissim[ , Taxon:= rownames(dissim_TPDs_plants$populations$dissimilarity) , ]
dimdissim[ , GDD := rowSums(dissim_TPDs_gdd$populations$dissimilarity)/(dim(dissim_TPDs_gdd$populations$dissimilarity)[2]-1)]
# dimdissim[ , gddTaxon:= rownames(dissim_TPDs_gdd$populations$dissimilarity) , ]
# dimdissim[!Taxon == gddTaxon, , ]
dimdissim[ , Light := rowSums(dissim_TPDs_light$populations$dissimilarity)/(dim(dissim_TPDs_light$populations$dissimilarity)[2]-1), ]
# dimdissim[ , lightTaxon:= rownames(dissim_TPDs_light$populations$dissimilarity) , ]
# dimdissim[!Taxon == lightTaxon, , ]
dimdissim[ , Depth := rowSums(dissim_TPDs_depth$populations$dissimilarity)/(dim(dissim_TPDs_depth$populations$dissimilarity)[2]-1), ]
# dimdissim[ , depthTaxon:= rownames(dissim_TPDs_depth$populations$dissimilarity) , ]
# dimdissim[!Taxon == depthTaxon, , ]
#order by total dissim
setcolorder(dimdissim, c("Taxon", "GDD", "Light","Depth", "Total"))
setorder(dimdissim, Total)
# Creating heatmap-style figure from table
dimdissim.m <- melt(dimdissim, id.vars = c("Taxon"))
dimdissim.m$Taxon <- factor(dimdissim.m$Taxon, levels = dimdissim$Taxon)Dimensional Dissimilarity
full_dimfig <- ggplot(data = dimdissim.m, aes(x=variable, y=Taxon, fill=value)) +
geom_tile() +
scale_fill_gradient(low = "white", high = "steelblue", name = "Dissimilarity") +
geom_text(aes(label = sprintf("%0.2f", dimdissim.m$value)), size=3) +
labs(x="Dimensional Dissimilarity") +
scale_x_discrete(expand = c(0, 0), position = "top") +
theme(axis.ticks = element_blank())+
theme(legend.position = "none")+
theme(axis.text.y = element_text(face=c(ifelse(dimdissim$Taxon == "Myriophyllum spicatum" |
dimdissim$Taxon == "Potamogeton crispus", "bold.italic", "italic"))))## Warning: Vectorized input to `element_text()` is not officially supported.
## Results may be unexpected or may change in future versions of ggplot2. ifelse(dimdissim$Taxon == "Myriophyllum spicatum" |
dimdissim$Taxon == "Potamogeton crispus", "bold", "plain")## [1] "plain" "plain" "plain" "plain" "plain" "plain" "plain" "plain" "plain"
## [10] "plain" "plain" "plain" "plain" "plain" "plain" "plain" "plain" "plain"
## [19] "plain" "bold" "plain" "plain" "plain" "plain" "plain" "plain" "plain"
## [28] "plain" "plain" "plain" "plain" "plain" "plain" "bold" simp_dim <- dimdissim.m[Taxon == "Myriophyllum spicatum" |
Taxon == "Potamogeton crispus", , ]
#Simplified heatmap for only the invaders
simp_dimfig <- ggplot(data = simp_dim, aes(x=variable, y=Taxon, fill=value)) +
geom_tile() +
scale_fill_gradient(low = "white", high = "steelblue", name = "Dissimilarity") +
geom_text(aes(label = sprintf("%0.2f", simp_dim$value)), size=3) +
labs(x="Dimensional Dissimilarity") +
scale_x_discrete(expand = c(0, 0), position = "top") +
theme(axis.ticks = element_blank())+
theme(legend.position = "none")+
theme(axis.text.y = element_text(face= "italic"))
# # Write to file:
tiff(file = "Fig_DimensFull.tif", width=8.5, height=10, units="in", res = 600, compression = "lzw")
full_dimfig
dev.off()## png
## 2 tiff(file = "Fig_DimensSimp.tif", width=8.5, height=2, units="in", res = 600, compression = "lzw")
simp_dimfig
dev.off()## png
## 2# plot all three axis of interest -----------------------------------------Niche visualizations
#custom plotting all species in a single plot GDD
gddTPD <- data.table(matrix(unlist(TPDs_plants_d1_a1$TPDs), ncol = length(TPDs_plants_d1_a1$TPDs), byrow = F))
setnames(gddTPD, names(TPDs_plants_d1_a1$TPDs))
gddTPD[, gddinterval := as.vector(TPDs_plants_d1_a1$data$evaluation_grid) , ]
gddTPD <- melt(gddTPD, id.vars = "gddinterval", variable.name = "taxon" )
levels(gddTPD$taxon)## [1] "Lemna trisulca" "Potamogeton zosteriformis"
## [3] "Schoenoplectus acutus" "Spirodela polyrhiza"
## [5] "Elodea canadensis" "Ceratophyllum demersum"
## [7] "Nuphar variegata" "Potamogeton amplifolius"
## [9] "Potamogeton crispus" "Potamogeton richardsonii"
## [11] "Myriophyllum sibiricum" "Stuckenia pectinata"
## [13] "Lemna minor" "Potamogeton natans"
## [15] "Vallisneria americana" "Potamogeton robbinsii"
## [17] "Nymphaea odorata" "Typha angustifolia"
## [19] "Zizania palustris" "Potamogeton praelongus"
## [21] "Brasenia schreberi" "Utricularia vulgaris"
## [23] "Najas flexilis" "Potamogeton illinoensis"
## [25] "Potamogeton gramineus" "Bidens beckii"
## [27] "Potamogeton friesii" "Heteranthera dubia"
## [29] "Potamogeton pusillus" "Najas guadalupensis"
## [31] "Myriophyllum spicatum" "Ranunculus aquatilis"
## [33] "Zannichellia palustris" "Potamogeton foliosus" gddTPD$taxon <- factor(gddTPD$taxon,levels(gddTPD$taxon)[c(1:8,10:30,32:34,9,31)])
gdd_inv <- ggplot(gddTPD, aes(gddinterval, value))+
geom_line(aes(color = taxon), size = 1)+
scale_color_manual( values = c(ifelse(levels(gddTPD$taxon) %in% c("Potamogeton crispus") ==T, "black" , ifelse(levels(gddTPD$taxon) %in% c("Myriophyllum spicatum") ==T, "red" , "gray")))) +
theme(legend.position = "none") +
xlab("Growing Degrees Days")+
xlim(c(0,3000))+
ylab("Density")
#custom plotting all species in a single plot Proplight
plTPD <- data.table(matrix(unlist(TPDs_plants_d2_a2$TPDs), ncol = length(TPDs_plants_d2_a2$TPDs), byrow = F))
setnames(plTPD, names(TPDs_plants_d2_a2$TPDs))
plTPD[, plinterval := as.vector(TPDs_plants_d2_a2$data$evaluation_grid) , ]
plTPD <- melt(plTPD, id.vars = "plinterval", variable.name = "taxon" )
#bound probability viz to force null at 0 > x > 1
plTPD[plinterval < 0 | plinterval > 1, value := 0]
plTPD[taxon=="Potamogeton crispus"]## plinterval taxon value
## 1: -0.147297530 Potamogeton crispus 0.000000000
## 2: -0.134402797 Potamogeton crispus 0.000000000
## 3: -0.121508063 Potamogeton crispus 0.000000000
## 4: -0.108613330 Potamogeton crispus 0.000000000
## 5: -0.095718597 Potamogeton crispus 0.000000000
## 6: -0.082823864 Potamogeton crispus 0.000000000
## 7: -0.069929130 Potamogeton crispus 0.000000000
## 8: -0.057034397 Potamogeton crispus 0.000000000
## 9: -0.044139664 Potamogeton crispus 0.000000000
## 10: -0.031244931 Potamogeton crispus 0.000000000
## 11: -0.018350197 Potamogeton crispus 0.000000000
## 12: -0.005455464 Potamogeton crispus 0.000000000
## 13: 0.007439269 Potamogeton crispus 0.263327250
## 14: 0.020334002 Potamogeton crispus 0.116536979
## 15: 0.033228736 Potamogeton crispus 0.079357602
## 16: 0.046123469 Potamogeton crispus 0.063974497
## 17: 0.059018202 Potamogeton crispus 0.054067249
## 18: 0.071912935 Potamogeton crispus 0.043225680
## 19: 0.084807669 Potamogeton crispus 0.035972688
## 20: 0.097702402 Potamogeton crispus 0.033450618
## 21: 0.110597135 Potamogeton crispus 0.030965574
## 22: 0.123491868 Potamogeton crispus 0.021356387
## 23: 0.136386602 Potamogeton crispus 0.016125257
## 24: 0.149281335 Potamogeton crispus 0.022887822
## 25: 0.162176068 Potamogeton crispus 0.018763972
## 26: 0.175070801 Potamogeton crispus 0.011921513
## 27: 0.187965535 Potamogeton crispus 0.013295858
## 28: 0.200860268 Potamogeton crispus 0.013560863
## 29: 0.213755001 Potamogeton crispus 0.010975580
## 30: 0.226649734 Potamogeton crispus 0.011367194
## 31: 0.239544468 Potamogeton crispus 0.007352632
## 32: 0.252439201 Potamogeton crispus 0.009275314
## 33: 0.265333934 Potamogeton crispus 0.011109828
## 34: 0.278228667 Potamogeton crispus 0.008434925
## 35: 0.291123401 Potamogeton crispus 0.004806268
## 36: 0.304018134 Potamogeton crispus 0.006230392
## 37: 0.316912867 Potamogeton crispus 0.008091421
## 38: 0.329807600 Potamogeton crispus 0.005772068
## 39: 0.342702334 Potamogeton crispus 0.000000000
## 40: 0.355597067 Potamogeton crispus 0.004579633
## 41: 0.368491800 Potamogeton crispus 0.000000000
## 42: 0.381386533 Potamogeton crispus 0.000000000
## 43: 0.394281267 Potamogeton crispus 0.004361070
## 44: 0.407176000 Potamogeton crispus 0.006075157
## 45: 0.420070733 Potamogeton crispus 0.000000000
## 46: 0.432965466 Potamogeton crispus 0.000000000
## 47: 0.445860200 Potamogeton crispus 0.000000000
## 48: 0.458754933 Potamogeton crispus 0.000000000
## 49: 0.471649666 Potamogeton crispus 0.000000000
## 50: 0.484544399 Potamogeton crispus 0.000000000
## 51: 0.497439133 Potamogeton crispus 0.000000000
## 52: 0.510333866 Potamogeton crispus 0.000000000
## 53: 0.523228599 Potamogeton crispus 0.000000000
## 54: 0.536123332 Potamogeton crispus 0.000000000
## 55: 0.549018066 Potamogeton crispus 0.000000000
## 56: 0.561912799 Potamogeton crispus 0.000000000
## 57: 0.574807532 Potamogeton crispus 0.000000000
## 58: 0.587702265 Potamogeton crispus 0.000000000
## 59: 0.600596999 Potamogeton crispus 0.000000000
## 60: 0.613491732 Potamogeton crispus 0.000000000
## 61: 0.626386465 Potamogeton crispus 0.000000000
## 62: 0.639281198 Potamogeton crispus 0.000000000
## 63: 0.652175932 Potamogeton crispus 0.000000000
## 64: 0.665070665 Potamogeton crispus 0.000000000
## 65: 0.677965398 Potamogeton crispus 0.000000000
## 66: 0.690860131 Potamogeton crispus 0.000000000
## 67: 0.703754865 Potamogeton crispus 0.000000000
## 68: 0.716649598 Potamogeton crispus 0.000000000
## 69: 0.729544331 Potamogeton crispus 0.000000000
## 70: 0.742439064 Potamogeton crispus 0.000000000
## 71: 0.755333798 Potamogeton crispus 0.000000000
## 72: 0.768228531 Potamogeton crispus 0.000000000
## 73: 0.781123264 Potamogeton crispus 0.000000000
## 74: 0.794017997 Potamogeton crispus 0.000000000
## 75: 0.806912731 Potamogeton crispus 0.000000000
## 76: 0.819807464 Potamogeton crispus 0.000000000
## 77: 0.832702197 Potamogeton crispus 0.000000000
## 78: 0.845596930 Potamogeton crispus 0.000000000
## 79: 0.858491664 Potamogeton crispus 0.000000000
## 80: 0.871386397 Potamogeton crispus 0.000000000
## 81: 0.884281130 Potamogeton crispus 0.000000000
## 82: 0.897175863 Potamogeton crispus 0.000000000
## 83: 0.910070597 Potamogeton crispus 0.000000000
## 84: 0.922965330 Potamogeton crispus 0.000000000
## 85: 0.935860063 Potamogeton crispus 0.000000000
## 86: 0.948754796 Potamogeton crispus 0.000000000
## 87: 0.961649530 Potamogeton crispus 0.000000000
## 88: 0.974544263 Potamogeton crispus 0.000000000
## 89: 0.987438996 Potamogeton crispus 0.000000000
## 90: 1.000333729 Potamogeton crispus 0.000000000
## 91: 1.013228463 Potamogeton crispus 0.000000000
## 92: 1.026123196 Potamogeton crispus 0.000000000
## 93: 1.039017929 Potamogeton crispus 0.000000000
## 94: 1.051912662 Potamogeton crispus 0.000000000
## 95: 1.064807396 Potamogeton crispus 0.000000000
## 96: 1.077702129 Potamogeton crispus 0.000000000
## 97: 1.090596862 Potamogeton crispus 0.000000000
## 98: 1.103491595 Potamogeton crispus 0.000000000
## 99: 1.116386329 Potamogeton crispus 0.000000000
## 100: 1.129281062 Potamogeton crispus 0.000000000
## plinterval taxon value levels(plTPD$taxon)## [1] "Lemna trisulca" "Potamogeton zosteriformis"
## [3] "Schoenoplectus acutus" "Spirodela polyrhiza"
## [5] "Elodea canadensis" "Ceratophyllum demersum"
## [7] "Nuphar variegata" "Potamogeton amplifolius"
## [9] "Potamogeton crispus" "Potamogeton richardsonii"
## [11] "Myriophyllum sibiricum" "Stuckenia pectinata"
## [13] "Lemna minor" "Potamogeton natans"
## [15] "Vallisneria americana" "Potamogeton robbinsii"
## [17] "Nymphaea odorata" "Typha angustifolia"
## [19] "Zizania palustris" "Potamogeton praelongus"
## [21] "Brasenia schreberi" "Utricularia vulgaris"
## [23] "Najas flexilis" "Potamogeton illinoensis"
## [25] "Potamogeton gramineus" "Bidens beckii"
## [27] "Potamogeton friesii" "Heteranthera dubia"
## [29] "Potamogeton pusillus" "Najas guadalupensis"
## [31] "Myriophyllum spicatum" "Ranunculus aquatilis"
## [33] "Zannichellia palustris" "Potamogeton foliosus" plTPD$taxon <- factor(plTPD$taxon,levels(plTPD$taxon)[c(1:8,10:30,32:34,9,31)])
#add a color selector for legend
plTPD[ , colorInv := as.factor(ifelse(taxon == "Potamogeton crispus", "black" , ifelse(taxon=="Myriophyllum spicatum", "red" , "gray"))) , ]
plTPD$colorInv <- factor(plTPD$colorInv, levels = c("black", "red", "gray"))
pl_inv <- ggplot(plTPD, aes(group = taxon, x = plinterval, y = value, color = colorInv))+
geom_line( size = 1)+
scale_color_manual(name = '', values=c('black','red', 'gray'),
labels = c( expression(italic("Potamogeton crispus")),
expression(italic("Myriophyllum spicatum")),
"Native species")) +
theme(legend.position = c(0.75,0.5),
legend.key = element_blank(),
legend.background = element_blank(),
legend.text=element_text(size=12.5),
legend.text.align = 0)+
xlab("Proportion of Surface Light Remaining at Sediments")+
xlim(c(-0.01,1.01))+
# coord_cartesian(ylim = c(0,0.20))+
ylab("Density")+
scale_y_continuous(breaks = c(0.00, 0.10, 0.20), labels = c("0.00", "0.10", "0.20"))
#custom plotting all species in a single plot Depth
depTPD <- data.table(matrix(unlist(TPDs_plants_d3_a3$TPDs), ncol = length(TPDs_plants_d3_a3$TPDs), byrow = F))
setnames(depTPD, names(TPDs_plants_d2_a2$TPDs))
depTPD[, depinterval := as.vector(TPDs_plants_d3_a3$data$evaluation_grid) , ]
depTPD <- melt(depTPD, id.vars = "depinterval", variable.name = "taxon" )
depTPD[depinterval < 0 , value := 0]
levels(depTPD$taxon)## [1] "Lemna trisulca" "Potamogeton zosteriformis"
## [3] "Schoenoplectus acutus" "Spirodela polyrhiza"
## [5] "Elodea canadensis" "Ceratophyllum demersum"
## [7] "Nuphar variegata" "Potamogeton amplifolius"
## [9] "Potamogeton crispus" "Potamogeton richardsonii"
## [11] "Myriophyllum sibiricum" "Stuckenia pectinata"
## [13] "Lemna minor" "Potamogeton natans"
## [15] "Vallisneria americana" "Potamogeton robbinsii"
## [17] "Nymphaea odorata" "Typha angustifolia"
## [19] "Zizania palustris" "Potamogeton praelongus"
## [21] "Brasenia schreberi" "Utricularia vulgaris"
## [23] "Najas flexilis" "Potamogeton illinoensis"
## [25] "Potamogeton gramineus" "Bidens beckii"
## [27] "Potamogeton friesii" "Heteranthera dubia"
## [29] "Potamogeton pusillus" "Najas guadalupensis"
## [31] "Myriophyllum spicatum" "Ranunculus aquatilis"
## [33] "Zannichellia palustris" "Potamogeton foliosus" depTPD$taxon <- factor(depTPD$taxon,levels(depTPD$taxon)[c(1:8,10:30,32:34,9,31)])
dep_inv <- ggplot(depTPD, aes(depinterval*0.3048, value))+
geom_line(aes(color = taxon), size = 1)+
scale_color_manual( values = c(ifelse(levels(depTPD$taxon) %in% c("Potamogeton crispus") ==T, "black" , ifelse(levels(depTPD$taxon) %in% c("Myriophyllum spicatum") ==T, "red" , "gray")))) +
theme(legend.position = "none")+
xlab("Depth (m)")+
xlim( c(-0.05,7))+
ylab("Density")+
scale_y_continuous(breaks = c(0.00, 0.20, 0.40), labels = c("0.00", "0.20", "0.40"))
# Make panel by niche dimension plot -----------------------------------------------------------
# ## One dimensional vizualizations
# # Figure with 3 axes broken out with overlapping species niches
fig_3axes <- grid.arrange(gdd_inv, pl_inv, dep_inv)## Warning: Removed 986 row(s) containing missing values (geom_path).## Warning: Removed 714 row(s) containing missing values (geom_path).## Warning: Removed 204 row(s) containing missing values (geom_path).
tiff("fig_3axes.tiff", res = 600, width = 9, height = 12, units = "in", compression = "lzw")
plot(fig_3axes) # Make plot
dev.off()## png
## 2# Model TPD 2 dimensions, 3 ways for all species ------------------------------------------
#gdd v dep
TPDs_gdddep <- TPDs(species = ptrm[, TAXON ], traits = ptrm[, .(gdd,DEPTH_FT*0.3048) ] )## -------Calculating densities for One population_Multiple species----------- #gdd v proplight
TPDs_gddlight <- TPDs(species = ptrm[, TAXON ], traits = ptrm[, .(gdd,proplight) ] )## -------Calculating densities for One population_Multiple species----------- #dep v proplight
TPDs_deplight <- TPDs(species = ptrm[, TAXON ], traits = ptrm[, .(DEPTH_FT*0.3048, proplight) ] )## -------Calculating densities for One population_Multiple species-----------# Library of niches --------------------------------------------------Great. That covers plotting the sampling distribution in the background. Now scale this to allow calc and paneling by crossed niche and species.
## Sampling “Niche” visualization This is a way to condier how the sampling was biased across our dataset
#TPD for full cross of gdd, depth, proplight for sample effort dist
locssamp[ , surveyed := "Samplous distributionous",]
TPDs_effort_gd <- TPDs(species = locssamp[, surveyed , ], traits = locssamp[, .(gdd, DEPTH_FT*0.3048) ] )## -------Calculating densities for One population_One species-----------## Be careful, the integral of the pdf of some cases differ from 1.
## They have been reescaled, but you should consider increasing
## 'trait_ranges'
## Samplous distributionous TPDs_effort_gl <- TPDs(species = locssamp[, surveyed , ], traits = locssamp[, .(gdd, proplight) ] )## -------Calculating densities for One population_One species----------- TPDs_effort_ld <- TPDs(species = locssamp[, surveyed , ], traits = locssamp[, .(proplight, DEPTH_FT*0.3048) ] )## -------Calculating densities for One population_One species-----------
##
## Be careful, the integral of the pdf of some cases differ from 1.
## They have been reescaled, but you should consider increasing
## 'trait_ranges'
## Samplous distributionous # effort gdd by dep
e_gd <- data.table(gdd = TPDs_effort_gd$data$evaluation_grid$gdd,
dep = TPDs_effort_gd$data$evaluation_grid$V2 ,
tpd = TPDs_effort_gd$TPDs$`Samplous distributionous`)
e_gd[tpd == 0 , tpd := NA , ]
gd_samps <- ggplot(e_gd, aes(x = gdd, y = dep))+
geom_tile(data = e_gd[is.na(tpd) == F, ], aes(fill = tpd ))
# scale_fill_gradient(low = "gray", high = "gray")+
# new_scale_fill()+
# geom_tile(data = gd[is.na(tpdp) == F, ], aes(fill = tpdp ))
# effort light by dep
e_gp <- data.table(gdd = TPDs_effort_gl$data$evaluation_grid$gdd,
pl = TPDs_effort_gl$data$evaluation_grid$proplight ,
tpd = TPDs_effort_gl$TPDs$`Samplous distributionous`)
e_gp[tpd == 0 , tpd := NA , ]
gl_samps <- ggplot(e_gp, aes(x = gdd, y = pl))+
geom_tile(data = e_gp[is.na(tpd) == F, ], aes(fill = tpd ))
# effort light by dep
e_pd <- data.table(dep = TPDs_effort_ld$data$evaluation_grid$V2,
pl = TPDs_effort_ld$data$evaluation_grid$proplight ,
tpd = TPDs_effort_ld$TPDs$`Samplous distributionous`)
e_pd[tpd == 0 , tpd := NA , ]
dl_samps <- ggplot(e_pd, aes(x = pl, y = dep))+
geom_tile(data = e_pd[is.na(tpd) == F, ], aes(fill = tpd ))
fig_sampleeffort <- grid.arrange(gd_samps, gl_samps, dl_samps, ncol = 3)
# write to file
# tiff("fig_sampledist.tiff", res = 600, width = 12, height = 4, units = "in", compression = "lzw")
# plot(fig_sampleeffort) # Make plot
# dev.off()Okay, so we have a way to plot what we want, now we have to build a table from which we can plot all of the species data. This need to be by crossed set
Library of Niches Visualizations
# library of niches gdd x depth -------------------------------------------
# gdd x depth
#effort data
e_gd <- data.table(gdd = TPDs_effort_gd$data$evaluation_grid$gdd,
dep = TPDs_effort_gd$data$evaluation_grid$V2 ,
tpd = TPDs_effort_gd$TPDs$`Samplous distributionous`)
#make species TPDs into data.table
setDT(TPDs_gdddep$TPDs)
speciesdat <- data.table(gdd = TPDs_gdddep$data$evaluation_grid$gdd,
dep = TPDs_gdddep$data$evaluation_grid$V2,
TPDs_gdddep$TPDs
)
#make long format
speciesdat.m <- melt(speciesdat, id.vars = c('gdd', 'dep' ))
#plot:
lib_GDDDepth <- ggplot(e_gd, aes(x = gdd, y = dep))+
geom_tile(data = e_gd[tpd > 0 , ], aes(fill = tpd ))+
scale_fill_gradient(low = "gray", high = "gray")+
theme(legend.position = NULL)+
new_scale_fill()+
geom_tile(data = speciesdat.m[value > 0, ], aes(fill = value, group = variable ))+
facet_wrap(~ variable, ncol = 5 )+
theme(legend.position = NULL)+
xlab("Growing degree days")+
ylab("Depth (m)")+
theme(strip.text = element_text(face = "italic"))
#plot only the ususal suspects
inv_GDDDepth <- ggplot(e_gd, aes(x = gdd, y = dep))+
geom_tile(data = e_gd[tpd > 0 , ], aes(fill = tpd ))+
scale_fill_gradient(low = "gray", high = "gray")+
theme(legend.position = NULL)+
new_scale_fill()+
geom_tile(data = speciesdat.m[value > 0 &
variable %in% c("Potamogeton crispus", "Myriophyllum spicatum"), ], aes(fill = value))+
facet_wrap(~ variable, ncol = 2 )+
# scale_fill_gradient(low = "lightblue" , high = "steelblue", name = "Dissimilarity") +
theme(legend.position = "none")+
xlab("Growing degree days")+
ylab("Depth (m)")+
theme(strip.text = element_text(face = "italic"))
# library of niches gdd x proplight -------------------------------------------
# gdd x light
#effort data
e_gl <- data.table(gdd = TPDs_effort_gl$data$evaluation_grid$gdd,
light = TPDs_effort_gl$data$evaluation_grid$proplight ,
tpd = TPDs_effort_gl$TPDs$`Samplous distributionous`)
#make species TPDs into data.table
setDT(TPDs_gddlight$TPDs)
speciesdat <- data.table(gdd = TPDs_gddlight$data$evaluation_grid$gdd,
light = TPDs_gddlight$data$evaluation_grid$proplight,
TPDs_gddlight$TPDs
)
#make long format
speciesdat.m <- melt(speciesdat, id.vars = c('gdd', 'light' ))
#plot:
lib_GDDLight <- ggplot(e_gl, aes(x = gdd, y = light))+
geom_tile(data = e_gl[tpd > 0 , ], aes(fill = tpd ))+
scale_fill_gradient(low = "gray", high = "gray")+
theme(legend.position = 'none')+
new_scale_fill()+
geom_tile(data = speciesdat.m[value > 0, ], aes(fill = value, group = variable ))+
facet_wrap(~ variable, ncol = 5 )+
theme(legend.position = NULL)+
xlab("Growing degree days")+
ylab("Proportion of surface irradiance\nremaining at substrate")+
theme(strip.text = element_text(face = "italic"))
#plot only the ususal suspects
inv_GDDLight <- ggplot(e_gl, aes(x = gdd, y = light))+
geom_tile(data = e_gl[tpd > 0 , ], aes(fill = tpd ))+
scale_fill_gradient(low = "gray", high = "gray")+
theme(legend.position = NULL)+
new_scale_fill()+
geom_tile(data = speciesdat.m[value > 0 &
variable %in% c("Potamogeton crispus", "Myriophyllum spicatum"), ], aes(fill = value))+
facet_wrap(~ variable, ncol = 2 )+
# scale_fill_gradient(low = "steelblue", high = "red", name = "Dissimilarity") +
theme(legend.position = "none")+
xlab("Growing degree days")+
ylab("Proportion of surface irradiance\nremaining at substrate")+
theme(strip.text = element_text(face = "italic"))
# library of niches depth x light -----------------------------------------
# depth x light
#effort data
e_dl <- data.table(dep = TPDs_effort_ld$data$evaluation_grid$V2,
light = TPDs_effort_ld$data$evaluation_grid$proplight ,
tpd = TPDs_effort_ld$TPDs$`Samplous distributionous`)
#make species TPDs into data.table
setDT(TPDs_deplight$TPDs)
speciesdat <- data.table(dep = TPDs_deplight$data$evaluation_grid$V1,
light = TPDs_deplight$data$evaluation_grid$proplight,
TPDs_deplight$TPDs
)
#make long format
speciesdat.m <- melt(speciesdat, id.vars = c('dep', 'light' ))
#plot:
lib_DepthLight <- ggplot(e_dl, aes(x = dep, y = light))+
geom_tile(data = e_dl[tpd > 0 , ], aes(fill = tpd ))+
scale_fill_gradient(low = "gray", high = "gray")+
theme(legend.position = 'none')+
new_scale_fill()+
geom_tile(data = speciesdat.m[value > 0, ], aes(fill = value, group = variable ))+
facet_wrap(~ variable, ncol = 5 )+
theme(legend.position = NULL)+
ylab("Proportion of surface irradiance\nremaining at substrate")+
xlab("Depth (m)")+
theme(strip.text = element_text(face = "italic"))
#plot only the ususal suspects
inv_DepthLight <- ggplot(e_dl, aes(x = dep, y = light))+
geom_tile(data = e_dl[tpd > 0 , ], aes(fill = tpd ))+
scale_fill_gradient(low = "gray", high = "gray")+
theme(legend.position = NULL)+
new_scale_fill()+
geom_tile(data = speciesdat.m[value > 0 &
variable %in% c("Potamogeton crispus", "Myriophyllum spicatum"), ], aes(fill = value))+
facet_wrap(~ variable, ncol = 2 )+
# scale_fill_gradient(low = "steelblue", high = "red", name = "Dissimilarity") +
theme(legend.position = "none")+
ylab("Proportion of surface irradiance\nremaining at substrate")+
xlab("Depth (m)")+
theme(strip.text = element_text(face = "italic"))
# invader niche plots -----------------------------------------------------
fig_invaderNiche <- grid.arrange(inv_GDDDepth, inv_GDDLight, inv_DepthLight) 
# Write to file:
tiff(file = "Fig_InvaderNiches.tif", width=8.5, height=11, units="in", res = 600, compression = "lzw")
grid.arrange(inv_GDDDepth, inv_GDDLight, inv_DepthLight)
dev.off()## png
## 2# library of niches dimensional -------------------------------------------------------
# lib_GDDDepth
# lib_GDDLight
# lib_DepthLight
# Write to file:
tiff(file = "Fig_LibraryGDDDepth.tif", width=10, height=13, units="in", res = 600, compression = "lzw")
lib_GDDDepth
dev.off()## png
## 2 tiff(file = "Fig_LibraryGDDLight.tif", width=10, height=13, units="in", res = 600, compression = "lzw")
lib_GDDLight
dev.off()## png
## 2 tiff(file = "Fig_LibraryDepthLight.tif", width=10, height=13, units="in", res = 600, compression = "lzw")
lib_DepthLight
dev.off()## png
## 2# stat diff in niche dimensional diffs between clp & ewm -------------------------------Statistical tests referenced in text
mean_traits <- plantocc[ ,.(gdd_hat = mean(gdd), proplight_hat = mean(proplight) , depthmeter_hat = mean(DEPTH_FT*0.3048) ) , by = TAXON ]
community_means <- plantocc[!TAXON %in% c('Myriophyllum spicatum', 'Potamogeton crispus'),.(gdd_hat = mean(gdd), proplight_hat = mean(proplight) , depthmeter_hat = mean(DEPTH_FT*0.3048) ) ]
invader_means <- plantocc[TAXON %in% c('Myriophyllum spicatum', 'Potamogeton crispus'),.(gdd_hat = mean(gdd), proplight_hat = mean(proplight) , depthmeter_hat = mean(DEPTH_FT*0.3048) ), TAXON ]
inv_occ <- plantocc[ TAXON %in% c('Myriophyllum spicatum', 'Potamogeton crispus'),.(TAXON,gdd, proplight, depm=DEPTH_FT*0.3048 )]
inv_occ[ , gdd_cs := scale(gdd) ,]
inv_occ[ , proplight_cs := scale(proplight) ,]
inv_occ[ , depm_cs := scale(depm) ,]
#' comparing TPDs by axis:
t.test(inv_occ[TAXON == 'Myriophyllum spicatum', gdd ], inv_occ[TAXON == 'Potamogeton crispus', gdd ])##
## Welch Two Sample t-test
##
## data: inv_occ[TAXON == "Myriophyllum spicatum", gdd] and inv_occ[TAXON == "Potamogeton crispus", gdd]
## t = 95.999, df = 21337, p-value < 2.2e-16
## alternative hypothesis: true difference in means is not equal to 0
## 95 percent confidence interval:
## 591.8439 616.5158
## sample estimates:
## mean of x mean of y
## 1350.9295 746.7496 t.test(inv_occ[TAXON == 'Myriophyllum spicatum', proplight ], inv_occ[TAXON == 'Potamogeton crispus', proplight ])##
## Welch Two Sample t-test
##
## data: inv_occ[TAXON == "Myriophyllum spicatum", proplight] and inv_occ[TAXON == "Potamogeton crispus", proplight]
## t = 61.633, df = 20926, p-value < 2.2e-16
## alternative hypothesis: true difference in means is not equal to 0
## 95 percent confidence interval:
## 0.09822972 0.10468287
## sample estimates:
## mean of x mean of y
## 0.18077178 0.07931548 t.test(inv_occ[TAXON == 'Myriophyllum spicatum', depm ], inv_occ[TAXON == 'Potamogeton crispus', depm ])##
## Welch Two Sample t-test
##
## data: inv_occ[TAXON == "Myriophyllum spicatum", depm] and inv_occ[TAXON == "Potamogeton crispus", depm]
## t = -36.437, df = 30426, p-value < 2.2e-16
## alternative hypothesis: true difference in means is not equal to 0
## 95 percent confidence interval:
## -0.3691724 -0.3314826
## sample estimates:
## mean of x mean of y
## 1.687721 2.038048 #with z-transformed vars
t.test(inv_occ[TAXON == 'Myriophyllum spicatum', gdd_cs ], inv_occ[TAXON == 'Potamogeton crispus', gdd_cs ])##
## Welch Two Sample t-test
##
## data: inv_occ[TAXON == "Myriophyllum spicatum", gdd_cs] and inv_occ[TAXON == "Potamogeton crispus", gdd_cs]
## t = 95.999, df = 21337, p-value < 2.2e-16
## alternative hypothesis: true difference in means is not equal to 0
## 95 percent confidence interval:
## 0.9477379 0.9872458
## sample estimates:
## mean of x mean of y
## 0.6800391 -0.2874527 t.test(inv_occ[TAXON == 'Myriophyllum spicatum', proplight_cs ], inv_occ[TAXON == 'Potamogeton crispus', proplight_cs ])##
## Welch Two Sample t-test
##
## data: inv_occ[TAXON == "Myriophyllum spicatum", proplight_cs] and inv_occ[TAXON == "Potamogeton crispus", proplight_cs]
## t = 61.633, df = 20926, p-value < 2.2e-16
## alternative hypothesis: true difference in means is not equal to 0
## 95 percent confidence interval:
## 0.6460729 0.6885163
## sample estimates:
## mean of x mean of y
## 0.4690339 -0.1982607 t.test(inv_occ[TAXON == 'Myriophyllum spicatum', depm_cs ], inv_occ[TAXON == 'Potamogeton crispus', depm_cs ])##
## Welch Two Sample t-test
##
## data: inv_occ[TAXON == "Myriophyllum spicatum", depm_cs] and inv_occ[TAXON == "Potamogeton crispus", depm_cs]
## t = -36.437, df = 30426, p-value < 2.2e-16
## alternative hypothesis: true difference in means is not equal to 0
## 95 percent confidence interval:
## -0.3576708 -0.3211552
## sample estimates:
## mean of x mean of y
## -0.2385696 0.1008434 #Kolmogorov–Smirnov test
ks.test(inv_occ[TAXON == 'Myriophyllum spicatum', gdd ], inv_occ[TAXON == 'Potamogeton crispus', gdd ])## Warning in ks.test(inv_occ[TAXON == "Myriophyllum spicatum", gdd], inv_occ[TAXON
## == : p-value will be approximate in the presence of ties##
## Two-sample Kolmogorov-Smirnov test
##
## data: inv_occ[TAXON == "Myriophyllum spicatum", gdd] and inv_occ[TAXON == "Potamogeton crispus", gdd]
## D = 0.41283, p-value < 2.2e-16
## alternative hypothesis: two-sided ks.test(inv_occ[TAXON == 'Myriophyllum spicatum', proplight ], inv_occ[TAXON == 'Potamogeton crispus', proplight ])## Warning in ks.test(inv_occ[TAXON == "Myriophyllum spicatum", proplight], : p-
## value will be approximate in the presence of ties##
## Two-sample Kolmogorov-Smirnov test
##
## data: inv_occ[TAXON == "Myriophyllum spicatum", proplight] and inv_occ[TAXON == "Potamogeton crispus", proplight]
## D = 0.31699, p-value < 2.2e-16
## alternative hypothesis: two-sided ks.test(inv_occ[TAXON == 'Myriophyllum spicatum', depm ], inv_occ[TAXON == 'Potamogeton crispus', depm ])## Warning in ks.test(inv_occ[TAXON == "Myriophyllum spicatum", depm],
## inv_occ[TAXON == : p-value will be approximate in the presence of ties##
## Two-sample Kolmogorov-Smirnov test
##
## data: inv_occ[TAXON == "Myriophyllum spicatum", depm] and inv_occ[TAXON == "Potamogeton crispus", depm]
## D = 0.15968, p-value < 2.2e-16
## alternative hypothesis: two-sided #with z-transformed vars
ks.test(inv_occ[TAXON == 'Myriophyllum spicatum', gdd_cs ], inv_occ[TAXON == 'Potamogeton crispus', gdd_cs ])## Warning in ks.test(inv_occ[TAXON == "Myriophyllum spicatum", gdd_cs],
## inv_occ[TAXON == : p-value will be approximate in the presence of ties##
## Two-sample Kolmogorov-Smirnov test
##
## data: inv_occ[TAXON == "Myriophyllum spicatum", gdd_cs] and inv_occ[TAXON == "Potamogeton crispus", gdd_cs]
## D = 0.41283, p-value < 2.2e-16
## alternative hypothesis: two-sided ks.test(inv_occ[TAXON == 'Myriophyllum spicatum', proplight_cs ], inv_occ[TAXON == 'Potamogeton crispus', proplight_cs ])## Warning in ks.test(inv_occ[TAXON == "Myriophyllum spicatum", proplight_cs], : p-
## value will be approximate in the presence of ties##
## Two-sample Kolmogorov-Smirnov test
##
## data: inv_occ[TAXON == "Myriophyllum spicatum", proplight_cs] and inv_occ[TAXON == "Potamogeton crispus", proplight_cs]
## D = 0.31699, p-value < 2.2e-16
## alternative hypothesis: two-sided ks.test(inv_occ[TAXON == 'Myriophyllum spicatum', depm_cs ], inv_occ[TAXON == 'Potamogeton crispus', depm_cs ])## Warning in ks.test(inv_occ[TAXON == "Myriophyllum spicatum", depm_cs],
## inv_occ[TAXON == : p-value will be approximate in the presence of ties##
## Two-sample Kolmogorov-Smirnov test
##
## data: inv_occ[TAXON == "Myriophyllum spicatum", depm_cs] and inv_occ[TAXON == "Potamogeton crispus", depm_cs]
## D = 0.15968, p-value < 2.2e-16
## alternative hypothesis: two-sidedNow testing for who has more cases of non-shared exceeding
# map of data -------------------------------------------------------------Map the data
Here’s all the data we started from:
#plant surveys
plants_map <- fread(file = "plants_DRUM.csv", drop = 1)
locs <- plants_map[is.na(TAXON) == F , .(.N, unique(center_utm), unique(center_u_1)), by = DOWLKNUM]
weatherlocs <- plants_map[is.na(TAXON) == F , .(.N, longitude = unique(Long_WX), latitude = unique(Lat_WX)), by = STATION.NAME]
#drop added plants layer
rm(plants_map)
summary(locs)## DOWLKNUM N V2 V3
## Length:1499 Min. : 1.0 Min. :208956 Min. :4817240
## Class :character 1st Qu.: 50.5 1st Qu.:299707 1st Qu.:4980617
## Mode :character Median : 129.0 Median :393954 Median :5062676
## Mean : 381.1 Mean :389070 Mean :5075833
## 3rd Qu.: 327.0 3rd Qu.:460782 3rd Qu.:5184050
## Max. :11082.0 Max. :721551 Max. :5429064
## NA's :39 NA's :39 #drop waterbodies w/o any location data (these are NWI numbered wetlands)
locs <- na.omit(locs, cols = c("V2","V3"))
locsutm <- SpatialPointsDataFrame(locs[,.(V2,V3),], locs, proj4string=CRS("+proj=utm +zone=15N +datum=WGS84"))
locsgeo <- spTransform(locsutm, CRS("+proj=longlat +datum=WGS84"))
str(locsgeo)## Formal class 'SpatialPointsDataFrame' [package "sp"] with 5 slots
## ..@ data :'data.frame': 1460 obs. of 4 variables:
## .. ..$ DOWLKNUM: chr [1:1460] "1000100" "1001600" "1003400" "1003500" ...
## .. ..$ N : int [1:1460] 285 40 82 153 147 111 127 103 133 99 ...
## .. ..$ V2 : num [1:1460] 494097 492017 483970 479140 482832 ...
## .. ..$ V3 : num [1:1460] 5114741 5184698 5173669 5170765 5190458 ...
## ..@ coords.nrs : num(0)
## ..@ coords : num [1:1460, 1:2] -93.1 -93.1 -93.2 -93.3 -93.2 ...
## .. ..- attr(*, "dimnames")=List of 2
## .. .. ..$ : NULL
## .. .. ..$ : chr [1:2] "V2" "V3"
## ..@ bbox : num [1:2, 1:2] -96.7 43.5 -90 49
## .. ..- attr(*, "dimnames")=List of 2
## .. .. ..$ : chr [1:2] "V2" "V3"
## .. .. ..$ : chr [1:2] "min" "max"
## ..@ proj4string:Formal class 'CRS' [package "sp"] with 1 slot
## .. .. ..@ projargs: chr "+proj=longlat +datum=WGS84 +no_defs" locs_sf <- st_as_sf(locsgeo)
summary(locs_sf$N)## Min. 1st Qu. Median Mean 3rd Qu. Max.
## 1 50 130 380 322 11082 locs_sf$bin <- cut(locs_sf$N, breaks = c(0,50, 100, 250, 500, 1000, 5000, 10000, 15000),
labels = c("50", "100", "250", "500", "1000", "5000", "10000", "12000") )
#drop WX station NAs:
weatherlocs <- na.omit(weatherlocs, cols = c("longitude","latitude"))
wxlocsutm <- SpatialPointsDataFrame(weatherlocs[,.(longitude, latitude),], weatherlocs, proj4string=CRS("+proj=longlat +zone=15N +datum=WGS84"))
wxlocsgeo <- spTransform(wxlocsutm, CRS("+proj=longlat +datum=WGS84"))
str(wxlocsgeo)## Formal class 'SpatialPointsDataFrame' [package "sp"] with 5 slots
## ..@ data :'data.frame': 16 obs. of 4 variables:
## .. ..$ STATION.NAME: chr [1:16] "DULUTH INTERNATIONAL AIRPORT" "CHISHOLM-HIBBING AIRPORT" "ST CLOUD REGIONAL AIRPORT" "PK RPDS MUNI-KONSHOK FD AP" ...
## .. ..$ N : int [1:16] 20791 25808 64577 140771 198868 7281 5908 62677 11960 3999 ...
## .. ..$ longitude : num [1:16] -92.2 -92.8 -94.1 -95.1 -93.2 ...
## .. ..$ latitude : num [1:16] 46.8 47.4 45.5 46.9 44.9 ...
## ..@ coords.nrs : num(0)
## ..@ coords : num [1:16, 1:2] -92.2 -92.8 -94.1 -95.1 -93.2 ...
## .. ..- attr(*, "dimnames")=List of 2
## .. .. ..$ : NULL
## .. .. ..$ : chr [1:2] "longitude" "latitude"
## ..@ bbox : num [1:2, 1:2] -97.5 43.2 -88.5 49
## .. ..- attr(*, "dimnames")=List of 2
## .. .. ..$ : chr [1:2] "longitude" "latitude"
## .. .. ..$ : chr [1:2] "min" "max"
## ..@ proj4string:Formal class 'CRS' [package "sp"] with 1 slot
## .. .. ..@ projargs: chr "+proj=longlat +datum=WGS84 +no_defs" wxlocs_sf <- st_as_sf(wxlocsgeo)
#other data for fig
usa <- map_data("usa")
canada <- map_data("world", region = "canada")
states <- map_data("state")
mn_df <- subset(states, region == c("minnesota"))
# MN DNR Ecological Classification System retrievd from: https://gisdata.mn.gov/dataset/geos-ecological-class-system
mn_er <- st_read("ecs_provinces_of_mn_v99a.shp")## Reading layer `ecs_provinces_of_mn_v99a' from data source `E:\My Drive\Documents\UMN\Grad School\Larkin Lab\R_projects\macropheno\UMN_DRUM_upload\ecs_provinces_of_mn_v99a.shp' using driver `ESRI Shapefile'
## Simple feature collection with 4 features and 9 fields
## geometry type: POLYGON
## dimension: XY
## bbox: xmin: 189775.3 ymin: 4816305 xmax: 761662.3 ymax: 5472414
## projected CRS: NAD83 / UTM zone 15N #TNC ecoregions from resilient sites analysis retrievd from: https://gisdata.mn.gov/dataset/env-resilient-sites-tnc
tnc_er<- st_read("Ecoregions.shp")## Reading layer `Ecoregions' from data source `E:\My Drive\Documents\UMN\Grad School\Larkin Lab\R_projects\macropheno\UMN_DRUM_upload\Ecoregions.shp' using driver `ESRI Shapefile'
## Simple feature collection with 8 features and 20 fields
## geometry type: MULTIPOLYGON
## dimension: XY
## bbox: xmin: -1015929 ymin: 4282384 xmax: 2054137 ymax: 6195680
## projected CRS: NAD83 / UTM zone 15N tnc_er <- tnc_er[c(2,3,4,8) , ]
#main
study_map <- ggplot(locs_sf)+
geom_sf(data = tnc_er, aes(fill = ECO_NAME), alpha = 0.2 )+
#geom_sf(data = mn_er, aes(fill = PROVNAME), alpha = 0.2 )+
scale_fill_viridis_d()+
geom_sf( size = 3, shape = 16, colour = "black", alpha = 0.3)+
geom_sf(data = wxlocs_sf, shape = 10 , colour = "red", alpha = 1.0, size = 5)+
geom_polygon(data = states, aes(x=long, y = lat, group = group), color = "black", alpha = .0, lwd = .75)+
coord_sf(xlim = c(-98, -87.5), ylim = c(43, 49.5))+
theme(text = element_text(size=15))+
ylab("Latitude")+
xlab("Longitude")+
theme(legend.position = c(.755,.29),
legend.background = element_rect(fill = "white"))+
guides(fill=guide_legend(title="Ecoregion"))
#inset
inset <- ggplotGrob(
ggplot() +
geom_polygon(data = usa,
aes(x=long, y = lat, group = group),
fill = "gray90", color = "gray50", size = 0.3) +
# geom_polygon(data = canada, aes(x=long, y = lat, group = group),
# fill = "gray90", color = "gray50", size = 0.3) +
coord_fixed(xlim = c(-110, -55), ylim = c(25, 60), ratio = 1.2) +
geom_polygon(data = mn_df, aes(x = long, y = lat), color = "black", size = 0.3,
fill = "gray40") +
geom_polygon(data = states, aes(x=long, y = lat, group = group), color = "black", alpha = .0, lwd = .75)+
theme_bw() +
theme(axis.title = element_blank(), axis.ticks = element_blank(), axis.text = element_blank())+
coord_map("polyconic")
) ## Coordinate system already present. Adding new coordinate system, which will replace the existing one. #combine
g3 <- study_map +
north(mn_df, anchor = c(x = -88.15, y = 44)) +
annotation_custom(grob = inset, xmin = -91.5, xmax = -87.75,
ymin = 42.7)+
scalebar(mn_df, dist = 50, dist_unit = "km",
transform = TRUE, model = "WGS84", st.size = 3)
# figure construction -----------------------------------------------------
#Map
tiff("Map.tiff", res = 600, width = 9, height = 12, units = "in", compression = "lzw", type = "cairo")
plot(g3) # Make plot
dev.off()## png
## 2 # footer ------------------------------------------------------------------
#' ## Document footer
#'
#' Session Information:
#+ sessionInfo
sessionInfo()## R version 3.6.1 (2019-07-05)
## Platform: x86_64-w64-mingw32/x64 (64-bit)
## Running under: Windows 10 x64 (build 18363)
##
## Matrix products: default
##
## locale:
## [1] LC_COLLATE=English_United States.1252
## [2] LC_CTYPE=English_United States.1252
## [3] LC_MONETARY=English_United States.1252
## [4] LC_NUMERIC=C
## [5] LC_TIME=English_United States.1252
##
## attached base packages:
## [1] grid stats graphics grDevices utils datasets methods
## [8] base
##
## other attached packages:
## [1] moments_0.14 ggsn_0.5.0 rgdal_1.5-18 sp_1.4-4
## [5] maps_3.3.0 sf_0.9-6 ngram_3.0.4 GGally_2.0.0
## [9] ggnewscale_0.4.3 gridExtra_2.3 TPD_1.1.0 ks_1.11.7
## [13] data.table_1.13.2 ezknitr_0.6 ggplot2_3.3.2 knitr_1.30
##
## loaded via a namespace (and not attached):
## [1] mclust_5.4.6 Rcpp_1.0.5 mvtnorm_1.1-1
## [4] lattice_0.20-41 tidyr_1.1.2 png_0.1-7
## [7] class_7.3-17 digest_0.6.26 R6_2.4.1
## [10] plyr_1.8.6 evaluate_0.14 e1071_1.7-4
## [13] httr_1.4.2 highr_0.8 pillar_1.4.6
## [16] RgoogleMaps_1.4.5.3 rlang_0.4.8 Matrix_1.2-18
## [19] rmarkdown_2.5 splines_3.6.1 labeling_0.4.2
## [22] stringr_1.4.0 foreign_0.8-71 munsell_0.5.0
## [25] compiler_3.6.1 xfun_0.18 pkgconfig_2.0.3
## [28] mgcv_1.8-33 htmltools_0.5.0 tidyselect_1.1.0
## [31] tibble_3.0.4 viridisLite_0.3.0 reshape_0.8.8
## [34] crayon_1.3.4 dplyr_1.0.2 withr_2.3.0
## [37] bitops_1.0-6 nlme_3.1-149 gtable_0.3.0
## [40] lifecycle_0.2.0 DBI_1.1.0 magrittr_1.5
## [43] units_0.6-7 scales_1.1.1 KernSmooth_2.23-17
## [46] stringi_1.5.3 mapproj_1.2.7 farver_2.0.3
## [49] ellipsis_0.3.1 generics_0.0.2 vctrs_0.3.4
## [52] rjson_0.2.20 RColorBrewer_1.1-2 tools_3.6.1
## [55] ggmap_3.0.0 glue_1.4.2 purrr_0.3.4
## [58] jpeg_0.1-8.1 yaml_2.2.1 colorspace_1.4-1
## [61] maptools_1.0-2 classInt_0.4-3