PublicationIBT
Katherine Duchesneau
2020-07-08
#Housekeeping
##Packages
library(readr)
library(tidyverse)
library(lme4)
library(MASS)
library(MuMIn)
library(vegan)
library(fitdistrplus)
library(car)
library(bbmle)
library(arm)
library(corrplot)
library(ade4)
library(lmtest)
library(RColorBrewer)
library(plyr)
library(ggplot2)
library(gridExtra)
library(caret)
library(glmnet)Session info
SI<-sessionInfo()
print(SI, locale = FALSE)## R version 4.0.2 (2020-06-22)
## Platform: x86_64-apple-darwin17.0 (64-bit)
## Running under: macOS Catalina 10.15.6
##
## Matrix products: default
## BLAS: /Library/Frameworks/R.framework/Versions/4.0/Resources/lib/libRblas.dylib
## LAPACK: /Library/Frameworks/R.framework/Versions/4.0/Resources/lib/libRlapack.dylib
##
## attached base packages:
## [1] stats4 stats graphics grDevices utils datasets methods
## [8] base
##
## other attached packages:
## [1] glmnet_4.0-2 caret_6.0-86 gridExtra_2.3 plyr_1.8.6
## [5] RColorBrewer_1.1-2 lmtest_0.9-37 zoo_1.8-8 ade4_1.7-15
## [9] corrplot_0.84 arm_1.11-2 bbmle_1.0.23.1 car_3.0-9
## [13] carData_3.0-4 fitdistrplus_1.1-1 survival_3.1-12 vegan_2.5-6
## [17] lattice_0.20-41 permute_0.9-5 MuMIn_1.43.17 MASS_7.3-51.6
## [21] lme4_1.1-23 Matrix_1.2-18 forcats_0.5.0 stringr_1.4.0
## [25] dplyr_1.0.2 purrr_0.3.4 tidyr_1.1.1 tibble_3.0.3
## [29] ggplot2_3.3.2 tidyverse_1.3.0 readr_1.3.1
##
## loaded via a namespace (and not attached):
## [1] minqa_1.2.4 colorspace_1.4-1 class_7.3-17
## [4] ellipsis_0.3.1 rio_0.5.16 htmlTable_2.0.1
## [7] base64enc_0.1-3 fs_1.5.0 rstudioapi_0.11
## [10] prodlim_2019.11.13 fansi_0.4.1 mvtnorm_1.1-1
## [13] lubridate_1.7.9 xml2_1.3.2 codetools_0.2-16
## [16] splines_4.0.2 knitr_1.29 Formula_1.2-3
## [19] jsonlite_1.7.0 nloptr_1.2.2.2 pROC_1.16.2
## [22] broom_0.7.0 cluster_2.1.0 dbplyr_1.4.4
## [25] png_0.1-7 compiler_4.0.2 httr_1.4.2
## [28] backports_1.1.9 assertthat_0.2.1 cli_2.0.2
## [31] htmltools_0.5.0 tools_4.0.2 coda_0.19-3
## [34] gtable_0.3.0 glue_1.4.1 reshape2_1.4.4
## [37] Rcpp_1.0.5 cellranger_1.1.0 vctrs_0.3.2
## [40] nlme_3.1-148 iterators_1.0.12 timeDate_3043.102
## [43] gower_0.2.2 xfun_0.16 openxlsx_4.1.5
## [46] rvest_0.3.6 lifecycle_0.2.0 statmod_1.4.34
## [49] scales_1.1.1 ipred_0.9-9 hms_0.5.3
## [52] parallel_4.0.2 yaml_2.2.1 curl_4.3
## [55] bdsmatrix_1.3-4 rpart_4.1-15 latticeExtra_0.6-29
## [58] stringi_1.4.6 foreach_1.5.0 checkmate_2.0.0
## [61] boot_1.3-25 zip_2.1.0 shape_1.4.4
## [64] lava_1.6.7 rlang_0.4.7 pkgconfig_2.0.3
## [67] evaluate_0.14 recipes_0.1.13 htmlwidgets_1.5.1
## [70] tidyselect_1.1.0 magrittr_1.5 R6_2.4.1
## [73] generics_0.0.2 Hmisc_4.4-1 DBI_1.1.0
## [76] pillar_1.4.6 haven_2.3.1 foreign_0.8-80
## [79] withr_2.2.0 mgcv_1.8-31 abind_1.4-5
## [82] nnet_7.3-14 modelr_0.1.8 crayon_1.3.4
## [85] rmarkdown_2.3 jpeg_0.1-8.1 grid_4.0.2
## [88] readxl_1.3.1 data.table_1.13.0 blob_1.2.1
## [91] ModelMetrics_1.2.2.2 reprex_0.3.0 digest_0.6.25
## [94] numDeriv_2016.8-1.1 munsell_0.5.0Import
Abundance and richness data
Combined <- read.csv("Combined.csv")
row.names(Combined)<-Combined$Site
Combined$Site<-NULL
Combined.species<-Combined[19,]
Combined<-Combined[-c(19),]
Combined<-mutate_all(Combined, function(x) as.numeric(as.character(x)))Distance matrix
# from previous work, calculated from plot to plot.
distmat <- matrix(c(0,5.17,1.53,3.05,1.44,3.27,0.73,0.93,0.37,3.85,20.43,2.36,3.02,6.04,7.06,8.44,16.98,8.41,5.17,0,5.67,3.42,4.27,2.94,5.81,5.40,5.33,3.22,17.43,4.74,2.14,3.11,2.03,5.31,14.07,6.57,1.53,5.67,0,2.69,1.41,3.09,1.19,0.60,1.16,3.44,19.56,1.41,3.58,7.16,7.34,9.63,16.12,9.79,3.05,3.42,2.69,0,1.63,0.48,3.38,2.71,2.97,0.8,17.43,1.43,1.89,5.79,4.78,8.27,13.98,8.99,1.44,4.27,1.41,1.63,0,1.90,1.78,1.18,1.34,2.43,19.01,1.15,2.17,5.86,6.00,8.35,15.55,8.68,3.27,2.94,3.09,0.48,1.90,0,3.67,3.04,3.23,0.71,17.34,1.89,1.56,5.39,4.32,7.85,13.88,8.64,0.73,5.81,1.19,3.38,1.78,3.67,0,0.73,0.48,4.18,20.60,2.39,3.66,6.76,7.66,9.16,17.16,9.09,0.93,5.40,0.6,2.71,1.18,3.04,0.73,0,0.57,3.51,19.87,1.66,3.26,6.68,7.17,9.14,16.42,9.24,0.37,5.33,1.16,2.97,1.34,3.23,0.48,0.57,0,3.77,20.28,2.12,3.17,6.36,7.18,8.77,16.83,8.77,3.85,3.22,3.44,0.8,2.43,0.71,4.18,3.51,3.77,0,16.66,2.09,2.21,5.91,4.27,8.32,13.21,9.22,20.43,17.43,19.56,17.43,19.01,17.34,20.60,19.87,20.28,16.66,0,18.21,18.40,20.22,15.84,21.46,3.45,23.39,2.36,4.74,1.41,1.43,1.15,1.89,2.39,1.66,2.12,2.09,18.21,0,2.87,6.76,6.21,9.27,14.77,9.74,3.02,2.14,3.58,1.89,2.17,1.56,3.66,3.26,3.17,2.21,18.40,2.87,0,3.95,4.05,6.43,14.99,7.12,6.04,3.11,7.16,5.79,5.86,5.39,6.76,6.68,6.36,5.91,20.22,6.76,3.95,0,4.36,2.51,16.92,3.46,7.06,2.03,7.34,4.78,6,4.32,7.66,7.17,7.18,4.27,15.84,6.21,4.05,4.36,0,5.92,12.57,7.63,8.44,5.31,9.63,8.27,8.35,7.85,9.16,9.14,8.77,8.32,21.46,9.27,6.43,2.51,5.92,0,18.29,2.13,16.98,14.07,16.12,13.98,15.55,13.88,17.16,16.42,16.83,13.21,3.45,14.77,14.99,16.92,12.57,18.29,0,20.16,8.41,6.57,9.79,8.99,8.68,8.64,9.09,9.24,8.77,9.22,23.39,9.74,7.12,3.46,7.63,2.13,20.16,0), nrow = 18, ncol = 18, byrow = TRUE)Native and non native assignment
non.native_USDA <- read.csv("non-native_USDA.csv")Explanatory variables
Make DB
Percent_area_Zone <- read.csv("Percent_area_Zone.csv")
totalA <- aggregate(area~descriptio,Percent_area_Zone,sum)
Percent_area_Zone<-merge(Percent_area_Zone,totalA, by="descriptio")
Percent_area_Zone$Per_area<-(Percent_area_Zone$area.x/Percent_area_Zone$area.y)*100
Percent_area_Zone$Name<-NULL
Percent_area_Zone$area.x<-NULL
Percent_area_Zone$area.y<-NULL
Percent_area_Zone$perimeter<-NULL
Zone_cover<-tidyr::spread(Percent_area_Zone,key="CLASS",value= "Per_area")
Zone_cover[is.na(Zone_cover)] <- 0
colnames(Zone_cover) <- c( "Site", "Airport_cover", "Commercial_cover", "Rural_cover", "Environment_cover", "Institution_cover","Minning_cover","Park_cover","Residential_cover")
Dist_to_Roads<- read.csv("Dist_to_roads.csv")
Dist_to_Roads<-tidyr::spread(Dist_to_Roads,key="id_2",value= "distance")
colnames(Dist_to_Roads)<-c("Site","Great_Northern_Rd","HW_17","HW_75","Meeting_Point","Trunk_Road")
IBT_data_1 <- read.csv("IBT_Data.csv")
IBT_data_1$Site<-as.factor(IBT_data_1$Site)
IBT_data_year <- read.csv("Island_surrounding_build_year_kat.csv")
colnames(IBT_data_year)[2]<-"Site"
IBT_data_year$Site<-as.factor(IBT_data_year$Site)
length(IBT_data_year$YearBuilt2)## [1] 3287IBT_data_year %>% tidyr::drop_na()->IBT_data_year
length(IBT_data_year$YearBuilt2)## [1] 2984IBT_data_year %>% dplyr::group_by(Site) %>% dplyr::summarise(Newest = max(YearBuilt2), Oldest= min(YearBuilt2),Average_building_age= mean(YearBuilt2), SD_building_age = sd(YearBuilt2) ,DiffYears=length(unique(YearBuilt2)))->IBT_data_events## `summarise()` ungrouping output (override with `.groups` argument)IBT_data_1[is.na(IBT_data_1$PopulationDensity2),]$PopulationDensity2<-0
IBT_data_1$Av_PopDens<-(IBT_data_1$PopulationDensity2+IBT_data_1$PopulationDensity1)/2
IBT_data_1$PopulationDensity1<-NULL
IBT_data_1$PopulationDensity2<-NULL
IBT_data<-IBT_data_1
IBT_data<-merge(IBT_data_1,Zone_cover,by="Site")
IBT_data<-merge(IBT_data,IBT_data_events,by="Site")Var selection
colnames(IBT_data)## [1] "Site" "Name" "Site_number"
## [4] "Area_m2" "Perimeter_m" "Distance_m"
## [7] "Closest_Road" "sp_rich" "Diversity"
## [10] "Evenness" "Exotic_Native_Ratio" "Native"
## [13] "Exotic" "per_exo" "sw_div"
## [16] "per_nat" "Av_PopDens" "Airport_cover"
## [19] "Commercial_cover" "Rural_cover" "Environment_cover"
## [22] "Institution_cover" "Minning_cover" "Park_cover"
## [25] "Residential_cover" "Newest" "Oldest"
## [28] "Average_building_age" "SD_building_age" "DiffYears"IBT_data$P_A_ratio<-IBT_data$Perimeter_m/IBT_data$Area_m2
IBT_data$Site<-rownames(IBT_data)
Y<-IBT_data[c(4,5:7,17,19,20,21,25,28,31)]
covar_cor = cor(Y)
write.csv(covar_cor,"covar_cor.csv")
mydata.rcorr = Hmisc::rcorr(as.matrix(Y))
covar_cor.p = mydata.rcorr$P
corrplot.mixed(covar_cor, upper = "ellipse", lower = "number")
IBT_PCA = dudi.pca(df = Y, scannf = FALSE, nf = 2, center = TRUE, scale = TRUE)
IBT_PCA$co## Comp1 Comp2
## Area_m2 -0.83055978 -0.23118749
## Perimeter_m -0.92049970 -0.10852220
## Distance_m 0.39930461 -0.77326984
## Closest_Road -0.50651299 -0.73943756
## Av_PopDens -0.62280650 0.29502854
## Commercial_cover 0.22103993 -0.20630908
## Rural_cover -0.91910930 0.10673178
## Environment_cover -0.58573112 0.24184569
## Residential_cover 0.85690525 0.18060601
## Average_building_age -0.07559943 0.82093034
## P_A_ratio 0.86303462 0.02170061s.corcircle(IBT_PCA$co)
IBT_PCA$eig/sum(IBT_PCA$eig)## [1] 0.4602335361 0.1924037486 0.1037752140 0.0771352020 0.0629644116
## [6] 0.0326285525 0.0274360411 0.0240211060 0.0107620722 0.0079182427
## [11] 0.0007218732# Notes about what is correlated.
# area and perimeter -> use area only
# Area and rural cover -> only use area because it is a critical component.
# important thing is that the area of an island is related to the rural cover, so more likely to be rural.
TRY<-IBT_dataH1
Size and diversity model
Large forest islands close to the “forest mainland” will have greater species richness as compared to small forest islands far from the "forest mainland
Linear Model Selection
Model
car::qqp(IBT_data$sp_rich, "norm")
## [1] 5 3plot(density(IBT_data$sp_rich))
shapiro.test(IBT_data$sp_rich)##
## Shapiro-Wilk normality test
##
## data: IBT_data$sp_rich
## W = 0.91968, p-value = 0.1276model_sprich <- lm(sp_rich~ Area_m2 +
Distance_m +
P_A_ratio+
Closest_Road+
Commercial_cover+
Residential_cover+
Av_PopDens,
data = TRY)
plot(model_sprich)
## Warning in sqrt(crit * p * (1 - hh)/hh): NaNs produced
## Warning in sqrt(crit * p * (1 - hh)/hh): NaNs produced
olsrr::ols_test_breusch_pagan(model_sprich)##
## Breusch Pagan Test for Heteroskedasticity
## -----------------------------------------
## Ho: the variance is constant
## Ha: the variance is not constant
##
## Data
## -----------------------------------
## Response : sp_rich
## Variables: fitted values of sp_rich
##
## Test Summary
## ----------------------------
## DF = 1
## Chi2 = 0.5169307
## Prob > Chi2 = 0.4721539vif(model_sprich)## Area_m2 Distance_m P_A_ratio Closest_Road
## 2.541593 2.812171 3.154475 4.540278
## Commercial_cover Residential_cover Av_PopDens
## 1.147406 4.693327 1.697393####MuMIn
options(na.action="na.fail")
modAll1<-lm(sp_rich~ Area_m2 +Distance_m +
Commercial_cover+
Residential_cover+
Closest_Road+
Av_PopDens+
Average_building_age,
data = TRY)
modAll.d<-dredge(modAll1)## Fixed term is "(Intercept)"modAll.d## Global model call: lm(formula = sp_rich ~ Area_m2 + Distance_m + Commercial_cover +
## Residential_cover + Closest_Road + Av_PopDens + Average_building_age,
## data = TRY)
## ---
## Model selection table
## (Int) Are_m2 Av_PpD Avr_bld_age Cls_Rod Cmm_cvr Dst_m
## 37 -1917.00 0.98640 0.0052780
## 1 53.78
## 33 40.40 0.0031560
## 17 54.97 -0.4833
## 2 56.45 -1.191e-05
## 49 41.26 -0.5138 0.0032520
## 65 48.22
## 35 26.83 0.027930 0.0037700
## 53 -1685.00 0.86970 -0.3119 0.0050850
## 10 53.90 -2.293e-05 2.580e-03
## 39 -1763.00 0.018130 0.90460 0.0055000
## 5 -571.80 0.31680
## 9 52.02 9.045e-04
## 34 43.66 -9.057e-06 0.0028670
## 38 -1824.00 -6.394e-06 0.94090 0.0049760
## 3 48.60 0.013170
## 45 -2023.00 1.04000 6.215e-04 0.0050280
## 73 37.82 2.489e-03
## 14 -1609.00 -2.871e-05 0.84140 4.052e-03
## 97 38.30 0.0028480
## 101 -1924.00 0.99020 0.0052970
## 18 58.26 -1.384e-05 -0.5599
## 41 40.35 -1.067e-04 0.0032180
## 67 28.03 0.039350
## 75 10.80 0.048790 2.965e-03
## 99 11.03 0.049180 0.0034370
## 81 49.10 -0.5091
## 4 47.71 -1.495e-05 0.023960
## 50 45.31 -1.098e-05 -0.5714 0.0029120
## 36 27.37 -1.300e-05 0.036410 0.0035420
## 26 55.71 -2.430e-05 2.477e-03 -0.5279
## 6 -543.10 -1.173e-05 0.30360
## 25 53.49 7.278e-04 -0.4568
## 13 -1001.00 0.53240 1.568e-03
## 77 -1240.00 0.64560 3.442e-03
## 66 52.76 -8.558e-06
## 21 -316.10 0.18780 -0.4359
## 51 30.97 0.020870 -0.4240 0.0036940
## 19 52.59 0.005891 -0.4568
## 69 -456.70 0.25590
## 74 43.82 -1.707e-05 3.224e-03
## 113 38.98 -0.5269 0.0029170
## 46 -2083.00 -1.892e-05 1.07500 2.574e-03 0.0033520
## 12 45.60 -2.566e-05 0.022860 2.541e-03
## 89 39.38 2.303e-03 -0.4477
## 57 41.08 -3.965e-04 -0.5304 0.0034880
## 11 48.14 0.010410 7.959e-04
## 43 23.90 0.033040 -8.095e-04 0.0043570
## 7 -475.70 0.008291 0.26650
## 42 47.58 -1.743e-05 1.670e-03 0.0016190
## 78 -1635.00 -2.278e-05 0.84910 4.723e-03
## 76 16.49 -1.760e-05 0.049710 3.731e-03
## 40 -1563.00 -9.534e-06 0.025440 0.80380 0.0051400
## 54 -1522.00 -8.077e-06 0.78940 -0.3729 0.0046660
## 105 36.43 1.526e-03 0.0015560
## 98 42.67 -7.993e-06 0.0028120
## 83 32.69 0.031580 -0.3832
## 55 -1599.00 0.014690 0.82300 -0.2595 0.0052980
## 103 -1394.00 0.030330 0.71400 0.0049770
## 20 51.99 -1.573e-05 0.016620 -0.4956
## 30 -1332.00 -2.861e-05 0.70180 3.743e-03 -0.3288
## 61 -1768.00 0.91170 3.751e-04 -0.2864 0.0049500
## 82 54.91 -1.079e-05 -0.5546
## 102 -1975.00 -9.583e-06 1.01900 0.0053240
## 68 32.60 -7.927e-06 0.038660
## 117 -1642.00 0.84820 -0.3189 0.0049870
## 22 -232.50 -1.362e-05 0.14720 -0.5216
## 47 -1807.00 0.016700 0.92730 1.873e-04 0.0054080
## 109 -1908.00 0.98120 1.413e-03 0.0040780
## 71 145.10 0.041300 -0.05982
## 115 15.65 0.041680 -0.3636 0.0033950
## 52 31.90 -1.377e-05 0.029260 -0.4602 0.0034450
## 16 -1514.00 -2.907e-05 0.005903 0.79190 3.955e-03
## 107 8.07 0.050580 1.775e-03 0.0019510
## 85 -160.80 0.10630 -0.4815
## 90 46.16 -1.874e-05 3.085e-03 -0.5057
## 79 -654.70 0.039050 0.33900 3.371e-03
## 91 14.97 0.042920 2.799e-03 -0.2633
## 100 15.31 -7.085e-06 0.048440 0.0033970
## 8 -301.30 -1.438e-05 0.020300 0.17740
## 29 -716.90 0.38940 1.259e-03 -0.3394
## 58 48.62 -1.813e-05 1.441e-03 -0.5485 0.0018330
## 70 -498.20 -9.133e-06 0.27940
## 28 49.70 -2.605e-05 0.015990 2.461e-03 -0.4662
## 27 52.05 0.003723 6.949e-04 -0.4413
## 15 -1047.00 -0.002676 0.55660 1.626e-03
## 44 31.08 -1.753e-05 0.033230 9.734e-04 0.0027550
## 23 -284.50 0.003252 0.17120 -0.4255
## 114 44.65 -1.027e-05 -0.5700 0.0028760
## 93 -992.80 0.52120 3.138e-03 -0.2898
## 59 27.83 0.026360 -9.069e-04 -0.4385 0.0043500
## 121 37.70 1.010e-03 -0.4947 0.0020570
## 106 43.65 -1.685e-05 3.148e-03 0.0001061
## 80 -1108.00 -2.128e-05 0.033410 0.57330 4.577e-03
## 62 -1819.00 -1.911e-05 0.94140 2.338e-03 -0.2975 0.0032530
## 48 -1873.00 -1.882e-05 0.016240 0.96460 2.141e-03 0.0037300
## 84 38.92 -9.811e-06 0.029740 -0.4318
## 110 -1987.00 -1.859e-05 1.02500 3.198e-03 0.0025910
## 92 21.91 -1.859e-05 0.042570 3.570e-03 -0.3224
## 94 -1381.00 -2.291e-05 0.72140 4.415e-03 -0.3001
## 56 -1345.00 -1.053e-05 0.022070 0.69480 -0.3128 0.0048580
## 87 313.50 0.035830 -0.14340 -0.4035
## 24 -70.26 -1.551e-05 0.015520 0.06208 -0.4837
## 108 15.40 -1.651e-05 0.050130 3.362e-03 0.0005267
## 86 -184.50 -1.092e-05 0.12130 -0.5235
## 111 -1340.00 0.032360 0.68530 1.606e-03 0.0035700
## 118 -1663.00 -1.078e-05 0.86160 -0.3608 0.0049770
## 104 -1476.00 -8.647e-06 0.028260 0.75850 0.0050230
## 72 83.38 -7.860e-06 0.039510 -0.02598
## 119 -1212.00 0.027210 0.62340 -0.2674 0.0047500
## 63 -1605.00 0.014540 0.82580 2.075e-05 -0.2586 0.0052880
## 125 -1663.00 0.85820 1.135e-03 -0.2803 0.0040450
## 116 21.58 -8.879e-06 0.039850 -0.4080 0.0033390
## 32 -1278.00 -2.883e-05 0.003682 0.67360 3.689e-03 -0.3226
## 123 12.60 0.043910 1.419e-03 -0.3096 0.0022120
## 31 -797.50 -0.004995 0.43120 1.360e-03 -0.3477
## 60 35.40 -1.809e-05 0.026270 9.290e-04 -0.4568 0.0026960
## 95 -516.20 0.036350 0.27010 3.155e-03 -0.2097
## 122 45.29 -1.761e-05 2.684e-03 -0.5151 0.0005629
## 88 256.10 -9.587e-06 0.033070 -0.11100 -0.4465
## 112 -1453.00 -1.791e-05 0.030240 0.74680 3.314e-03 0.0021700
## 96 -959.40 -2.152e-05 0.030360 0.49960 4.352e-03 -0.2326
## 64 -1662.00 -1.901e-05 0.013980 0.85870 1.987e-03 -0.2708 0.0035880
## 126 -1733.00 -1.880e-05 0.89710 2.928e-03 -0.2923 0.0025400
## 124 20.52 -1.705e-05 0.042940 3.031e-03 -0.3335 0.0007621
## 120 -1276.00 -9.810e-06 0.024360 0.65910 -0.3109 0.0047660
## 127 -1200.00 0.029540 0.61610 1.375e-03 -0.2162 0.0035890
## 128 -1304.00 -1.815e-05 0.027170 0.67310 3.087e-03 -0.2330 0.0021730
## Rsd_cvr df logLik AICc delta weight
## 37 4 -73.917 158.9 0.00 0.116
## 1 2 -77.482 159.8 0.85 0.075
## 33 3 -76.037 159.8 0.88 0.075
## 17 3 -76.837 161.4 2.48 0.034
## 2 3 -76.852 161.4 2.51 0.033
## 49 4 -75.171 161.4 2.51 0.033
## 65 0.124700 3 -76.974 161.7 2.75 0.029
## 35 4 -75.384 161.8 2.93 0.027
## 53 5 -73.559 162.1 3.21 0.023
## 10 4 -75.527 162.1 3.22 0.023
## 39 5 -73.592 162.2 3.27 0.023
## 5 3 -77.251 162.2 3.31 0.022
## 9 3 -77.275 162.3 3.35 0.022
## 34 4 -75.628 162.3 3.42 0.021
## 38 5 -73.667 162.3 3.42 0.021
## 3 3 -77.351 162.4 3.50 0.020
## 45 5 -73.807 162.6 3.70 0.018
## 73 0.249300 4 -75.777 162.6 3.72 0.018
## 14 5 -73.844 162.7 3.78 0.017
## 97 0.076330 4 -75.832 162.7 3.83 0.017
## 101 -0.002674 5 -73.917 162.8 3.92 0.016
## 18 4 -75.925 162.9 4.02 0.016
## 41 4 -76.034 163.1 4.23 0.014
## 67 0.230400 4 -76.073 163.2 4.31 0.013
## 75 0.404300 5 -74.186 163.4 4.46 0.012
## 99 0.198500 5 -74.230 163.5 4.55 0.012
## 81 0.133000 4 -76.213 163.5 4.59 0.012
## 4 4 -76.420 163.9 5.01 0.009
## 50 5 -74.511 164.0 5.11 0.009
## 36 5 -74.520 164.0 5.13 0.009
## 26 5 -74.574 164.1 5.24 0.008
## 6 4 -76.624 164.3 5.41 0.008
## 25 4 -76.695 164.5 5.55 0.007
## 13 4 -76.711 164.5 5.59 0.007
## 77 0.272600 5 -74.790 164.6 5.67 0.007
## 66 0.065750 4 -76.755 164.6 5.68 0.007
## 21 4 -76.756 164.6 5.68 0.007
## 51 5 -74.804 164.6 5.70 0.007
## 19 4 -76.810 164.7 5.79 0.006
## 69 0.115000 4 -76.818 164.7 5.80 0.006
## 74 0.168500 5 -74.863 164.7 5.81 0.006
## 113 0.083760 5 -74.899 164.8 5.89 0.006
## 46 6 -72.643 164.9 6.01 0.006
## 12 5 -75.072 165.1 6.23 0.005
## 89 0.247300 5 -75.116 165.2 6.32 0.005
## 57 5 -75.130 165.3 6.35 0.005
## 11 4 -77.194 165.5 6.55 0.004
## 43 5 -75.236 165.5 6.56 0.004
## 7 4 -77.204 165.5 6.57 0.004
## 42 5 -75.279 165.6 6.65 0.004
## 78 0.172200 6 -73.001 165.6 6.73 0.004
## 76 0.323900 6 -73.011 165.7 6.75 0.004
## 40 6 -73.061 165.8 6.85 0.004
## 54 6 -73.153 165.9 7.03 0.003
## 105 0.174700 5 -75.602 166.2 7.29 0.003
## 98 0.021920 5 -75.616 166.2 7.32 0.003
## 83 0.215800 5 -75.643 166.3 7.37 0.003
## 55 6 -73.348 166.3 7.42 0.003
## 103 0.094720 6 -73.377 166.4 7.48 0.003
## 20 5 -75.711 166.4 7.51 0.003
## 30 6 -73.464 166.6 7.65 0.003
## 61 6 -73.519 166.7 7.76 0.002
## 82 0.059420 5 -75.838 166.7 7.76 0.002
## 102 -0.070190 6 -73.531 166.7 7.79 0.002
## 68 0.174000 5 -75.866 166.7 7.82 0.002
## 117 0.013150 6 -73.552 166.7 7.83 0.002
## 22 5 -75.871 166.7 7.83 0.002
## 47 6 -73.584 166.8 7.89 0.002
## 109 0.089120 6 -73.673 167.0 8.07 0.002
## 71 0.237900 5 -76.066 167.1 8.22 0.002
## 115 0.185000 6 -73.754 167.1 8.23 0.002
## 52 6 -73.764 167.2 8.25 0.002
## 16 6 -73.815 167.3 8.35 0.002
## 107 0.316400 6 -73.856 167.3 8.44 0.002
## 85 0.128500 5 -76.186 167.4 8.46 0.002
## 90 0.158400 6 -73.925 167.5 8.57 0.002
## 79 0.385500 6 -73.944 167.5 8.61 0.002
## 91 0.384400 6 -73.945 167.5 8.61 0.002
## 100 0.148500 6 -74.028 167.7 8.78 0.001
## 8 5 -76.350 167.7 8.79 0.001
## 29 5 -76.418 167.8 8.92 0.001
## 58 6 -74.220 168.1 9.16 0.001
## 70 0.051230 5 -76.566 168.1 9.22 0.001
## 28 6 -74.343 168.3 9.41 0.001
## 27 5 -76.684 168.4 9.46 0.001
## 15 5 -76.707 168.4 9.50 0.001
## 44 6 -74.397 168.4 9.52 0.001
## 23 5 -76.749 168.5 9.59 0.001
## 114 0.014420 6 -74.505 168.6 9.74 0.001
## 93 0.266800 6 -74.527 168.7 9.78 0.001
## 59 6 -74.606 168.8 9.94 0.001
## 121 0.148400 6 -74.791 169.2 10.31 0.001
## 106 0.164500 6 -74.862 169.4 10.45 0.001
## 80 0.275400 7 -72.259 169.7 10.81 0.001
## 62 7 -72.289 169.8 10.87 0.001
## 48 7 -72.403 170.0 11.09 0.000
## 84 0.144100 6 -75.316 170.3 11.36 0.000
## 110 0.074060 7 -72.538 170.3 11.37 0.000
## 92 0.295100 7 -72.599 170.4 11.49 0.000
## 94 0.165600 7 -72.655 170.5 11.60 0.000
## 56 7 -72.688 170.6 11.66 0.000
## 87 0.233000 6 -75.602 170.8 11.93 0.000
## 24 6 -75.702 171.0 12.13 0.000
## 108 0.305100 7 -72.988 171.2 12.26 0.000
## 86 0.053470 6 -75.802 171.2 12.33 0.000
## 111 0.205600 7 -73.042 171.3 12.37 0.000
## 118 -0.060730 7 -73.046 171.3 12.38 0.000
## 104 0.027170 7 -73.047 171.3 12.38 0.000
## 72 0.177700 6 -75.865 171.4 12.45 0.000
## 119 0.097970 7 -73.111 171.4 12.51 0.000
## 63 7 -73.348 171.9 12.98 0.000
## 125 0.084970 7 -73.394 172.0 13.08 0.000
## 116 0.120600 7 -73.424 172.0 13.14 0.000
## 32 7 -73.452 172.1 13.19 0.000
## 123 0.281300 7 -73.515 172.2 13.32 0.000
## 31 6 -76.402 172.4 13.53 0.000
## 60 7 -73.642 172.5 13.57 0.000
## 95 0.373600 7 -73.798 172.8 13.88 0.000
## 122 0.136700 7 -73.902 173.0 14.09 0.000
## 88 0.159000 7 -75.291 175.8 16.87 0.000
## 112 0.183400 8 -71.916 175.8 16.92 0.000
## 96 0.260900 8 -72.041 176.1 17.17 0.000
## 64 8 -72.108 176.2 17.30 0.000
## 126 0.069560 8 -72.193 176.4 17.48 0.000
## 124 0.267000 8 -72.551 177.1 18.19 0.000
## 120 0.021860 8 -72.679 177.4 18.45 0.000
## 127 0.192200 8 -72.870 177.7 18.83 0.000
## 128 0.168800 9 -71.689 183.9 24.97 0.000
## Models ranked by AICc(x)dim(modAll.d)## [1] 128 13dim(subset(modAll.d, cumsum(weight) <= .95) )## [1] 63 13avgmod1.95p <- model.avg(modAll.d, cumsum(weight) <= .95)
confint(avgmod1.95p)## 2.5 % 97.5 %
## (Intercept) -2.585186e+03 1.548056e+03
## Average_building_age -3.599484e-01 2.027279e+00
## Distance_m -6.624450e-04 8.989557e-03
## Commercial_cover -1.410364e+00 4.747541e-01
## Area_m2 -4.399938e-05 1.349979e-05
## Residential_cover -2.087168e-01 5.346398e-01
## Av_PopDens -3.778130e-02 9.203323e-02
## Closest_Road -2.510064e-03 6.465953e-03DF<-data.frame(TRY$Distance_m, TRY$Commercial_cover, TRY$Residential_cover, TRY$Closest_Road, TRY$Average_building_age)
hier.part::hier.part(TRY$sp_rich,DF,"gaussian")
## $gfs
## [1] 0.00000000 0.14836012 0.06925920 0.05496087 0.02278603 0.02535045
## [7] 0.22650228 0.16755712 0.14861869 0.32706677 0.13158517 0.08380366
## [13] 0.07750842 0.17261738 0.07117181 0.08210893 0.24956694 0.22999090
## [19] 0.35335397 0.18859200 0.32708764 0.33530767 0.23120226 0.13415844
## [25] 0.11151621 0.25856709 0.25848024 0.35384530 0.35618003 0.34511032
## [31] 0.27992569 0.36508136
##
## $IJ
## I J Total
## TRY.Distance_m 0.14598463 0.002375490 0.14836012
## TRY.Commercial_cover 0.04954193 0.019717270 0.06925920
## TRY.Residential_cover 0.05177638 0.003184484 0.05496087
## TRY.Closest_Road 0.03683461 -0.014048573 0.02278603
## TRY.Average_building_age 0.08094382 -0.055593371 0.02535045
##
## $I.perc
## ind.exp.var
## TRY.Distance_m 39.98687
## TRY.Commercial_cover 13.57011
## TRY.Residential_cover 14.18215
## TRY.Closest_Road 10.08942
## TRY.Average_building_age 22.17145
##
## $params
## $params$full.model
## [1] "y ~ TRY.Distance_m + TRY.Commercial_cover + TRY.Residential_cover + TRY.Closest_Road + TRY.Average_building_age"
##
## $params$family
## [1] "gaussian"
##
## $params$link
## [1] "default"
##
## $params$gof
## [1] "Rsqu"model_sprichmBest<-lm(sp_rich~ Area_m2 +Distance_m +Average_building_age,data = TRY)Plots
model_sprich<-update(model_sprichmBest,.~.+Area_m2)
(Area_AND_richness_plot<-jtools::effect_plot(model_sprich, pred=Area_m2, interval = TRUE,
int.type = "confidence",data = IBT_data)+
geom_smooth(color=c('red'))+
geom_point(aes(Area_m2,sp_rich),TRY)+
labs(x="\nArea (m^2)", y="Species richness")+
theme_classic() +
theme(
axis.text = element_text(family="Times",colour = "black"),
axis.title.x = element_text(family="Times",size=22),
axis.text.x = element_text(family="Times",size=18,angle=90),
axis.title.y = element_text(family="Times",size=22, angle=90),
axis.text.y = element_text(family="Times",size=18),
axis.ticks = element_blank(),
panel.background = element_rect(fill="white"),
panel.border = element_rect(colour = "black", fill=NA, size=1),
plot.title=element_text(face="bold", size=24)
))## `geom_smooth()` using method = 'loess' and formula 'y ~ x'
(Dist_AND_richness_plot<-jtools::effect_plot(model_sprich, pred=Distance_m, interval = TRUE,
int.type = "confidence",data = IBT_data)+
geom_smooth(color=c('red'))+
geom_point(aes(Distance_m,sp_rich),TRY)+
labs(x="\nDistance (m)", y="Species richness")+
theme_classic() +
theme(
axis.text = element_text(family="Times",colour = "black"),
axis.title.x = element_text(family="Times",size=22),
axis.text.x = element_text(family="Times",size=18,angle=90),
axis.title.y = element_text(family="Times",size=22, angle=90),
axis.text.y = element_text(family="Times",size=18),
axis.ticks = element_blank(),
panel.background = element_rect(fill="white"),
panel.border = element_rect(colour = "black", fill=NA, size=1),
plot.title=element_text(face="bold", size=24)
))## `geom_smooth()` using method = 'loess' and formula 'y ~ x'
(Controlled_richness_plots<-grid.arrange(Area_AND_richness_plot,Dist_AND_richness_plot, nrow = 1))## `geom_smooth()` using method = 'loess' and formula 'y ~ x'
## `geom_smooth()` using method = 'loess' and formula 'y ~ x'
## TableGrob (1 x 2) "arrange": 2 grobs
## z cells name grob
## 1 1 (1-1,1-1) arrange gtable[layout]
## 2 2 (1-1,2-2) arrange gtable[layout]ggsave("Controlled_richness_plots.pdf", Controlled_richness_plots,width = 10, height = 5)H2
Urban forest islands that are closer together will be more similar in species composition than those farther apart.
Mantel Corr
Overall
simmat<- vegdist(Combined, "bray",diag=T)
mantel(simmat, distmat)##
## Mantel statistic based on Pearson's product-moment correlation
##
## Call:
## mantel(xdis = simmat, ydis = distmat)
##
## Mantel statistic r: 0.0316
## Significance: 0.37
##
## Upper quantiles of permutations (null model):
## 90% 95% 97.5% 99%
## 0.127 0.164 0.189 0.221
## Permutation: free
## Number of permutations: 999dis.mancor <- mantel.correlog(D.eco = simmat, D.geo = distmat, nperm = 1000)
(mancor_plot<-plot(dis.mancor))
## NULLinvasive
Combined_invasive <-data.frame(c(1:18))
Combined_native <-data.frame(c(1:18))
V <- vector()
c<-1
d<-1
non.native_USDA$Scientific.Name <- as.character(non.native_USDA$Scientific.Name)
for (i in 1:length(non.native_USDA$Scientific.Name))for(j in 1:length(names(Combined.species))) if(non.native_USDA$Scientific.Name[i]==Combined.species[1,j]){
Combined_invasive[,c] <- c(1:18)
Combined_invasive[,c] <- Combined[,j]
V[c] <- j
names(Combined_invasive)[c]<-non.native_USDA$Scientific.Name[i]
c<-c+1
}
Combined_native <- Combined[,-c(unique(V))]
rowSums(Combined_invasive)## [1] 27 47 289 104 40 317 240 40 11 129 0 40 57 4 48 28 48 0Combined_invasive2<-Combined_invasive[-c(11,18),]
distmat.inv<-distmat[-c(11,18), -c(11,18)]
simmat.inv<- vegdist(Combined_invasive2, "bray",diag=T)
mantel(simmat.inv,distmat.inv)##
## Mantel statistic based on Pearson's product-moment correlation
##
## Call:
## mantel(xdis = simmat.inv, ydis = distmat.inv)
##
## Mantel statistic r: 0.3133
## Significance: 0.003
##
## Upper quantiles of permutations (null model):
## 90% 95% 97.5% 99%
## 0.156 0.198 0.223 0.264
## Permutation: free
## Number of permutations: 999inv.mancor <- mantel.correlog(D.eco = simmat.inv, D.geo = distmat.inv, nperm = 1000)
summary(inv.mancor)## Length Class Mode
## mantel.res 40 -none- numeric
## n.class 1 -none- numeric
## break.pts 9 -none- numeric
## mult 1 -none- character
## n.tests 1 -none- numeric
## call 4 -none- callinv.mancor##
## Mantel Correlogram Analysis
##
## Call:
##
## mantel.correlog(D.eco = simmat.inv, D.geo = distmat.inv, nperm = 1000)
##
## class.index n.dist Mantel.cor Pr(Mantel) Pr(corrected)
## D.cl.1 1.490000 74.000000 0.393146 0.0020 0.001998 **
## D.cl.2 3.730000 66.000000 -0.018273 0.4665 0.466533
## D.cl.3 5.970000 40.000000 -0.249692 0.0020 0.005994 **
## D.cl.4 8.210000 28.000000 -0.148346 0.0709 0.141858
## D.cl.5 10.450000 2.000000 NA NA NA
## D.cl.6 12.690000 4.000000 NA NA NA
## D.cl.7 14.930000 12.000000 NA NA NA
## D.cl.8 17.170000 14.000000 NA NA NA
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1inv.mancor_plot <-plot(inv.mancor)
native
Combined_native2<-Combined_native[-c(11,18),]
simmat.nat<- vegdist(Combined_native2, "bray",diag=T)
mantel(simmat.nat, distmat.inv)##
## Mantel statistic based on Pearson's product-moment correlation
##
## Call:
## mantel(xdis = simmat.nat, ydis = distmat.inv)
##
## Mantel statistic r: 0.001443
## Significance: 0.45
##
## Upper quantiles of permutations (null model):
## 90% 95% 97.5% 99%
## 0.129 0.155 0.186 0.209
## Permutation: free
## Number of permutations: 999nat.mancor <- mantel.correlog(D.eco = simmat.nat, D.geo = distmat.inv, nperm = 1000)
nat.mancor##
## Mantel Correlogram Analysis
##
## Call:
##
## mantel.correlog(D.eco = simmat.nat, D.geo = distmat.inv, nperm = 1000)
##
## class.index n.dist Mantel.cor Pr(Mantel) Pr(corrected)
## D.cl.1 1.490000 74.000000 0.135157 0.0769 0.07692 .
## D.cl.2 3.730000 66.000000 -0.159627 0.0420 0.08392 .
## D.cl.3 5.970000 40.000000 0.019671 0.3816 0.38162
## D.cl.4 8.210000 28.000000 -0.072381 0.2268 0.45355
## D.cl.5 10.450000 2.000000 NA NA NA
## D.cl.6 12.690000 4.000000 NA NA NA
## D.cl.7 14.930000 12.000000 NA NA NA
## D.cl.8 17.170000 14.000000 NA NA NA
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1nat.mancor_plot <-plot(nat.mancor)
Overall subsamp
Combined_sub<-Combined[-c(11,18),]
simmat<- vegdist(Combined_sub, "bray",diag=T)
mantel(simmat, distmat.inv)##
## Mantel statistic based on Pearson's product-moment correlation
##
## Call:
## mantel(xdis = simmat, ydis = distmat.inv)
##
## Mantel statistic r: 0.1024
## Significance: 0.158
##
## Upper quantiles of permutations (null model):
## 90% 95% 97.5% 99%
## 0.132 0.165 0.196 0.230
## Permutation: free
## Number of permutations: 999dis.mancor <- mantel.correlog(D.eco = simmat, D.geo = distmat.inv, nperm = 1000)
summary(dis.mancor)## Length Class Mode
## mantel.res 40 -none- numeric
## n.class 1 -none- numeric
## break.pts 9 -none- numeric
## mult 1 -none- character
## n.tests 1 -none- numeric
## call 4 -none- calldis.mancor_plot <-plot(dis.mancor)
###Plot
#create a dataframe for the correlograms - all
bs.df <- data.frame(dis.mancor$mantel.res)
bs.df$Method <- "Bray"
bs.df$Species <- "All species"
#create a dataframe for the correlograms - invasive
ls.df <- data.frame(inv.mancor$mantel.res)
ls.df$Method <- "Bray"
ls.df$Species <- "Invasive species"
#create a dataframe for the correlograms - native
rs.df <- data.frame(nat.mancor$mantel.res)
rs.df$Method <- "Bray"
rs.df$Species <- "Native species"
#make a single data frame
man.df <- rbind(bs.df, ls.df,rs.df)
#make a significance factor
man.df$Significance <- ifelse(man.df$Pr.corrected. < 0.05, "Significant", "Insignificant")
#make a distance class factor
man.df$Distance.Class <- rownames(man.df)
#plot the correlograms
(mantel.cor.plot<-ggplot(man.df, aes(class.index, Mantel.cor, group=Species, shape=Species, col = Species)) +
geom_hline(yintercept=0.00, color="grey") +
geom_line() +
geom_point(aes(fill=Significance, size=3)) +
scale_shape_manual(values = c(21, 22, 24)) +
scale_fill_manual(values=c("white", "black"))+
scale_color_manual(values=c("darkorchid3", "red", "royalblue3"))+
ylab("Mantel Correlation") +
xlab("Distance Class") +
xlim(0, 10)+
theme_classic() +
theme(
axis.text = element_text(family="Times",colour = "black"),
axis.title.x = element_text(family="Times",size=22),
axis.text.x = element_text(family="Times",size=18,angle=90),
axis.title.y = element_text(family="Times",size=22, angle=90),
axis.text.y = element_text(family="Times",size=18),
axis.ticks = element_blank(),
panel.background = element_rect(fill="white"),
panel.border = element_rect(colour = "black", fill=NA, size=1),
plot.title=element_text(face="bold", size=24)
))## Warning: Removed 12 row(s) containing missing values (geom_path).## Warning: Removed 12 rows containing missing values (geom_point).
ggsave("mantel_cor_plot.pdf", plot=mantel.cor.plot, height = 5, width =7)## Warning: Removed 12 row(s) containing missing values (geom_path).
## Warning: Removed 12 rows containing missing values (geom_point).H3
Small forest islands will have a greater number of exotic species as compared to large forest islands in part due to their less species rich native communities and the higher level of disturbances in small fragments.
Model Exotic:Native to space
FIT
bin <- cbind(TRY$Exotic, TRY$Native)
modelExo<- glm(bin ~ Area_m2 +
Distance_m +
P_A_ratio+
Closest_Road+
Commercial_cover+
Residential_cover+
Rural_cover,
data = TRY, family = "binomial")
plot(modelExo)
## Warning in sqrt(crit * p * (1 - hh)/hh): NaNs produced
## Warning in sqrt(crit * p * (1 - hh)/hh): NaNs produced
lmtest::bptest(modelExo)##
## studentized Breusch-Pagan test
##
## data: modelExo
## BP = 5.0857, df = 7, p-value = 0.6495vif(modelExo)## Area_m2 Distance_m P_A_ratio Closest_Road
## 4.929727 4.118133 2.288017 7.096241
## Commercial_cover Residential_cover Rural_cover
## 1.020422 4.779306 6.871009###MuMI
modAll2<-lm(Exotic_Native_Ratio~ Area_m2 +Distance_m +
Commercial_cover+
Residential_cover+
Closest_Road+
Average_building_age,
data = TRY)
modAll2.d<-dredge(modAll2)## Fixed term is "(Intercept)"dim(modAll2.d)## [1] 64 12dim(subset(modAll2.d, cumsum(weight) <= .95) )## [1] 9 12avgmod2.95p <- model.avg(modAll2.d, cumsum(weight) <= .95)
confint(avgmod2.95p)## 2.5 % 97.5 %
## (Intercept) -2.662714e+01 -5.513492e+00
## Average_building_age 3.111891e-03 1.331171e-02
## Closest_Road 1.231977e-05 5.104410e-05
## Residential_cover 1.480596e-03 5.129473e-03
## Distance_m -1.164063e-05 6.122029e-05
## Area_m2 -1.886306e-07 7.479477e-08
## Commercial_cover -2.852716e-03 5.548209e-03DF<-data.frame( TRY$Residential_cover, TRY$Closest_Road, TRY$Average_building_age)
hier.part::hier.part(TRY$Exotic_Native_Ratio,DF,"gaussian")
## $gfs
## [1] 0.000000000 0.341997886 0.000143118 0.124155367 0.513458976 0.416164966
## [7] 0.154220493 0.747625544
##
## $IJ
## I J Total
## TRY.Residential_cover 0.4460219 -0.10402400 0.341997886
## TRY.Closest_Road 0.1441223 -0.14397915 0.000143118
## TRY.Average_building_age 0.1574814 -0.03332602 0.124155367
##
## $I.perc
## ind.exp.var
## TRY.Residential_cover 59.65846
## TRY.Closest_Road 19.27733
## TRY.Average_building_age 21.06421
##
## $params
## $params$full.model
## [1] "y ~ TRY.Residential_cover + TRY.Closest_Road + TRY.Average_building_age"
##
## $params$family
## [1] "gaussian"
##
## $params$link
## [1] "default"
##
## $params$gof
## [1] "Rsqu"modAll.lm<-lm(Exotic_Native_Ratio~Residential_cover+
Closest_Road+
Average_building_age,
data = TRY)
summary(modAll.lm)$adj.r.squared ## [1] 0.6935453Plot Exotic:Native to space
modelExo.signif<-update(modAll.lm,.~.+Area_m2+Distance_m)
summary(modelExo.signif)##
## Call:
## lm(formula = Exotic_Native_Ratio ~ Residential_cover + Closest_Road +
## Average_building_age + Area_m2 + Distance_m, data = TRY)
##
## Residuals:
## Min 1Q Median 3Q Max
## -0.095916 -0.055393 0.001226 0.031492 0.125740
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -1.800e+01 4.771e+00 -3.773 0.00266 **
## Residential_cover 2.814e-03 8.738e-04 3.221 0.00735 **
## Closest_Road 2.822e-05 1.194e-05 2.363 0.03586 *
## Average_building_age 9.099e-03 2.408e-03 3.778 0.00263 **
## Area_m2 -3.733e-08 6.469e-08 -0.577 0.57452
## Distance_m 1.355e-05 1.439e-05 0.941 0.36508
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.07366 on 12 degrees of freedom
## Multiple R-squared: 0.7813, Adjusted R-squared: 0.6902
## F-statistic: 8.575 on 5 and 12 DF, p-value: 0.001175summary.lm(modelExo.signif)$sigma## [1] 0.07365543confint(modelExo.signif)## 2.5 % 97.5 %
## (Intercept) -2.839191e+01 -7.603313e+00
## Residential_cover 9.103469e-04 4.718064e-03
## Closest_Road 2.198953e-06 5.424194e-05
## Average_building_age 3.851541e-03 1.434554e-02
## Area_m2 -1.782785e-07 1.036112e-07
## Distance_m -1.781131e-05 4.491043e-05(Area_AND_propEXO_plot<-jtools::effect_plot(modelExo.signif, pred=Area_m2, interval = TRUE,
int.type = "confidence",data = TRY)+
geom_smooth(color=c('red'), formula = y~x)+
geom_point(aes(Area_m2,Exotic_Native_Ratio),TRY)+
labs(x="\nArea (m^2)", y="Proportion exotic to native")+
ylim(0,0.6)+
theme_classic() +
theme(
axis.text = element_text(family="Times",colour = "black"),
axis.title.x = element_text(family="Times",size=22),
axis.text.x = element_text(family="Times",size=18,angle=90),
axis.title.y = element_text(family="Times",size=22, angle=90),
axis.text.y = element_text(family="Times",size=18),
axis.ticks = element_blank(),
panel.background = element_rect(fill="white"),
panel.border = element_rect(colour = "black", fill=NA, size=1),
plot.title=element_text(face="bold", size=24)
))## `geom_smooth()` using method = 'loess'
(Dist_AND_propEXO_plot<-jtools::effect_plot(modelExo.signif, pred=Distance_m, interval = TRUE,
int.type = "confidence",data = TRY)+
geom_smooth(color=c('red'), formula = y~x)+
geom_point(aes(Distance_m,Exotic_Native_Ratio),TRY)+
labs(x="\nDistance (m)", y="")+
ylim(0,0.6)+
theme_classic() +
theme(
axis.text = element_text(family="Times",colour = "black"),
axis.title.x = element_text(family="Times",size=22),
axis.text.x = element_text(family="Times",size=18,angle=90),
axis.title.y = element_text(family="Times",size=22, angle=90),
axis.text.y = element_text(family="Times",size=18),
axis.ticks = element_blank(),
panel.background = element_rect(fill="white"),
panel.border = element_rect(colour = "black", fill=NA, size=1),
plot.title=element_text(face="bold", size=24)
))## `geom_smooth()` using method = 'loess'
(Building_AND_propEXO_plot<-jtools::effect_plot(modelExo.signif, pred=Average_building_age , interval = TRUE,
int.type = "confidence",data = TRY)+
geom_smooth(color=c('red'))+
geom_point(aes(Average_building_age ,Exotic_Native_Ratio),TRY)+
labs(x="\nAvg. Year of Constuction", y="")+
ylim(0,0.6)+
theme(
axis.text = element_text(family="Times",colour = "black"),
axis.title.x = element_text(family="Times",size=22),
axis.text.x = element_text(family="Times",size=18, angle=90),
axis.title.y = element_text(family="Times",size=22, angle=90),
axis.text.y = element_text(family="Times",size=18),
axis.ticks = element_blank(),
panel.background = element_rect(fill="white"),
panel.border = element_rect(colour = "black", fill=NA, size=1),
plot.title=element_text(face="bold", size=24)
))## `geom_smooth()` using method = 'loess' and formula 'y ~ x'
(Road_AND_propEXO_plot<-jtools::effect_plot(modelExo, pred=Closest_Road, interval = TRUE,
int.type = "confidence",data = TRY)+
geom_smooth(color=c('red'))+
geom_point(aes(Closest_Road,Exotic_Native_Ratio),TRY)+
labs(x="\nClosest Road (m)", y="")+
ylim(0,0.6)+
theme_classic() +
theme(
axis.text = element_text(family="Times",colour = "black"),
axis.title.x = element_text(family="Times",size=22),
axis.text.x = element_text(family="Times",size=18, angle=90),
axis.title.y = element_text(family="Times",size=22, angle=90),
axis.text.y = element_text(family="Times",size=18),
axis.ticks = element_blank(),
panel.background = element_rect(fill="white"),
panel.border = element_rect(colour = "black", fill=NA, size=1),
plot.title=element_text(face="bold", size=24)
))## `geom_smooth()` using method = 'loess' and formula 'y ~ x'
(Controlled_propEXO_plots<-grid.arrange(Area_AND_propEXO_plot,Dist_AND_propEXO_plot,Building_AND_propEXO_plot, Road_AND_propEXO_plot, nrow = 1))## `geom_smooth()` using method = 'loess'
## `geom_smooth()` using method = 'loess'## `geom_smooth()` using method = 'loess' and formula 'y ~ x'
## `geom_smooth()` using method = 'loess' and formula 'y ~ x'
## TableGrob (1 x 4) "arrange": 4 grobs
## z cells name grob
## 1 1 (1-1,1-1) arrange gtable[layout]
## 2 2 (1-1,2-2) arrange gtable[layout]
## 3 3 (1-1,3-3) arrange gtable[layout]
## 4 4 (1-1,4-4) arrange gtable[layout]ggsave("Controlled_propEXO_plots.pdf", Controlled_propEXO_plots,width = 19, height = 5)Model Exotic and Native independently
###LRT #### Exotic
modAll3<-lm(Exotic~ Area_m2 +Distance_m +
Commercial_cover+
Residential_cover+
Closest_Road+
Average_building_age,
data = TRY)
modAll3.d<-dredge(modAll3)## Fixed term is "(Intercept)"dim(modAll3.d)## [1] 64 12dim(subset(modAll3.d, cumsum(weight) <= .95) )## [1] 19 12avgmod3.95p <- model.avg(modAll3.d, cumsum(weight) <= .95)
confint(avgmod3.95p)## 2.5 % 97.5 %
## (Intercept) -6.971203e+02 3.902370e+02
## Closest_Road 3.971701e-04 2.700223e-03
## Residential_cover 7.322410e-02 2.689672e-01
## Commercial_cover -4.373788e-01 5.728961e-02
## Average_building_age -6.945090e-02 5.280615e-01
## Area_m2 -1.257737e-05 3.792542e-06
## Distance_m -2.358516e-03 2.186925e-03DF<-data.frame( TRY$Residential_cover, TRY$Closest_Road)
hier.part::hier.part(TRY$Exotic,DF,"gaussian")
## $gfs
## [1] 0.00000000 0.28558769 0.01589601 0.55375769
##
## $IJ
## I J Total
## TRY.Residential_cover 0.4117247 -0.126137 0.28558769
## TRY.Closest_Road 0.1420330 -0.126137 0.01589601
##
## $I.perc
## ind.exp.var
## TRY.Residential_cover 74.35106
## TRY.Closest_Road 25.64894
##
## $params
## $params$full.model
## [1] "y ~ TRY.Residential_cover + TRY.Closest_Road"
##
## $params$family
## [1] "gaussian"
##
## $params$link
## [1] "default"
##
## $params$gof
## [1] "Rsqu"modAll3.lm<-lm(Exotic~Residential_cover+
Closest_Road,
data = TRY)
summary(modAll3.lm)$adj.r.squared ## [1] 0.4942587Native
modAll4<-lm(Native~ Area_m2 +
Distance_m +
Commercial_cover+
Residential_cover+
Closest_Road+
Average_building_age,
data = TRY)
plot(modAll4)
## Warning in sqrt(crit * p * (1 - hh)/hh): NaNs produced
## Warning in sqrt(crit * p * (1 - hh)/hh): NaNs produced
lmtest::bptest(modAll4)##
## studentized Breusch-Pagan test
##
## data: modAll4
## BP = 2.8373, df = 6, p-value = 0.829modAll4.d<-dredge(modAll4)## Fixed term is "(Intercept)"dim(modAll4.d)## [1] 64 12dim(subset(modAll4.d, cumsum(weight) <= .95) )## [1] 31 12avgmod4.95p <- model.avg(modAll4.d, cumsum(weight) <= .95)
confint(avgmod4.95p)## 2.5 % 97.5 %
## (Intercept) -1.452386e+03 1.049575e+03
## Distance_m -1.304899e-03 6.108806e-03
## Average_building_age -4.717770e-01 1.420527e+00
## Area_m2 -3.006897e-05 1.214050e-05
## Commercial_cover -1.011542e+00 4.653252e-01
## Closest_Road -2.052519e-03 4.005600e-03
## Residential_cover -2.388505e-01 2.429632e-01DF<-data.frame( TRY$Residential_cover, TRY$Closest_Road)
hier.part::hier.part(TRY$Exotic,DF,"gaussian")
## $gfs
## [1] 0.00000000 0.28558769 0.01589601 0.55375769
##
## $IJ
## I J Total
## TRY.Residential_cover 0.4117247 -0.126137 0.28558769
## TRY.Closest_Road 0.1420330 -0.126137 0.01589601
##
## $I.perc
## ind.exp.var
## TRY.Residential_cover 74.35106
## TRY.Closest_Road 25.64894
##
## $params
## $params$full.model
## [1] "y ~ TRY.Residential_cover + TRY.Closest_Road"
##
## $params$family
## [1] "gaussian"
##
## $params$link
## [1] "default"
##
## $params$gof
## [1] "Rsqu"modAll4.lm<-lm(Exotic~Residential_cover+
Closest_Road,
data = TRY)
summary(modAll4.lm)$adj.r.squared ## [1] 0.4942587PLOT
####Exotic:Native#####
modelExoALONE.best<-update(modAll3.lm,.~.+Area_m2+Distance_m)
logit2prob <- function(logit){
odds <- exp(logit)
prob <- odds / (1 + odds)
return(prob)
}
logit2prob(coef(modelExoALONE.best))## (Intercept) Residential_cover Closest_Road Area_m2
## 0.8738109 0.5492324 0.5005612 0.4999987
## Distance_m
## 0.4997464logit2prob(confint(modelExoALONE.best))## 2.5 % 97.5 %
## (Intercept) 0.005314348 0.9998886
## Residential_cover 0.519435594 0.5786802
## Closest_Road 0.500139271 0.5009832
## Area_m2 0.499996380 0.5000009
## Distance_m 0.499281584 0.5002113(Area_AND_EXO_plot<-jtools::effect_plot(modelExoALONE.best, pred=Area_m2, interval = TRUE,
int.type = "confidence",data = IBT_data)+
geom_smooth(color=c('red'), formula = y~x)+
geom_point(aes(Area_m2,Exotic),IBT_data)+
labs(x="\nArea (m^2)", y="Number of exotic species")+
theme_classic() +
theme(
axis.text = element_text(family="Times",colour = "black"),
axis.title.x = element_text(family="Times",size=22),
axis.text.x = element_text(family="Times",size=18,angle=90),
axis.title.y = element_text(family="Times",size=22, angle=90),
axis.text.y = element_text(family="Times",size=18),
axis.ticks = element_blank(),
panel.background = element_rect(fill="white"),
panel.border = element_rect(colour = "black", fill=NA, size=1),
plot.title=element_text(face="bold", size=24)
))## `geom_smooth()` using method = 'loess'
(Dist_AND_EXO_plot<-jtools::effect_plot(modelExoALONE.best, pred=Distance_m, interval = TRUE,
int.type = "confidence",data = IBT_data)+
geom_smooth(color=c('red'), formula = y~x)+
geom_point(aes(Distance_m,Exotic),IBT_data)+
labs(x="\nDistance (m)", y="Number of exotic species")+
theme_classic() +
theme(
axis.text = element_text(family="Times",colour = "black"),
axis.title.x = element_text(family="Times",size=22),
axis.text.x = element_text(family="Times",size=18,angle=90),
axis.title.y = element_text(family="Times",size=22, angle=90),
axis.text.y = element_text(family="Times",size=18),
axis.ticks = element_blank(),
panel.background = element_rect(fill="white"),
panel.border = element_rect(colour = "black", fill=NA, size=1),
plot.title=element_text(face="bold", size=24)
))## `geom_smooth()` using method = 'loess'