Final

Data

These data are from mpg ranch located in SW Montana. The ranch was previously used for agriculture and livestock, but is now dedicated to conservation. Investigating invasive species, such as Potentilla recta, is part of their conservation goals. The data I have are a species matrix of ~360 species in ~ 430 plots and an environmental variable set extracted from a digital elevation model.

Objectives

1.Is there spatial structure to the occurrence of Potentilla recta “POTREC”

2.Can P. recta be predicted by its relatives “POTGLA” & “POTGRA”?

3.Can P. recta be predicted by variables from a DEM?

4.How do the DEM variables relate to the species composition of MPG?

Make a map of survey sites for next week’s survey.

Plan

View the spatial structure of the potentillas by using variograms, Moran’s I, and spline correlogs.
View linear models of P. recta as a function of its relatives.
View linear models of P. recta as a function of the environmental variables.
Ordinate the species matrix.
Plot co-occurrence sites and the sites from a 2013 survey.

Results

In the 2012 survey P. recta is not incredibly abundant.The variograms, Moran’s I, and spline correlogs suggest that there isn’t spatial structure to the distribution of P.recta

library(gstat)
library(sp)
library(spatstat)

## Loading required package: nlme

## 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

setwd("E:/mpg")
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: C:/Users/Faythe/Documents/R/win-library/3.3/rgdal/gdal
##  Loaded PROJ.4 runtime: Rel. 4.9.2, 08 September 2015, [PJ_VERSION: 492]
##  Path to PROJ.4 shared files: C:/Users/Faythe/Documents/R/win-library/3.3/rgdal/proj
##  Linking to sp version: 1.2-3

library(raster)

## 
## Attaching package: 'raster'

## The following objects are masked from 'package:spatstat':
## 
##     area, rotate, shift

## The following object is masked from 'package:nlme':
## 
##     getData

library(vegan)

## Loading required package: permute

## Loading required package: lattice

## 
## Attaching package: 'lattice'

## The following object is masked from 'package:spatstat':
## 
##     panel.histogram

## This is vegan 2.3-5

load("matrices/threeMatrices.Rdata")
spp.pot <- spp.df[,277:279]
spp.pot$x <- xy.df$x
spp.pot$y <- xy.df$y
coordinates(spp.pot) <- ~x+y
spp.pot$elev<-env.df$elev
spp.pot$twi<-env.df$twi
spp.pot$slope<-env.df$slope
spp.pot$aspect<-env.df$aspect
spp.pot$southness<-env.df$southness
spp.pot$westness<-env.df$westness
spp.pot$tri<-env.df$tri
spp.pot$tpi<-env.df$tpi
#data plots
plot(spp.pot)

bubble(spp.pot,"POTREC")

bubble(spp.pot,"POTGLA")

bubble(spp.pot,"POTGRA")

spplot(spp.pot,"POTREC",colorkey=T)

p.recta.Var <- variogram(POTREC~1, spp.pot, cloud = F)
plot(p.recta.Var,pch=20,cex=1.5,col="black",
     ylab=expression("Semivariance ("*gamma*")"),
     xlab="Distance (m)",main = "Potentilla recta abundance")

library(ncf)

## 
## Attaching package: 'ncf'

## The following object is masked from 'package:vegan':
## 
##     mantel.correlog

library(spdep)

## Loading required package: Matrix

w <- knn2nb(knearneigh(spp.pot,k=8))
p.rec.I=moran.test(spp.pot$POTREC,nb2listw(w))
p.rec.I

## 
##  Moran I test under randomisation
## 
## data:  spp.pot$POTREC  
## weights: nb2listw(w)  
## 
## Moran I statistic standard deviate = 0.68634, p-value = 0.2462
## alternative hypothesis: greater
## sample estimates:
## Moran I statistic       Expectation          Variance 
##      0.0123726928     -0.0023419204      0.0004596387

recta <- spline.correlog(x=coordinates(spp.pot)[,1], y=coordinates(spp.pot)[,2],
                         z=spp.pot$POTREC, resamp=100, quiet=TRUE)
plot(recta)

The other Potentilla spp. did not seem to be a good predictors of Potentilla recta.

lm.rectaxgla=lm(POTREC~POTGLA, spp.pot)
summary(lm.rectaxgla)

## 
## Call:
## lm(formula = POTREC ~ POTGLA, data = spp.pot)
## 
## Residuals:
##    Min     1Q Median     3Q    Max 
## -2.174 -0.720 -0.720 -0.720 38.780 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept)   0.7202     0.1825   3.946 9.29e-05 ***
## POTGLA        0.3231     0.5676   0.569     0.57    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 3.725 on 426 degrees of freedom
## Multiple R-squared:  0.00076,    Adjusted R-squared:  -0.001586 
## F-statistic: 0.324 on 1 and 426 DF,  p-value: 0.5695

spp.pot$lmrecXglaResids <- residuals(lm.rectaxgla)
bubble(spp.pot, zcol = 'lmrecXglaResids')

w <- knn2nb(knearneigh(spp.pot, k=8))
recXgla.resids=moran.test(spp.pot$lmrecXglaResids, nb2listw(w))
recXgla.resids

## 
##  Moran I test under randomisation
## 
## data:  spp.pot$lmrecXglaResids  
## weights: nb2listw(w)  
## 
## Moran I statistic standard deviate = 0.72045, p-value = 0.2356
## alternative hypothesis: greater
## sample estimates:
## Moran I statistic       Expectation          Variance 
##      0.0130995659     -0.0023419204      0.0004593841

lm.rectaxgra=lm(POTREC~POTGRA, spp.pot)
summary(lm.rectaxgra)

## 
## Call:
## lm(formula = POTREC ~ POTGRA, data = spp.pot)
## 
## Residuals:
##    Min     1Q Median     3Q    Max 
## -2.344 -0.723 -0.723 -0.723 38.777 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept)   0.7225     0.1813   3.986 7.91e-05 ***
## POTGRA        0.4633     0.6763   0.685    0.494    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 3.724 on 426 degrees of freedom
## Multiple R-squared:  0.0011, Adjusted R-squared:  -0.001245 
## F-statistic: 0.4692 on 1 and 426 DF,  p-value: 0.4937

lm.both=lm(POTREC~POTGRA+POTGLA, spp.pot)
summary(lm.both)

## 
## Call:
## lm(formula = POTREC ~ POTGRA + POTGLA, data = spp.pot)
## 
## Residuals:
##    Min     1Q Median     3Q    Max 
## -2.549 -0.719 -0.719 -0.719 38.781 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept)   0.7192     0.1827   3.937 9.65e-05 ***
## POTGRA        0.3703     0.8995   0.412    0.681    
## POTGLA        0.1185     0.7548   0.157    0.875    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 3.729 on 425 degrees of freedom
## Multiple R-squared:  0.001158,   Adjusted R-squared:  -0.003542 
## F-statistic: 0.2464 on 2 and 425 DF,  p-value: 0.7817

The DEM variables did not make good predictors of

lm.rectaxgra=lm(POTREC~POTGRA, spp.pot)
summary(lm.rectaxgra)

## 
## Call:
## lm(formula = POTREC ~ POTGRA, data = spp.pot)
## 
## Residuals:
##    Min     1Q Median     3Q    Max 
## -2.344 -0.723 -0.723 -0.723 38.777 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept)   0.7225     0.1813   3.986 7.91e-05 ***
## POTGRA        0.4633     0.6763   0.685    0.494    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 3.724 on 426 degrees of freedom
## Multiple R-squared:  0.0011, Adjusted R-squared:  -0.001245 
## F-statistic: 0.4692 on 1 and 426 DF,  p-value: 0.4937

lm.rectaxelev=lm(POTREC~elev,spp.pot)
summary(lm.rectaxelev)

## 
## Call:
## lm(formula = POTREC ~ elev, data = spp.pot)
## 
## Residuals:
##    Min     1Q Median     3Q    Max 
## -0.998 -0.850 -0.701 -0.497 38.641 
## 
## Coefficients:
##               Estimate Std. Error t value Pr(>|t|)
## (Intercept)  1.7668244  1.1411402   1.548    0.122
## elev        -0.0007972  0.0008724  -0.914    0.361
## 
## Residual standard error: 3.723 on 426 degrees of freedom
## Multiple R-squared:  0.001956,   Adjusted R-squared:  -0.0003867 
## F-statistic: 0.8349 on 1 and 426 DF,  p-value: 0.3614

lm.rectaxtwi=lm(POTREC~twi,spp.pot)
summary(lm.rectaxtwi)

## 
## Call:
## lm(formula = POTREC ~ twi, data = spp.pot)
## 
## Residuals:
##    Min     1Q Median     3Q    Max 
## -0.971 -0.824 -0.746 -0.580 38.776 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)  
## (Intercept)  1.37235    0.80205   1.711   0.0878 .
## twi         -0.09343    0.11496  -0.813   0.4168  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 3.723 on 426 degrees of freedom
## Multiple R-squared:  0.001548,   Adjusted R-squared:  -0.0007958 
## F-statistic: 0.6605 on 1 and 426 DF,  p-value: 0.4168

lm.rectaxslope=lm(POTREC~slope,spp.pot)
summary(lm.rectaxslope)

## 
## Call:
## lm(formula = POTREC ~ slope, data = spp.pot)
## 
## Residuals:
##    Min     1Q Median     3Q    Max 
## -1.526 -0.902 -0.640 -0.446 38.650 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)
## (Intercept)  0.28299    0.36202   0.782    0.435
## slope        0.02822    0.01953   1.445    0.149
## 
## Residual standard error: 3.717 on 426 degrees of freedom
## Multiple R-squared:  0.004878,   Adjusted R-squared:  0.002542 
## F-statistic: 2.088 on 1 and 426 DF,  p-value: 0.1492

lm.rectaxaspect=lm(POTREC~aspect,spp.pot)
summary(lm.rectaxaspect)

## 
## Call:
## lm(formula = POTREC ~ aspect, data = spp.pot)
## 
## Residuals:
##    Min     1Q Median     3Q    Max 
## -0.762 -0.741 -0.733 -0.727 38.775 
## 
## Coefficients:
##               Estimate Std. Error t value Pr(>|t|)
## (Intercept)  0.7619550  0.5210514   1.462    0.144
## aspect      -0.0001087  0.0021426  -0.051    0.960
## 
## Residual standard error: 3.726 on 426 degrees of freedom
## Multiple R-squared:  6.042e-06,  Adjusted R-squared:  -0.002341 
## F-statistic: 0.002574 on 1 and 426 DF,  p-value: 0.9596

lm.rectaxS=lm(POTREC~southness,spp.pot)
summary(lm.rectaxS)

## 
## Call:
## lm(formula = POTREC ~ southness, data = spp.pot)
## 
## Residuals:
##    Min     1Q Median     3Q    Max 
## -0.782 -0.763 -0.730 -0.691 38.813 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept)   0.7336     0.1813   4.047 6.16e-05 ***
## southness     0.0483     0.2743   0.176     0.86    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 3.726 on 426 degrees of freedom
## Multiple R-squared:  7.276e-05,  Adjusted R-squared:  -0.002274 
## F-statistic: 0.031 on 1 and 426 DF,  p-value: 0.8603

The ordination of the plant species matrix is pretty neat. The DEM variable that fits the whole community structure was elevation. This makes sense because elevation corresponds to other environmental variables such as temperature. POTGLA and POTGRA are near each other, and POTREC is on the opposite side. To me this means that the two natives (POTGLA and POTGRA) are more common in the same type of community, and that the community where POTREC is found is different. This is interesting because although we did not find any direct relationship between POTREC and its congeners they appear in different communities. It will be awesome to get the soil microbial data and see if we get the same or similar findings.

library(rgdal) 
library(raster) 
library(vegan)
load("matrices/threeMatrices.Rdata")
pairs(env.df[,-4])

load("matrices/ord.Rdata")
(fit <- envfit(ord, env.df[,-4], perm = 999))

## 
## ***VECTORS
## 
##              NMDS1    NMDS2     r2 Pr(>r)   
## elev       0.87399  0.48595 0.0259  0.005 **
## twi       -0.76017 -0.64973 0.0020  0.666   
## slope      0.80399  0.59464 0.0044  0.381   
## southness -0.83329 -0.55283 0.0024  0.587   
## westness   0.25956  0.96573 0.0059  0.295   
## tri        0.81512  0.57929 0.0042  0.399   
## tpi       -0.72006  0.69391 0.0022  0.643   
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## Permutation: free
## Number of permutations: 999

plot(ord,type="t")
plot(fit)

fit

## 
## ***VECTORS
## 
##              NMDS1    NMDS2     r2 Pr(>r)   
## elev       0.87399  0.48595 0.0259  0.005 **
## twi       -0.76017 -0.64973 0.0020  0.666   
## slope      0.80399  0.59464 0.0044  0.381   
## southness -0.83329 -0.55283 0.0024  0.587   
## westness   0.25956  0.96573 0.0059  0.295   
## tri        0.81512  0.57929 0.0042  0.399   
## tpi       -0.72006  0.69391 0.0022  0.643   
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## Permutation: free
## Number of permutations: 999

spp=ord$species
spp.1=spp[277:279,]
plot(spp.1)

Below is a map of locations, form the 2012 survey, where P.recta was found with a congener as well as the sites from the 2013 survey.

library(ggmap)

## Loading required package: ggplot2

library(ggplot2)
mpg<- c(lon = -114.0048532,lat = 46.6856164)
myMPG <- get_map(location=mpg,source="stamen", maptype="terrain", crop=FALSE, zoom=13)

## Map from URL : http://maps.googleapis.com/maps/api/staticmap?center=46.685616,-114.004853&zoom=13&size=640x640&scale=2&maptype=terrain&sensor=false

## Map from URL : http://tile.stamen.com/terrain/13/1500/2890.jpg

## Map from URL : http://tile.stamen.com/terrain/13/1501/2890.jpg

## Map from URL : http://tile.stamen.com/terrain/13/1502/2890.jpg

## Map from URL : http://tile.stamen.com/terrain/13/1503/2890.jpg

## Map from URL : http://tile.stamen.com/terrain/13/1500/2891.jpg

## Map from URL : http://tile.stamen.com/terrain/13/1501/2891.jpg

## Map from URL : http://tile.stamen.com/terrain/13/1502/2891.jpg

## Map from URL : http://tile.stamen.com/terrain/13/1503/2891.jpg

## Map from URL : http://tile.stamen.com/terrain/13/1500/2892.jpg

## Map from URL : http://tile.stamen.com/terrain/13/1501/2892.jpg

## Map from URL : http://tile.stamen.com/terrain/13/1502/2892.jpg

## Map from URL : http://tile.stamen.com/terrain/13/1503/2892.jpg

## Map from URL : http://tile.stamen.com/terrain/13/1500/2893.jpg

## Map from URL : http://tile.stamen.com/terrain/13/1501/2893.jpg

## Map from URL : http://tile.stamen.com/terrain/13/1502/2893.jpg

## Map from URL : http://tile.stamen.com/terrain/13/1503/2893.jpg

ggmap(myMPG)

sites=read.csv("sites.csv",header = TRUE)
ggmap(myMPG)+
  geom_point(aes(x =long, y =lat, show_guide=T,colour=type), data = sites,
             alpha = .5, na.rm = T, size = 4) +
  scale_colour_brewer(sites$type,type="qual",palette="Dark2",direction = 1) +
  xlab("Longitude") +
  ylab("Latitdude") +
  ggtitle("Potentilla Map")

Conclusions

There is no spatial structure in the distribution of P.recta from the 2012 survey. The congeners and DEM environmental variables were not good predictors of P.recta. Elevation was the only DEM variable to fit the ordinated species matrix.The map at the end will be a helpful visual this week and next as we plan for the Potentilla survey. Interestingly, the two natives are clumped together and apart from the exotic. This analysis is not complete; MPG ranch was a ranch and used for agriculture and raising livestock for over 100 years. Future analysis would investigate the land use history effect of current plant populations on MPG.

spatial autocorrelation of the whole community

# Detrend the species data by regression on the site coordinates
spp.hel <- decostand(spp.df, "hellinger")
spp.bcd <- vegdist(spp.df,method="bray")


# model by space
spp.hel.resid <- resid(lm(as.matrix(spp.hel) ~ ., data=data.frame(xy.df)))
spp.hel.D <- dist(spp.hel.resid)
summary(spp.hel.D)

##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
##  0.1353  1.0100  1.1630  1.1250  1.2710  1.4620

The community is auto correlated up to a distance of just over 1000m?. Filled in black = alpha < 0.05.