Week 3 Home range and resource selection
MD
13/02/2023
In today’s practical we will look at two spatially oriented and closely aligned areas of species ecology: range and resource selection. We will explore some of the ideas covered in the lecture using a telemetry dataset of Canis lupus (grey wolf) from a 2006 monitoring campaign in the Alberta region of Canada. We will use these data to explore options for calculating and visualizing animal home ranges and look at method for computing estimates of resource selection. The results will provide the input to the Week 8 practical of this course where we will carry out a least cost path analysis.The home ranges estimated this week will provide the beginning and end points of a potential habitat corridor and the resource preferences calculated later in today’s practical will be used to compute a least cost path between these points.
Before we begin, copy the data from the “Canis” folder in the GEOG70922 dropbox link provided and paste it to a folder under your GEOG70922 directory called “Week3”. Open R studio and save a new script to this same folder - call it “Prac3.R”.
First, load the packages that we will use in today’s practical:
install.packages("terra") # for spatial (vector/raster) objects
install.packages("dismo") # for generating background data
install.packages("adehabitatHR") # for computing home range
install.packages("adehabitatHS") # for computing resource selection
install.packages("rgdal") # for projections - this can take a while to load so please be patient
install.packages("maptools") # for mapping reference data
install.packages("reshape2")
install.packages("rgeos")
install.packages("raster")Now, let’s load in and have a look at the telemetry data for today’s practical.
setwd("P:/Week3")Canislupus<-read.csv("Wolves.csv")View the first six rows of the dataset:
head(Canislupus)## event.id visible timestamp Lon Lat Taxon
## 1 11952917594 TRUE 02:52.0 -113.3277 55.44190 Canis lupus
## 2 11952917595 TRUE 02:53.0 -113.3163 55.44307 Canis lupus
## 3 11952917596 TRUE 03:08.0 -113.2985 55.44496 Canis lupus
## 4 11952917597 TRUE 00:54.0 -113.3148 55.43769 Canis lupus
## 5 11952917598 TRUE 01:52.0 -113.3262 55.42736 Canis lupus
## 6 11952917599 TRUE 00:54.0 -113.3264 55.42768 Canis lupus
## tag.local.identifier animal_ID study.name
## 1 3 3 Latham Alberta Wolves
## 2 3 3 Latham Alberta Wolves
## 3 3 3 Latham Alberta Wolves
## 4 3 3 Latham Alberta Wolves
## 5 3 3 Latham Alberta Wolves
## 6 3 3 Latham Alberta WolvesThe dataset contains information on the taxon, individual animal ID and coordinates of each location detected. The key difference between this dataset and the one from last week is the “event.id” field. Note that here each row stands for a location of an individual animal, the ID of which is contained in the “animal_ID” field. These locations will be the basis of our definition of resource use by these animals.
Now let’s create a SpatialPointsDataFrame from these data and plot them. We need to create a SpatialPointsDataFrame object rather than simply a SpatialPoints object (as we used last week) as the functions we will be using require access to the data table associated with these points, specifically the “animal_ID” field, so that we can perform analyses on individual animals.
library(rgdal)## Loading required package: sp## Please note that rgdal will be retired by the end of 2023,
## plan transition to sf/stars/terra functions using GDAL and PROJ
## at your earliest convenience.
##
## rgdal: version: 1.5-32, (SVN revision 1176)
## Geospatial Data Abstraction Library extensions to R successfully loaded
## Loaded GDAL runtime: GDAL 3.4.3, released 2022/04/22
## Path to GDAL shared files: C:/Users/mlibxmda/AppData/Local/R/win-library/4.2/rgdal/gdal
## GDAL binary built with GEOS: TRUE
## Loaded PROJ runtime: Rel. 7.2.1, January 1st, 2021, [PJ_VERSION: 721]
## Path to PROJ shared files: C:/Users/mlibxmda/AppData/Local/R/win-library/4.2/rgdal/proj
## PROJ CDN enabled: FALSE
## Linking to sp version:1.5-0
## To mute warnings of possible GDAL/OSR exportToProj4() degradation,
## use options("rgdal_show_exportToProj4_warnings"="none") before loading sp or rgdal.library(raster)
library(terra)## terra 1.6.7##
## Attaching package: 'terra'## The following object is masked from 'package:rgdal':
##
## projectCanis.latlong<-data.frame(x=Canislupus$Lon,y=Canislupus$Lat) #create a data frame containing coordinates
head(Canis.latlong) #inspect## x y
## 1 -113.3277 55.44190
## 2 -113.3163 55.44307
## 3 -113.2985 55.44496
## 4 -113.3148 55.43769
## 5 -113.3262 55.42736
## 6 -113.3264 55.42768crs.latlong<-crs("+init=epsg:4326") # object for the desired projection (WGS84)
Canis_attr<-Canislupus[,4:9] #object containing the data we want (columns 6 to 9 of the .csv file)
Canis<-vect(Canis.latlong,geom=c("x","y"))
Canis<-cbind(Canis,Canis_attr)
#inspect
head(Canis)## Lon Lat Taxon tag.local.identifier animal_ID
## 1 -113.3277 55.44190 Canis lupus 3 3
## 2 -113.3163 55.44307 Canis lupus 3 3
## 3 -113.2985 55.44496 Canis lupus 3 3
## 4 -113.3148 55.43769 Canis lupus 3 3
## 5 -113.3262 55.42736 Canis lupus 3 3
## 6 -113.3264 55.42768 Canis lupus 3 3
## study.name
## 1 Latham Alberta Wolves
## 2 Latham Alberta Wolves
## 3 Latham Alberta Wolves
## 4 Latham Alberta Wolves
## 5 Latham Alberta Wolves
## 6 Latham Alberta Wolvescrs(Canis)<-"epsg:4326"
plot(Canis)
Let’s visualize each animal by creating a vector of colours to assign to each
#check unique animal IDs
unique(Canis$animal_ID)## [1] 3 6 10 13 23 27 29 30 31 12#object for animal ID as a factor (categorical variable so we can us this to set colours later)
canisCol <- factor(Canis$animal_ID)
#set colours to number of levels (IDs) in the animalID field (now a categorical variable/factor)- "pch" means 'point character', 16 is the code for a round point symbol
plot(Canis, pch = 16,col=canisCol)
This doesn’t tell us much about the geographical context of the study area so let’s call in some data from R using the maptools library:
library(rgeos)## rgeos version: 0.5-9, (SVN revision 684)
## GEOS runtime version: 3.9.1-CAPI-1.14.2
## Please note that rgeos will be retired by the end of 2023,
## plan transition to sf functions using GEOS at your earliest convenience.
## GEOS using OverlayNG
## Linking to sp version: 1.5-0
## Polygon checking: TRUElibrary(maptools)## Checking rgeos availability: TRUE
## Please note that 'maptools' will be retired by the end of 2023,
## plan transition at your earliest convenience;
## some functionality will be moved to 'sp'.data(wrld_simpl)
plot(wrld_simpl,axes=T)
plot(Canis,col='red', add=T)
Okay, we’re in North America, let’s zoom in a little more:
plot(wrld_simpl, xlim=c(-110,-100), ylim=c(53,70), axes=TRUE)
plot(Canis,col='red', cex=0.75,add=T)
OK, now we are familar with the locations of these data, we’ll plot these against a raster of habitat types for this part of the world (the raster is already clipped to the extent of the points data to save space). Consult the metadata file on the “Canis” folder for details on this dataset. The habitat categories are as follows (have a look at the metadata file also):
11”Unclassified”, 21”Settlement”, 24”Roads”, 31”Water”, 41”Forest”, 42”Forest Wetland”, 45”Trees”, 46”Treed Wetland”, 51”Cropland”, 61”Grassland Managed” 62”Grassland Unmanaged”, 71”Wetland”, 73”Wetland Shrub”, 74”Wetland Herb”, 91”Other land”
CA<-rast("AlbertaCanis.tif")Now, our point data are in WGS84 but the habitat raster CRS is set to the UTM zone 22N for Canada (EPSG3979). We prefer the CRS of the raster object as this is a projected CRS with units in metres rather than degrees (so with this we can calculate the area for our animal home ranges). So, we need to reproject our points to fit on our habitat raster. We will do this in much the same way as we did for the Sciurus vulgaris data last week.
canisProj<-crs("+init=epsg:3979")
Canis<-project(Canis,canisProj)
crs(Canis)## [1] "PROJCRS[\"unknown\",\n BASEGEOGCRS[\"unknown\",\n DATUM[\"Unknown based on GRS80 ellipsoid\",\n ELLIPSOID[\"GRS 1980\",6378137,298.257222101,\n LENGTHUNIT[\"metre\",1],\n ID[\"EPSG\",7019]]],\n PRIMEM[\"Greenwich\",0,\n ANGLEUNIT[\"degree\",0.0174532925199433],\n ID[\"EPSG\",8901]]],\n CONVERSION[\"unknown\",\n METHOD[\"Lambert Conic Conformal (2SP)\",\n ID[\"EPSG\",9802]],\n PARAMETER[\"Latitude of false origin\",49,\n ANGLEUNIT[\"degree\",0.0174532925199433],\n ID[\"EPSG\",8821]],\n PARAMETER[\"Longitude of false origin\",-95,\n ANGLEUNIT[\"degree\",0.0174532925199433],\n ID[\"EPSG\",8822]],\n PARAMETER[\"Latitude of 1st standard parallel\",49,\n ANGLEUNIT[\"degree\",0.0174532925199433],\n ID[\"EPSG\",8823]],\n PARAMETER[\"Latitude of 2nd standard parallel\",77,\n ANGLEUNIT[\"degree\",0.0174532925199433],\n ID[\"EPSG\",8824]],\n PARAMETER[\"Easting at false origin\",0,\n LENGTHUNIT[\"metre\",1],\n ID[\"EPSG\",8826]],\n PARAMETER[\"Northing at false origin\",0,\n LENGTHUNIT[\"metre\",1],\n ID[\"EPSG\",8827]]],\n CS[Cartesian,2],\n AXIS[\"(E)\",east,\n ORDER[1],\n LENGTHUNIT[\"metre\",1,\n ID[\"EPSG\",9001]]],\n AXIS[\"(N)\",north,\n ORDER[2],\n LENGTHUNIT[\"metre\",1,\n ID[\"EPSG\",9001]]]]"plot(CA)
plot(Canis, pch = 16,col=canisCol,add=TRUE)
The first thing we will do in terms of analysis is to calculate home ranges for each animal (each unique value under the animal_ID field). For this we use the mcp() function from the AdehabitatHR library. ‘mcp’ stands for “minimum convex polygon” and is conceptually identical to the convex hull tool that some of you will have used in ArcGIS. To compute separate polygons for each animal we need to provide the field in the object that contains the animal IDs. Then we set the ‘percent’ argument to 95% to remove the 5 percent of the locations that are farthest away from the centroid of the polygon.
library(adehabitatHR)## Loading required package: deldir## deldir 1.0-6 Nickname: "Mendacious Cosmonaut"##
## The syntax of deldir() has had an important change.
## The arguments have been re-ordered (the first three
## are now "x, y, z") and some arguments have been
## eliminated. The handling of the z ("tags")
## argument has been improved.
##
## The "dummy points" facility has been removed.
## This facility was a historical artefact, was really
## of no use to anyone, and had hung around much too
## long. Since there are no longer any "dummy points",
## the structure of the value returned by deldir() has
## changed slightly. The arguments of plot.deldir()
## have been adjusted accordingly; e.g. the character
## string "wpoints" ("which points") has been
## replaced by the logical scalar "showpoints".
## The user should consult the help files.## Loading required package: ade4## Loading required package: adehabitatMA## Registered S3 methods overwritten by 'adehabitatMA':
## method from
## print.SpatialPixelsDataFrame sp
## print.SpatialPixels sp##
## Attaching package: 'adehabitatMA'## The following object is masked from 'package:terra':
##
## buffer## The following object is masked from 'package:raster':
##
## buffer## Loading required package: adehabitatLT## Loading required package: CircStats## Loading required package: MASS##
## Attaching package: 'MASS'## The following object is masked from 'package:terra':
##
## area## The following objects are masked from 'package:raster':
##
## area, select## Loading required package: boot#check the column number of animal IDs
head(Canis) #it's col 5## Lon Lat Taxon tag.local.identifier animal_ID
## 1 -113.3277 55.44190 Canis lupus 3 3
## 2 -113.3163 55.44307 Canis lupus 3 3
## 3 -113.2985 55.44496 Canis lupus 3 3
## 4 -113.3148 55.43769 Canis lupus 3 3
## 5 -113.3262 55.42736 Canis lupus 3 3
## 6 -113.3264 55.42768 Canis lupus 3 3
## study.name
## 1 Latham Alberta Wolves
## 2 Latham Alberta Wolves
## 3 Latham Alberta Wolves
## 4 Latham Alberta Wolves
## 5 Latham Alberta Wolves
## 6 Latham Alberta Wolves#delineate home range of each animal based on minimum bounding polygon. Note that the AdeHabitta library only accepts objects of class SpatialPoints from the SP library (this is a pre-cursor to the terra library). So here we convert the object to a "Spatial" (SP) object using as(object, "Spatial") from the raster library.
Canis<-as(Canis,"Spatial")
WCH<-mcp(Canis[,5],percent=95)
#add to map!
plot(CA)
plot(WCH,add=T)
Suppose we want to calculate the home range of the entire population (rather than indivduals). For this we need to use the first column of the dataset (this column only contains one unique value - the name of the species - and so will treat the data as one individual).
#minimum convex hull for all animals:
allCP<-mcp(Canis[,3],percent=99)
plot(CA)
plot(WCH,add=TRUE)
plot(allCP, add=TRUE)
Now let’s use one of the other options in the ‘adehabitatHR’ library to calcualte the size (area) of home ranges of each animal for different ‘percent’ values. Here we use the mcp.area() function and programme it to plot range sizes for increments of 5% between a min of 20 and max of 100. We set the ‘in’ units to metres (the units of the projection) and resulting area units to km2.
RangeCanis<-mcp.area(Canis[,5],percent = seq(20,100, by = 5),
unin = "m",
unout = "km2")
We can also use a kernel approach to map home ranges. A common method is to use a normal distribution to establish kernel weights (i.e. a Gaussian kernel) and a bivariate approach (i.e. x and y directions considered). Below we use such an approach to compute kernel density estimations (KDE) based on our telemetry data. The kernelUD (‘UD’ = Use Density) requires a column of animal IDs, a value for the bandwidth ‘h’ to use (we will set it to “href” to use the default bandwidth) and a kernel approach (here we use a bivariate normal kernel - “bivnorm”). We will initially set the ID column to column 5 of the Canis object - “animal_ID” (Canis[,5]) to create separate estimations for each animal and then to column 3 (‘Taxon’) to recreate for at the population level.
Caniskernel <- kernelUD(Canis[,5], h = "href", kern ="bivnorm")
kernelCanisAll<-kernelUD(Canis[,3], h = "href", kern ="bivnorm")
class(Caniskernel)## [1] "estUDm"image(Caniskernel)
class(kernelCanisAll)## [1] "estUDm"image(kernelCanisAll)
The kernalCanisAll is returned as an estUD (estimated UseDensity) object and, for Caniskernel (multiple animals), an estUDm file is returned. We can convert these to a more familiar raster format, first by converting to a spatial pixels data frame (a kind of raster model from the sp package) and then to raster object with the raster() function. We will remove zero values from our map by setting all values <=0 to NA (‘right=TRUE’ here means that values below and including zero become NA).
r<-estUDm2spixdf(kernelCanisAll)
rCanis<-raster(r)
#remove zeros values from our map by setting all values <=0 to NA.
rCanis <- reclassify(rCanis, cbind(-Inf, 0, NA), right=TRUE)We can also create a polygon boundary based on the kernel estimates. Here we use the getverticeshr() function fom the adehabitatHR library () and select ‘percent’ = 95 so that 95% of the estimated distribution is within the contour of our polygon:
polyCanisAll<-getverticeshr(kernelCanisAll, percent = 95)
#plot the raster and 95% contour together.
plot(rCanis)
plot(polyCanisAll,add=TRUE)
Instead of manually reproducng the above for all animal IDs let’s use a for loop to iterate through each unique ID and create our density estimation for each, plotting our 95% contour as we go:
plot(CA,axes=T)
for (i in unique(Canis$animal_ID)){
Canis.i<-Canis[which(Canis$animal_ID==i),]
CanisID.i<-data.frame(Canis.i$animal_ID)
xy.i<-cbind(Canis.i$Lon,Canis.i$Lat)
coordinates(CanisID.i)<-xy.i
kud.i <- kernelUD(Canis.i[,3], h="href", kern ="bivnorm")
r.i<-estUDm2spixdf(kud.i)
#plot(r.i)
Ras.i <-raster(r.i)
rCanis.i <- reclassify(Ras.i, cbind(-Inf, 0, NA), right=TRUE)
polyCanis.i<-getverticeshr(kud.i, percent = 95)
clipCanis.i<-mask(rCanis.i,polyCanis.i)
plot(clipCanis.i,add=TRUE,axes=T,box=F,legend=F)
plot(polyCanis.i,add=TRUE,main=i)
}
Part 2. Resource Selection
Now we know something about the options for establishing the range of our animals using telemetry data, next we’ll look at resource selection. Conceptually, resource selection is closely linked to scale response (Week 2) and to species distribution modelling (Weeks 4 and 6) in that we try to establish the preference of habitat types and model the likelihood of their use by an animal or species of interest. We will be using the ‘adehabitatHS’ library to compute selection ratios (Wi) for Design I scenarios (whole population, no differentiation between animals in terms of use and availability) and Design II scenarios (same availability but use calculated individually for each animal). Before we get to that we need to create data by extracting values from our habitat layer for use (our Canis SpatialPointsDataFrame) and for overall availability (by sampling a large number of background points).
#make sure results are replicable
set.seed(5)
back.xy <- spatSample(CA,size=20000,as.points=TRUE) #sample a regular grid of background points. Returns a points spatVector
#now we need to convert back to a spatVector to use the extract function
pres.cov <- extract(CA, vect(Canis)) #extracts values from layers at pres locations
pres.cov<-data.frame(pres.cov,pres=1)
#back.xy is a spatVector so no need to convert
back.cov<-extract(CA,back.xy)
back.cov<-data.frame(back.cov,pres=0)
#inspect
head(pres.cov)## ID AlbertaCanis pres
## 1 1 41 1
## 2 2 42 1
## 3 3 42 1
## 4 4 74 1
## 5 5 45 1
## 6 6 41 1head(back.cov)## ID AlbertaCanis pres
## 1 1 41 0
## 2 2 74 0
## 3 3 41 0
## 4 4 51 0
## 5 5 42 0
## 6 6 41 0pres.cov<-pres.cov[,2:3]
back.cov<-back.cov[,2:3]
colnames(pres.cov)<-c("Cover","Pres")
colnames(back.cov)<-c("Cover","Pres")Now we use the tapply() function which works here in a similar way as the table() function that we used last week to calculate background counts for each habitat type.This will provide our estimate of resource availability.
avail <- tapply(back.cov$Cover, back.cov$Cover, length) # here 'length' returns the length of the vector for each habitat (i.e. the count)
avail## 21 25 31 41 42 45 46 51 62 71 73 74 91
## 86 55 946 10699 4647 44 160 782 7 40 1714 811 9avail<-as.numeric(avail)Now we have our background points categorized, next we will carry out a similar cross-tabulation, calculating the number of times each animal is recorded within each habitat. We will call the dcast() function from the ‘reshape2’ library to re-format our data frame ‘pres.cov’ so that each animal is assigned to a row, columns denote habitat type and values give the frequency of occurence. You can think of this as performing a similar function to a regular pivot table (with count as the cell values). We create two data objects representing use. One with individual animal IDs (for Design II analyses comparing use to availability ratios for each animal individually) and one combining all data together (in order to model overall use to availability ratios at the population level). See this week’s lecture and suggested reading for more in-depth explanations.
#create a column in pres.cov combining all animals into one group (for Design I analysis) and another for individual animal IDs (for Design II analysis)
pres.cov$WID<-Canis$animal_ID
pres.cov$ID<-as.numeric((Canis$Taxon))
library(reshape2)
#use dcast function to format data
USECanis<-dcast(pres.cov, WID ~ Cover, length, value.var = "WID")
USEAll<-dcast(pres.cov, ID ~ Cover, length, value.var = "ID")
USECanis## WID 21 25 31 41 42 45 46 51 62 71 73 74
## 1 3 0 3 49 896 896 68 31 7 0 2 373 237
## 2 6 0 0 135 567 297 0 0 0 0 0 27 16
## 3 10 0 0 180 713 551 2 20 0 1 11 102 163
## 4 12 2 1 6 57 82 0 2 0 0 0 20 15
## 5 13 0 1 88 296 357 1 19 0 0 7 218 88
## 6 23 0 0 205 253 26 0 0 0 0 0 0 12
## 7 27 0 0 68 559 111 0 1 0 0 0 16 17
## 8 29 0 0 8 138 110 0 1 0 0 0 14 25
## 9 30 1 0 45 843 595 1 10 0 0 0 98 149
## 10 31 1 6 10 927 224 7 7 0 0 9 129 110USEAll## ID 21 25 31 41 42 45 46 51 62 71 73 74
## 1 NA 4 11 794 5249 3249 79 91 7 1 29 997 832To enter our data into the adehabitatHR wides() set of functions it is important that both our use and availability data have matching numbers of habitat types. In this case we need to remove type 91 (barren land) from ‘avail’ and we also remove NA values from both our use data frames.
#remove value 13 from the availability data as this does not appear in our use data
avail## [1] 86 55 946 10699 4647 44 160 782 7 40 1714 811
## [13] 9avail<-avail[c(-13)]
#remove NAs (col 14)
USECanis<-USECanis[-14]
USEAll<-USEAll[-14]
#check length of availability and use data to see that they now match.
length(avail)## [1] 12length(USECanis[,c(2:ncol(USECanis))])## [1] 12So we have all our data in the correct format for use with the adehabitatHS funtions.
In the next part of the practical, the widesI() function will be used to explore resource selection by animals with the Design I approach where habitat use USEAll and availability avail are measured at the population level and animals are not identified. This function tests habitat selection with the Chi-squared method (‘Khi2’ in the model results) where the Pearson and log-likelihood statistics are returned. The Manly selectivity measure (selection ratio = used/available) is computed, the preference / avoidance is tested for each habitat, and the differences between selection ratios are computed.
For widesII() we pass our Design II data (USECanis) and the same availability (avail) for all animals, where use is measured for each one. Tests of identical habitat use for all animals, and of habitat selection are also provided.
Have another look at the function descriptions, with ?function(), for widesI and widesII functions before entering the code below.
USEAll<-as.matrix(USEAll)
library(adehabitatHS)
#use selection for all animals
sel.ratioI <- widesI(u = USEAll[,c(2:ncol(USEAll))],a = avail, avknown = T, alpha = 0.05)
summary(sel.ratioI)#this shows the list of statistics available to view## Length Class Mode
## used.prop 12 -none- numeric
## se.used 12 -none- numeric
## avail.prop 12 -none- numeric
## se.avail 12 -none- numeric
## wi 12 -none- numeric
## se.wi 12 -none- numeric
## chisquwi 2 data.frame list
## Bi 12 -none- numeric
## Khi2P 3 -none- numeric
## Khi2L 3 -none- numeric
## avknown 1 -none- logical
## comparisons 4 -none- list
## profile 13 -none- character
## alpha 1 -none- numericsel.ratioI # returns a matrix of Wi values and significance##
##
## ************** Manly's Selection ratios for design I ********
##
## Significance of habitat selection:
## Khi2L df pvalue
## 1579.721 11.000 0.000
##
##
## Table of ratios (p-values should be
## compared with Bonferroni level= 0.004166667 )
## used avail Wi SE.Wi P Bi
## 21 0.000 0.004 0.082 0.041 0.000 0.007
## 25 0.001 0.003 0.352 0.106 0.000 0.028
## 31 0.070 0.047 1.479 0.051 0.000 0.118
## 41 0.463 0.535 0.865 0.009 0.000 0.069
## 42 0.286 0.232 1.232 0.018 0.000 0.098
## 45 0.007 0.002 3.164 0.355 0.000 0.252
## 46 0.008 0.008 1.002 0.105 0.982 0.080
## 51 0.001 0.039 0.016 0.006 0.000 0.001
## 62 0.000 0.000 0.252 0.252 0.003 0.020
## 71 0.003 0.002 1.278 0.237 0.241 0.102
## 73 0.088 0.086 1.025 0.031 0.417 0.082
## 74 0.073 0.041 1.808 0.060 0.000 0.144
##
##
## Bonferroni classement
## Based on 95 % confidence intervals on the differences of Wi :
##
## habitat 45 74 31 71 42 73 46 41 25 62 21 51
## 45 ----
## 74 ---- ----
## 31 --------
## 71 ---------------------------- ----
## 42 -------- ----
## 73 ---- -------- ----
## 46 -------------------- ----
## 41 ---- -------- ----
## 25 ----------------
## 62 ---- ----------------------------
## 21 ----------------
## 51 ----------------The results provided by sel.ratio I suggest that, overall, resource selection is taking place at a significant level according to the Chi-squared statistic (Khi2L = 1561.421, p<0.001) and that for the majority of our habitat types selection/availability ratios are significant as seen in the Table of Ratios output. Here values for Wi > 1 indicate preference and Wi < 1 indicates aversion, where P is significant below the Bonferroni level given and Bi gives the proportion of use by each habitat type: Bi = Wi/sum(Wi))
Let’s plot these values in the graphics pane.
plot(sel.ratioI)
The wide I() function has provided us with a breakdown of selection ratios for the population considered together. Now we’ll repeat the analysis but with plots for all animal IDs.
#use test for each wolf
sel.ratioII <- widesII(u = USECanis[,c(2:ncol(USECanis))],a = avail, avknown = T, alpha = 0.05)
row.names(sel.ratioII$wij)<-as.character(c(1:nrow(sel.ratioII$wij)))# this is to give the plot() legend a label for each animal
sel.ratioII$Khi2L1 # test of identical use of resource by animals## Khi2L1 df pvalue
## 2539.232 99.000 0.000sel.ratioII$Khi2L2 #test of overall habitat selection## Khi2L2 df pvalue
## 4118.953 110.000 0.000sel.ratioII$wij #selection ratios for each animal and resource category## 21 25 31 41 42 45 46
## 1 0.0000000 0.4256121 0.4041665 0.6534618 1.5044949 12.0590093 1.5118096
## 2 0.0000000 0.0000000 2.7378486 1.0167324 1.2261697 0.0000000 0.0000000
## 3 0.0000000 0.0000000 2.1823203 0.7643345 1.3599284 0.5213321 1.4336632
## 4 2.5130107 1.9647174 0.6853666 0.5756977 1.9067949 0.0000000 1.3507432
## 5 0.0000000 0.3381142 1.7298864 0.5144872 1.4286360 0.4226427 2.2083081
## 6 0.0000000 0.0000000 8.7340479 0.9530818 0.2255036 0.0000000 0.0000000
## 7 0.0000000 0.0000000 1.8613798 1.3529627 0.6185396 0.0000000 0.1618442
## 8 0.0000000 0.0000000 0.5711388 0.8711214 1.5986847 0.0000000 0.4221073
## 9 0.1334406 0.0000000 0.5458933 0.9042130 1.4693683 0.2608157 0.7172431
## 10 0.1625549 1.5250604 0.1477772 1.2112535 0.6738664 2.2240464 0.6116128
## 51 62 71 73 74
## 1 0.06984683 0.000000 0.3901444 1.6980615 2.2802522
## 2 0.00000000 0.000000 0.0000000 0.3022176 0.3785000
## 3 0.00000000 1.638472 3.1540591 0.6825375 2.3051749
## 4 0.00000000 0.000000 0.0000000 1.2609038 1.9986337
## 5 0.00000000 0.000000 3.2543488 2.3652210 2.0178453
## 6 0.00000000 0.000000 0.0000000 0.0000000 0.5963665
## 7 0.00000000 0.000000 0.0000000 0.2417277 0.5428068
## 8 0.00000000 0.000000 0.0000000 0.5516454 2.0819101
## 9 0.00000000 0.000000 0.0000000 0.6561477 2.1083941
## 10 0.00000000 0.000000 3.1454371 1.0521493 1.8961396sel.ratioII$wi #selection ratios overall per resource category ## 21 25 31 41 42 45 46
## 0.08197249 0.35248171 1.47923084 0.86464926 1.23220687 3.16432441 1.00236985
## 51 62 71 73 74
## 0.01577604 0.25177265 1.27774619 1.02515829 1.80804427sel.ratioII$se.wi #selection ratio SEs## 21 25 31 41 42 45 46
## 0.05169156 0.18447732 0.51134978 0.08754273 0.12775984 2.27081985 0.23326791
## 51 62 71 73 74
## 0.01378622 0.24202965 0.55614340 0.24642861 0.22945818plot(sel.ratioII)
Note that in the widesII() function, selection ratios are plotted for each animal.
Hypotheses tests performed by the wides I() and wides II() functions suggest that our animals do not exhibit identical resource use (DesignII), although resource selection is significant overall (Design I). Given the data, we have reasonable justification to assume that, on the whole, our population have a preference for patches of woodland, herbacious wetland, and locations in close proximity to water.However, selection ratios related to land-cover type 45 and 31 appear to be subject to the influence of outliers so we might consider removing these points or constraining a subsequent analysis to only home ranges. The most consistent resource use appears to be associated with wet forest and herbacious wetlands. We can also see an overall aversion to more anthropogenc land-uses. We will use these assumptions later in the course to calculate landscape resistance and least cost paths for these animals. For now continue to the following exercises:
Week 3 exercises
Use the kernel Use Density approach to map the core home ranges (polygon) of animals 1 to 3 in our dataset (animal IDs 3,6 and 10). Produce spatial plots (polygons of core ranges) and provide code that will calculate the core range in km2 of each animal and all three animals combined.
Write a 300 word overview of resource selection analysis, summarizing its role in movement ecology from a mammal conservation perspective.
Here are some links to relevant reading on this topic:
[https://royalsocietypublishing.org/doi/full/10.1098/rstb.2010.0087]
[https://link.springer.com/chapter/10.1007/978-3-319-51272-3_1]