title: “Final Project” author: “Shannon” date: “May 23, 2016” output: html_document — #Final Project: All Is Not Loss: Plant Biodiversity in the Anthropocene. Looking at Mediterranean, Temperate Grasslands, and Tropical dry broadleaf forest biomes and correlation between Native Species Richness and land cover type (pasture, crop, and urban) using two regression models: ordinary least squares (OLS) and generalized least squares (GLM).
library(sp)
library(maptools)## Checking rgeos availability: FALSE
## Note: when rgeos is not available, polygon geometry computations in maptools depend on gpclib,
## which has a restricted licence. It is disabled by default;
## to enable gpclib, type gpclibPermit()library(rgdal)## rgdal: version: 1.1-10, (SVN revision 622)
## Geospatial Data Abstraction Library extensions to R successfully loaded
## Loaded GDAL runtime: GDAL 2.0.1, released 2015/09/15
## Path to GDAL shared files: U:/R/win-library/3.2/rgdal/gdal
## GDAL does not use iconv for recoding strings.
## Loaded PROJ.4 runtime: Rel. 4.9.1, 04 March 2015, [PJ_VERSION: 491]
## Path to PROJ.4 shared files: U:/R/win-library/3.2/rgdal/proj
## Linking to sp version: 1.2-3library(raster)
library(gstat)
library(ggplot2)
library(spatstat)## Loading required package: nlme##
## Attaching package: 'nlme'## The following object is masked from 'package:raster':
##
## getData## Loading required package: rpart##
## spatstat 1.45-2 (nickname: 'Caretaker Mode')
## For an introduction to spatstat, type 'beginner'##
## Attaching package: 'spatstat'## The following object is masked from 'package:gstat':
##
## idw## The following objects are masked from 'package:raster':
##
## area, rotate, shiftlibrary(ncf)
library(nlme)
library(spdep)## Loading required package: Matrixlibrary(automap)Set the working directory, set headers to TRUE and get a summary to see what variable we are working with. This is the entire dataset, which is an incredible amount of information. I’ve plotted it here to give you an idea of the scale.
setwd("C:/Shannon")
ellis<- read.csv(file.choose())
med <- ellis[3628:4038,] #mediterranean forests, woodlands & shrub
trop <-ellis[13885:14250,] #tropical & subtropical dry broadleaf forests
grass <-ellis[7036:9680,] #temperate grasslands, savannas & shrublands
#####Check the first and last part of each variable to make sure we are using all of the right data set.
class(grass)## [1] "data.frame"summary(grass)## X.1 X Y hex_area
## Min. : 1777 Min. :-97.367 Min. :-34.769 Min. :7789
## 1st Qu.:27125 1st Qu.: -6.666 1st Qu.:-16.779 1st Qu.:7792
## Median :35883 Median : 20.663 Median : -5.685 Median :7792
## Mean :34892 Mean : 18.103 Mean : -4.056 Mean :7792
## 3rd Qu.:42093 3rd Qu.: 33.710 3rd Qu.: 9.703 3rd Qu.:7792
## Max. :51314 Max. :151.804 Max. : 70.836 Max. :7793
## land_km2 pland Olsoncl Olsonrlm
## Min. :4025 Min. :0.5166 Min. :6 Afrotropic :1840
## 1st Qu.:7568 1st Qu.:0.9713 1st Qu.:7 Australasia: 283
## Median :7719 Median :0.9907 Median :7 Indo-Malay : 1
## Mean :7643 Mean :0.9809 Mean :7 Nearctic : 10
## 3rd Qu.:7812 3rd Qu.:1.0026 3rd Qu.:7 Neotropic : 510
## Max. :8295 Max. :1.0646 Max. :7 Palearctic : 1
## anthrov2 totpdens pcult ppast
## Min. :11.00 Min. : 0.000 Min. :0.000000 Min. :0.0000
## 1st Qu.:41.00 1st Qu.: 1.769 1st Qu.:0.005453 1st Qu.:0.2207
## Median :42.00 Median : 7.495 Median :0.044729 Median :0.4402
## Mean :42.95 Mean : 22.144 Mean :0.089721 Mean :0.4187
## 3rd Qu.:43.00 3rd Qu.: 20.964 3rd Qu.:0.120019 3rd Qu.:0.5830
## Max. :62.00 Max. :892.969 Max. :0.976126 Max. :1.0000
## purbn HL N ASL
## Min. :0.0000000 Min. :0.0000 Min. : 226.7 Min. : 0.00
## 1st Qu.:0.0000000 1st Qu.:0.2205 1st Qu.: 792.1 1st Qu.: 46.77
## Median :0.0001213 Median :0.3809 Median :1040.7 Median : 80.49
## Mean :0.0013714 Mean :0.3702 Mean :1082.9 Mean : 92.33
## 3rd Qu.:0.0007514 3rd Qu.:0.5087 3rd Qu.:1327.6 3rd Qu.:126.42
## Max. :0.2237110 Max. :0.9903 Max. :3422.5 Max. :761.08
## IS CS OS ASI
## Min. : 13.54 Min. : 0.00 Min. : 0.00 Min. : 13.54
## 1st Qu.: 37.76 1st Qu.: 7.00 1st Qu.: 0.00 1st Qu.: 48.69
## Median : 47.23 Median :14.00 Median : 0.00 Median : 64.78
## Mean : 48.33 Mean :15.28 Mean : 11.95 Mean : 75.55
## 3rd Qu.: 57.67 3rd Qu.:22.00 3rd Qu.: 0.00 3rd Qu.: 80.63
## Max. :125.36 Max. :54.00 Max. :400.00 Max. :482.63
## ASR ASR_N ASL_N ASI_N
## Min. : 240.0 Min. :0.5190 Min. :0.00000 Min. :0.03780
## 1st Qu.: 754.6 1st Qu.:0.9440 1st Qu.:0.04385 1st Qu.:0.05184
## Median :1016.4 Median :0.9808 Median :0.08269 Median :0.05805
## Mean :1066.1 Mean :0.9821 Mean :0.08777 Mean :0.06986
## 3rd Qu.:1334.1 3rd Qu.:1.0238 3rd Qu.:0.12007 3rd Qu.:0.06731
## Max. :3520.6 Max. :1.5723 Max. :0.56573 Max. :0.58782
## IS_N ASI_ASL_N
## Min. :0.03663 Min. :0.04694
## 1st Qu.:0.04344 1st Qu.:0.10299
## Median :0.04538 Median :0.14227
## Mean :0.04573 Mean :0.15764
## 3rd Qu.:0.04767 3rd Qu.:0.18769
## Max. :0.07966 Max. :0.65494class(med)## [1] "data.frame"summary(med)## X.1 X Y hex_area
## Min. : 9705 Min. :-122.364 Min. :-38.32 Min. :7789
## 1st Qu.:12280 1st Qu.: -3.901 1st Qu.:-30.98 1st Qu.:7792
## Median :14176 Median : 10.688 Median : 34.43 Median :7792
## Mean :25190 Mean : 28.145 Mean : 14.27 Mean :7792
## 3rd Qu.:49497 3rd Qu.: 39.667 3rd Qu.: 38.57 3rd Qu.:7792
## Max. :52954 Max. : 144.790 Max. : 44.55 Max. :7794
## land_km2 pland Olsoncl Olsonrlm
## Min. :3942 Min. :0.5059 Min. :12 Afrotropic : 12
## 1st Qu.:6940 1st Qu.:0.8908 1st Qu.:12 Australasia:100
## Median :7652 Median :0.9821 Median :12 Indo-Malay : 0
## Mean :7091 Mean :0.9101 Mean :12 Nearctic : 16
## 3rd Qu.:7739 3rd Qu.:0.9933 3rd Qu.:12 Neotropic : 23
## Max. :8093 Max. :1.0386 Max. :12 Palearctic :260
## anthrov2 totpdens pcult ppast
## Min. :11.0 Min. : 0.000 Min. :0.00000 Min. :0.001089
## 1st Qu.:32.0 1st Qu.: 2.871 1st Qu.:0.04851 1st Qu.:0.172404
## Median :34.0 Median : 33.489 Median :0.25226 Median :0.289222
## Mean :36.8 Mean : 89.192 Mean :0.26639 Mean :0.369447
## 3rd Qu.:42.0 3rd Qu.: 92.511 3rd Qu.:0.40145 3rd Qu.:0.511609
## Max. :62.0 Max. :2167.880 Max. :0.98547 Max. :1.000000
## purbn HL N ASL
## Min. :0.0000000 Min. :0.08858 Min. : 473.2 Min. : 16.40
## 1st Qu.:0.0001052 1st Qu.:0.40839 1st Qu.: 863.4 1st Qu.: 99.79
## Median :0.0028499 Median :0.52303 Median :1217.4 Median :172.50
## Mean :0.0148183 Mean :0.52750 Mean :1284.3 Mean :202.87
## 3rd Qu.:0.0085479 3rd Qu.:0.63697 3rd Qu.:1558.9 3rd Qu.:263.07
## Max. :0.8257600 Max. :0.99393 Max. :3417.0 Max. :923.20
## IS CS OS ASI
## Min. : 48.25 Min. : 0.00 Min. : 0.00 Min. : 53.25
## 1st Qu.: 79.01 1st Qu.:16.00 1st Qu.: 0.00 1st Qu.: 99.64
## Median :104.73 Median :40.00 Median : 0.00 Median :183.99
## Mean :108.17 Mean :35.52 Mean : 89.54 Mean :233.23
## 3rd Qu.:128.27 3rd Qu.:56.00 3rd Qu.:200.00 3rd Qu.:369.21
## Max. :244.11 Max. :78.00 Max. :400.00 Max. :651.34
## ASR ASR_N ASL_N ASI_N
## Min. : 500.4 Min. :0.6190 Min. :0.01838 Min. :0.07584
## 1st Qu.: 883.8 1st Qu.:0.9596 1st Qu.:0.09966 1st Qu.:0.10633
## Median :1207.7 Median :1.0044 Median :0.13762 Median :0.12601
## Mean :1314.6 Mean :1.0291 Mean :0.15439 Mean :0.18346
## 3rd Qu.:1644.8 3rd Qu.:1.0787 3rd Qu.:0.18344 3rd Qu.:0.24818
## Max. :3391.9 Max. :1.8181 Max. :0.63969 Max. :0.91529
## IS_N ASI_ASL_N
## Min. :0.07144 Min. :0.1086
## 1st Qu.:0.08228 1st Qu.:0.2196
## Median :0.08602 Median :0.2915
## Mean :0.08662 Mean :0.3379
## 3rd Qu.:0.09151 3rd Qu.:0.4086
## Max. :0.10197 Max. :1.0125class(trop)## [1] "data.frame"summary(trop)## X.1 X Y hex_area
## Min. :15977 Min. :-110.19 Min. :-21.03 Min. :7790
## 1st Qu.:20735 1st Qu.: -65.35 1st Qu.: -4.73 1st Qu.:7792
## Median :23920 Median : 74.88 Median : 15.58 Median :7792
## Mean :27221 Mean : 21.19 Mean : 9.78 Mean :7792
## 3rd Qu.:35380 3rd Qu.: 81.98 3rd Qu.: 21.37 3rd Qu.:7792
## Max. :44416 Max. : 125.38 Max. : 30.62 Max. :7792
## land_km2 pland Olsoncl Olsonrlm
## Min. :3916 Min. :0.5026 Min. :1.000 Afrotropic : 27
## 1st Qu.:7436 1st Qu.:0.9543 1st Qu.:2.000 Australasia: 8
## Median :7636 Median :0.9800 Median :2.000 Indo-Malay :189
## Mean :7353 Mean :0.9437 Mean :1.997 Nearctic : 9
## 3rd Qu.:7790 3rd Qu.:0.9998 3rd Qu.:2.000 Neotropic :132
## Max. :8134 Max. :1.0439 Max. :2.000 Palearctic : 1
## anthrov2 totpdens pcult ppast
## Min. :21.00 Min. : 0.00 Min. :0.00000 Min. :0.0001579
## 1st Qu.:23.00 1st Qu.: 16.03 1st Qu.:0.08626 1st Qu.:0.0317180
## Median :32.00 Median : 86.71 Median :0.28876 Median :0.0717338
## Mean :34.08 Mean : 145.86 Mean :0.32136 Mean :0.1897127
## 3rd Qu.:42.00 3rd Qu.: 204.26 3rd Qu.:0.50887 3rd Qu.:0.3182775
## Max. :61.00 Max. :1249.28 Max. :0.93095 Max. :0.8460300
## purbn HL N
## Min. :0.0000000 Min. :0.002707 Min. : 546.4
## 1st Qu.:0.0002384 1st Qu.:0.312356 1st Qu.:1124.3
## Median :0.0016723 Median :0.467280 Median :1539.5
## Mean :0.0047885 Mean :0.452619 Mean :1571.5
## 3rd Qu.:0.0052725 3rd Qu.:0.574748 3rd Qu.:1925.6
## Max. :0.0732057 Max. :0.933265 Max. :3796.1
## ASL IS CS OS
## Min. : 0.7485 Min. : 77.34 Min. : 4.00 Min. : 0.00
## 1st Qu.:110.6149 1st Qu.:162.40 1st Qu.:24.00 1st Qu.: 0.00
## Median :161.1935 Median :209.04 Median :34.00 Median : 0.00
## Mean :189.7042 Mean :211.79 Mean :32.06 Mean : 54.64
## 3rd Qu.:256.3926 3rd Qu.:253.13 3rd Qu.:41.00 3rd Qu.:200.00
## Max. :680.1391 Max. :441.63 Max. :60.00 Max. :400.00
## ASI ASR ASR_N ASL_N
## Min. :114.1 Min. : 651.1 Min. :0.7553 Min. :0.000569
## 1st Qu.:207.1 1st Qu.:1173.0 1st Qu.:1.0309 1st Qu.:0.075629
## Median :260.5 Median :1670.7 Median :1.0685 Median :0.124555
## Mean :298.5 Mean :1680.3 Mean :1.0695 Mean :0.129935
## 3rd Qu.:373.1 3rd Qu.:2080.1 3rd Qu.:1.1064 3rd Qu.:0.164367
## Max. :761.8 Max. :4070.6 Max. :1.3274 Max. :0.433611
## ASI_N IS_N ASI_ASL_N
## Min. :0.07232 Min. :0.04072 Min. :0.1490
## 1st Qu.:0.14730 1st Qu.:0.13143 1st Qu.:0.2283
## Median :0.17509 Median :0.13686 Median :0.3095
## Mean :0.19941 Mean :0.13801 Mean :0.3293
## 3rd Qu.:0.20856 3rd Qu.:0.14482 3rd Qu.:0.3917
## Max. :0.51369 Max. :0.16491 Max. :0.7683plot(ellis$X, ellis$Y) #plot the datapoints
coordinates(ellis)<-c("X","Y")
class(ellis)## [1] "SpatialPointsDataFrame"
## attr(,"package")
## [1] "sp"str(ellis)## Formal class 'SpatialPointsDataFrame' [package "sp"] with 5 slots
## ..@ data :'data.frame': 16805 obs. of 24 variables:
## .. ..$ X.1 : int [1:16805] 7287 7347 7426 7430 7454 7470 7486 7558 7566 7574 ...
## .. ..$ hex_area : num [1:16805] 7792 7792 7792 7792 7792 ...
## .. ..$ land_km2 : num [1:16805] 7358 7584 7564 7786 7576 ...
## .. ..$ pland : num [1:16805] 0.944 0.973 0.971 0.999 0.972 ...
## .. ..$ Olsoncl : int [1:16805] 13 13 13 13 13 13 13 13 13 13 ...
## .. ..$ Olsonrlm : Factor w/ 6 levels "Afrotropic","Australasia",..: 6 6 6 6 6 6 6 6 6 6 ...
## .. ..$ anthrov2 : int [1:16805] 42 42 61 43 54 43 42 43 42 42 ...
## .. ..$ totpdens : num [1:16805] 1.519 1.502 1.438 0.213 0.957 ...
## .. ..$ pcult : num [1:16805] 0.0975 0.0668 0.0259 0.0373 0.0216 ...
## .. ..$ ppast : num [1:16805] 0.238 0.501 0.125 0.855 0.226 ...
## .. ..$ purbn : num [1:16805] 1.29e-04 2.48e-05 1.85e-04 0.00 1.26e-05 ...
## .. ..$ HL : num [1:16805] 0.257 0.401 0.109 0.607 0.172 ...
## .. ..$ N : num [1:16805] 724 580 677 535 795 ...
## .. ..$ ASL : num [1:16805] 23.24 31.77 8.58 52.26 16.35 ...
## .. ..$ IS : num [1:16805] 42.2 35.2 39.9 32.9 45.5 ...
## .. ..$ CS : int [1:16805] 12 11 11 10 9 17 15 11 19 2 ...
## .. ..$ OS : int [1:16805] 0 0 0 0 0 200 0 0 0 0 ...
## .. ..$ ASI : num [1:16805] 54.2 46.2 50.9 42.9 54.5 ...
## .. ..$ ASR : num [1:16805] 755 594 720 525 833 ...
## .. ..$ ASR_N : num [1:16805] 1.043 1.025 1.063 0.982 1.048 ...
## .. ..$ ASL_N : num [1:16805] 0.0321 0.0548 0.0127 0.0977 0.0206 ...
## .. ..$ ASI_N : num [1:16805] 0.0748 0.0796 0.0752 0.0802 0.0686 ...
## .. ..$ IS_N : num [1:16805] 0.0582 0.0606 0.0589 0.0615 0.0573 ...
## .. ..$ ASI_ASL_N: num [1:16805] 0.1069 0.1344 0.0878 0.178 0.0892 ...
## ..@ coords.nrs : int [1:2] 2 3
## ..@ coords : num [1:16805, 1:2] 93.3 91.9 95.5 64 90.5 ...
## .. ..- attr(*, "dimnames")=List of 2
## .. .. ..$ : NULL
## .. .. ..$ : chr [1:2] "X" "Y"
## ..@ bbox : num [1:2, 1:2] -177.5 -54.9 179.1 82.3
## .. ..- attr(*, "dimnames")=List of 2
## .. .. ..$ : chr [1:2] "X" "Y"
## .. .. ..$ : chr [1:2] "min" "max"
## ..@ proj4string:Formal class 'CRS' [package "sp"] with 1 slot
## .. .. ..@ projargs: chr NAproj4string(ellis)<- CRS("+proj=longlat + ellps=WGS84 +datum=WGS84 +no_defs +towgs84=0,0,0")
ellis@proj4string## CRS arguments:
## +proj=longlat +datum=WGS84 +no_defs +towgs84=0,0,0 +ellps=WGS84plot(ellis, axes=TRUE)
Mediterranean Biome
plot(med$X, med$Y) #plot the datapoints
coordinates(med)<-c("X","Y")
class(med)## [1] "SpatialPointsDataFrame"
## attr(,"package")
## [1] "sp"str(med)## Formal class 'SpatialPointsDataFrame' [package "sp"] with 5 slots
## ..@ data :'data.frame': 411 obs. of 24 variables:
## .. ..$ X.1 : int [1:411] 9705 9779 9833 9849 9864 9887 9935 9936 10028 10103 ...
## .. ..$ hex_area : num [1:411] 7792 7792 7792 7792 7794 ...
## .. ..$ land_km2 : num [1:411] 6010 7105 7624 7625 7807 ...
## .. ..$ pland : num [1:411] 0.771 0.912 0.978 0.979 1.002 ...
## .. ..$ Olsoncl : int [1:411] 12 12 12 12 12 12 12 12 12 12 ...
## .. ..$ Olsonrlm : Factor w/ 6 levels "Afrotropic","Australasia",..: 6 6 6 6 6 6 6 6 6 6 ...
## .. ..$ anthrov2 : int [1:411] 51 32 32 32 32 32 42 32 32 32 ...
## .. ..$ totpdens : num [1:411] 494.3 128.3 26.3 113.1 266.3 ...
## .. ..$ pcult : num [1:411] 0.333 0.256 0.292 0.317 0.419 ...
## .. ..$ ppast : num [1:411] 0.1121 0.0873 0.2129 0.1032 0.1214 ...
## .. ..$ purbn : num [1:411] 0.07052 0.02353 0.00503 0.02348 0.05616 ...
## .. ..$ HL : num [1:411] 0.478 0.338 0.439 0.409 0.556 ...
## .. ..$ N : num [1:411] 1417 1667 1181 1505 1614 ...
## .. ..$ ASL : num [1:411] 173 132 129 151 242 ...
## .. ..$ IS : num [1:411] 119 136 102 125 132 ...
## .. ..$ CS : int [1:411] 63 60 52 50 61 43 40 45 48 51 ...
## .. ..$ OS : int [1:411] 400 200 200 200 400 200 200 200 400 200 ...
## .. ..$ ASI : num [1:411] 582 396 354 375 593 ...
## .. ..$ ASR : num [1:411] 1825 1930 1406 1729 1965 ...
## .. ..$ ASR_N : num [1:411] 1.29 1.16 1.19 1.15 1.22 ...
## .. ..$ ASL_N : num [1:411] 0.1219 0.0792 0.1091 0.1 0.1497 ...
## .. ..$ ASI_N : num [1:411] 0.411 0.237 0.3 0.249 0.367 ...
## .. ..$ IS_N : num [1:411] 0.0837 0.0813 0.0865 0.0828 0.0818 ...
## .. ..$ ASI_ASL_N: num [1:411] 0.532 0.317 0.409 0.349 0.517 ...
## ..@ coords.nrs : int [1:2] 2 3
## ..@ coords : num [1:411, 1:2] 8.71 9.98 3.06 4.32 11.25 ...
## .. ..- attr(*, "dimnames")=List of 2
## .. .. ..$ : chr [1:411] "3628" "3629" "3630" "3631" ...
## .. .. ..$ : chr [1:2] "X" "Y"
## ..@ bbox : num [1:2, 1:2] -122.4 -38.3 144.8 44.6
## .. ..- attr(*, "dimnames")=List of 2
## .. .. ..$ : chr [1:2] "X" "Y"
## .. .. ..$ : chr [1:2] "min" "max"
## ..@ proj4string:Formal class 'CRS' [package "sp"] with 1 slot
## .. .. ..@ projargs: chr NAproj4string(med)<- CRS("+proj=longlat + ellps=WGS84 +datum=WGS84 +no_defs +towgs84=0,0,0")
med@proj4string## CRS arguments:
## +proj=longlat +datum=WGS84 +no_defs +towgs84=0,0,0 +ellps=WGS84plot(med, axes=TRUE) #Plot the spatial data with WGS84 projection
Grassland Biome
plot(grass$X, grass$Y) #plot the datapoints
coordinates(grass)<-c("X","Y")
class(grass)## [1] "SpatialPointsDataFrame"
## attr(,"package")
## [1] "sp"str(grass)## Formal class 'SpatialPointsDataFrame' [package "sp"] with 5 slots
## ..@ data :'data.frame': 2645 obs. of 24 variables:
## .. ..$ X.1 : int [1:2645] 16191 16257 16369 16425 16537 16600 16698 16875 16987 17064 ...
## .. ..$ hex_area : num [1:2645] 7792 7792 7792 7792 7792 ...
## .. ..$ land_km2 : num [1:2645] 6417 7767 5808 4025 4903 ...
## .. ..$ pland : num [1:2645] 0.824 0.997 0.745 0.517 0.629 ...
## .. ..$ Olsoncl : int [1:2645] 7 7 7 7 7 7 7 7 7 7 ...
## .. ..$ Olsonrlm : Factor w/ 6 levels "Afrotropic","Australasia",..: 4 4 4 4 4 4 4 4 3 4 ...
## .. ..$ anthrov2 : int [1:2645] 32 11 33 32 11 32 33 33 21 33 ...
## .. ..$ totpdens : num [1:2645] 33.7 308.5 65.8 260.6 351.3 ...
## .. ..$ pcult : num [1:2645] 0.272 0.414 0.49 0.263 0.447 ...
## .. ..$ ppast : num [1:2645] 0.141 0.284 0.251 0.178 0.339 ...
## .. ..$ purbn : num [1:2645] 0.0183 0.1889 0.0421 0.1611 0.2237 ...
## .. ..$ HL : num [1:2645] 0.385 0.792 0.699 0.543 0.897 ...
## .. ..$ N : num [1:2645] 1344 1132 1288 1360 1524 ...
## .. ..$ ASL : num [1:2645] 113 279 251 179 512 ...
## .. ..$ IS : num [1:2645] 58.2 50.6 56.2 58.8 64.6 ...
## .. ..$ CS : int [1:2645] 10 12 7 3 10 6 10 10 43 11 ...
## .. ..$ OS : int [1:2645] 200 400 200 400 400 200 200 200 0 0 ...
## .. ..$ ASI : num [1:2645] 268 463 263 462 475 ...
## .. ..$ ASR : num [1:2645] 1499 1316 1300 1643 1487 ...
## .. ..$ ASR_N : num [1:2645] 1.116 1.162 1.01 1.208 0.975 ...
## .. ..$ ASL_N : num [1:2645] 0.0837 0.2463 0.1946 0.1314 0.3358 ...
## .. ..$ ASI_N : num [1:2645] 0.2 0.409 0.204 0.34 0.311 ...
## .. ..$ IS_N : num [1:2645] 0.0433 0.0447 0.0437 0.0433 0.0424 ...
## .. ..$ ASI_ASL_N: num [1:2645] 0.283 0.655 0.399 0.471 0.647 ...
## ..@ coords.nrs : int [1:2] 2 3
## ..@ coords : num [1:2645, 1:2] -93.4 -95.8 -94.2 -90.2 -95 ...
## .. ..- attr(*, "dimnames")=List of 2
## .. .. ..$ : chr [1:2645] "7036" "7037" "7038" "7039" ...
## .. .. ..$ : chr [1:2] "X" "Y"
## ..@ bbox : num [1:2, 1:2] -97.4 -34.8 151.8 70.8
## .. ..- attr(*, "dimnames")=List of 2
## .. .. ..$ : chr [1:2] "X" "Y"
## .. .. ..$ : chr [1:2] "min" "max"
## ..@ proj4string:Formal class 'CRS' [package "sp"] with 1 slot
## .. .. ..@ projargs: chr NAproj4string(grass)<- CRS("+proj=longlat + ellps=WGS84 +datum=WGS84 +no_defs +towgs84=0,0,0")
grass@proj4string## CRS arguments:
## +proj=longlat +datum=WGS84 +no_defs +towgs84=0,0,0 +ellps=WGS84plot(grass, axes=TRUE) 
Tropical Biome
plot(trop$X, trop$Y) #plot the datapoints
coordinates(trop)<-c("X","Y")
class(trop)## [1] "SpatialPointsDataFrame"
## attr(,"package")
## [1] "sp"str(trop)## Formal class 'SpatialPointsDataFrame' [package "sp"] with 5 slots
## ..@ data :'data.frame': 366 obs. of 24 variables:
## .. ..$ X.1 : int [1:366] 16706 16938 17132 17385 17591 17711 17836 17893 17961 18067 ...
## .. ..$ hex_area : num [1:366] 7792 7792 7792 7792 7792 ...
## .. ..$ land_km2 : num [1:366] 7860 7253 7514 7891 7059 ...
## .. ..$ pland : num [1:366] 1.009 0.931 0.964 1.013 0.906 ...
## .. ..$ Olsoncl : int [1:366] 2 2 2 2 2 2 2 2 2 2 ...
## .. ..$ Olsonrlm : Factor w/ 6 levels "Afrotropic","Australasia",..: 4 4 4 4 4 3 4 3 3 5 ...
## .. ..$ anthrov2 : int [1:366] 42 43 42 42 42 22 42 22 22 42 ...
## .. ..$ totpdens : num [1:366] 2.26 1.52 2.91 1.31 62.64 ...
## .. ..$ pcult : num [1:366] 0.02879 0.17579 0.00189 0.08294 0.34277 ...
## .. ..$ ppast : num [1:366] 0.57 0.319 0.568 0.394 0.478 ...
## .. ..$ purbn : num [1:366] 0 0 0 0 0.00666 ...
## .. ..$ HL : num [1:366] 0.409 0.388 0.381 0.345 0.668 ...
## .. ..$ N : num [1:366] 849 1216 750 1181 865 ...
## .. ..$ ASL : num [1:366] 88.8 119.2 71.8 100.5 178.9 ...
## .. ..$ IS : num [1:366] 129 174 117 169 131 ...
## .. ..$ CS : int [1:366] 15 10 17 20 29 27 26 27 40 20 ...
## .. ..$ OS : int [1:366] 0 0 0 0 200 0 0 200 0 0 ...
## .. ..$ ASI : num [1:366] 144 184 134 189 360 ...
## .. ..$ ASR : num [1:366] 905 1280 813 1270 1047 ...
## .. ..$ ASR_N : num [1:366] 1.07 1.05 1.08 1.08 1.21 ...
## .. ..$ ASL_N : num [1:366] 0.1046 0.0981 0.0957 0.0851 0.2068 ...
## .. ..$ ASI_N : num [1:366] 0.17 0.151 0.178 0.16 0.417 ...
## .. ..$ IS_N : num [1:366] 0.152 0.143 0.156 0.144 0.152 ...
## .. ..$ ASI_ASL_N: num [1:366] 0.275 0.249 0.274 0.246 0.623 ...
## ..@ coords.nrs : int [1:2] 2 3
## ..@ coords : num [1:366, 1:2] -110 -109 -110 -109 -110 ...
## .. ..- attr(*, "dimnames")=List of 2
## .. .. ..$ : chr [1:366] "13885" "13886" "13887" "13888" ...
## .. .. ..$ : chr [1:2] "X" "Y"
## ..@ bbox : num [1:2, 1:2] -110.2 -21 125.4 30.6
## .. ..- attr(*, "dimnames")=List of 2
## .. .. ..$ : chr [1:2] "X" "Y"
## .. .. ..$ : chr [1:2] "min" "max"
## ..@ proj4string:Formal class 'CRS' [package "sp"] with 1 slot
## .. .. ..@ projargs: chr NAproj4string(trop)<- CRS("+proj=longlat + ellps=WGS84 +datum=WGS84 +no_defs +towgs84=0,0,0")
trop@proj4string## CRS arguments:
## +proj=longlat +datum=WGS84 +no_defs +towgs84=0,0,0 +ellps=WGS84plot(trop, axes=TRUE) 
Create variograms of the Native Species Richness (N) and land cover types: pasture, crop, and urban for the Mediterranean biome.
NmedVarCloud<-variogram(N~1, med, cloud = FALSE)
plot(NmedVarCloud, pch=20,cex=1.5, col="dodgerblue3",
ylab=expression("Semivariance ("*gamma*")"),
xlab="Distance (km)", main="Native Species Richness")
pcultmedVarCloud<-variogram(pcult~1, med, cloud = FALSE, cutoff=1200)
plot(pcultmedVarCloud, pch=20,cex=1.5, col="dodgerblue3",
ylab=expression("Semivariance ("*gamma*")"),
xlab="Distance (km)", main="Crops")
ppastmedVarCloud<-variogram(ppast~1, med, cloud = FALSE, cutoff=1500)
plot(ppastmedVarCloud, pch=20,cex=1.5, col="dodgerblue3",
ylab=expression("Semivariance ("*gamma*")"),
xlab="Distance (km)", main="Pasture")
purbanmedVarCloud<-variogram(purbn~1, med, cloud = FALSE, cutoff=1000)
plot(purbanmedVarCloud, pch=20,cex=1.5, col="dodgerblue3",
ylab=expression("Semivariance ("*gamma*")"),
xlab="Distance (km)", main="Native Species Richness")
Create variograms of the Native Species Richness (N) and land cover types: pasture, crop, and urban for the Tropical biome.
NtropVarCloud<-variogram(N~1, trop, cloud = FALSE)
plot(NtropVarCloud, pch=20,cex=1.5, col="dodgerblue3",
ylab=expression("Semivariance ("*gamma*")"),
xlab="Distance (km)", main="Native Species Richness")
pculttropVarCloud<-variogram(pcult~1, trop, cloud = FALSE, cutoff=1200)
plot(pculttropVarCloud, pch=20,cex=1.5, col="dodgerblue3",
ylab=expression("Semivariance ("*gamma*")"),
xlab="Distance (km)", main="Crops")
ppasttropVarCloud<-variogram(ppast~1, trop, cloud = FALSE, cutoff=1000)
plot(ppasttropVarCloud, pch=20,cex=1.5, col="dodgerblue3",
ylab=expression("Semivariance ("*gamma*")"),
xlab="Distance (km)", main="Pasture")
purbantropVarCloud<-variogram(purbn~1, trop, cloud = FALSE, cutoff=2000)
plot(purbantropVarCloud, pch=20,cex=1.5, col="dodgerblue3",
ylab=expression("Semivariance ("*gamma*")"),
xlab="Distance (km)", main="Urban")
Create variograms of the Native Species Richness (N) and land cover types: pasture, crop, and urban for the Grassland biome.
NgVarCloud<-variogram(N~1, grass, cloud = FALSE)
plot(NgVarCloud, pch=20,cex=1.5, col="dodgerblue3",
ylab=expression("Semivariance ("*gamma*")"),
xlab="Distance (km)", main="Native Species Richness")
pcultgVarCloud<-variogram(pcult~1, grass, cloud = FALSE, cutoff=1200)
plot(pcultgVarCloud, pch=20,cex=1.5, col="dodgerblue3",
ylab=expression("Semivariance ("*gamma*")"),
xlab="Distance (km)", main="Crops")
ppastgVarCloud<-variogram(ppast~1, grass, cloud = FALSE, cutoff=1000)
plot(ppastgVarCloud, pch=20,cex=1.5, col="dodgerblue3",
ylab=expression("Semivariance ("*gamma*")"),
xlab="Distance (km)", main="Pasture")
purbangVarCloud<-variogram(purbn~1, grass, cloud = FALSE, cutoff=1000)
plot(purbangVarCloud, pch=20,cex=1.5, col="dodgerblue3",
ylab=expression("Semivariance ("*gamma*")"),
xlab="Distance (km)", main="Urban")
Bubble Plots:
understanding the magnitude of correlation of population density and land cover types.
#*Mediterranean*
bubble(med, xcol="X", ycol="Y", zcol="totpdens", maxsize=1, main="Population Density")
bubble(med, xcol="X", ycol="Y", zcol="pcult", maxsize=1, main="Crop")
bubble(med, xcol="X", ycol="Y", zcol="ppast", maxsize=1, main="Pasture")
bubble(med, xcol="X", ycol="Y", zcol="purbn", maxsize=1, main="Urban")
#*Tropical*
bubble(trop, xcol="X", ycol="Y", zcol="totpdens", maxsize=1, main="Population Density")
bubble(trop, xcol="X", ycol="Y", zcol="pcult", maxsize=1, main="Crop")
bubble(trop, xcol="X", ycol="Y", zcol="ppast", maxsize=1, main="Pasture")
bubble(trop, xcol="X", ycol="Y", zcol="purbn", maxsize=1, main="Urban")
#*Grassland*
bubble(grass, xcol="X", ycol="Y", zcol="totpdens", maxsize=1, main="Population Density")
bubble(grass, xcol="X", ycol="Y", zcol="pcult", maxsize=1, main="Crop")
bubble(grass, xcol="X", ycol="Y", zcol="ppast", maxsize=1, main="Pasture")
bubble(grass, xcol="X", ycol="Y", zcol="purbn", maxsize=1, main="Urban")
Land Cover: Crops, pasture, and urban areas
#*Mediterranean*
hscat(N~1,med,breaks=seq(0,1000,by=100))
hscat(pcult~1,med,breaks=seq(0,1000,by=100))
hscat(ppast~1,med,breaks=seq(0,1000,by=100))
hscat(purbn~1,med,breaks=seq(0,1000,by=100))
#*Tropical*
hscat(N~1,trop,breaks=seq(0,1000,by=100))
hscat(pcult~1,trop,breaks=seq(0,1000,by=100))
hscat(ppast~1,trop,breaks=seq(0,1000,by=100))
hscat(purbn~1,trop,breaks=seq(0,1000,by=100))
#*Grassland*
hscat(N~1,grass,breaks=seq(0,1000,by=100))
hscat(pcult~1,grass,breaks=seq(0,1000,by=100))
hscat(ppast~1,grass,breaks=seq(0,1000,by=100))
hscat(purbn~1,grass,breaks=seq(0,1000,by=100))
Moran’s I and Splines for the Mediterranean biome
#*Native Species Richness*
Nmed.8<-knn2nb(knearneigh(med,k=8))
moran.test(med$N, nb2listw(Nmed.8)) #Moran's I = 0.83; highly positive autocorrelation.##
## Moran I test under randomisation
##
## data: med$N
## weights: nb2listw(Nmed.8)
##
## Moran I statistic standard deviate = 36.045, p-value < 2.2e-16
## alternative hypothesis: greater
## sample estimates:
## Moran I statistic Expectation Variance
## 0.833926655 -0.002439024 0.000538388Nmed.I<-spline.correlog(x=med$X, y=med$Y, z=med$N, resamp=5, latlon=TRUE, quiet=TRUE, xmax=3000)
plot(Nmed.I) #Autocorrelation up to around 1 km.
#*Crops*
cropmed.8<-knn2nb(knearneigh(med,k=8))
moran.test(med$pcult, nb2listw(cropmed.8)) #Moran's I; 0.56 indicating positive autocorrelation##
## Moran I test under randomisation
##
## data: med$pcult
## weights: nb2listw(cropmed.8)
##
## Moran I statistic standard deviate = 24.332, p-value < 2.2e-16
## alternative hypothesis: greater
## sample estimates:
## Moran I statistic Expectation Variance
## 0.5630721332 -0.0024390244 0.0005401726cropmed.I<-spline.correlog(x=med$X, y=med$Y, z=med$pcult, resamp=5, latlon=TRUE, quiet=TRUE, xmax=4000)
plot(cropmed.I) #Autocorrelation up to 1/2 a km
#*Pasture*
pastmed.8<-knn2nb(knearneigh(med,k=8))
moran.test(med$ppast, nb2listw(pastmed.8)) #Moran's I; 0.65 highly positive autocorrelation##
## Moran I test under randomisation
##
## data: med$ppast
## weights: nb2listw(pastmed.8)
##
## Moran I statistic standard deviate = 28.08, p-value < 2.2e-16
## alternative hypothesis: greater
## sample estimates:
## Moran I statistic Expectation Variance
## 0.650546202 -0.002439024 0.000540783pastmed.I<-spline.correlog(x=med$X, y=med$Y, z=med$ppast, resamp=5, latlon=TRUE, quiet=TRUE, xmax=1500)
plot(pastmed.I) #Autocorrelation up to ~3/4 of a km.
#*Urban*
urbanmed.8<-knn2nb(knearneigh(med,k=8))
moran.test(med$purbn, nb2listw(urbanmed.8)) #Not highly autocorrelated, Moran's I; 0.17##
## Moran I test under randomisation
##
## data: med$purbn
## weights: nb2listw(urbanmed.8)
##
## Moran I statistic standard deviate = 8.777, p-value < 2.2e-16
## alternative hypothesis: greater
## sample estimates:
## Moran I statistic Expectation Variance
## 0.1657291941 -0.0024390244 0.0003671066urbanmed.I<-spline.correlog(x=med$X, y=med$Y, z=med$purbn, resamp=5, latlon=TRUE, quiet=TRUE, xmax=1000)
plot(urbanmed.I) #Not really much autocorrelation, possibly some at very short distances <200m...
Moran’s I and Splines for the Tropical biome
#*Native Species Richness*
Ntrop.8<-knn2nb(knearneigh(trop,k=8))
moran.test(trop$N, nb2listw(Ntrop.8)) #Moran's I; 0.72 high positive spatial autocorrelation##
## Moran I test under randomisation
##
## data: trop$N
## weights: nb2listw(Ntrop.8)
##
## Moran I statistic standard deviate = 29.728, p-value < 2.2e-16
## alternative hypothesis: greater
## sample estimates:
## Moran I statistic Expectation Variance
## 0.7252692609 -0.0027397260 0.0005997065Ntrop.I<-spline.correlog(x=trop$X, y=trop$Y, z=trop$N, resamp=5, latlon=TRUE, quiet=TRUE, xmax=3000)
plot(Ntrop.I) #Spatial autocorrelation up to around 2 km
#*Crops*
croptrop.8<-knn2nb(knearneigh(trop,k=8))
moran.test(trop$pcult, nb2listw(croptrop.8)) #High autocorrelation; Moran's I; 0.87##
## Moran I test under randomisation
##
## data: trop$pcult
## weights: nb2listw(croptrop.8)
##
## Moran I statistic standard deviate = 35.64, p-value < 2.2e-16
## alternative hypothesis: greater
## sample estimates:
## Moran I statistic Expectation Variance
## 0.8721964736 -0.0027397260 0.0006026656croptrop.I<-spline.correlog(x=trop$X, y=trop$Y, z=trop$pcult, resamp=5, latlon=TRUE, quiet=TRUE, xmax=4000)
plot(croptrop.I) #Whoa! Autocorrelation up to 3 km! Crops seem like crops affecting native species in the tropics is a big deal
#*Pasture*
pasttrop.8<-knn2nb(knearneigh(trop,k=8))
moran.test(trop$ppast, nb2listw(pasttrop.8)) # Again, high positive autocorrelation.##
## Moran I test under randomisation
##
## data: trop$ppast
## weights: nb2listw(pasttrop.8)
##
## Moran I statistic standard deviate = 35.793, p-value < 2.2e-16
## alternative hypothesis: greater
## sample estimates:
## Moran I statistic Expectation Variance
## 0.8745684579 -0.0027397260 0.0006007678pasttrop.I<-spline.correlog(x=trop$X, y=trop$Y, z=trop$ppast, resamp=5, latlon=TRUE, quiet=TRUE, xmax=5000)
plot(pasttrop.I) #autocorrelation up to 4.5 km
#*Urban*
urbantrop.8<-knn2nb(knearneigh(trop,k=8))
moran.test(trop$purbn, nb2listw(urbantrop.8)) #Low Moran's I.##
## Moran I test under randomisation
##
## data: trop$purbn
## weights: nb2listw(urbantrop.8)
##
## Moran I statistic standard deviate = 4.3142, p-value = 8.01e-06
## alternative hypothesis: greater
## sample estimates:
## Moran I statistic Expectation Variance
## 0.0999670752 -0.0027397260 0.0005667629urbantrop.I<-spline.correlog(x=trop$X, y=trop$Y, z=trop$purbn, resamp=5, latlon=TRUE, quiet=TRUE, xmax=1000)
plot(urbantrop.I) #low spatial autocorrelation; Moran's I 0.1
Fitting a variogram model
#*Mediterranean*
medNVar <- variogram(N~1, med)
summary(medNVar)## np dist gamma dir.hor dir.ver
## Min. : 322 Min. : 182.4 Min. : 46008 Min. :0 Min. :0
## 1st Qu.:1273 1st Qu.:1153.1 1st Qu.:173354 1st Qu.:0 1st Qu.:0
## Median :2719 Median :2174.8 Median :197587 Median :0 Median :0
## Mean :2594 Mean :2171.2 Mean :181244 Mean :0 Mean :0
## 3rd Qu.:3473 3rd Qu.:3183.0 3rd Qu.:215909 3rd Qu.:0 3rd Qu.:0
## Max. :4905 Max. :4183.9 Max. :244358 Max. :0 Max. :0
## id
## var1:15
##
##
##
##
## plot(medNVar, pch=20, cex=1.5,col="magenta",
ylab=expression("Semivariance ("*gamma*")"),
xlab="Distance (km)", main = "Native Species Richness")
sph.model <- vgm(psill=2e05, model="Sph", range=1800, nugget=5e04)
sph.fit <-fit.variogram(object = medNVar, model = sph.model)
plot(medNVar, pch=20, cex=1.25, col="magenta",
ylab=expression("Semivariance ("*gamma*")"),
xlab="Distance (km)", main = "Native Species Richness", model=sph.fit)
med.NVar <- autofitVariogram(N~1,med)
summary(med.NVar)## Experimental variogram:
## np dist gamma dir.hor dir.ver id
## 1 274 87.51469 30452.21 0 0 var1
## 2 1344 139.09962 39257.55 0 0 var1
## 3 1429 236.17040 54754.35 0 0 var1
## 4 2217 350.74833 64935.59 0 0 var1
## 5 2023 481.77031 83662.64 0 0 var1
## 6 2348 613.02305 96873.68 0 0 var1
## 7 7663 908.97277 141260.52 0 0 var1
## 8 5711 1356.99636 193314.44 0 0 var1
## 9 6830 1927.71159 202549.90 0 0 var1
## 10 4856 2587.57792 209028.30 0 0 var1
## 11 2792 3288.06516 234103.68 0 0 var1
## 12 1512 3942.12332 205616.78 0 0 var1
##
## Fitted variogram model:
## model psill range
## 1 Nug 18814.84 0.000
## 2 Sph 192658.90 2031.973
## Sums of squares betw. var. model and sample var.[1] 1336128plot(med.NVar)
med.cultVar <- autofitVariogram(pcult~1,med)
summary(med.cultVar)## Experimental variogram:
## np dist gamma dir.hor dir.ver id
## 1 274 87.51469 0.01523781 0 0 var1
## 2 1344 139.09962 0.02399912 0 0 var1
## 3 1429 236.17040 0.04091604 0 0 var1
## 4 2217 350.74833 0.04254972 0 0 var1
## 5 2023 481.77031 0.04801703 0 0 var1
## 6 2348 613.02305 0.05594117 0 0 var1
## 7 7663 908.97277 0.06026914 0 0 var1
## 8 5711 1356.99636 0.05752098 0 0 var1
## 9 6830 1927.71159 0.04867414 0 0 var1
## 10 4856 2587.57792 0.04718771 0 0 var1
## 11 2792 3288.06516 0.04994671 0 0 var1
## 12 1512 3942.12332 0.06248566 0 0 var1
##
## Fitted variogram model:
## model psill range kappa
## 1 Nug 6.204688e-05 0.0000 0.0
## 2 Ste 5.601899e-02 275.3207 0.8
## Sums of squares betw. var. model and sample var.[1] 1.102465e-06plot(med.cultVar)
med.pastVar <- autofitVariogram(ppast~1,med)
summary(med.pastVar)## Experimental variogram:
## np dist gamma dir.hor dir.ver id
## 1 274 87.51469 0.01694635 0 0 var1
## 2 1344 139.09962 0.02497053 0 0 var1
## 3 1429 236.17040 0.03905366 0 0 var1
## 4 2217 350.74833 0.03850254 0 0 var1
## 5 2023 481.77031 0.04850897 0 0 var1
## 6 2348 613.02305 0.06286127 0 0 var1
## 7 7663 908.97277 0.10022936 0 0 var1
## 8 5711 1356.99636 0.09487425 0 0 var1
## 9 6830 1927.71159 0.07181500 0 0 var1
## 10 4856 2587.57792 0.07856138 0 0 var1
## 11 2792 3288.06516 0.08154839 0 0 var1
## 12 1512 3942.12332 0.07680450 0 0 var1
##
## Fitted variogram model:
## model psill range
## 1 Nug 0.009415596 0.000
## 2 Sph 0.083309543 1190.851
## Sums of squares betw. var. model and sample var.[1] 5.309736e-06plot(med.pastVar)
med.urbnVar <- autofitVariogram(purbn~1,med)
summary(med.urbnVar)## Experimental variogram:
## np dist gamma dir.hor dir.ver id
## 1 274 87.51469 0.0004516839 0 0 var1
## 2 1344 139.09962 0.0012522535 0 0 var1
## 3 1429 236.17040 0.0017103363 0 0 var1
## 4 2217 350.74833 0.0009054361 0 0 var1
## 5 2023 481.77031 0.0007542836 0 0 var1
## 6 2348 613.02305 0.0005974064 0 0 var1
## 7 7663 908.97277 0.0004998116 0 0 var1
## 8 5711 1356.99636 0.0005852270 0 0 var1
## 9 6830 1927.71159 0.0005452296 0 0 var1
## 10 4856 2587.57792 0.0003778553 0 0 var1
## 11 2792 3288.06516 0.0003835389 0 0 var1
## 12 1512 3942.12332 0.0004286593 0 0 var1
##
## Fitted variogram model:
## model psill range
## 1 Nug 0.000000000 0.0000
## 2 Sph 0.001085034 152.4056
## Sums of squares betw. var. model and sample var.[1] 2.545932e-08plot(med.urbnVar)
#*Tropical*
tropNVar <- variogram(N~1, trop)
summary(tropNVar)## np dist gamma dir.hor dir.ver
## Min. : 496 Min. : 200.8 Min. : 48704 Min. :0 Min. :0
## 1st Qu.:1112 1st Qu.:1271.4 1st Qu.:157321 1st Qu.:0 1st Qu.:0
## Median :1779 Median :2413.3 Median :328181 Median :0 Median :0
## Mean :1760 Mean :2401.5 Mean :346814 Mean :0 Mean :0
## 3rd Qu.:2550 3rd Qu.:3500.9 3rd Qu.:484807 3rd Qu.:0 3rd Qu.:0
## Max. :3716 Max. :4645.5 Max. :732273 Max. :0 Max. :0
## id
## var1:15
##
##
##
##
## plot(tropNVar, pch=20, cex=1.5,col="gold2",
ylab=expression("Semivariance ("*gamma*")"),
xlab="Distance (km)", main = "Native Species Richness")
Gau.model.trop.n <- vgm(psill=6e05, model="Gau", range=3700, nugget=2e05)
Gau.fit.trop.n <-fit.variogram(object = tropNVar, model = Gau.model.trop.n)
plot(tropNVar, pch=20, cex=1.25, col="gold2",
ylab=expression("Semivariance ("*gamma*")"),
xlab="Distance (km)", main = "Native Species Richness", model=Gau.model.trop.n)
trop.NVar <- autofitVariogram(N~1,trop) #Gaussian model
summary(trop.NVar)## Experimental variogram:
## np dist gamma dir.hor dir.ver id
## 1 563 94.40304 40123.05 0 0 var1
## 2 855 172.84714 40895.60 0 0 var1
## 3 1006 260.18125 57704.64 0 0 var1
## 4 1876 377.84163 61102.59 0 0 var1
## 5 1617 523.98659 87792.51 0 0 var1
## 6 1427 675.31568 100138.94 0 0 var1
## 7 3355 988.84732 113765.19 0 0 var1
## 8 2026 1504.84064 218548.08 0 0 var1
## 9 3558 2163.26130 317374.53 0 0 var1
## 10 6152 2904.51558 462516.80 0 0 var1
## 11 2539 3548.52640 687036.03 0 0 var1
## 12 2121 4600.60958 422581.86 0 0 var1
##
## Fitted variogram model:
## model psill range kappa
## 1 Nug 39987.12 0.000 0.0
## 2 Ste 764601.72 4030.748 1.4
## Sums of squares betw. var. model and sample var.[1] 14650696plot(trop.NVar)
trop.cultVar <- autofitVariogram(pcult~1,trop)
summary(trop.cultVar)## Experimental variogram:
## np dist gamma dir.hor dir.ver id
## 1 563 94.40304 0.005520187 0 0 var1
## 2 855 172.84714 0.009219947 0 0 var1
## 3 1006 260.18125 0.013957716 0 0 var1
## 4 1876 377.84163 0.019267134 0 0 var1
## 5 1617 523.98659 0.024464709 0 0 var1
## 6 1427 675.31568 0.028584906 0 0 var1
## 7 3355 988.84732 0.032754399 0 0 var1
## 8 2026 1504.84064 0.033871676 0 0 var1
## 9 3558 2163.26130 0.045135233 0 0 var1
## 10 6152 2904.51558 0.044986124 0 0 var1
## 11 2539 3548.52640 0.050057785 0 0 var1
## 12 2121 4600.60958 0.068066260 0 0 var1
##
## Fitted variogram model:
## model psill range kappa
## 1 Nug 0.00000000 0.000 0.0
## 2 Ste 0.04629835 1040.433 0.5
## Sums of squares betw. var. model and sample var.[1] 1.165952e-07plot(trop.cultVar)
trop.pastVar <- autofitVariogram(ppast~1,trop)
summary(trop.pastVar)## Experimental variogram:
## np dist gamma dir.hor dir.ver id
## 1 563 94.40304 0.002160417 0 0 var1
## 2 855 172.84714 0.003356253 0 0 var1
## 3 1006 260.18125 0.003676607 0 0 var1
## 4 1876 377.84163 0.003992843 0 0 var1
## 5 1617 523.98659 0.005166124 0 0 var1
## 6 1427 675.31568 0.004535636 0 0 var1
## 7 3355 988.84732 0.003234269 0 0 var1
## 8 2026 1504.84064 0.007418553 0 0 var1
## 9 3558 2163.26130 0.011670356 0 0 var1
## 10 6152 2904.51558 0.005783572 0 0 var1
## 11 2539 3548.52640 0.009077099 0 0 var1
## 12 2121 4600.60958 0.072506344 0 0 var1
##
## Fitted variogram model:
## model psill range kappa
## 1 Nug 0.002903069 0.00 0.0
## 2 Ste 4.491055464 97555.49 1.6
## Sums of squares betw. var. model and sample var.[1] 3.493877e-07plot(trop.pastVar)
trop.urbnVar <- autofitVariogram(purbn~1,trop)
summary(trop.urbnVar)## Experimental variogram:
## np dist gamma dir.hor dir.ver id
## 1 563 94.40304 6.414682e-05 0 0 var1
## 2 855 172.84714 4.383328e-05 0 0 var1
## 3 1006 260.18125 5.620040e-05 0 0 var1
## 4 1876 377.84163 5.929421e-05 0 0 var1
## 5 1617 523.98659 6.007356e-05 0 0 var1
## 6 1427 675.31568 6.447135e-05 0 0 var1
## 7 3355 988.84732 9.359233e-05 0 0 var1
## 8 2026 1504.84064 1.031828e-04 0 0 var1
## 9 3558 2163.26130 4.825537e-05 0 0 var1
## 10 6152 2904.51558 7.303613e-05 0 0 var1
## 11 2539 3548.52640 1.129056e-04 0 0 var1
## 12 2121 4600.60958 8.477464e-05 0 0 var1
##
## Fitted variogram model:
## model psill range kappa
## 1 Nug 5.708341e-05 0.000 0
## 2 Ste 2.901931e-05 1024.694 10
## Sums of squares betw. var. model and sample var.[1] 1.182197e-11plot(trop.urbnVar)
#*Grassland*
grassNVar <- variogram(N~1, grass)
summary(grassNVar)## np dist gamma dir.hor
## Min. : 41122 Min. : 212.5 Min. : 18370 Min. :0
## 1st Qu.: 95860 1st Qu.:1276.3 1st Qu.: 93151 1st Qu.:0
## Median :125219 Median :2391.1 Median :115775 Median :0
## Mean :113920 Mean :2395.0 Mean :103930 Mean :0
## 3rd Qu.:139085 3rd Qu.:3505.6 3rd Qu.:127562 3rd Qu.:0
## Max. :146365 Max. :4620.1 Max. :156265 Max. :0
## dir.ver id
## Min. :0 var1:15
## 1st Qu.:0
## Median :0
## Mean :0
## 3rd Qu.:0
## Max. :0plot(grassNVar, pch=20, cex=1.5,col="firebrick", cutoff=2500,
ylab=expression("Semivariance ("*gamma*")"),
xlab="Distance (km)", main = "Native Species Richness")
Mat.model.grass.n <- vgm(psill=126144, model="Sph", range=2347, nugget=1153)
Mat.fit.grass.n <-fit.variogram(object = grassNVar, model = Mat.model.grass.n)
plot(grassNVar, pch=20, cex=1.25, col="firebrick",
ylab=expression("Semivariance ("*gamma*")"),
xlab="Distance (km)", main = "Native Species Richness", model=Mat.model.grass.n)
grass.NVar <- autofitVariogram(N~1,grass)
summary(grass.NVar)## Experimental variogram:
## np dist gamma dir.hor dir.ver id
## 1 5804 93.64138 9175.603 0 0 var1
## 2 12817 168.26411 14663.213 0 0 var1
## 3 17631 256.56728 20784.604 0 0 var1
## 4 40539 377.96649 30139.800 0 0 var1
## 5 48730 527.13390 42324.066 0 0 var1
## 6 56877 678.49745 53653.644 0 0 var1
## 7 211549 1010.22379 77354.396 0 0 var1
## 8 227826 1508.43676 105219.729 0 0 var1
## 9 337446 2130.74662 126410.245 0 0 var1
## 10 311651 2882.49709 124316.891 0 0 var1
## 11 260525 3623.17991 113090.890 0 0 var1
## 12 219770 4476.02031 145884.998 0 0 var1
##
## Fitted variogram model:
## model psill range
## 1 Nug 1152.23 0.000
## 2 Sph 124984.99 2346.725
## Sums of squares betw. var. model and sample var.[1] 8708604plot(grass.NVar)
grass.cultVar <- autofitVariogram(pcult~1,grass)
summary(grass.cultVar)## Experimental variogram:
## np dist gamma dir.hor dir.ver id
## 1 5804 93.64138 0.003148945 0 0 var1
## 2 12817 168.26411 0.005656262 0 0 var1
## 3 17631 256.56728 0.008272418 0 0 var1
## 4 40539 377.96649 0.010843040 0 0 var1
## 5 48730 527.13390 0.012810019 0 0 var1
## 6 56877 678.49745 0.013180807 0 0 var1
## 7 211549 1010.22379 0.014549152 0 0 var1
## 8 227826 1508.43676 0.016847358 0 0 var1
## 9 337446 2130.74662 0.017755577 0 0 var1
## 10 311651 2882.49709 0.017054447 0 0 var1
## 11 260525 3623.17991 0.019931636 0 0 var1
## 12 219770 4476.02031 0.020183370 0 0 var1
##
## Fitted variogram model:
## model psill range
## 1 Nug 0.00000000 0.0000
## 2 Exp 0.01722531 417.0169
## Sums of squares betw. var. model and sample var.[1] 8.321684e-07plot(grass.cultVar)
grass.pastVar <- autofitVariogram(ppast~1,grass)
summary(grass.pastVar)## Experimental variogram:
## np dist gamma dir.hor dir.ver id
## 1 5804 93.64138 0.01088045 0 0 var1
## 2 12817 168.26411 0.02156288 0 0 var1
## 3 17631 256.56728 0.03154042 0 0 var1
## 4 40539 377.96649 0.04003607 0 0 var1
## 5 48730 527.13390 0.04485294 0 0 var1
## 6 56877 678.49745 0.05104765 0 0 var1
## 7 211549 1010.22379 0.05925459 0 0 var1
## 8 227826 1508.43676 0.06195969 0 0 var1
## 9 337446 2130.74662 0.07155919 0 0 var1
## 10 311651 2882.49709 0.06865372 0 0 var1
## 11 260525 3623.17991 0.05993260 0 0 var1
## 12 219770 4476.02031 0.05713138 0 0 var1
##
## Fitted variogram model:
## model psill range kappa
## 1 Nug 0.0003143646 0.0000 0.0
## 2 Ste 0.0644507537 546.6965 0.6
## Sums of squares betw. var. model and sample var.[1] 9.926226e-06plot(grass.pastVar)
grass.urbnVar <- autofitVariogram(purbn~1,grass) #What in the world is happening?
summary(grass.urbnVar)## Experimental variogram:
## np dist gamma dir.hor dir.ver id
## 1 5804 93.64138 1.829923e-05 0 0 var1
## 2 12817 168.26411 1.399609e-05 0 0 var1
## 3 17631 256.56728 1.228600e-05 0 0 var1
## 4 40539 377.96649 1.001973e-05 0 0 var1
## 5 48730 527.13390 1.024027e-05 0 0 var1
## 6 56877 678.49745 9.564109e-06 0 0 var1
## 7 211549 1010.22379 9.914540e-06 0 0 var1
## 8 227826 1508.43676 1.105793e-05 0 0 var1
## 9 337446 2130.74662 1.021372e-05 0 0 var1
## 10 311651 2882.49709 8.950448e-06 0 0 var1
## 11 260525 3623.17991 1.950155e-05 0 0 var1
## 12 219770 4476.02031 1.405007e-05 0 0 var1
##
## Fitted variogram model:
## model psill range kappa
## 1 Nug 9.277551e-06 0.000 0.00
## 2 Ste 5.979734e-06 2350.574 0.05
## Sums of squares betw. var. model and sample var.[1] 3.496303e-11plot(grass.urbnVar)
OLS and GLM for Mediterranean Sub-Set
lm(N~ppast+pcult+purbn, data=med)##
## Call:
## lm(formula = N ~ ppast + pcult + purbn, data = med)
##
## Coefficients:
## (Intercept) ppast pcult purbn
## 1238.6 -224.0 490.0 -141.4med.lm<-lm(N~ppast+pcult+purbn, data=med)
summary(med.lm) #Looks like crops have the biggest influence on native species in mediterranean climates. However, let's look at the residuals since they were all auto correlated based on previous Moran's I statistics and variograms.##
## Call:
## lm(formula = N ~ ppast + pcult + purbn, data = med)
##
## Residuals:
## Min 1Q Median 3Q Max
## -815.06 -335.84 -88.43 210.16 2231.33
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 1238.59 83.29 14.870 < 2e-16 ***
## ppast -223.99 128.88 -1.738 0.082971 .
## pcult 490.03 143.00 3.427 0.000673 ***
## purbn -141.38 475.20 -0.298 0.766214
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 503.5 on 407 degrees of freedom
## Multiple R-squared: 0.08301, Adjusted R-squared: 0.07625
## F-statistic: 12.28 on 3 and 407 DF, p-value: 1.046e-07med$lmResidsmed<-residuals(med.lm)
bubble(med, zcol='lmResidsmed') #Autocorrelated residuals! The linear model is not clean, and using a GLS may be more appropriate.
med.gls <- knn2nb(knearneigh(med,k=8))
moran.test(med$lmResidsmed,nb2listw(med.gls)) #Residuals are autocorrelated; OLS may not be the model best suited for this analysis. ##
## Moran I test under randomisation
##
## data: med$lmResidsmed
## weights: nb2listw(med.gls)
##
## Moran I statistic standard deviate = 34.819, p-value < 2.2e-16
## alternative hypothesis: greater
## sample estimates:
## Moran I statistic Expectation Variance
## 0.8034568373 -0.0024390244 0.0005356955med.resids <- spline.correlog(x=coordinates(med)[,1], y=coordinates(med)[,2],z=med$lmResidsmed, resamp=50, quiet=TRUE, xmax=200)
plot(med.resids)
lmResidsmedVar <- variogram(lmResidsmed~1, data=med)
plot(lmResidsmedVar)
sph.model.med <- vgm(psill=1.5e5, model="Sph", range=1800, nugget=50000)
sph.fit.med <- fit.variogram(object = lmResidsmedVar, model = sph.model.med)
plot(lmResidsmedVar,model=sph.fit.med) # Residuals are autocorrelated up to around 2.2km
gls1.med <- gls(N~ppast+pcult+purbn,data=med)
summary(gls1.med)## Generalized least squares fit by REML
## Model: N ~ ppast + pcult + purbn
## Data: med
## AIC BIC logLik
## 6241.405 6261.449 -3115.702
##
## Coefficients:
## Value Std.Error t-value p-value
## (Intercept) 1238.5916 83.2940 14.870123 0.0000
## ppast -223.9855 128.8769 -1.737980 0.0830
## pcult 490.0254 143.0042 3.426650 0.0007
## purbn -141.3847 475.1961 -0.297529 0.7662
##
## Correlation:
## (Intr) ppast pcult
## ppast -0.881
## pcult -0.832 0.629
## purbn -0.315 0.258 0.182
##
## Standardized residuals:
## Min Q1 Med Q3 Max
## -1.6186516 -0.6669605 -0.1756071 0.4173561 4.4312690
##
## Residual standard error: 503.5423
## Degrees of freedom: 411 total; 407 residualcs1.med <-corSpher(c(2030,40000/175000),form=~Y+X,nugget=TRUE) #Range (2032, Nugget:Sill ratio (40000/175000))
gls2.med <- update(gls1.med,correlation=cs1.med) #Update the model with proper correlation structure
summary(gls1.med) #AIC 6241## Generalized least squares fit by REML
## Model: N ~ ppast + pcult + purbn
## Data: med
## AIC BIC logLik
## 6241.405 6261.449 -3115.702
##
## Coefficients:
## Value Std.Error t-value p-value
## (Intercept) 1238.5916 83.2940 14.870123 0.0000
## ppast -223.9855 128.8769 -1.737980 0.0830
## pcult 490.0254 143.0042 3.426650 0.0007
## purbn -141.3847 475.1961 -0.297529 0.7662
##
## Correlation:
## (Intr) ppast pcult
## ppast -0.881
## pcult -0.832 0.629
## purbn -0.315 0.258 0.182
##
## Standardized residuals:
## Min Q1 Med Q3 Max
## -1.6186516 -0.6669605 -0.1756071 0.4173561 4.4312690
##
## Residual standard error: 503.5423
## Degrees of freedom: 411 total; 407 residualsummary(gls2.med) #AIC 5549## Generalized least squares fit by REML
## Model: N ~ ppast + pcult + purbn
## Data: med
## AIC BIC logLik
## 5549.008 5577.069 -2767.504
##
## Correlation Structure: Spherical spatial correlation
## Formula: ~Y + X
## Parameter estimate(s):
## range nugget
## 136.85981123 0.00147364
##
## Coefficients:
## Value Std.Error t-value p-value
## (Intercept) 1674.0475 768.6473 2.177914 0.0300
## ppast -428.2755 95.2263 -4.497450 0.0000
## pcult -431.7521 97.4145 -4.432114 0.0000
## purbn -273.7530 210.5963 -1.299894 0.1944
##
## Correlation:
## (Intr) ppast pcult
## ppast -0.051
## pcult -0.053 0.651
## purbn -0.017 0.186 0.158
##
## Standardized residuals:
## Min Q1 Med Q3 Max
## -0.64203466 -0.33999642 -0.08415708 0.11371311 1.27896474
##
## Residual standard error: 1623.974
## Degrees of freedom: 411 total; 407 residualmed$glsResids.med <-residuals(gls2.med,type="normalized")
bubble(med,zcol='glsResids.med')
moran.test(med$glsResids.med,nb2listw(med.gls)) #Woot! Moran's I; 0.03. ##
## Moran I test under randomisation
##
## data: med$glsResids.med
## weights: nb2listw(med.gls)
##
## Moran I statistic standard deviate = 1.5009, p-value = 0.06669
## alternative hypothesis: greater
## sample estimates:
## Moran I statistic Expectation Variance
## 0.0323591888 -0.0024390244 0.0005375588residsI.pmed <-spline.correlog(x=coordinates(med)[,1], y=coordinates(med)[,2],z=med$glsResids.med, resamp=50, quiet=TRUE)
plot(residsI.pmed)
Looks like crop usage and pastures are significant predictors of native species richness in the Mediterranean.
Tropical OLS and GLM
lm(N~ppast+pcult+purbn, data=trop)##
## Call:
## lm(formula = N ~ ppast + pcult + purbn, data = trop)
##
## Coefficients:
## (Intercept) ppast pcult purbn
## 1844.9 118.3 -1054.6 8986.6trop.lm<-lm(N~ppast+pcult+purbn, data=trop)
summary(trop.lm) #Looks like crops and urbanization are big predictors for native species richness, however, I'm assuming these data are autocorrelated.##
## Call:
## lm(formula = N ~ ppast + pcult + purbn, data = trop)
##
## Residuals:
## Min 1Q Median 3Q Max
## -1290.98 -420.40 -42.66 305.21 2047.64
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 1844.92 75.08 24.572 < 2e-16 ***
## ppast 118.29 169.19 0.699 0.48493
## pcult -1054.65 143.41 -7.354 1.29e-12 ***
## purbn 8986.57 2990.99 3.005 0.00284 **
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 504.6 on 362 degrees of freedom
## Multiple R-squared: 0.2238, Adjusted R-squared: 0.2174
## F-statistic: 34.79 on 3 and 362 DF, p-value: < 2.2e-16trop$lmResidstrop<-residuals(trop.lm)
bubble(trop, zcol='lmResidstrop') #Autocorrelated residuals! The linear models is not clean, and using a GLS may be more appropriate.
trop.gls <- knn2nb(knearneigh(trop,k=8))
moran.test(trop$lmResidstrop,nb2listw(trop.gls)) #Residuals are positively autocorrelated: Moran's I is 0.63##
## Moran I test under randomisation
##
## data: trop$lmResidstrop
## weights: nb2listw(trop.gls)
##
## Moran I statistic standard deviate = 25.837, p-value < 2.2e-16
## alternative hypothesis: greater
## sample estimates:
## Moran I statistic Expectation Variance
## 0.6295687233 -0.0027397260 0.0005989311trop.resids <- spline.correlog(x=coordinates(trop)[,1], y=coordinates(trop)[,2],z=trop$lmResidstrop, resamp=50, quiet=TRUE, xmax=300)
plot(trop.resids) #Oh yikes, that's a mess.... 
lmResidstropVar <- variogram(lmResidstrop~1, data=trop)
plot(lmResidstropVar)
sph.model.trop <- vgm(psill=1.5e5, model="Sph", range=1800, nugget=50000)
sph.fit.trop <- fit.variogram(object = lmResidstropVar, model = sph.model.trop)
plot(lmResidstropVar,model=sph.fit.trop)
gls1.trop <- gls(N~ppast+pcult+purbn,data=trop)
summary(gls1.trop) #The GLS model indicates crops and urbanization are significant predictors for native species richness## Generalized least squares fit by REML
## Model: N ~ ppast + pcult + purbn
## Data: trop
## AIC BIC logLik
## 5550.983 5570.442 -2770.492
##
## Coefficients:
## Value Std.Error t-value p-value
## (Intercept) 1844.924 75.0836 24.571590 0.0000
## ppast 118.285 169.1938 0.699113 0.4849
## pcult -1054.650 143.4079 -7.354196 0.0000
## purbn 8986.571 2990.9886 3.004549 0.0028
##
## Correlation:
## (Intr) ppast pcult
## ppast -0.809
## pcult -0.850 0.639
## purbn -0.047 -0.056 -0.195
##
## Standardized residuals:
## Min Q1 Med Q3 Max
## -2.55826405 -0.83308822 -0.08453566 0.60481777 4.05768681
##
## Residual standard error: 504.6325
## Degrees of freedom: 366 total; 362 residualcs1.trop <-corSpher(c(3847,40237/780235),form=~Y+X,nugget=TRUE)
gls2.trop <- update(gls1.trop,correlation=cs1.trop)
summary(gls1.trop) #AIC 5550## Generalized least squares fit by REML
## Model: N ~ ppast + pcult + purbn
## Data: trop
## AIC BIC logLik
## 5550.983 5570.442 -2770.492
##
## Coefficients:
## Value Std.Error t-value p-value
## (Intercept) 1844.924 75.0836 24.571590 0.0000
## ppast 118.285 169.1938 0.699113 0.4849
## pcult -1054.650 143.4079 -7.354196 0.0000
## purbn 8986.571 2990.9886 3.004549 0.0028
##
## Correlation:
## (Intr) ppast pcult
## ppast -0.809
## pcult -0.850 0.639
## purbn -0.047 -0.056 -0.195
##
## Standardized residuals:
## Min Q1 Med Q3 Max
## -2.55826405 -0.83308822 -0.08453566 0.60481777 4.05768681
##
## Residual standard error: 504.6325
## Degrees of freedom: 366 total; 362 residualsummary(gls2.trop) #AIC 5125## Generalized least squares fit by REML
## Model: N ~ ppast + pcult + purbn
## Data: trop
## AIC BIC logLik
## 5125.888 5153.13 -2555.944
##
## Correlation Structure: Spherical spatial correlation
## Formula: ~Y + X
## Parameter estimate(s):
## range nugget
## 1.801209e+02 7.197538e-03
##
## Coefficients:
## Value Std.Error t-value p-value
## (Intercept) 1309.393 1077.1048 1.215659 0.2249
## ppast -37.609 216.4767 -0.173734 0.8622
## pcult -342.223 151.9021 -2.252918 0.0249
## purbn 3664.410 1509.8242 2.427044 0.0157
##
## Correlation:
## (Intr) ppast pcult
## ppast -0.091
## pcult -0.030 0.071
## purbn -0.006 -0.057 -0.101
##
## Standardized residuals:
## Min Q1 Med Q3 Max
## -0.41247498 -0.01488615 0.18239818 0.36396032 1.39786999
##
## Residual standard error: 1826.823
## Degrees of freedom: 366 total; 362 residualtrop$glsResids.trop <-residuals(gls2.trop,type="normalized")
bubble(trop,zcol='glsResids.trop')
moran.test(trop$glsResids.trop,nb2listw(trop.gls)) #Slight negative autocorrelation; Moran's I: -0.01##
## Moran I test under randomisation
##
## data: trop$glsResids.trop
## weights: nb2listw(trop.gls)
##
## Moran I statistic standard deviate = -0.25146, p-value = 0.5993
## alternative hypothesis: greater
## sample estimates:
## Moran I statistic Expectation Variance
## -0.0088571305 -0.0027397260 0.0005918302residsI.trop <-spline.correlog(x=coordinates(trop)[,1], y=coordinates(trop)[,2],z=trop$glsResids.trop, resamp=50, quiet=TRUE)
plot(residsI.trop) #Much better...
Grassland
lm(N~ppast+pcult+purbn, data=grass)##
## Call:
## lm(formula = N ~ ppast + pcult + purbn, data = grass)
##
## Coefficients:
## (Intercept) ppast pcult purbn
## 1201.5 -352.7 272.3 3355.1grass.lm<-lm(N~ppast+pcult+purbn, data=grass)
summary(grass.lm) #Looks like crops, pasture, and urbanization have an impact on native species richness... but wait! Autocorrelation!##
## Call:
## lm(formula = N ~ ppast + pcult + purbn, data = grass)
##
## Residuals:
## Min 1Q Median 3Q Max
## -971.18 -264.66 -40.74 204.01 2238.49
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 1201.51 16.14 74.466 < 2e-16 ***
## ppast -352.70 31.44 -11.218 < 2e-16 ***
## pcult 272.27 61.67 4.415 1.05e-05 ***
## purbn 3355.07 1015.89 3.303 0.000971 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 381 on 2641 degrees of freedom
## Multiple R-squared: 0.05975, Adjusted R-squared: 0.05868
## F-statistic: 55.94 on 3 and 2641 DF, p-value: < 2.2e-16grass$lmResidsgrass<-residuals(grass.lm)
bubble(grass, zcol='lmResidsgrass') #Autocorrelated residuals! The linear model is probably not clean, and using a GLS may be more appropriate.
grass.gls <- knn2nb(knearneigh(grass,k=8))
moran.test(grass$lmResidsgrass,nb2listw(grass.gls)) #Moran's I shows negative autocorrelation of -0.88. Dissimilar values are clustered together. ##
## Moran I test under randomisation
##
## data: grass$lmResidsgrass
## weights: nb2listw(grass.gls)
##
## Moran I statistic standard deviate = 92.504, p-value < 2.2e-16
## alternative hypothesis: greater
## sample estimates:
## Moran I statistic Expectation Variance
## 8.769947e-01 -3.782148e-04 8.995864e-05grass.resids <- spline.correlog(x=coordinates(grass)[,1], y=coordinates(grass)[,2],z=grass$lmResidsgrass, resamp=2, quiet=TRUE, xmax=300)
plot(grass.resids)
lmResidsgrassVar <- variogram(lmResidsgrass~1, data=grass)
plot(lmResidsgrassVar)
sph.model.grass <- vgm(psill=1.2e5, model="Sph", range=2500, nugget=1000)
sph.fit.grass <- fit.variogram(object = lmResidsgrassVar, model = sph.model.grass)
plot(lmResidsgrassVar,model=sph.fit.grass)
gls1.grass <- gls(N~ppast+pcult+purbn,data=grass)
summary(gls1.grass) #The GLS model indicates crops and urbanization are significant predictors for native species richness## Generalized least squares fit by REML
## Model: N ~ ppast + pcult + purbn
## Data: grass
## AIC BIC logLik
## 38909.95 38939.34 -19449.97
##
## Coefficients:
## Value Std.Error t-value p-value
## (Intercept) 1201.512 16.1350 74.46615 0.000
## ppast -352.697 31.4406 -11.21789 0.000
## pcult 272.273 61.6749 4.41465 0.000
## purbn 3355.069 1015.8909 3.30259 0.001
##
## Correlation:
## (Intr) ppast pcult
## ppast -0.822
## pcult -0.335 0.011
## purbn -0.043 0.030 -0.198
##
## Standardized residuals:
## Min Q1 Med Q3 Max
## -2.5487014 -0.6945547 -0.1069146 0.5354006 5.8745620
##
## Residual standard error: 381.0474
## Degrees of freedom: 2645 total; 2641 residualcs1.grass <-corSpher(c(2500,1000/1.2e5),form=~Y+X,nugget=TRUE)
gls2.grass <- update(gls1.grass,correlation=cs1.grass)
summary(gls1.grass) #AIC 38,910## Generalized least squares fit by REML
## Model: N ~ ppast + pcult + purbn
## Data: grass
## AIC BIC logLik
## 38909.95 38939.34 -19449.97
##
## Coefficients:
## Value Std.Error t-value p-value
## (Intercept) 1201.512 16.1350 74.46615 0.000
## ppast -352.697 31.4406 -11.21789 0.000
## pcult 272.273 61.6749 4.41465 0.000
## purbn 3355.069 1015.8909 3.30259 0.001
##
## Correlation:
## (Intr) ppast pcult
## ppast -0.822
## pcult -0.335 0.011
## purbn -0.043 0.030 -0.198
##
## Standardized residuals:
## Min Q1 Med Q3 Max
## -2.5487014 -0.6945547 -0.1069146 0.5354006 5.8745620
##
## Residual standard error: 381.0474
## Degrees of freedom: 2645 total; 2641 residualsummary(gls2.grass) #AIC 5,125## Generalized least squares fit by REML
## Model: N ~ ppast + pcult + purbn
## Data: grass
## AIC BIC logLik
## 32155.84 32196.99 -16070.92
##
## Correlation Structure: Spherical spatial correlation
## Formula: ~Y + X
## Parameter estimate(s):
## range nugget
## 7.147312e+01 3.235081e-09
##
## Coefficients:
## Value Std.Error t-value p-value
## (Intercept) 1040.3906 253.4676 4.104629 0.0000
## ppast -47.8544 19.6979 -2.429422 0.0152
## pcult 93.3392 34.3748 2.715339 0.0067
## purbn 1321.0738 348.7618 3.787897 0.0002
##
## Correlation:
## (Intr) ppast pcult
## ppast -0.018
## pcult -0.015 0.203
## purbn -0.026 0.025 -0.024
##
## Standardized residuals:
## Min Q1 Med Q3 Max
## -1.10550405 -0.32176820 0.01751551 0.39401160 3.24070621
##
## Residual standard error: 735.5896
## Degrees of freedom: 2645 total; 2641 residualgrass$glsResids.grass <-residuals(gls2.grass,type="normalized")
bubble(grass,zcol='glsResids.grass')
moran.test(grass$glsResids.grass,nb2listw(grass.gls)) ##
## Moran I test under randomisation
##
## data: grass$glsResids.grass
## weights: nb2listw(grass.gls)
##
## Moran I statistic standard deviate = -2.1282, p-value = 0.9833
## alternative hypothesis: greater
## sample estimates:
## Moran I statistic Expectation Variance
## -2.052530e-02 -3.782148e-04 8.962113e-05residsI.grass <-spline.correlog(x=coordinates(grass)[,1], y=coordinates(grass)[,2],z=grass$glsResids.grass, resamp=2, quiet=TRUE)
plot(residsI.grass) #Much better..