Final Project: All Is Not Loss: Plant Biodiversity in the Anthropocene.

Looking at three subsets of the dataset: Mediterranean and Tropical and Subtropical Dry Broadleaf Forest and Temperate Grassland biomes. The overall objectve is to assess whether total population (persons per km^2) and land cover type (crops, pasture, and urban) affect native species (N).

Objective 1

Assess whether data between each land cover type and biome are autocorrelated using variograms and calculating Moran’s I.

Objective 2

Compare the ordinary least squares (OLS) and generalalized least squares (GLM) regression models for spatially autocorrelated data. I’m assuming this combination of variables will be autocorrelated based on an initial exploratory analysis of the Mediterranean biome

Objective 3

Evaluate wheter the (spatial autoregressive) SAR model for each biome and cover is appropriate given the data.

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.0-4, (SVN revision 548)
##  Geospatial Data Abstraction Library extensions to R successfully loaded
##  Loaded GDAL runtime: GDAL 1.11.2, released 2015/02/10
##  Path to GDAL shared files: /Users/shannoncall/Library/R/3.2/library/rgdal/gdal
##  Loaded PROJ.4 runtime: Rel. 4.9.1, 04 March 2015, [PJ_VERSION: 491]
##  Path to PROJ.4 shared files: /Users/shannoncall/Library/R/3.2/library/rgdal/proj
##  Linking to sp version: 1.1-1
library(raster)
library(gstat)
library(ggplot2)
library(spatstat)
## Loading required package: nlme
## 
## Attaching package: 'nlme'
## The following object is masked from 'package:raster':
## 
##     getData
## 
## spatstat 1.45-0       (nickname: 'One Lakh') 
## For an introduction to spatstat, type 'beginner'
## 
## Note: spatstat version 1.45-0 is out of date by more than 12 weeks; a newer version should be available.
## 
## Attaching package: 'spatstat'
## The following object is masked from 'package:gstat':
## 
##     idw
## The following objects are masked from 'package:raster':
## 
##     area, rotate, shift
library(ncf)
library(nlme)
library(spdep)
## Loading required package: Matrix
library(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("/Users/shannoncall/Desktop/ESCI597A/FinalProject")
ellis<- read.csv("/Users/shannoncall/Desktop/ESCI597A/FinalProject/ellis_2012_data.csv")
med <- ellis[3628:4038,] #mediterranean forests, woodlands & shrub
trop <-ellis[13885:14250,] #tropical & subtropical dry broadleaf forests
grass <-ellis[7036:9680,] #temperate grasslands, savannas & shurblands
#####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)
##       L8ID             X                 Y              hexarea    
##  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  
##     landkm2         pland           Olsoncl         Olsonrlm   
##  Min.   :4025   Min.   :0.5170   Min.   :6   Afrotropic :1840  
##  1st Qu.:7568   1st Qu.:0.9710   1st Qu.:7   Australasia: 283  
##  Median :7719   Median :0.9910   Median :7   Indo-Malay :   1  
##  Mean   :7643   Mean   :0.9809   Mean   :7   Nearctic   :  10  
##  3rd Qu.:7812   3rd Qu.:1.0030   3rd Qu.:7   Neotropic  : 510  
##  Max.   :8295   Max.   :1.0650   Max.   :7   Palearctic :   1  
##     anthrov2        totpdens          pcult             ppast       
##  Min.   :11.00   Min.   :  0.00   Min.   :0.00000   Min.   :0.0000  
##  1st Qu.:41.00   1st Qu.:  1.80   1st Qu.:0.00500   1st Qu.:0.2210  
##  Median :42.00   Median :  7.50   Median :0.04500   Median :0.4400  
##  Mean   :42.95   Mean   : 22.14   Mean   :0.08971   Mean   :0.4187  
##  3rd Qu.:43.00   3rd Qu.: 21.00   3rd Qu.:0.12000   3rd Qu.:0.5830  
##  Max.   :62.00   Max.   :893.00   Max.   :0.97600   Max.   :1.0000  
##      purbn                HL               N             ASL        
##  Min.   :0.000000   Min.   :0.0000   Min.   : 227   Min.   :  0.00  
##  1st Qu.:0.000000   1st Qu.:0.2200   1st Qu.: 792   1st Qu.: 47.00  
##  Median :0.000000   Median :0.3810   Median :1041   Median : 80.00  
##  Mean   :0.001338   Mean   :0.3702   Mean   :1083   Mean   : 92.32  
##  3rd Qu.:0.001000   3rd Qu.:0.5090   3rd Qu.:1328   3rd Qu.:126.00  
##  Max.   :0.224000   Max.   :0.9900   Max.   :3423   Max.   :761.00  
##        IS               CS              OS              ASI        
##  Min.   : 14.00   Min.   : 0.00   Min.   :  0.00   Min.   : 14.00  
##  1st Qu.: 38.00   1st Qu.: 7.00   1st Qu.:  0.00   1st Qu.: 49.00  
##  Median : 47.00   Median :14.00   Median :  0.00   Median : 65.00  
##  Mean   : 48.33   Mean   :15.28   Mean   : 11.95   Mean   : 75.56  
##  3rd Qu.: 58.00   3rd Qu.:22.00   3rd Qu.:  0.00   3rd Qu.: 81.00  
##  Max.   :125.00   Max.   :54.00   Max.   :400.00   Max.   :483.00  
##       ASR            ASRN             ASLN              ASIN        
##  Min.   : 240   Min.   :0.5190   Min.   :0.00000   Min.   :0.03800  
##  1st Qu.: 755   1st Qu.:0.9440   1st Qu.:0.04400   1st Qu.:0.05200  
##  Median :1016   Median :0.9810   Median :0.08300   Median :0.05800  
##  Mean   :1066   Mean   :0.9821   Mean   :0.08777   Mean   :0.06986  
##  3rd Qu.:1334   3rd Qu.:1.0240   3rd Qu.:0.12000   3rd Qu.:0.06700  
##  Max.   :3521   Max.   :1.5720   Max.   :0.56600   Max.   :0.58800  
##       ISN            ASI.ASL.N      
##  Min.   :0.03700   Min.   :0.04694  
##  1st Qu.:0.04300   1st Qu.:0.10299  
##  Median :0.04500   Median :0.14227  
##  Mean   :0.04573   Mean   :0.15764  
##  3rd Qu.:0.04800   3rd Qu.:0.18769  
##  Max.   :0.08000   Max.   :0.65494
class(med)
## [1] "data.frame"
summary(med)
##       L8ID             X                  Y             hexarea    
##  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  
##     landkm2         pland           Olsoncl          Olsonrlm  
##  Min.   :3942   Min.   :0.5060   Min.   :12   Afrotropic : 12  
##  1st Qu.:6940   1st Qu.:0.8905   1st Qu.:12   Australasia:100  
##  Median :7652   Median :0.9820   Median :12   Indo-Malay :  0  
##  Mean   :7091   Mean   :0.9101   Mean   :12   Nearctic   : 16  
##  3rd Qu.:7739   3rd Qu.:0.9930   3rd Qu.:12   Neotropic  : 23  
##  Max.   :8093   Max.   :1.0390   Max.   :12   Palearctic :260  
##     anthrov2       totpdens           pcult            ppast       
##  Min.   :11.0   Min.   :   0.00   Min.   :0.0000   Min.   :0.0010  
##  1st Qu.:32.0   1st Qu.:   2.90   1st Qu.:0.0485   1st Qu.:0.1725  
##  Median :34.0   Median :  33.50   Median :0.2520   Median :0.2890  
##  Mean   :36.8   Mean   :  89.19   Mean   :0.2664   Mean   :0.3695  
##  3rd Qu.:42.0   3rd Qu.:  92.50   3rd Qu.:0.4015   3rd Qu.:0.5115  
##  Max.   :62.0   Max.   :2167.90   Max.   :0.9850   Max.   :1.0000  
##      purbn               HL               N               ASL       
##  Min.   :0.00000   Min.   :0.0890   Min.   : 473.0   Min.   : 16.0  
##  1st Qu.:0.00000   1st Qu.:0.4085   1st Qu.: 863.5   1st Qu.:100.0  
##  Median :0.00300   Median :0.5230   Median :1217.0   Median :173.0  
##  Mean   :0.01482   Mean   :0.5275   Mean   :1284.3   Mean   :202.9  
##  3rd Qu.:0.00850   3rd Qu.:0.6370   3rd Qu.:1559.0   3rd Qu.:263.0  
##  Max.   :0.82600   Max.   :0.9940   Max.   :3417.0   Max.   :923.0  
##        IS              CS              OS              ASI       
##  Min.   : 48.0   Min.   : 0.00   Min.   :  0.00   Min.   : 53.0  
##  1st Qu.: 79.0   1st Qu.:16.00   1st Qu.:  0.00   1st Qu.: 99.5  
##  Median :105.0   Median :40.00   Median :  0.00   Median :184.0  
##  Mean   :108.2   Mean   :35.52   Mean   : 89.54   Mean   :233.2  
##  3rd Qu.:128.0   3rd Qu.:56.00   3rd Qu.:200.00   3rd Qu.:369.5  
##  Max.   :244.0   Max.   :78.00   Max.   :400.00   Max.   :651.0  
##       ASR            ASRN             ASLN             ASIN       
##  Min.   : 500   Min.   :0.6190   Min.   :0.0180   Min.   :0.0760  
##  1st Qu.: 884   1st Qu.:0.9595   1st Qu.:0.0995   1st Qu.:0.1060  
##  Median :1208   Median :1.0040   Median :0.1380   Median :0.1260  
##  Mean   :1315   Mean   :1.0291   Mean   :0.1544   Mean   :0.1835  
##  3rd Qu.:1644   3rd Qu.:1.0790   3rd Qu.:0.1835   3rd Qu.:0.2480  
##  Max.   :3392   Max.   :1.8180   Max.   :0.6400   Max.   :0.9150  
##       ISN           ASI.ASL.N     
##  Min.   :0.0710   Min.   :0.1086  
##  1st Qu.:0.0820   1st Qu.:0.2196  
##  Median :0.0860   Median :0.2915  
##  Mean   :0.0866   Mean   :0.3379  
##  3rd Qu.:0.0915   3rd Qu.:0.4086  
##  Max.   :0.1020   Max.   :1.0125
class(trop)
## [1] "data.frame"
summary(trop)
##       L8ID             X                 Y             hexarea    
##  Min.   :15977   Min.   :-110.19   Min.   :-21.03   Min.   :7790  
##  1st Qu.:20734   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  
##     landkm2         pland           Olsoncl             Olsonrlm  
##  Min.   :3916   Min.   :0.5030   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.:1.0000   3rd Qu.:2.000   Neotropic  :132  
##  Max.   :8134   Max.   :1.0440   Max.   :2.000   Palearctic :  1  
##     anthrov2        totpdens           pcult             ppast       
##  Min.   :21.00   Min.   :   0.00   Min.   :0.00000   Min.   :0.0000  
##  1st Qu.:23.00   1st Qu.:  16.05   1st Qu.:0.08625   1st Qu.:0.0320  
##  Median :32.00   Median :  86.75   Median :0.28850   Median :0.0715  
##  Mean   :34.08   Mean   : 145.86   Mean   :0.32137   Mean   :0.1897  
##  3rd Qu.:42.00   3rd Qu.: 204.25   3rd Qu.:0.50900   3rd Qu.:0.3180  
##  Max.   :61.00   Max.   :1249.30   Max.   :0.93100   Max.   :0.8460  
##      purbn               HL               N             ASL       
##  Min.   :0.00000   Min.   :0.0030   Min.   : 546   Min.   :  1.0  
##  1st Qu.:0.00000   1st Qu.:0.3123   1st Qu.:1124   1st Qu.:110.2  
##  Median :0.00200   Median :0.4670   Median :1540   Median :161.5  
##  Mean   :0.00479   Mean   :0.4526   Mean   :1572   Mean   :189.7  
##  3rd Qu.:0.00500   3rd Qu.:0.5747   3rd Qu.:1926   3rd Qu.:256.8  
##  Max.   :0.07300   Max.   :0.9330   Max.   :3796   Max.   :680.0  
##        IS              CS              OS              ASI       
##  Min.   : 77.0   Min.   : 4.00   Min.   :  0.00   Min.   :114.0  
##  1st Qu.:162.2   1st Qu.:24.00   1st Qu.:  0.00   1st Qu.:207.0  
##  Median :209.0   Median :34.00   Median :  0.00   Median :260.5  
##  Mean   :211.8   Mean   :32.06   Mean   : 54.64   Mean   :298.5  
##  3rd Qu.:253.0   3rd Qu.:41.00   3rd Qu.:200.00   3rd Qu.:372.8  
##  Max.   :442.0   Max.   :60.00   Max.   :400.00   Max.   :762.0  
##       ASR            ASRN            ASLN            ASIN       
##  Min.   : 651   Min.   :0.755   Min.   :0.001   Min.   :0.0720  
##  1st Qu.:1173   1st Qu.:1.031   1st Qu.:0.076   1st Qu.:0.1470  
##  Median :1671   Median :1.069   Median :0.125   Median :0.1750  
##  Mean   :1680   Mean   :1.069   Mean   :0.130   Mean   :0.1994  
##  3rd Qu.:2080   3rd Qu.:1.106   3rd Qu.:0.164   3rd Qu.:0.2087  
##  Max.   :4071   Max.   :1.327   Max.   :0.434   Max.   :0.5140  
##       ISN          ASI.ASL.N     
##  Min.   :0.041   Min.   :0.1490  
##  1st Qu.:0.131   1st Qu.:0.2283  
##  Median :0.137   Median :0.3095  
##  Mean   :0.138   Mean   :0.3293  
##  3rd Qu.:0.145   3rd Qu.:0.3917  
##  Max.   :0.165   Max.   :0.7683
plot(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:
##   .. ..$ L8ID     : int [1:16805] 7287 7347 7426 7430 7454 7470 7486 7558 7566 7574 ...
##   .. ..$ hexarea  : num [1:16805] 7792 7792 7792 7792 7792 ...
##   .. ..$ landkm2  : 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.5 1.5 1.4 0.2 1 24 2 0.6 2 1.6 ...
##   .. ..$ pcult    : num [1:16805] 0.097 0.067 0.026 0.037 0.022 0.112 0.012 0.028 0.047 0.013 ...
##   .. ..$ ppast    : num [1:16805] 0.238 0.501 0.125 0.855 0.226 0.632 0.763 0.866 0.848 0.916 ...
##   .. ..$ purbn    : num [1:16805] 0 0 0 0 0 0.006 0 0 0 0 ...
##   .. ..$ HL       : num [1:16805] 0.257 0.401 0.109 0.607 0.172 0.539 0.52 0.606 0.612 0.624 ...
##   .. ..$ N        : int [1:16805] 724 580 677 535 795 712 451 591 534 468 ...
##   .. ..$ ASL      : int [1:16805] 23 32 9 52 16 58 35 58 53 48 ...
##   .. ..$ IS       : int [1:16805] 42 35 40 33 46 42 29 36 33 30 ...
##   .. ..$ 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      : int [1:16805] 54 46 51 43 55 259 44 47 52 32 ...
##   .. ..$ ASR      : int [1:16805] 755 594 720 525 833 913 460 580 533 452 ...
##   .. ..$ ASRN     : num [1:16805] 1.043 1.025 1.063 0.982 1.048 ...
##   .. ..$ ASLN     : num [1:16805] 0.032 0.055 0.013 0.098 0.021 0.082 0.078 0.097 0.099 0.102 ...
##   .. ..$ ASIN     : num [1:16805] 0.075 0.08 0.075 0.08 0.069 0.363 0.097 0.079 0.097 0.067 ...
##   .. ..$ ISN      : num [1:16805] 0.058 0.061 0.059 0.062 0.057 0.058 0.063 0.06 0.062 0.063 ...
##   .. ..$ ASI.ASL.N: num [1:16805] 0.1069 0.1343 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 NA
proj4string(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=WGS84
plot(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:
##   .. ..$ L8ID     : int [1:411] 9705 9779 9833 9849 9864 9887 9935 9936 10028 10103 ...
##   .. ..$ hexarea  : num [1:411] 7792 7792 7792 7792 7794 ...
##   .. ..$ landkm2  : 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 0.393 0.159 0.363 0.611 0.306 ...
##   .. ..$ ppast    : num [1:411] 0.112 0.087 0.213 0.103 0.121 0.145 0.173 0.372 0.15 0.129 ...
##   .. ..$ purbn    : num [1:411] 0.071 0.024 0.005 0.023 0.056 0.007 0.01 0.011 0.061 0.007 ...
##   .. ..$ HL       : num [1:411] 0.478 0.338 0.439 0.409 0.556 0.497 0.284 0.621 0.772 0.399 ...
##   .. ..$ N        : int [1:411] 1417 1667 1181 1505 1614 1394 2592 2117 1442 1283 ...
##   .. ..$ ASL      : int [1:411] 173 132 129 151 242 179 168 374 369 124 ...
##   .. ..$ IS       : int [1:411] 119 136 102 125 132 117 195 165 120 109 ...
##   .. ..$ 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      : int [1:411] 582 396 354 375 593 360 435 410 568 360 ...
##   .. ..$ ASR      : int [1:411] 1825 1930 1406 1729 1965 1575 2859 2154 1641 1519 ...
##   .. ..$ ASRN     : num [1:411] 1.29 1.16 1.19 1.15 1.22 ...
##   .. ..$ ASLN     : num [1:411] 0.122 0.079 0.109 0.1 0.15 0.128 0.065 0.176 0.256 0.097 ...
##   .. ..$ ASIN     : num [1:411] 0.411 0.237 0.3 0.249 0.367 0.258 0.168 0.194 0.394 0.281 ...
##   .. ..$ ISN      : num [1:411] 0.084 0.081 0.086 0.083 0.082 0.084 0.075 0.078 0.083 0.085 ...
##   .. ..$ 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 NA
proj4string(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=WGS84
plot(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:
##   .. ..$ L8ID     : int [1:2645] 16191 16257 16369 16425 16537 16600 16698 16875 16987 17064 ...
##   .. ..$ hexarea  : num [1:2645] 7792 7792 7792 7792 7792 ...
##   .. ..$ landkm2  : num [1:2645] 6417 7767 5808 4025 4903 ...
##   .. ..$ pland    : num [1:2645] 0.824 0.997 0.745 0.517 0.629 0.523 0.943 0.821 0.983 0.952 ...
##   .. ..$ 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 0.277 0.533 0.499 0.298 0.522 ...
##   .. ..$ ppast    : num [1:2645] 0.141 0.284 0.251 0.178 0.339 0.232 0.428 0.497 0.029 0.484 ...
##   .. ..$ purbn    : num [1:2645] 0.018 0.189 0.042 0.161 0.224 0.026 0.022 0.011 0.001 0.002 ...
##   .. ..$ HL       : num [1:2645] 0.385 0.792 0.699 0.543 0.897 0.458 0.84 0.842 0.319 0.846 ...
##   .. ..$ N        : int [1:2645] 1344 1132 1288 1360 1525 1519 1320 1209 1164 1267 ...
##   .. ..$ ASL      : int [1:2645] 113 279 251 179 512 158 371 342 78 362 ...
##   .. ..$ IS       : int [1:2645] 58 51 56 59 65 64 57 53 52 56 ...
##   .. ..$ 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      : int [1:2645] 268 463 263 462 475 270 267 263 95 67 ...
##   .. ..$ ASR      : int [1:2645] 1499 1316 1300 1643 1487 1631 1216 1131 1181 971 ...
##   .. ..$ ASRN     : num [1:2645] 1.116 1.162 1.01 1.208 0.975 ...
##   .. ..$ ASLN     : num [1:2645] 0.084 0.246 0.195 0.131 0.336 0.104 0.281 0.283 0.067 0.286 ...
##   .. ..$ ASIN     : num [1:2645] 0.2 0.409 0.204 0.34 0.311 0.178 0.203 0.218 0.081 0.052 ...
##   .. ..$ ISN      : num [1:2645] 0.043 0.045 0.044 0.043 0.042 0.042 0.043 0.044 0.044 0.044 ...
##   .. ..$ 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 NA
proj4string(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=WGS84
plot(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:
##   .. ..$ L8ID     : int [1:366] 16706 16938 17132 17385 17591 17711 17836 17893 17961 18067 ...
##   .. ..$ hexarea  : num [1:366] 7792 7792 7792 7792 7792 ...
##   .. ..$ landkm2  : 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.3 1.5 2.9 1.3 62.6 ...
##   .. ..$ pcult    : num [1:366] 0.029 0.176 0.002 0.083 0.343 0.659 0.124 0.633 0.646 0.047 ...
##   .. ..$ ppast    : num [1:366] 0.57 0.319 0.568 0.394 0.478 0.031 0.513 0.034 0.034 0.263 ...
##   .. ..$ purbn    : num [1:366] 0 0 0 0 0.007 0.003 0.003 0.025 0.003 0 ...
##   .. ..$ HL       : num [1:366] 0.409 0.388 0.381 0.345 0.668 0.683 0.469 0.681 0.672 0.223 ...
##   .. ..$ N        : int [1:366] 849 1216 751 1181 865 751 1104 770 749 1674 ...
##   .. ..$ ASL      : int [1:366] 89 119 72 100 179 161 138 164 156 86 ...
##   .. ..$ IS       : int [1:366] 129 174 117 169 131 117 160 119 117 226 ...
##   .. ..$ 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      : int [1:366] 144 184 134 189 360 144 186 346 157 246 ...
##   .. ..$ ASR      : int [1:366] 905 1280 813 1270 1047 734 1153 952 749 1833 ...
##   .. ..$ ASRN     : num [1:366] 1.06 1.05 1.08 1.07 1.21 ...
##   .. ..$ ASLN     : num [1:366] 0.105 0.098 0.096 0.085 0.207 0.214 0.125 0.213 0.209 0.052 ...
##   .. ..$ ASIN     : num [1:366] 0.17 0.151 0.178 0.16 0.417 0.192 0.169 0.45 0.209 0.147 ...
##   .. ..$ ISN      : num [1:366] 0.152 0.143 0.156 0.144 0.152 0.156 0.145 0.155 0.156 0.135 ...
##   .. ..$ 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 NA
proj4string(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=WGS84
plot(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 across the globe. This helps me visualize the scope of data, where the biomes are located and helps provide another source of information for underlying mechanisms.

#*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")

Lagged Scatter plots per biome and land cover types: 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.046, p-value < 2.2e-16
## alternative hypothesis: greater
## sample estimates:
## Moran I statistic       Expectation          Variance 
##      0.8339409151     -0.0024390244      0.0005383879
Nmed.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.328, p-value < 2.2e-16
## alternative hypothesis: greater
## sample estimates:
## Moran I statistic       Expectation          Variance 
##      0.5629734237     -0.0024390244      0.0005401727
cropmed.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.081, p-value < 2.2e-16
## alternative hypothesis: greater
## sample estimates:
## Moran I statistic       Expectation          Variance 
##      0.6505699770     -0.0024390244      0.0005407834
pastmed.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.7973, p-value < 2.2e-16
## alternative hypothesis: greater
## sample estimates:
## Moran I statistic       Expectation          Variance 
##      0.1661412314     -0.0024390244      0.0003672107
urbanmed.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.729, p-value < 2.2e-16
## alternative hypothesis: greater
## sample estimates:
## Moran I statistic       Expectation          Variance 
##      0.7252913909     -0.0027397260      0.0005997074
Ntrop.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.643, p-value < 2.2e-16
## alternative hypothesis: greater
## sample estimates:
## Moran I statistic       Expectation          Variance 
##      0.8722565393     -0.0027397260      0.0006026645
croptrop.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.791, p-value < 2.2e-16
## alternative hypothesis: greater
## sample estimates:
## Moran I statistic       Expectation          Variance 
##      0.8745221079     -0.0027397260      0.0006007681
pasttrop.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.3752, p-value = 6.065e-06
## alternative hypothesis: greater
## sample estimates:
## Moran I statistic       Expectation          Variance 
##      0.1014304761     -0.0027397260      0.0005668741
urbantrop.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

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.   : 46009   Min.   :0   Min.   :0  
##  1st Qu.:1273   1st Qu.:1153.1   1st Qu.:173356   1st Qu.:0   1st Qu.:0  
##  Median :2719   Median :2174.8   Median :197598   Median :0   Median :0  
##  Mean   :2594   Mean   :2171.2   Mean   :181253   Mean   :0   Mean   :0  
##  3rd Qu.:3473   3rd Qu.:3183.0   3rd Qu.:215900   3rd Qu.:0   3rd Qu.:0  
##  Max.   :4905   Max.   :4183.9   Max.   :244328   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  30454.95       0       0 var1
## 2  1344  139.09961  39258.65       0       0 var1
## 3  1429  236.17040  54754.92       0       0 var1
## 4  2217  350.74833  64935.16       0       0 var1
## 5  2023  481.77032  83665.96       0       0 var1
## 6  2348  613.02305  96891.25       0       0 var1
## 7  7663  908.97277 141286.26       0       0 var1
## 8  5711 1356.99636 193347.15       0       0 var1
## 9  6830 1927.71159 202587.99       0       0 var1
## 10 4856 2587.57793 209040.19       0       0 var1
## 11 2792 3288.06517 234097.31       0       0 var1
## 12 1512 3942.12332 205591.89       0       0 var1
## 
## Fitted variogram model:
##   model     psill    range
## 1   Nug  18811.45    0.000
## 2   Sph 192689.13 2031.883
## Sums of squares betw. var. model and sample var.[1] 1336337
plot(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.01524141       0       0 var1
## 2  1344  139.09961 0.02400387       0       0 var1
## 3  1429  236.17040 0.04092103       0       0 var1
## 4  2217  350.74833 0.04255403       0       0 var1
## 5  2023  481.77032 0.04801990       0       0 var1
## 6  2348  613.02305 0.05594666       0       0 var1
## 7  7663  908.97277 0.06027060       0       0 var1
## 8  5711 1356.99636 0.05752174       0       0 var1
## 9  6830 1927.71159 0.04867014       0       0 var1
## 10 4856 2587.57793 0.04718487       0       0 var1
## 11 2792 3288.06517 0.04994567       0       0 var1
## 12 1512 3942.12332 0.06248239       0       0 var1
## 
## Fitted variogram model:
##   model        psill    range kappa
## 1   Nug 6.199426e-05   0.0000   0.0
## 2   Ste 5.601977e-02 275.2665   0.8
## Sums of squares betw. var. model and sample var.[1] 1.102736e-06
plot(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.01694860       0       0 var1
## 2  1344  139.09961 0.02497957       0       0 var1
## 3  1429  236.17040 0.03906286       0       0 var1
## 4  2217  350.74833 0.03851386       0       0 var1
## 5  2023  481.77032 0.04851853       0       0 var1
## 6  2348  613.02305 0.06286870       0       0 var1
## 7  7663  908.97277 0.10022954       0       0 var1
## 8  5711 1356.99636 0.09487867       0       0 var1
## 9  6830 1927.71159 0.07182112       0       0 var1
## 10 4856 2587.57793 0.07856313       0       0 var1
## 11 2792 3288.06517 0.08155669       0       0 var1
## 12 1512 3942.12332 0.07681496       0       0 var1
## 
## Fitted variogram model:
##   model       psill    range
## 1   Nug 0.009422413    0.000
## 2   Sph 0.083305471 1190.746
## Sums of squares betw. var. model and sample var.[1] 5.308069e-06
plot(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.0004516770       0       0 var1
## 2  1344  139.09961 0.0012530632       0       0 var1
## 3  1429  236.17040 0.0017116872       0       0 var1
## 4  2217  350.74833 0.0009065203       0       0 var1
## 5  2023  481.77032 0.0007547635       0       0 var1
## 6  2348  613.02305 0.0005986767       0       0 var1
## 7  7663  908.97277 0.0005005286       0       0 var1
## 8  5711 1356.99636 0.0005852467       0       0 var1
## 9  6830 1927.71159 0.0005455258       0       0 var1
## 10 4856 2587.57793 0.0003783851       0       0 var1
## 11 2792 3288.06517 0.0003839431       0       0 var1
## 12 1512 3942.12332 0.0004290341       0       0 var1
## 
## Fitted variogram model:
##   model       psill    range
## 1   Nug 0.000000000   0.0000
## 2   Sph 0.001177381 217.7267
## Sums of squares betw. var. model and sample var.[1] 2.614558e-08
plot(med.urbnVar)

#*Tropical*
tropNVar <- variogram(N~1, trop)
summary(tropNVar)
##        np            dist            gamma           dir.hor     dir.ver 
##  Min.   : 496   Min.   : 200.8   Min.   : 48697   Min.   :0   Min.   :0  
##  1st Qu.:1112   1st Qu.:1271.4   1st Qu.:157336   1st Qu.:0   1st Qu.:0  
##  Median :1779   Median :2413.3   Median :328166   Median :0   Median :0  
##  Mean   :1760   Mean   :2401.5   Mean   :346824   Mean   :0   Mean   :0  
##  3rd Qu.:2550   3rd Qu.:3500.9   3rd Qu.:484849   3rd Qu.:0   3rd Qu.:0  
##  Max.   :3716   Max.   :4645.5   Max.   :732320   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  40120.19       0       0 var1
## 2   855  172.84714  40888.89       0       0 var1
## 3  1006  260.18125  57701.40       0       0 var1
## 4  1876  377.84163  61090.80       0       0 var1
## 5  1617  523.98659  87784.43       0       0 var1
## 6  1427  675.31568 100139.23       0       0 var1
## 7  3355  988.84731 113771.70       0       0 var1
## 8  2026 1504.84064 218557.98       0       0 var1
## 9  3558 2163.26132 317372.61       0       0 var1
## 10 6152 2904.51559 462516.18       0       0 var1
## 11 2539 3548.52640 687070.38       0       0 var1
## 12 2121 4600.60957 422589.62       0       0 var1
## 
## Fitted variogram model:
##   model     psill    range kappa
## 1   Nug  40236.74    0.000   0.0
## 2   Ste 739997.96 3846.651   1.5
## Sums of squares betw. var. model and sample var.[1] 14650349
plot(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.005514824       0       0 var1
## 2   855  172.84714 0.009212895       0       0 var1
## 3  1006  260.18125 0.013948206       0       0 var1
## 4  1876  377.84163 0.019254794       0       0 var1
## 5  1617  523.98659 0.024448382       0       0 var1
## 6  1427  675.31568 0.028569645       0       0 var1
## 7  3355  988.84731 0.032753031       0       0 var1
## 8  2026 1504.84064 0.033862958       0       0 var1
## 9  3558 2163.26132 0.045132409       0       0 var1
## 10 6152 2904.51559 0.044958775       0       0 var1
## 11 2539 3548.52640 0.050016620       0       0 var1
## 12 2121 4600.60957 0.068065910       0       0 var1
## 
## Fitted variogram model:
##   model      psill    range kappa
## 1   Nug 0.00000000    0.000   0.0
## 2   Ste 0.04629478 1041.202   0.5
## Sums of squares betw. var. model and sample var.[1] 1.163732e-07
plot(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.002160974       0       0 var1
## 2   855  172.84714 0.003355720       0       0 var1
## 3  1006  260.18125 0.003675748       0       0 var1
## 4  1876  377.84163 0.003991965       0       0 var1
## 5  1617  523.98659 0.005165840       0       0 var1
## 6  1427  675.31568 0.004535222       0       0 var1
## 7  3355  988.84731 0.003234505       0       0 var1
## 8  2026 1504.84064 0.007414587       0       0 var1
## 9  3558 2163.26132 0.011671621       0       0 var1
## 10 6152 2904.51559 0.005784703       0       0 var1
## 11 2539 3548.52640 0.009073648       0       0 var1
## 12 2121 4600.60957 0.072476617       0       0 var1
## 
## Fitted variogram model:
##   model       psill    range kappa
## 1   Nug 0.002903048     0.00   0.0
## 2   Ste 4.493031686 97590.07   1.6
## Sums of squares betw. var. model and sample var.[1] 3.49028e-07
plot(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.398401e-05       0       0 var1
## 2   855  172.84714 4.384211e-05       0       0 var1
## 3  1006  260.18125 5.630616e-05       0       0 var1
## 4  1876  377.84163 5.925986e-05       0       0 var1
## 5  1617  523.98659 5.983952e-05       0       0 var1
## 6  1427  675.31568 6.445760e-05       0       0 var1
## 7  3355  988.84731 9.335782e-05       0       0 var1
## 8  2026 1504.84064 1.031093e-04       0       0 var1
## 9  3558 2163.26132 4.831872e-05       0       0 var1
## 10 6152 2904.51559 7.313874e-05       0       0 var1
## 11 2539 3548.52640 1.126363e-04       0       0 var1
## 12 2121 4600.60957 8.495545e-05       0       0 var1
## 
## Fitted variogram model:
##   model        psill    range kappa
## 1   Nug 5.675377e-05   0.0000     0
## 2   Ste 2.445177e-05 795.6194    10
## Sums of squares betw. var. model and sample var.[1] 1.181187e-11
plot(trop.urbnVar)

#It seems like there is some different mechanism influencing the variograms in the tropical biome. Which makes sense, since the nutrient inputs in the system are so drastically different than mediterranean and grass land biomes.

#*Grassland*
grassNVar <- variogram(N~1, grass)
summary(grassNVar)
##        np              dist            gamma           dir.hor 
##  Min.   : 41122   Min.   : 212.5   Min.   : 18372   Min.   :0  
##  1st Qu.: 95860   1st Qu.:1276.3   1st Qu.: 93154   1st Qu.:0  
##  Median :125219   Median :2391.1   Median :115782   Median :0  
##  Mean   :113920   Mean   :2395.0   Mean   :103935   Mean   :0  
##  3rd Qu.:139084   3rd Qu.:3505.6   3rd Qu.:127569   3rd Qu.:0  
##  Max.   :146365   Max.   :4620.1   Max.   :156272   Max.   :0  
##     dir.ver     id    
##  Min.   :0   var1:15  
##  1st Qu.:0            
##  Median :0            
##  Mean   :0            
##  3rd Qu.:0            
##  Max.   :0
plot(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   9176.217       0       0 var1
## 2   12817  168.26411  14665.972       0       0 var1
## 3   17631  256.56728  20786.560       0       0 var1
## 4   40539  377.96649  30142.105       0       0 var1
## 5   48730  527.13390  42325.899       0       0 var1
## 6   56877  678.49745  53655.744       0       0 var1
## 7  211549 1010.22379  77357.390       0       0 var1
## 8  227826 1508.43676 105223.563       0       0 var1
## 9  337446 2130.74662 126418.368       0       0 var1
## 10 311651 2882.49709 124322.117       0       0 var1
## 11 260525 3623.17991 113096.582       0       0 var1
## 12 219770 4476.02031 145891.950       0       0 var1
## 
## Fitted variogram model:
##   model      psill    range
## 1   Nug   1153.479    0.000
## 2   Sph 124990.593 2346.776
## Sums of squares betw. var. model and sample var.[1] 8709228
plot(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.003149632       0       0 var1
## 2   12817  168.26411 0.005656960       0       0 var1
## 3   17631  256.56728 0.008273131       0       0 var1
## 4   40539  377.96649 0.010843567       0       0 var1
## 5   48730  527.13390 0.012810510       0       0 var1
## 6   56877  678.49745 0.013180889       0       0 var1
## 7  211549 1010.22379 0.014548992       0       0 var1
## 8  227826 1508.43676 0.016847260       0       0 var1
## 9  337446 2130.74662 0.017754809       0       0 var1
## 10 311651 2882.49709 0.017055363       0       0 var1
## 11 260525 3623.17991 0.019932955       0       0 var1
## 12 219770 4476.02031 0.020186006       0       0 var1
## 
## Fitted variogram model:
##   model      psill    range
## 1   Nug 0.00000000   0.0000
## 2   Exp 0.01722486 416.9518
## Sums of squares betw. var. model and sample var.[1] 8.3252e-07
plot(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.01087973       0       0 var1
## 2   12817  168.26411 0.02156178       0       0 var1
## 3   17631  256.56728 0.03153914       0       0 var1
## 4   40539  377.96649 0.04003507       0       0 var1
## 5   48730  527.13390 0.04485106       0       0 var1
## 6   56877  678.49745 0.05104508       0       0 var1
## 7  211549 1010.22379 0.05925443       0       0 var1
## 8  227826 1508.43676 0.06195645       0       0 var1
## 9  337446 2130.74662 0.07155589       0       0 var1
## 10 311651 2882.49709 0.06865125       0       0 var1
## 11 260525 3623.17991 0.05993092       0       0 var1
## 12 219770 4476.02031 0.05713202       0       0 var1
## 
## Fitted variogram model:
##   model        psill    range kappa
## 1   Nug 0.0003136473   0.0000   0.0
## 2   Ste 0.0644493800 546.6906   0.6
## Sums of squares betw. var. model and sample var.[1] 9.92466e-06
plot(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.840911e-05       0       0 var1
## 2   12817  168.26411 1.413213e-05       0       0 var1
## 3   17631  256.56728 1.240806e-05       0       0 var1
## 4   40539  377.96649 1.015735e-05       0       0 var1
## 5   48730  527.13390 1.038306e-05       0       0 var1
## 6   56877  678.49745 9.699720e-06       0       0 var1
## 7  211549 1010.22379 1.004936e-05       0       0 var1
## 8  227826 1508.43676 1.121150e-05       0       0 var1
## 9  337446 2130.74662 1.035195e-05       0       0 var1
## 10 311651 2882.49709 9.080749e-06       0       0 var1
## 11 260525 3623.17991 1.964251e-05       0       0 var1
## 12 219770 4476.02031 1.416946e-05       0       0 var1
## 
## Fitted variogram model:
##   model        psill    range kappa
## 1   Nug 9.409422e-06    0.000  0.00
## 2   Ste 5.993546e-06 2354.979  0.05
## Sums of squares betw. var. model and sample var.[1] 3.480182e-11
plot(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.2       -223.4        490.7       -140.7
med.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.11 -335.82  -88.35  210.12 2230.89 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept)  1238.24      83.29  14.867  < 2e-16 ***
## ppast        -223.38     128.86  -1.734 0.083762 .  
## pcult         490.65     143.00   3.431 0.000663 ***
## purbn        -140.69     474.98  -0.296 0.767235    
## ---
## 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.08302,    Adjusted R-squared:  0.07626 
## F-statistic: 12.28 on 3 and 407 DF,  p-value: 1.044e-07
med$lmResidsmed<-residuals(med.lm)
bubble(med, zcol='lmResidsmed') #Autocorrelated residuals! The linear models 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.818, p-value < 2.2e-16
## alternative hypothesis: greater
## sample estimates:
## Moran I statistic       Expectation          Variance 
##       0.803425555      -0.002439024       0.000535697
med.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.412 6261.456 -3115.706
## 
## Coefficients:
##                 Value Std.Error   t-value p-value
## (Intercept) 1238.2390   83.2889 14.866789  0.0000
## ppast       -223.3779  128.8584 -1.733514  0.0838
## pcult        490.6510  143.0015  3.431090  0.0007
## purbn       -140.6861  474.9839 -0.296191  0.7672
## 
##  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.6187423 -0.6669121 -0.1754527  0.4172844  4.4303611 
## 
## Residual standard error: 503.5458 
## Degrees of freedom: 411 total; 407 residual
cs1.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.412 6261.456 -3115.706
## 
## Coefficients:
##                 Value Std.Error   t-value p-value
## (Intercept) 1238.2390   83.2889 14.866789  0.0000
## ppast       -223.3779  128.8584 -1.733514  0.0838
## pcult        490.6510  143.0015  3.431090  0.0007
## purbn       -140.6861  474.9839 -0.296191  0.7672
## 
##  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.6187423 -0.6669121 -0.1754527  0.4172844  4.4303611 
## 
## Residual standard error: 503.5458 
## Degrees of freedom: 411 total; 407 residual
summary(gls2.med) #AIC 5549
## Generalized least squares fit by REML
##   Model: N ~ ppast + pcult + purbn 
##   Data: med 
##       AIC      BIC   logLik
##   5549.04 5577.102 -2767.52
## 
## Correlation Structure: Spherical spatial correlation
##  Formula: ~Y + X 
##  Parameter estimate(s):
##        range       nugget 
## 1.368596e+02 1.470961e-03 
## 
## Coefficients:
##                 Value Std.Error   t-value p-value
## (Intercept) 1673.7083  768.7691  2.177127  0.0300
## ppast       -427.4284   95.2171 -4.488986  0.0000
## pcult       -431.2618   97.4214 -4.426765  0.0000
## purbn       -274.9230  210.5842 -1.305526  0.1925
## 
##  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.64214233 -0.33994820 -0.08396528  0.11346874  1.27855416 
## 
## Residual standard error: 1624.231 
## Degrees of freedom: 411 total; 407 residual
med$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.499, p-value = 0.06694
## alternative hypothesis: greater
## sample estimates:
## Moran I statistic       Expectation          Variance 
##      0.0323149767     -0.0024390244      0.0005375577
residsI.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  
##      1845.4        117.8      -1055.4       8961.8
trop.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 
## -1291.24  -418.28   -43.71   304.52  2051.36 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept)  1845.39      75.09  24.576  < 2e-16 ***
## ppast         117.84     169.23   0.696   0.4867    
## pcult       -1055.35     143.46  -7.357 1.27e-12 ***
## purbn        8961.80    2988.99   2.998   0.0029 ** 
## ---
## 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.8 on 3 and 362 DF,  p-value: < 2.2e-16
trop$lmResidstrop<-residuals(trop.lm)
bubble(trop, zcol='lmResidstrop') #Autocorrelated residuals! The linear model 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.835, p-value < 2.2e-16
## alternative hypothesis: greater
## sample estimates:
## Moran I statistic       Expectation          Variance 
##       0.629508160      -0.002739726       0.000598925
trop.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.972 5570.43 -2770.486
## 
## Coefficients:
##                 Value Std.Error   t-value p-value
## (Intercept)  1845.393   75.0894 24.575934  0.0000
## ppast         117.842  169.2291  0.696349  0.4867
## pcult       -1055.355  143.4564 -7.356621  0.0000
## purbn        8961.797 2988.9925  2.998267  0.0029
## 
##  Correlation: 
##       (Intr) ppast  pcult 
## ppast -0.809              
## pcult -0.850  0.640       
## purbn -0.046 -0.057 -0.196
## 
## Standardized residuals:
##         Min          Q1         Med          Q3         Max 
## -2.55881277 -0.82889269 -0.08661567  0.60346634  4.06511716 
## 
## Residual standard error: 504.6239 
## Degrees of freedom: 366 total; 362 residual
cs1.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.972 5570.43 -2770.486
## 
## Coefficients:
##                 Value Std.Error   t-value p-value
## (Intercept)  1845.393   75.0894 24.575934  0.0000
## ppast         117.842  169.2291  0.696349  0.4867
## pcult       -1055.355  143.4564 -7.356621  0.0000
## purbn        8961.797 2988.9925  2.998267  0.0029
## 
##  Correlation: 
##       (Intr) ppast  pcult 
## ppast -0.809              
## pcult -0.850  0.640       
## purbn -0.046 -0.057 -0.196
## 
## Standardized residuals:
##         Min          Q1         Med          Q3         Max 
## -2.55881277 -0.82889269 -0.08661567  0.60346634  4.06511716 
## 
## Residual standard error: 504.6239 
## Degrees of freedom: 366 total; 362 residual
summary(gls2.trop) #AIC 5125
## Generalized least squares fit by REML
##   Model: N ~ ppast + pcult + purbn 
##   Data: trop 
##        AIC      BIC    logLik
##   5125.775 5153.016 -2555.887
## 
## Correlation Structure: Spherical spatial correlation
##  Formula: ~Y + X 
##  Parameter estimate(s):
##        range       nugget 
## 1.673412e+02 7.734042e-03 
## 
## Coefficients:
##                Value Std.Error    t-value p-value
## (Intercept) 1323.603 1011.3421  1.3087585  0.1914
## ppast        -37.096  216.5424 -0.1713111  0.8641
## pcult       -342.579  152.0266 -2.2534156  0.0248
## purbn       3673.761 1510.4997  2.4321498  0.0155
## 
##  Correlation: 
##       (Intr) ppast  pcult 
## ppast -0.095              
## pcult -0.031  0.070       
## purbn -0.006 -0.058 -0.103
## 
## Standardized residuals:
##         Min          Q1         Med          Q3         Max 
## -0.43597313 -0.02316476  0.18114759  0.36930016  1.44230388 
## 
## Residual standard error: 1761.63 
## Degrees of freedom: 366 total; 362 residual
trop$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.26526, p-value = 0.6046
## alternative hypothesis: greater
## sample estimates:
## Moran I statistic       Expectation          Variance 
##     -0.0091928791     -0.0027397260      0.0005918221
residsI.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.8       -352.8        272.5       3309.8
grass.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.2 -265.2  -40.3  204.0 2238.7 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept)  1201.75      16.13  74.486  < 2e-16 ***
## ppast        -352.79      31.44 -11.220  < 2e-16 ***
## pcult         272.55      61.68   4.419 1.03e-05 ***
## purbn        3309.83    1014.33   3.263  0.00112 ** 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 381.1 on 2641 degrees of freedom
## Multiple R-squared:  0.05966,    Adjusted R-squared:  0.05859 
## F-statistic: 55.85 on 3 and 2641 DF,  p-value: < 2.2e-16
grass$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.509, p-value < 2.2e-16
## alternative hypothesis: greater
## sample estimates:
## Moran I statistic       Expectation          Variance 
##      8.770397e-01     -3.782148e-04      8.995861e-05
grass.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. AIC 38,910
## Generalized least squares fit by REML
##   Model: N ~ ppast + pcult + purbn 
##   Data: grass 
##        AIC      BIC    logLik
##   38910.26 38939.65 -19450.13
## 
## Coefficients:
##                Value Std.Error   t-value p-value
## (Intercept) 1201.752   16.1340  74.48556  0.0000
## ppast       -352.791   31.4426 -11.22015  0.0000
## pcult        272.545   61.6800   4.41870  0.0000
## purbn       3309.825 1014.3333   3.26305  0.0011
## 
##  Correlation: 
##       (Intr) ppast  pcult 
## ppast -0.822              
## pcult -0.335  0.011       
## purbn -0.040  0.030 -0.198
## 
## Standardized residuals:
##        Min         Q1        Med         Q3        Max 
## -2.5486788 -0.6960177 -0.1057454  0.5353624  5.8748106 
## 
## Residual standard error: 381.0696 
## Degrees of freedom: 2645 total; 2641 residual
cs1.grass <-corSpher(c(2500,1000/1.2e5),form=~Y+X,nugget=TRUE)
#gls2.grass <- update(gls1.grass,correlation=cs1.grass)
#summary(gls2.grass) 
#grass$glsResids.grass<-residuals(gls2.grass,type="normalized")
#bubble(grass,zcol='glsResids.grass')
#moran.test(grass$glsResids.grass,nb2listw(grass.gls)) 
#residsI.grass <-spline.correlog(x=coordinates(grass)[,1], y=coordinates(grass)[,2],z=grass$glsResids.grass, resamp=1, quiet=TRUE)
#plot(residsI.grass)
Pastures, crops, and urbanization all seem to have an effect on native species in grasslands. This may have to do with the fact that the authors did not distinguish between which species were invasive, only that native species loss opened niches for new plants to colonize. This availability of niche space seems plausible in grasslands where a lot of crops are being grown.