Birds Survey

Birds Survey

Detection and Occupancy prac 4. Instructions and code

#package probably already installed. If not:
# click in the toolbar Tools -> Install Package -> type in the package name and click install, or try:

install.packages(“unmarked”)
# Once you have installed a package onto your computer, you won’t have to do it again. Unless there’s an updated version.

# load the library into your current analysis session (needs to be reloaded each session)

library(unmarked)
##### Part 1: THE DATA #####

#download the data from cloud deakin, and put it somewhere where you can access it and know the directory name.

### You need to tell R which folder has the file you will use, and where to save output by setting the working directory (it might be different from the one given
below)

## Step 1: Set your Working Directory – Tell R where to find that data on your computer
#you can use the menus to do this: File>change.dir> and then select the directory you have your data file in.

# or using code: for example (you will need to modify that to point to your directory:
setwd(“E:\AA Deakin\Landscape Ecology Teaching SLE322\2017 Pracs\prac 1. neph unmarked prac”)
## Step 2: load data – Is the data you want to load in saved as a .csv file?
#if you have set your directory correctly:

birdata<-read.csv(“neph_unmarked_sle322.csv”)
# change to your “file name”, in your case “siteandsurveycovars4analysis.csv”

# If you can’t remember the exact pathway your data is saved on. Instead use a function file.choose() to open up a window and select where the file is saved.

birdata<-read.csv(file.choose())
## Step 3: Save R script
#STOP!!!! Before you move on…
#You will want to save this modified version of your R script as a different file than R_code_landscape_ecology. So..
#Step 1. Click File
#Step 2. Click Save As
#Step 3. Change the file name (e.g. Landscape_project_script)
#Step 4. Click Save
#
#Anytime that the text on the file name of your script is red that means R has not saved your changes to the script. When you move between different sections of this
script I would highly recommend you to remember to click save.
## Step 4: check that your data were loaded into R correctly
#check the dimensions of the data
dim(birdata)

#check that the top of the dataset looks ok
head(birdata)

#check that the bottom of the dataset looks ok; no extra lines or messy bits
tail(birdata)

#to see a list of column headings. Identify if the data file: the three observation columns,
# the two site-level covariates, and the two observation-level covariates.
names(birdata)

# other ways to double check that your data has been read in correctly
summary(birdata) # will give you a summary of each column. I find this helpful for looking to make sure there are
# no erroneous data entries (e.g. typed in year as 20017) or to see if you have any missing data (e.g Na= 1)
str(birdata) # will give you a read out about the str of your data. This is helpful when running statistics and
# you want to double check to see if R has read in your factor levels
##### Part 2: UNMARKED ANALYSIS OF DETECTION AND OCCUPANCY WITH SITE AND OBSERVATION COVARIATES ###
#library(unmarked) #just in case not yet opened

## Step 1: Create our dataframe
#make a data frame for the unmarked analysis using covariates for sites and for observations.
#This is where unmarked takes the data from to do the analysis
#you must first choose a species from birdata (there are three to choose from, selected because the models would fit; much more work is needed to find models to fit
other species)
#you need to update the column.number in the code below to refer to columns in the birdata (there are 8 columns for each species representing your presence/absence
survey data)

frame2<-unmarkedFrameOccu(
y=birdata[start.column.number:end.column.number],
siteCovs= as.data.frame(birdata[,c(4,5,6,7)]),
obsCovs=list(
time=as.data.frame(birdata[,c(grep(“time”,colnames(birdata)))]),
wind=as.data.frame(birdata[,c(grep(“wind”,colnames(birdata)))]),
weather=as.data.frame(birdata[,c(grep(“weather”,colnames(birdata)))]))
)
## Step 2: Model selection process for detection part of our model
# We need to find a model, or set of models, that best describes the factors influencing detection and occupancy.
# We set up a range of alternative detection models first, keeping occupancy as just an intercept model (~1) for this first step.
# This minimises the number of model combinations that need to be assessed by not having to use every combination of detection model with every combination of
occupancy model.
# the detection covariates we will use are weather, time and wind, only as main effects. We will not test for interactions to keep things simple

#use this code to evaluate detection covariates

model.null<-occu(~1~1,data=frame2) #no covariates fitted
model.1.1<-occu(~weather ~1,data=frame2)
model.1.2<-occu(~time ~1,data=frame2)
model.1.3<-occu(~wind ~1,data=frame2)
model.2.1<-occu(~weather+time ~1,data=frame2)
model.2.2<-occu(~weather+wind ~1,data=frame2)
model.2.3<-occu(~wind+time ~1,data=frame2)
model.3<-occu(~weather+wind+time ~1,data=frame2)

model.list<-fitList(null=model.null, model.1.1, model.1.2,model.1.3, model.2.1, model.2.2, model.2.3, model.3)
#fitList is the unmarked function to organise models for selection

model.comparison<-as(modSel(model.list,nullmod=”null”),”data.frame”)

# and we can print to a csv file to open in excel etc.
write.csv(model.comparison, “model_selection_detection.csv”)#the part in “” will be the file name in the directory you specified in setwd() or with the drop down menu
previously
### Step 3: compare models and select the model with the best fit (i.e. lowest AICc)
#You need to determine which model of your survey covariates has the best fit. Use the AIC to determine the best fit; the lowest value = best model.

#the best detection covariates are those in the top ranked model. (there may be other similarly-ranked models, but in this study we will just use the top ranked model)
#find your top ranked model as the first line in model.comparison

model.comparison

#or you can use this code to get it (just paste it in to R for the prac, you can try to understand it in your spare time if you are aiming to become a better R user)
best.detection.covars<-gsub(” ~ 1″, “”, model.comparison$formula[1])
if(length(best.detection.covars)==0) best.detection.covars <-“~1″
best.detection.covars
##################################################
### Step 4: Model selection process for occupancy part of our model
#OCCUPANCY AND DETECTION COVARIATE MODELS
#You will need copy the formula for the survey covariates from the best survey covariates models above into each model formula below
#to do that you paste into the green bits below the formula, should be something like ~wind + weather for example
#the models allow us to explore:
# two landscape variables:
# Roads (length of roads within 1km)
# Native (area of native vegetation within 1km)
#and we will use the site variables pc_exoticcanopy and site_type
#these covariates were selected to minimise correlations among explanatory variables (eg; mown grass highly correlated with site_type and native canopy cover)
#to keep things manageable, we will not look at interactions in these models

#note: if you are using eastern yellow robin, some of the models don’t work properly. Delete model 3.3 and 4.1 from the list below, and from the model.comps line below
that.

model.null<-occu(put.your.detect.covars.here~1,data=frame2)
model.1.1<-occu(put.your.detect.covars.here ~site_type,data=frame2)
model.1.2<-occu(put.your.detect.covars.here ~pc_exoticcanopy,data=frame2)
model.1.3<-occu(put.your.detect.covars.here ~tot_road_lengths_meters,data=frame2)
model.1.4<-occu(put.your.detect.covars.here ~tot_native_wood_area_meters,data=frame2)
model.2.1<-occu(put.your.detect.covars.here ~site_type+pc_exoticcanopy,data=frame2)
model.2.2<-occu(put.your.detect.covars.here ~site_type+tot_road_lengths_meters,data=frame2)
model.2.3<-occu(put.your.detect.covars.here ~site_type+tot_native_wood_area_meters,data=frame2)
model.2.4<-occu(put.your.detect.covars.here ~pc_exoticcanopy+tot_road_lengths_meters,data=frame2)
model.2.5<-occu(put.your.detect.covars.here ~pc_exoticcanopy+tot_native_wood_area_meters,data=frame2)
model.2.6<-occu(put.your.detect.covars.here ~tot_road_lengths_meters+tot_native_wood_area_meters,data=frame2)
model.3.1<-occu(put.your.detect.covars.here ~site_type+pc_exoticcanopy+tot_road_lengths_meters,data=frame2)
model.3.2<-occu(put.your.detect.covars.here ~site_type+pc_exoticcanopy+tot_native_wood_area_meters,data=frame2)
model.3.3<-occu(put.your.detect.covars.here ~site_type+tot_road_lengths_meters+tot_native_wood_area_meters,data=frame2)
model.3.4<-occu(put.your.detect.covars.here ~pc_exoticcanopy+tot_road_lengths_meters+tot_native_wood_area_meters,data=frame2)
model.4.1<-occu(put.your.detect.covars.here ~site_type+pc_exoticcanopy+tot_road_lengths_meters+tot_native_wood_area_meters,data=frame2)

#note: if you are using eastern yellow robin, some of the models don’t work properly. Delete model 3.3 and 4.1 from the model.comps line below.

#Put the models we can fit into the list
model.comps<-fitList(null=model.null,model.1.1, model.1.2, model.1.3, model.1.4, model.2.1, model.2.2, model.2.3, model.2.4, model.2.5, model.2.6, model.3.1,
model.3.2, model.3.3, model.3.4, model.4.1)

#do the model ranking based on AIC
model.comps2<-as(modSel(model.comps,nullmod=”null”),”data.frame”)

### Step 5: compare models and select the model with the best fit (i.e. lowest AIC)
model.comps2

#save our results to csv for later.
write.csv(model.comps2, “model_selection_sitespres_abs.csv”)

######################################
##### Part 3: PLOT THE RESULTS
#now we’ve got the best model, we want to see how detection and occupancy depend on our covariates.

#include just the models you want to model average over (ie: all those within 2 AIC of the best).
# In your data, some species have more than one model within 2 AIC, so you can average over those
#The function will automatically name your models, with a warning
# (ie: don’t worry about warning).

fmList<-fitList(model.tc_nands) # Change the arbitrary name of your best fit model

#for eg, it might look something like:
fmList<-fitList(model.1.3, model.2.6, model.2.4)

#once you know which models you are plotting, you know what variables you need to plot (the covariates in the model). You need to make a plot of each covariate
separately.
#To create the file called”newdat” below, you must include values for every variable that was in the best model for that component of the model (where component means
either the detection component or the occupancy component; in other words, include either all the detection, or all the occupancy covariates).

#Specify the range of values you want to get predictions for (that will be for one covariate only). For covariates that were in the model but you are not plotting at
the moment, you need to specify what value they should take. Use the mean for continuous variables (when values are numbers and not categories). Use all levels for
variables that are factors (categorical) (where levels means the different values a factor could take
#eg for the factor weather
#> levels(birdata$weather.1)
#[1] “cloud” “rain” “sun”
#so it has 3 levels; cloud, rain sun

#if you are plotting a detection covariate, there are multiple columns of data, and we want to explore the whole range, so specify the 8 columns for that variable. A
hypothetical example given here

#if the best models included ~weather as the detection covariate and tot_road_lengths_meters, tot_native_wood_area_meters, pc_exoticcanopy were in the best (or best
set) model(s) as site covariates (ie: they influenced occupancy)

# then first we want to plot the detection covariate; weather

newdat<-expand.grid(
weather= c(“cloud”, “rain”, “sun”),
name.of.other.covariate.if.there.was.one=c(“Put the levels of that covariate here, using quotes, with commas separating each levelor the mean of a continuous
variable”)) #see further examples below if you need to use the mean of a continuous variable

#in this hypothetical example, there isn’t an extra detection covariate, so newdat looks like:
newdat<-expand.grid(
weather= c(“cloud”, “rain”, “sun”))
#the predict function takes the models in fmList and predicts values, given the values of the explanatory variables that were in the model and specified in newdat. It
also calculates confidence intervals as a measure of variability in the data

det.predict<-predict(fmList, type=”det”, newdata=newdat, appendData=TRUE)
#note type= “det” or “state” (“state” for sites, “det” for detection)

#name the upper and lower confidence intervals for use in plotting and get the mean values
upperci<-det.predict$upper
lowerci<- det.predict$lower
meanvalues<-det.predict$Predicted

#if you are plotting a factor with specific levels, like weather, then you need to use the next code. See further down for if you are plotting a continuous variable.
ylabel<-“Probability of detection”
maxy<- 1.05
barnames<- c(“cloud”, “rain”, “sun”)#update levels to be the ones for the factor you are plotting. in this case, it is set up for weather. but you might also have
wind, which has levels “light”, “moderate”, “strong”, or other categorical covariates.

par(mar= c(10, 6, 2, 2) + 0.1) # this sets parameters for your plot. no need to change them most likely c(bottom, left, top, right) number of lines of margin; adjust
bottom to suit length of barnames
mp<- barplot(meanvalues, space = 0.5,
names.arg = barnames, beside = TRUE,
horiz = FALSE, prcol = NULL, border = par(“fg”),
ylim = c(0, maxy),ylab = ylabel, xpd = FALSE, las=1,lwd=2,xaxt=”s”,
cex.axis = 2.0, cex.names = 2.0, cex.lab=2.0, mgp=c(3.5,1,0))

#add error bars to plot
arrows(mp, lowerci ,mp, upperci, code=3, angle=90, length=.1, lwd=2)

#add bottom line
abline(0,0,lwd=3)
#now you’ve got your graph you need to save it. Two ways to do this.

savePlot(filename=”My detection bar graph.wmf”,type=”wmf”)#rename as you prefer

#or you can right-click on the graph and select save as metafile (or copy as metafile and paste straight into a word document).

#on the other hand, if you needed to plot a continuous variable like time, then you need to update the code below. Note that time is in minutes, and values in the
afternoon had minutes subtracted after 13:00. So, for eg. 16:00 has the same value as 10:00 (10 x 60 = 600). This is because 1pm is typically the worst time see
birds, but it is good in the morning and improves later in the afternoon

newdat<-expand.grid(
time= seq(min(birdata[,grep(“time”, colnames(birdata))], na.rm=TRUE), max(birdata[,grep(“time”, colnames(birdata))],na.rm=TRUE),length.out=100),
name.of.other.covariate.if.there.was.one=c(“if a factor, put the levels of that covariate here, using quotes, with commas separating each level, or put mean
value if a continuous variable”))

#in this hypothetical example where time is the only detection covariate, newdat looks like:
newdat<-expand.grid(
time= seq(min(birdata[,grep(“time”, colnames(birdata))], na.rm=TRUE), max(birdata[,grep(“time”, colnames(birdata))],na.rm=TRUE),length.out=100))
#the predict function takes the models in fmList and predicts values, given the values of the explanatory variables that were in the model and specified in newdat. It
also calculates confidence intervals as a measure of variability in the data

det.predict<-predict(fmList, type=”det”, newdata=newdat, appendData=TRUE)
#note type= “det” or “state” (“state” for sites, “det” for detection)
#name the upper and lower confidence intervals for use in plotting and get the mean values
upperci<-det.predict$upper
lowerci<- det.predict$lower
meanvalues<-det.predict$Predicted

#You will need to change the variables in green below to be the covariate that you want plotted on the x-axis. You will also need to modify xmin and xinterval to
match the spread of your x-axis data. EG. if you are plotting % tree cover, your values might span from 5 to 95%, so it might be sensible to set xmin to 0 and
xinterval to 20.
#specify this list of minimums, maximums, and intervals so that legend is put in the right place.
ymax<-1
ymin<-0
xmax<-ceiling(max(birdata[,grep(“time”, colnames(birdata))],na.rm=TRUE))#set maximum x value. change “time” to be the text of the variable you are plotting
xmin<-200#min time was 250
yinterval<-0.2
xinterval<-100#intervals of 100 ok for time which spans up to 833

#creates the plot
plot(Predicted~time,data=det.predict[1:100,],
type=”n”,ann=F,axes=F,ylim=c(0,ymax)) #sets up the space to do the plot, without drawing the plot
axis(1,at=seq(xmin, xmax,xinterval) ,cex.axis=1.3) #add axis
axis(2,at=seq(ymin,ymax,yinterval),las=2,cex.axis=1.3) #add axis
box(bty=’l’) #what kind of box to draw around plot
title(xlab=”time (minutes)”,ylab=”Probability of detection”,cex.lab=1.5) #use title() function to add x and y axis labels. Change the text here to be what you want you
x label to be

# cex “a value controlling the size of texts and symbols with respect to the default…”
# las ” a integer which controls the orientation of the axis labels …”

# adds the predicted values and confidence intervals for clear/rain level across the range of temperatures we specified in newdat
points(Predicted~time,data=det.predict[1:100,],type=’l’,lty=1,lwd=2) #draw on the predicted line
xs<- c(det.predict$time[1:100],sort(det.predict$time[1:100],decreasing=TRUE)) #it assumes that it joins first to last point
ys<- c(upperci[1:100],lowerci[1:100][order(det.predict$time[1:100],decreasing=TRUE)])
polygon(xs,ys, col=rgb(0,0,0,0.1), border=NA) #draw on the confidence bands as a polygon
savePlot(filename=”My detection continuous graph.wmf”,type=”wmf”) #update graph name
#Next plotting occupancy

#To create the file called”newdat” below, you must include values for every variable that was in the best model for that component of the model (where component means
either the detection component or the occupancy component; in other words, include either all the detection, or all the occupancy covariates).

#Specify the range of values you want to get predictions for (that will be for one covariate only). For covariates that were in the model but you are not plotting at
the moment, you need to specify what value they should take. Use the mean for continuous variables (when values are numbers and not categories). Use all levels for
variables that are factors (categorical) (where levels means the different values a factor could take).
#if the best models included ~weather as the detection covariate and site_type+tot_road_lengths_meters were in the best (or best set) of model(s) as site covariates
(ie: they influenced occupancy)

# as a hypothetical example, we now want to plot site_type, which is a categorical covariate (a factor).
#> levels(birdata$site_type)
#[1] “natural” “park” “restoration” #has three levels

newdat<-expand.grid(
site_type= c(“natural”, “park”, “restoration”),
name.of.other.covariate.if.there.was.one=c(“Put the levels of that covariate here, using quotes, with commas separating each level or the mean of a continuous
variable”))

#in this hypothetical example, there is also tot_road_lengths_meters. So newdat looks like:
newdat<-expand.grid(
site_type= c(“natural”, “park”, “restoration”),
tot_road_lengths_meters=(mean(birdata$tot_road_lengths_meters)))
#the predict function takes the models in fmList and predicts values, given the values of the explanatory variables that were in the model and specified in newdat. It
also calculates confidence intervals as a measure of variability in the data

det.predict<-predict(fmList, type=”state”, newdata=newdat, appendData=TRUE)
#note type= “det” or “state” (“state” for sites, “det” for detection)

#name the upper and lower confidence intervals for use in plotting and get the mean values
upperci<-det.predict$upper
lowerci<- det.predict$lower
meanvalues<-det.predict$Predicted

#if you are plotting a factor with specific levels, like weather, then you need to use the next code. See further down for if you are plotting a continuous variable.
ylabel<-“Proportion of sites occupied”
maxy<- 1.05
barnames<- c(“natural”, “park”, “restoration”)#update levels to be the ones for the factor you are plotting. in this case, it is set up for site_type.

par(mar= c(10, 6, 2, 2) + 0.1) # this sets parameters for your plot. no need to change them most likely c(bottom, left, top, right) number of lines of margin; adjust
bottom to suit length of barnames
mp<- barplot(meanvalues, space = 0.5,
names.arg = barnames, beside = TRUE,
horiz = FALSE, prcol = NULL, border = par(“fg”),
ylim = c(0, maxy),ylab = ylabel, xpd = FALSE, las=1,lwd=2,xaxt=”s”,
cex.axis = 2.0, cex.names = 2.0, cex.lab=2.0, mgp=c(3.5,1,0))

#add error bars to plot
arrows(mp, lowerci ,mp, upperci, code=3, angle=90, length=.1, lwd=2)

#add bottom line
abline(0,0,lwd=3)
#now you’ve got your graph you need to save it. Two ways to do this.

savePlot(filename=”My occurence bar graph.wmf”,type=”wmf”)#rename as you prefer

#or you can right-click on the graph and select save as metafile (or copy as metafile and paste straight into a word document).

#on the other hand, if you needed to plot a continuous variable like tot_road_lengths_meters, then you need to update the code below.

newdat<-expand.grid(
tot_road_lengths_meters = seq(min(birdata[,grep(“road”, colnames(birdata))], na.rm=TRUE), max(birdata[,grep(“road”, colnames(birdata))],
na.rm=TRUE),length.out=100),
name.of.other.covariate.if.there.was.one=c(“if a factor, put the levels of that covariate here, using quotes, with commas separating each level, or put mean
value if a continuous variable”))

#in this hypothetical example where tot_road_lengths_metersis in the model with site_type, , newdat looks like:
newdat<-expand.grid(
tot_road_lengths_meters = seq(min(birdata[,grep(“road”, colnames(birdata))], na.rm=TRUE), max(birdata[,grep(“road”, colnames(birdata))])),
site_type= c(“natural”, “park”, “restoration”))
#the predict function takes the models in fmList and predicts values, given the values of the explanatory variables that were in the model and specified in newdat. It
also calculates confidence intervals as a measure of variability in the data

det.predict<-predict(fmList, type=”state”, newdata=newdat, appendData=TRUE)
#note type= “det” or “state” (“state” for sites, “det” for detection)
#name the upper and lower confidence intervals for use in plotting and get the mean values
upperci<-det.predict$upper
lowerci<- det.predict$lower
meanvalues<-det.predict$Predicted
#You will need to change the variables in green below to be the covariate that you want plotted on the x-axis. You will also need to modify xmin and xinterval to
match the spread of your x-axis data. EG. if you are plotting % tree cover, your values might span from 5 to 95%, so it might be sensible to set xmin to 0 and
xinterval to 20.
#specify this list of minimums, maximums, and intervals so that legend is put in the right place.
ymax<-1
ymin<-0
xmax<-ceiling(max(birdata[,grep(“road”, colnames(birdata))],na.rm=TRUE))#set maximum x value. change “time” to be the text of the variable you are plotting
xmin<-0#min road was 1
yinterval<-0.2
xinterval<-20
#intervals of 20 ok for road which spans up to 120 (the scale is not in Meters though, unlike name of variable implies. it is km).

#creates the plot., note use of 1:100 in next code line: this refers to the first 100 lines of det.predict table. That means it will be plotting at values where
site_type is “natural” because “natural” is the first category listed in det.predict (the others are 101:200, and 201:300).

plot(Predicted~tot_road_lengths_meters,data=det.predict[1:100,],
type=”n”,ann=F,axes=F,ylim=c(0,ymax)) #sets up the space to do the plot, without drawing the plot
axis(1,at=seq(xmin, xmax,xinterval) ,cex.axis=1.3) #add axis
axis(2,at=seq(ymin,ymax,yinterval),las=2,cex.axis=1.3) #add axis
box(bty=’l’) #what kind of box to draw around plot
title(xlab=”Road length (km)”,ylab=”Probability of detection”,cex.lab=1.5) #use title() function to add x and y axis labels. Change the text here to be what you want
you x label to be

# cex “a value controlling the size of texts and symbols with respect to the default…”
# las ” a integer which controls the orientation of the axis labels …”

# adds the predicted values and confidence intervals for clear/rain level across the range of temperatures we specified in newdat
points(Predicted~tot_road_lengths_meters,data=det.predict[1:100,],type=’l’,lty=1,lwd=2) #draw on the predicted line
xs<- c(det.predict$tot_road_lengths_meters[1:100],sort(det.predict$tot_road_lengths_meters[1:100],decreasing=TRUE)) #it assumes that it joins first to last point
ys<- c(upperci[1:100],lowerci[1:100][order(det.predict$tot_road_lengths_meters[1:100],decreasing=TRUE)])
polygon(xs,ys, col=rgb(0,0,0,0.1), border=NA) #draw on the confidence bands as a polygon
savePlot(filename=”My detection continuous graph of occurrence.wmf”,type=”wmf”) #update graph name

# Landscape Ecology Code
# 2017

##### Prac 4 Bird Occupancy and Detection modelling
# Part 1. Revision
### Vectors, Matrices and Data Frames

# What is a vector?
# a sequence of data elements of the same basic type.
numbers <- c(1, 2 ,3)
letters <- c(“a”, “b”, “c”)
statements <- c(TRUE, FALSE, TRUE)

# What is a matrix?
# matrix(c(vector of numbers in column order),# of rows, # of columns)
our.first.matrix <- matrix(c(1, 2, 3, 4, 5, 6, 7, 8, 9, 10), 2, 5)

# What is a data frame?
# A bunch of vectors
# Just like a matrix except you can mix columns of different kinds of data. One can be numbers,
# while another can be text.

sex <- c(“M”, “F”, “F”, “M”)
thorax <- c(3.2, 15.4, 14.3, 4.1)
better.solution<- data.frame(sex, thorax)
rownames(better.solution) <- c(“Spider 1”, “Spider 2”, “Spider 3”, “Spider 4″)
better.solution

# what is a function?
# A function is a defined arguement. For instance you want the sum of a vector
# Like excel
numbers <- c(1, 2, 3)
sum(numbers)
help(sum)
# Subsetting data: e.g. select data from rows 2 and 3 (i.e. female spiders)
spider2_3 <- better.solution[2:3, ]
spider2_3

spiderF <- better.solution[2:3,]
spiderF

spiderF1 <- better.solution[-c(1,4),]
spiderF1

spiderFemale <- better.solution[better.solution$sex==”F”,]
spiderFemale

######
#
# Other helpful resources
#
##########

# google
# stackoverflow discussion boards
# github
# coursera
# codeacademy
##################################################################################################
# Part 2. Occupancy Modelling

install.packages(“unmarked”)

library(unmarked)

#############

setwd(“~/Desktop/Demo – Sem 1/Don/4) Assignment prac”) # Change for your working directory (e.g. your desktop)

birds <- read.csv(“siteandsurveycovars4analysis.csv”)

# Check that data were loaded correctly
dim(birds)
head(birds)
tail(birds)
names(birds)
summary(birds)
str(birds)

### Data ###

# Presence/absence of birds:
# Brown thornbill: columns 32:39
# Pied burrawong: columns 40:79
# Eastern yellow robin: columns 48:55

# 120 sites
# 10-11 surveys at each site

# Aim:
# 1) what variables affect the likeihood of detecting birds?
# 2) what varuables affect the likeihood that a site will be occupied?

# Observations affecting detection:
# 1) Weather: cloud, rain, sun
# 2) Time
# 3) Wind: light, moderate, strong

# Factors affecting occupancy:
# 1) Site type: largely natural, restored site, or manicured park
# 2) Exotic canopy cover (%)
# 3) Total road lengths (m) in 1 km buffer
# 4) Total native wood area
######### Example for Eastern Yellow Robin ##############

# UNMARKED ANALYSIS OF DETECTION AND OCCUPANCY (no covariates)

## Unmarked package ####

## Inputs – in general ###

# 1) Detection histories
# 2) Site-specific covariates: occupancy model
# 3) Sample-specific covariates: detection model

#create a simple unmarked data frame with just the occurrence data in it (i.e. columns 2, 3, 4)
frame1 <- unmarkedFrameOccu(y = birds[, 48:55])

## Simple occupancy model:
analysis1 <- occu (~1~1, frame1)
analysis1

# Backtransform
backTransform(analysis1,’det’) # detection estimates

backTransform(analysis1, ‘state’) # Occupancy estimates

# Data are binomial so they are automatically transformed to logit scale. To be able to interpret, transform
# back to an untransformed scale

##################################################

#UNMARKED ANALYSIS OF DETECTION AND OCCUPANCY WITH SITE AND OBSERVATION COVARIATES

# Original #
frame2 <- unmarkedFrameOccu(
y = birds[, 48:55], # define data/observations
siteCovs = as.data.frame(birds[,c(1:7)]), # set site covariates: in columns 7 and 8 (nearest occupied patch, spinifex)
obsCovs = list( # Set detection covariates: i.e. our survey covariates
weather = as.data.frame(birds[, 16:23]), # e.g. weather observations
wind = as.data.frame(birds[, 24:31]), # e.g. wind observations
time = as.data.frame(birds[, 8:15]))) # e.g. time of survey

### Model selection ###

## DETECTION COVARIATE MODELS – Use null model

model.null<-occu(~1~1, data = frame2) #no covariates fitted
model.1.1<-occu(~weather ~1, data=frame2)
model.1.2<-occu(~time ~1, data=frame2)
model.1.3<-occu(~wind ~1, data=frame2)
model.2.1<-occu(~weather+time ~1, data=frame2)
model.2.2<-occu(~weather+wind ~1, data=frame2)
model.2.3<-occu(~wind+time ~1, data=frame2)
model.3<-occu(~weather+wind+time ~1, data=frame2)

model.list<-fitList(null=model.null, model.1.1, model.1.2, model.1.3, model.2.1, model.2.2, model.2.3, model.3)
#fitList is the unmarked function to organise models for selection

#do the model ranking and model selection process
model.comparison<- as(modSel(model.list, nullmod = “null”), “data.frame”)
model.comparison

# Print to a csv file to open in excel etc.
write.csv(model.comparison, “model_selection_detection_robin.csv”) #the part in “” will be the file name in the directory you specified in setwd() or with the drop
down menu previously

##################################################

### OCCUPANCY AND DETECTION COVARIATE MODELS
# altering the occupancy part of the model
# For this assignment, we are only using main effects, no interactions between covariates

model.null<-occu(~weather ~1, data = frame2)
model.1.1<-occu(~weather ~site_type, data = frame2)
model.1.2<-occu(~weather ~pc_exoticcanopy, data = frame2)
model.1.3<-occu(~weather ~tot_road_lengths_meters, data = frame2)
model.1.4<-occu(~weather ~tot_native_wood_area_meters, data = frame2)
model.2.1<-occu(~weather ~site_type+pc_exoticcanopy, data = frame2)
model.2.2<-occu(~weather ~site_type+tot_road_lengths_meters, data = frame2)
model.2.3<-occu(~weather ~site_type+tot_native_wood_area_meters, data = frame2)
model.2.4<-occu(~weather ~pc_exoticcanopy+tot_road_lengths_meters, data = frame2)
model.2.5<-occu(~weather ~pc_exoticcanopy+tot_native_wood_area_meters, data = frame2)
model.2.6<-occu(~weather ~tot_road_lengths_meters+tot_native_wood_area_meters, data = frame2)
model.3.1<-occu(~weather ~site_type+pc_exoticcanopy+tot_road_lengths_meters, data = frame2)
model.3.2<-occu(~weather ~site_type+pc_exoticcanopy+tot_native_wood_area_meters, data = frame2)
model.3.3<-occu(~weather ~site_type+tot_road_lengths_meters+tot_native_wood_area_meters, data = frame2)
model.3.4<-occu(~weather ~pc_exoticcanopy+tot_road_lengths_meters+tot_native_wood_area_meters, data = frame2)
model.4.1<-occu(~weather ~site_type+ pc_exoticcanopy+tot_road_lengths_meters+tot_native_wood_area_meters, data = frame2)

#Put the models we can fit into the list
model.comps<- fitList(null = model.null, model.1.1, model.1.2, model.1.3, model.1.4,
model.2.1, model.2.2, model.2.3, model.2.4, model.2.5, model.2.6,
model.3.1, model.3.2, model.3.3, model.3.4, model.4.1)

#do the model ranking and model selection process
model.comps2 <- as(modSel(model.comps, nullmod = “null”), “data.frame”)

AICc <- model.comps2$AIC * (model.comps2$n / (model.comps2$n – model.comps2$nPars-1))
delta.small<- AICc – min(AICc)
likli.small<- exp(-0.5 * delta.small)
AICwt.small <- likli.small / (sum(likli.small))

model.comps3 <- data.frame(model.comps2, AICc, delta.small, likli.small, AICwt.small)[order(AICc),] #rank results by AICc

#ditch the model estimates from this table because we will extract them in a later step.
model.comps3a <- model.comps3[,c(1:2, which(colnames(model.comps3) == “negLogLike”) : length(model.comps3))]
model.comps3a

#save our results to csv for later.
write.csv(model.comps3a, “model_selection_bird_surveys_robin.csv”)

######################################

#### PLOT THE RESULTS

fmList <- fitList(model.1.3) # change to your best model
##### Plot 1: # How was detection influenced by weather?

newdat<-expand.grid(
weather= c(“cloud”, “rain”, “sun”))
det.predict<-predict(fmList, type=”det”, newdata=newdat, appendData=TRUE)
upperci<-det.predict$upper
lowerci<- det.predict$lower
meanvalues<-det.predict$Predicted
ylabel<-“Probability of detection”
maxy<- 1.05
barnames<- c(“cloud”, “rain”, “sun”) #update levels to be the ones for the factor you are plotting. in this case, it is set up for weather. but you might also have
wind, which has levels “light”, “moderate”, “strong”, or other categorical covariates.
par(mar= c(10, 6, 2, 2) + 0.1) # this sets parameters for your plot. no need to change them most likely c(bottom, left, top, right) number of lines of margin; adjust
bottom to suit length of barnames
mp<- barplot(meanvalues, space = 0.5,
names.arg = barnames, beside = TRUE,
horiz = FALSE, prcol = NULL, border = par(“fg”),
ylim = c(0, maxy),ylab = ylabel, xpd = FALSE, las=1,lwd=2,xaxt=”s”,
cex.axis = 2.0, cex.names = 2.0, cex.lab=2.0, mgp=c(3.5,1,0))

arrows(mp, lowerci ,mp, upperci, code=3, angle=90, length=.1, lwd=2)
abline(0,0,lwd=3)

#Save plot:
savePlot(filename=”Detection_plot_weather.wmf”,type=”wmf”) #rename as you prefer
### Plot 2: How occupancy is influenced by road length?

newdat<-expand.grid(
tot_road_lengths_meters = seq(min(birds[,grep(“road”, colnames(birds))], na.rm=TRUE), max(birds[,grep(“road”, colnames(birds))]) ,length.out=100),
site_type= c(“natural”, “park”, “restoration”),
tot_native_wood_area_meters=mean(birds$tot_native_wood_area_meters))
det.predict<-predict(fmList, type=”state”, newdata=newdat, appendData=TRUE)
upperci<-det.predict$upper
lowerci<- det.predict$lower
meanvalues<-det.predict$Predicted
ymax<-1
ymin<-0
xmax<-ceiling(max(birds[,grep(“road”, colnames(birds))], na.rm=TRUE)) #set maximum x value. change “time” to be the text of the variable you are plotting
xmin<-0 #min road was 1
yinterval<-0.2
xinterval<-20
plot(Predicted~ tot_road_lengths_meters,data=det.predict[1:100,],
type=”n”,ann=F,axes=F,ylim=c(0,ymax)) #sets up the space to do the plot, without drawing the plot
axis(1,at=seq(xmin, xmax,xinterval) ,cex.axis=1.3) #add axis
axis(2,at=seq(ymin,ymax,yinterval),las=2,cex.axis=1.3) #add axis
box(bty=’l’) #what kind of box to draw around plot
title(xlab=”Road length (km)”,ylab=”Probability of occupancy”,cex.lab=1.5) #use title() function to add x and y axis labels. Change the text here to be what you want
you x label to be
points(Predicted~ tot_road_lengths_meters, data=det.predict[1:100,], type=’l’,lty=1, lwd=2) #draw on the predicted line
xs<- c(det.predict$tot_road_lengths_meters[1:100], sort(det.predict$tot_road_lengths_meters [1:100],decreasing=TRUE)) #it assumes that it joins first to last point
ys<- c(upperci[1:100],lowerci[1:100][order(det.predict$tot_road_lengths_meters [1:100],decreasing=TRUE)])
polygon(xs,ys, col=rgb(0,0,0,0.1), border=NA) #draw on the confidence bands as a polygon

# Save plot
savePlot(filename=”Occurrence_plot_roads.wmf”,type=”wmf”) #update graph name
################################################################################

##### Plots for other species ######

fmList <- fitList(model.1.3) # change to your best model

### Plot for occupancy: Site Type

newdat<-expand.grid(
site_type= c(“natural”, “park”, “restoration”),
tot_road_lengths_meters=(mean(birds$tot_road_lengths_meters)),
tot_native_wood_area_meters=mean(birds$tot_native_wood_area_meters))
det.predict<-predict(fmList, type=”state”, newdata=newdat, appendData=TRUE)
upperci<-det.predict$upper
lowerci<- det.predict$lower
meanvalues<-det.predict$Predicted
ylabel<-“Proportion of sites occupied”
maxy<- 1.05
barnames<- c(“natural”, “park”, “restoration”) #update levels to be the ones for the factor you are plotting. in this case, it is set up for site_type.
par(mar= c(10, 6, 2, 2) + 0.1) # this sets parameters for your plot. no need to change them most likely c(bottom, left, top, right) number of lines of margin; adjust
bottom to suit length of barnames
mp<- barplot(meanvalues, space = 0.5,
names.arg = barnames, beside = TRUE,
horiz = FALSE, prcol = NULL, border = par(“fg”),
ylim = c(0, maxy),ylab = ylabel, xpd = FALSE, las=1,lwd=2,xaxt=”s”,
cex.axis = 2.0, cex.names = 2.0, cex.lab=2.0, mgp=c(3.5,1,0))
arrows(mp, lowerci ,mp, upperci, code=3, angle=90, length=.1, lwd=2)
abline(0,0,lwd=3)

# Save plot
savePlot(filename=”Occurence_plot_site_type.wmf”,type=”wmf”) #rename as you prefer
### Plot for occupancy: Native wood area

newdat<-expand.grid(
tot_native_wood_area_meters = seq(min(birds[,grep(“wood_area”, colnames(birds))], na.rm=TRUE), max(birds[,grep(“wood_area”, colnames(birds))]) ,length.out=100),
site_type= c(“natural”, “park”, “restoration”),
tot_road_lengths_meters=mean(birds$ tot_road_lengths_meters))
det.predict<-predict(fmList, type=”state”, newdata=newdat, appendData=TRUE)
upperci<-det.predict$upper
lowerci<- det.predict$lower
meanvalues<-det.predict$Predicted
ymax<-1
ymin<-0
xmax<-ceiling(max(birds[,grep(“wood_area”, colnames(birds))], na.rm=TRUE)) #set maximum x value. change “time” to be the text of the variable you are plotting
xmin<-0 #min road was 1
yinterval<-0.2
xinterval<-20
plot(Predicted~ tot_native_wood_area_meters,data=det.predict[1:100,],
type=”n”,ann=F,axes=F,ylim=c(0,ymax)) #sets up the space to do the plot, without drawing the plot
axis(1,at=seq(xmin, xmax,xinterval) ,cex.axis=1.3) #add axis
axis(2,at=seq(ymin,ymax,yinterval),las=2,cex.axis=1.3) #add axis
box(bty=’l’) #what kind of box to draw around plot
title(xlab=”Total Native Woody Veg (ha)”,ylab=”Probability of detection”,cex.lab=1.5) #use title() function to add x and y axis labels. Change the text here to be what
you want you x label to be
points(Predicted~tot_native_wood_area_meters, data=det.predict[1:100,], type=’l’,lty=1, lwd=2) #draw on the predicted line
xs<- c(det.predict$tot_native_wood_area_meters [1:100], sort(det.predict$tot_native_wood_area_meters [1:100],decreasing=TRUE)) #it assumes that it joins first to last
point
ys<- c(upperci[1:100],lowerci[1:100][order(det.predict$tot_native_wood_area_meters [1:100],decreasing=TRUE)])
polygon(xs,ys, col=rgb(0,0,0,0.1), border=NA) #draw on the confidence bands as a polygon

# Save plot
savePlot(filename=”Occurrence_plot_wood.wmf”,type=”wmf”) #update graph name

# Landscape Ecology Code
# 2017

##### Prac 4 Bird Occupancy and Detection modelling
# Part 1. Revision
### Vectors, Matrices and Data Frames

# What is a vector?
# a sequence of data elements of the same basic type.
numbers <- c(1, 2 ,3)
letters <- c(“a”, “b”, “c”)
statements <- c(TRUE, FALSE, TRUE)

# What is a matrix?
# matrix(c(vector of numbers in column order),# of rows, # of columns)
our.first.matrix <- matrix(c(1, 2, 3, 4, 5, 6, 7, 8, 9, 10), 2, 5)

# What is a data frame?
# A bunch of vectors
# Just like a matrix except you can mix columns of different kinds of data. One can be numbers,
# while another can be text.

sex <- c(“M”, “F”, “F”, “M”)
thorax <- c(3.2, 15.4, 14.3, 4.1)
better.solution<- data.frame(sex, thorax)
rownames(better.solution) <- c(“Spider 1”, “Spider 2”, “Spider 3”, “Spider 4″)
better.solution

# what is a function?
# A function is a defined arguement. For instance you want the sum of a vector
# Like excel
numbers <- c(1, 2, 3)
sum(numbers)
help(sum)
# Subsetting data: e.g. select data from rows 2 and 3 (i.e. female spiders)
spider2_3 <- better.solution[2:3, ]
spider2_3

spiderF <- better.solution[2:3,]
spiderF

spiderF1 <- better.solution[-c(1,4),]
spiderF1

spiderFemale <- better.solution[better.solution$sex==”F”,]
spiderFemale

######
#
# Other helpful resources
#
##########

# google
# stackoverflow discussion boards
# github
# coursera
# codeacademy
##################################################################################################
# Part 2. Occupancy Modelling

install.packages(“unmarked”)

library(unmarked)

#############

setwd(“~/Desktop/Demo – Sem 1/Don/4) Assignment prac”) # Change for your working directory (e.g. your desktop)

birds <- read.csv(“siteandsurveycovars4analysis.csv”)

# Check that data were loaded correctly
dim(birds)
head(birds)
tail(birds)
names(birds)
summary(birds)
str(birds)

### Data ###

# Presence/absence of birds:
# Brown thornbill: columns 32:39
# Pied burrawong: columns 40:79
# Eastern yellow robin: columns 48:55

# 120 sites
# 10-11 surveys at each site

# Aim:
# 1) what variables affect the likeihood of detecting birds?
# 2) what varuables affect the likeihood that a site will be occupied?

# Observations affecting detection:
# 1) Weather: cloud, rain, sun
# 2) Time
# 3) Wind: light, moderate, strong

# Factors affecting occupancy:
# 1) Site type: largely natural, restored site, or manicured park
# 2) Exotic canopy cover (%)
# 3) Total road lengths (m) in 1 km buffer
# 4) Total native wood area
######### Example for Eastern Yellow Robin ##############

# UNMARKED ANALYSIS OF DETECTION AND OCCUPANCY (no covariates)

## Unmarked package ####

## Inputs – in general ###

# 1) Detection histories
# 2) Site-specific covariates: occupancy model
# 3) Sample-specific covariates: detection model

#create a simple unmarked data frame with just the occurrence data in it (i.e. columns 2, 3, 4)
frame1 <- unmarkedFrameOccu(y = birds[, 48:55])

## Simple occupancy model:
analysis1 <- occu (~1~1, frame1)
analysis1

# Backtransform
backTransform(analysis1,’det’) # detection estimates

backTransform(analysis1, ‘state’) # Occupancy estimates

# Data are binomial so they are automatically transformed to logit scale. To be able to interpret, transform
# back to an untransformed scale

##################################################

#UNMARKED ANALYSIS OF DETECTION AND OCCUPANCY WITH SITE AND OBSERVATION COVARIATES

# Original #
frame2 <- unmarkedFrameOccu(
y = birds[, 48:55], # define data/observations
siteCovs = as.data.frame(birds[,c(1:7)]), # set site covariates: in columns 7 and 8 (nearest occupied patch, spinifex)
obsCovs = list( # Set detection covariates: i.e. our survey covariates
weather = as.data.frame(birds[, 16:23]), # e.g. weather observations
wind = as.data.frame(birds[, 24:31]), # e.g. wind observations
time = as.data.frame(birds[, 8:15]))) # e.g. time of survey

### Model selection ###

## DETECTION COVARIATE MODELS – Use null model

model.null<-occu(~1~1, data = frame2) #no covariates fitted
model.1.1<-occu(~weather ~1, data=frame2)
model.1.2<-occu(~time ~1, data=frame2)
model.1.3<-occu(~wind ~1, data=frame2)
model.2.1<-occu(~weather+time ~1, data=frame2)
model.2.2<-occu(~weather+wind ~1, data=frame2)
model.2.3<-occu(~wind+time ~1, data=frame2)
model.3<-occu(~weather+wind+time ~1, data=frame2)

model.list<-fitList(null=model.null, model.1.1, model.1.2, model.1.3, model.2.1, model.2.2, model.2.3, model.3)
#fitList is the unmarked function to organise models for selection

#do the model ranking and model selection process
model.comparison<- as(modSel(model.list, nullmod = “null”), “data.frame”)
model.comparison

# Print to a csv file to open in excel etc.
write.csv(model.comparison, “model_selection_detection_robin.csv”) #the part in “” will be the file name in the directory you specified in setwd() or with the drop
down menu previously

##################################################

### OCCUPANCY AND DETECTION COVARIATE MODELS
# altering the occupancy part of the model
# For this assignment, we are only using main effects, no interactions between covariates

model.null<-occu(~weather ~1, data = frame2)
model.1.1<-occu(~weather ~site_type, data = frame2)
model.1.2<-occu(~weather ~pc_exoticcanopy, data = frame2)
model.1.3<-occu(~weather ~tot_road_lengths_meters, data = frame2)
model.1.4<-occu(~weather ~tot_native_wood_area_meters, data = frame2)
model.2.1<-occu(~weather ~site_type+pc_exoticcanopy, data = frame2)
model.2.2<-occu(~weather ~site_type+tot_road_lengths_meters, data = frame2)
model.2.3<-occu(~weather ~site_type+tot_native_wood_area_meters, data = frame2)
model.2.4<-occu(~weather ~pc_exoticcanopy+tot_road_lengths_meters, data = frame2)
model.2.5<-occu(~weather ~pc_exoticcanopy+tot_native_wood_area_meters, data = frame2)
model.2.6<-occu(~weather ~tot_road_lengths_meters+tot_native_wood_area_meters, data = frame2)
model.3.1<-occu(~weather ~site_type+pc_exoticcanopy+tot_road_lengths_meters, data = frame2)
model.3.2<-occu(~weather ~site_type+pc_exoticcanopy+tot_native_wood_area_meters, data = frame2)
model.3.3<-occu(~weather ~site_type+tot_road_lengths_meters+tot_native_wood_area_meters, data = frame2)
model.3.4<-occu(~weather ~pc_exoticcanopy+tot_road_lengths_meters+tot_native_wood_area_meters, data = frame2)
model.4.1<-occu(~weather ~site_type+ pc_exoticcanopy+tot_road_lengths_meters+tot_native_wood_area_meters, data = frame2)

#Put the models we can fit into the list
model.comps<- fitList(null = model.null, model.1.1, model.1.2, model.1.3, model.1.4,
model.2.1, model.2.2, model.2.3, model.2.4, model.2.5, model.2.6,
model.3.1, model.3.2, model.3.3, model.3.4, model.4.1)

#do the model ranking and model selection process
model.comps2 <- as(modSel(model.comps, nullmod = “null”), “data.frame”)

AICc <- model.comps2$AIC * (model.comps2$n / (model.comps2$n – model.comps2$nPars-1))
delta.small<- AICc – min(AICc)
likli.small<- exp(-0.5 * delta.small)
AICwt.small <- likli.small / (sum(likli.small))

model.comps3 <- data.frame(model.comps2, AICc, delta.small, likli.small, AICwt.small)[order(AICc),] #rank results by AICc

#ditch the model estimates from this table because we will extract them in a later step.
model.comps3a <- model.comps3[,c(1:2, which(colnames(model.comps3) == “negLogLike”) : length(model.comps3))]
model.comps3a

#save our results to csv for later.
write.csv(model.comps3a, “model_selection_bird_surveys_robin.csv”)

######################################

#### PLOT THE RESULTS

fmList <- fitList(model.1.3) # change to your best model
##### Plot 1: # How was detection influenced by weather?

newdat<-expand.grid(
weather= c(“cloud”, “rain”, “sun”))
det.predict<-predict(fmList, type=”det”, newdata=newdat, appendData=TRUE)
upperci<-det.predict$upper
lowerci<- det.predict$lower
meanvalues<-det.predict$Predicted
ylabel<-“Probability of detection”
maxy<- 1.05
barnames<- c(“cloud”, “rain”, “sun”) #update levels to be the ones for the factor you are plotting. in this case, it is set up for weather. but you might also have
wind, which has levels “light”, “moderate”, “strong”, or other categorical covariates.
par(mar= c(10, 6, 2, 2) + 0.1) # this sets parameters for your plot. no need to change them most likely c(bottom, left, top, right) number of lines of margin; adjust
bottom to suit length of barnames
mp<- barplot(meanvalues, space = 0.5,
names.arg = barnames, beside = TRUE,
horiz = FALSE, prcol = NULL, border = par(“fg”),
ylim = c(0, maxy),ylab = ylabel, xpd = FALSE, las=1,lwd=2,xaxt=”s”,
cex.axis = 2.0, cex.names = 2.0, cex.lab=2.0, mgp=c(3.5,1,0))

arrows(mp, lowerci ,mp, upperci, code=3, angle=90, length=.1, lwd=2)
abline(0,0,lwd=3)

#Save plot:
savePlot(filename=”Detection_plot_weather.wmf”,type=”wmf”) #rename as you prefer
### Plot 2: How occupancy is influenced by road length?

newdat<-expand.grid(
tot_road_lengths_meters = seq(min(birds[,grep(“road”, colnames(birds))], na.rm=TRUE), max(birds[,grep(“road”, colnames(birds))]) ,length.out=100),
site_type= c(“natural”, “park”, “restoration”),
tot_native_wood_area_meters=mean(birds$tot_native_wood_area_meters))
det.predict<-predict(fmList, type=”state”, newdata=newdat, appendData=TRUE)
upperci<-det.predict$upper
lowerci<- det.predict$lower
meanvalues<-det.predict$Predicted
ymax<-1
ymin<-0
xmax<-ceiling(max(birds[,grep(“road”, colnames(birds))], na.rm=TRUE)) #set maximum x value. change “time” to be the text of the variable you are plotting
xmin<-0 #min road was 1
yinterval<-0.2
xinterval<-20
plot(Predicted~ tot_road_lengths_meters,data=det.predict[1:100,],
type=”n”,ann=F,axes=F,ylim=c(0,ymax)) #sets up the space to do the plot, without drawing the plot
axis(1,at=seq(xmin, xmax,xinterval) ,cex.axis=1.3) #add axis
axis(2,at=seq(ymin,ymax,yinterval),las=2,cex.axis=1.3) #add axis
box(bty=’l’) #what kind of box to draw around plot
title(xlab=”Road length (km)”,ylab=”Probability of occupancy”,cex.lab=1.5) #use title() function to add x and y axis labels. Change the text here to be what you want
you x label to be
points(Predicted~ tot_road_lengths_meters, data=det.predict[1:100,], type=’l’,lty=1, lwd=2) #draw on the predicted line
xs<- c(det.predict$tot_road_lengths_meters[1:100], sort(det.predict$tot_road_lengths_meters [1:100],decreasing=TRUE)) #it assumes that it joins first to last point
ys<- c(upperci[1:100],lowerci[1:100][order(det.predict$tot_road_lengths_meters [1:100],decreasing=TRUE)])
polygon(xs,ys, col=rgb(0,0,0,0.1), border=NA) #draw on the confidence bands as a polygon

# Save plot
savePlot(filename=”Occurrence_plot_roads.wmf”,type=”wmf”) #update graph name
################################################################################

##### Plots for other species ######

fmList <- fitList(model.1.3) # change to your best model

### Plot for occupancy: Site Type

newdat<-expand.grid(
site_type= c(“natural”, “park”, “restoration”),
tot_road_lengths_meters=(mean(birds$tot_road_lengths_meters)),
tot_native_wood_area_meters=mean(birds$tot_native_wood_area_meters))
det.predict<-predict(fmList, type=”state”, newdata=newdat, appendData=TRUE)
upperci<-det.predict$upper
lowerci<- det.predict$lower
meanvalues<-det.predict$Predicted
ylabel<-“Proportion of sites occupied”
maxy<- 1.05
barnames<- c(“natural”, “park”, “restoration”) #update levels to be the ones for the factor you are plotting. in this case, it is set up for site_type.
par(mar= c(10, 6, 2, 2) + 0.1) # this sets parameters for your plot. no need to change them most likely c(bottom, left, top, right) number of lines of margin; adjust
bottom to suit length of barnames
mp<- barplot(meanvalues, space = 0.5,
names.arg = barnames, beside = TRUE,
horiz = FALSE, prcol = NULL, border = par(“fg”),
ylim = c(0, maxy),ylab = ylabel, xpd = FALSE, las=1,lwd=2,xaxt=”s”,
cex.axis = 2.0, cex.names = 2.0, cex.lab=2.0, mgp=c(3.5,1,0))
arrows(mp, lowerci ,mp, upperci, code=3, angle=90, length=.1, lwd=2)
abline(0,0,lwd=3)

# Save plot
savePlot(filename=”Occurence_plot_site_type.wmf”,type=”wmf”) #rename as you prefer
### Plot for occupancy: Native wood area

newdat<-expand.grid(
tot_native_wood_area_meters = seq(min(birds[,grep(“wood_area”, colnames(birds))], na.rm=TRUE), max(birds[,grep(“wood_area”, colnames(birds))]) ,length.out=100),
site_type= c(“natural”, “park”, “restoration”),
tot_road_lengths_meters=mean(birds$ tot_road_lengths_meters))
det.predict<-predict(fmList, type=”state”, newdata=newdat, appendData=TRUE)
upperci<-det.predict$upper
lowerci<- det.predict$lower
meanvalues<-det.predict$Predicted
ymax<-1
ymin<-0
xmax<-ceiling(max(birds[,grep(“wood_area”, colnames(birds))], na.rm=TRUE)) #set maximum x value. change “time” to be the text of the variable you are plotting
xmin<-0 #min road was 1
yinterval<-0.2
xinterval<-20
plot(Predicted~ tot_native_wood_area_meters,data=det.predict[1:100,],
type=”n”,ann=F,axes=F,ylim=c(0,ymax)) #sets up the space to do the plot, without drawing the plot
axis(1,at=seq(xmin, xmax,xinterval) ,cex.axis=1.3) #add axis
axis(2,at=seq(ymin,ymax,yinterval),las=2,cex.axis=1.3) #add axis
box(bty=’l’) #what kind of box to draw around plot
title(xlab=”Total Native Woody Veg (ha)”,ylab=”Probability of detection”,cex.lab=1.5) #use title() function to add x and y axis labels. Change the text here to be what
you want you x label to be
points(Predicted~tot_native_wood_area_meters, data=det.predict[1:100,], type=’l’,lty=1, lwd=2) #draw on the predicted line
xs<- c(det.predict$tot_native_wood_area_meters [1:100], sort(det.predict$tot_native_wood_area_meters [1:100],decreasing=TRUE)) #it assumes that it joins first to last
point
ys<- c(upperci[1:100],lowerci[1:100][order(det.predict$tot_native_wood_area_meters [1:100],decreasing=TRUE)])
polygon(xs,ys, col=rgb(0,0,0,0.1), border=NA) #draw on the confidence bands as a polygon

# Save plot
savePlot(filename=”Occurrence_plot_wood.wmf”,type=”wmf”) #update graph name

# Landscape Ecology Code
# 2017

##### Prac 4 Bird Occupancy and Detection modelling
# Part 1. Revision
### Vectors, Matrices and Data Frames

# What is a vector?
# a sequence of data elements of the same basic type.
numbers <- c(1, 2 ,3)
letters <- c(“a”, “b”, “c”)
statements <- c(TRUE, FALSE, TRUE)

# What is a matrix?
# matrix(c(vector of numbers in column order),# of rows, # of columns)
our.first.matrix <- matrix(c(1, 2, 3, 4, 5, 6, 7, 8, 9, 10), 2, 5)

# What is a data frame?
# A bunch of vectors
# Just like a matrix except you can mix columns of different kinds of data. One can be numbers,
# while another can be text.

sex <- c(“M”, “F”, “F”, “M”)
thorax <- c(3.2, 15.4, 14.3, 4.1)
better.solution<- data.frame(sex, thorax)
rownames(better.solution) <- c(“Spider 1”, “Spider 2”, “Spider 3”, “Spider 4″)
better.solution

# what is a function?
# A function is a defined arguement. For instance you want the sum of a vector
# Like excel
numbers <- c(1, 2, 3)
sum(numbers)
help(sum)
# Subsetting data: e.g. select data from rows 2 and 3 (i.e. female spiders)
spider2_3 <- better.solution[2:3, ]
spider2_3

spiderF <- better.solution[2:3,]
spiderF

spiderF1 <- better.solution[-c(1,4),]
spiderF1

spiderFemale <- better.solution[better.solution$sex==”F”,]
spiderFemale

######
#
# Other helpful resources
#
##########

# google
# stackoverflow discussion boards
# github
# coursera
# codeacademy
##################################################################################################
# Part 2. Occupancy Modelling

install.packages(“unmarked”)

library(unmarked)

#############

setwd(“~/Desktop/Demo – Sem 1/Don/4) Assignment prac”) # Change for your working directory (e.g. your desktop)

birds <- read.csv(“siteandsurveycovars4analysis.csv”)

# Check that data were loaded correctly
dim(birds)
head(birds)
tail(birds)
names(birds)
summary(birds)
str(birds)

### Data ###

# Presence/absence of birds:
# Brown thornbill: columns 32:39
# Pied burrawong: columns 40:79
# Eastern yellow robin: columns 48:55

# 120 sites
# 10-11 surveys at each site

# Aim:
# 1) what variables affect the likeihood of detecting birds?
# 2) what varuables affect the likeihood that a site will be occupied?

# Observations affecting detection:
# 1) Weather: cloud, rain, sun
# 2) Time
# 3) Wind: light, moderate, strong

# Factors affecting occupancy:
# 1) Site type: largely natural, restored site, or manicured park
# 2) Exotic canopy cover (%)
# 3) Total road lengths (m) in 1 km buffer
# 4) Total native wood area
######### Example for Eastern Yellow Robin ##############

# UNMARKED ANALYSIS OF DETECTION AND OCCUPANCY (no covariates)

## Unmarked package ####

## Inputs – in general ###

# 1) Detection histories
# 2) Site-specific covariates: occupancy model
# 3) Sample-specific covariates: detection model

#create a simple unmarked data frame with just the occurrence data in it (i.e. columns 2, 3, 4)
frame1 <- unmarkedFrameOccu(y = birds[, 48:55])

## Simple occupancy model:
analysis1 <- occu (~1~1, frame1)
analysis1

# Backtransform
backTransform(analysis1,’det’) # detection estimates

backTransform(analysis1, ‘state’) # Occupancy estimates

# Data are binomial so they are automatically transformed to logit scale. To be able to interpret, transform
# back to an untransformed scale

##################################################

#UNMARKED ANALYSIS OF DETECTION AND OCCUPANCY WITH SITE AND OBSERVATION COVARIATES

# Original #
frame2 <- unmarkedFrameOccu(
y = birds[, 48:55], # define data/observations
siteCovs = as.data.frame(birds[,c(1:7)]), # set site covariates: in columns 7 and 8 (nearest occupied patch, spinifex)
obsCovs = list( # Set detection covariates: i.e. our survey covariates
weather = as.data.frame(birds[, 16:23]), # e.g. weather observations
wind = as.data.frame(birds[, 24:31]), # e.g. wind observations
time = as.data.frame(birds[, 8:15]))) # e.g. time of survey

### Model selection ###

## DETECTION COVARIATE MODELS – Use null model

model.null<-occu(~1~1, data = frame2) #no covariates fitted
model.1.1<-occu(~weather ~1, data=frame2)
model.1.2<-occu(~time ~1, data=frame2)
model.1.3<-occu(~wind ~1, data=frame2)
model.2.1<-occu(~weather+time ~1, data=frame2)
model.2.2<-occu(~weather+wind ~1, data=frame2)
model.2.3<-occu(~wind+time ~1, data=frame2)
model.3<-occu(~weather+wind+time ~1, data=frame2)

model.list<-fitList(null=model.null, model.1.1, model.1.2, model.1.3, model.2.1, model.2.2, model.2.3, model.3)
#fitList is the unmarked function to organise models for selection

#do the model ranking and model selection process
model.comparison<- as(modSel(model.list, nullmod = “null”), “data.frame”)
model.comparison

# Print to a csv file to open in excel etc.
write.csv(model.comparison, “model_selection_detection_robin.csv”) #the part in “” will be the file name in the directory you specified in setwd() or with the drop
down menu previously

##################################################

### OCCUPANCY AND DETECTION COVARIATE MODELS
# altering the occupancy part of the model
# For this assignment, we are only using main effects, no interactions between covariates

model.null<-occu(~weather ~1, data = frame2)
model.1.1<-occu(~weather ~site_type, data = frame2)
model.1.2<-occu(~weather ~pc_exoticcanopy, data = frame2)
model.1.3<-occu(~weather ~tot_road_lengths_meters, data = frame2)
model.1.4<-occu(~weather ~tot_native_wood_area_meters, data = frame2)
model.2.1<-occu(~weather ~site_type+pc_exoticcanopy, data = frame2)
model.2.2<-occu(~weather ~site_type+tot_road_lengths_meters, data = frame2)
model.2.3<-occu(~weather ~site_type+tot_native_wood_area_meters, data = frame2)
model.2.4<-occu(~weather ~pc_exoticcanopy+tot_road_lengths_meters, data = frame2)
model.2.5<-occu(~weather ~pc_exoticcanopy+tot_native_wood_area_meters, data = frame2)
model.2.6<-occu(~weather ~tot_road_lengths_meters+tot_native_wood_area_meters, data = frame2)
model.3.1<-occu(~weather ~site_type+pc_exoticcanopy+tot_road_lengths_meters, data = frame2)
model.3.2<-occu(~weather ~site_type+pc_exoticcanopy+tot_native_wood_area_meters, data = frame2)
model.3.3<-occu(~weather ~site_type+tot_road_lengths_meters+tot_native_wood_area_meters, data = frame2)
model.3.4<-occu(~weather ~pc_exoticcanopy+tot_road_lengths_meters+tot_native_wood_area_meters, data = frame2)
model.4.1<-occu(~weather ~site_type+ pc_exoticcanopy+tot_road_lengths_meters+tot_native_wood_area_meters, data = frame2)

#Put the models we can fit into the list
model.comps<- fitList(null = model.null, model.1.1, model.1.2, model.1.3, model.1.4,
model.2.1, model.2.2, model.2.3, model.2.4, model.2.5, model.2.6,
model.3.1, model.3.2, model.3.3, model.3.4, model.4.1)

#do the model ranking and model selection process
model.comps2 <- as(modSel(model.comps, nullmod = “null”), “data.frame”)

AICc <- model.comps2$AIC * (model.comps2$n / (model.comps2$n – model.comps2$nPars-1))
delta.small<- AICc – min(AICc)
likli.small<- exp(-0.5 * delta.small)
AICwt.small <- likli.small / (sum(likli.small))

model.comps3 <- data.frame(model.comps2, AICc, delta.small, likli.small, AICwt.small)[order(AICc),] #rank results by AICc

#ditch the model estimates from this table because we will extract them in a later step.
model.comps3a <- model.comps3[,c(1:2, which(colnames(model.comps3) == “negLogLike”) : length(model.comps3))]
model.comps3a

#save our results to csv for later.
write.csv(model.comps3a, “model_selection_bird_surveys_robin.csv”)

######################################

#### PLOT THE RESULTS

fmList <- fitList(model.1.3) # change to your best model
##### Plot 1: # How was detection influenced by weather?

newdat<-expand.grid(
weather= c(“cloud”, “rain”, “sun”))
det.predict<-predict(fmList, type=”det”, newdata=newdat, appendData=TRUE)
upperci<-det.predict$upper
lowerci<- det.predict$lower
meanvalues<-det.predict$Predicted
ylabel<-“Probability of detection”
maxy<- 1.05
barnames<- c(“cloud”, “rain”, “sun”) #update levels to be the ones for the factor you are plotting. in this case, it is set up for weather. but you might also have
wind, which has levels “light”, “moderate”, “strong”, or other categorical covariates.
par(mar= c(10, 6, 2, 2) + 0.1) # this sets parameters for your plot. no need to change them most likely c(bottom, left, top, right) number of lines of margin; adjust
bottom to suit length of barnames
mp<- barplot(meanvalues, space = 0.5,
names.arg = barnames, beside = TRUE,
horiz = FALSE, prcol = NULL, border = par(“fg”),
ylim = c(0, maxy),ylab = ylabel, xpd = FALSE, las=1,lwd=2,xaxt=”s”,
cex.axis = 2.0, cex.names = 2.0, cex.lab=2.0, mgp=c(3.5,1,0))

arrows(mp, lowerci ,mp, upperci, code=3, angle=90, length=.1, lwd=2)
abline(0,0,lwd=3)

#Save plot:
savePlot(filename=”Detection_plot_weather.wmf”,type=”wmf”) #rename as you prefer
### Plot 2: How occupancy is influenced by road length?

newdat<-expand.grid(
tot_road_lengths_meters = seq(min(birds[,grep(“road”, colnames(birds))], na.rm=TRUE), max(birds[,grep(“road”, colnames(birds))]) ,length.out=100),
site_type= c(“natural”, “park”, “restoration”),
tot_native_wood_area_meters=mean(birds$tot_native_wood_area_meters))
det.predict<-predict(fmList, type=”state”, newdata=newdat, appendData=TRUE)
upperci<-det.predict$upper
lowerci<- det.predict$lower
meanvalues<-det.predict$Predicted
ymax<-1
ymin<-0
xmax<-ceiling(max(birds[,grep(“road”, colnames(birds))], na.rm=TRUE)) #set maximum x value. change “time” to be the text of the variable you are plotting
xmin<-0 #min road was 1
yinterval<-0.2
xinterval<-20
plot(Predicted~ tot_road_lengths_meters,data=det.predict[1:100,],
type=”n”,ann=F,axes=F,ylim=c(0,ymax)) #sets up the space to do the plot, without drawing the plot
axis(1,at=seq(xmin, xmax,xinterval) ,cex.axis=1.3) #add axis
axis(2,at=seq(ymin,ymax,yinterval),las=2,cex.axis=1.3) #add axis
box(bty=’l’) #what kind of box to draw around plot
title(xlab=”Road length (km)”,ylab=”Probability of occupancy”,cex.lab=1.5) #use title() function to add x and y axis labels. Change the text here to be what you want
you x label to be
points(Predicted~ tot_road_lengths_meters, data=det.predict[1:100,], type=’l’,lty=1, lwd=2) #draw on the predicted line
xs<- c(det.predict$tot_road_lengths_meters[1:100], sort(det.predict$tot_road_lengths_meters [1:100],decreasing=TRUE)) #it assumes that it joins first to last point
ys<- c(upperci[1:100],lowerci[1:100][order(det.predict$tot_road_lengths_meters [1:100],decreasing=TRUE)])
polygon(xs,ys, col=rgb(0,0,0,0.1), border=NA) #draw on the confidence bands as a polygon

# Save plot
savePlot(filename=”Occurrence_plot_roads.wmf”,type=”wmf”) #update graph name
################################################################################

##### Plots for other species ######

fmList <- fitList(model.1.3) # change to your best model

### Plot for occupancy: Site Type

newdat<-expand.grid(
site_type= c(“natural”, “park”, “restoration”),
tot_road_lengths_meters=(mean(birds$tot_road_lengths_meters)),
tot_native_wood_area_meters=mean(birds$tot_native_wood_area_meters))
det.predict<-predict(fmList, type=”state”, newdata=newdat, appendData=TRUE)
upperci<-det.predict$upper
lowerci<- det.predict$lower
meanvalues<-det.predict$Predicted
ylabel<-“Proportion of sites occupied”
maxy<- 1.05
barnames<- c(“natural”, “park”, “restoration”) #update levels to be the ones for the factor you are plotting. in this case, it is set up for site_type.
par(mar= c(10, 6, 2, 2) + 0.1) # this sets parameters for your plot. no need to change them most likely c(bottom, left, top, right) number of lines of margin; adjust
bottom to suit length of barnames
mp<- barplot(meanvalues, space = 0.5,
names.arg = barnames, beside = TRUE,
horiz = FALSE, prcol = NULL, border = par(“fg”),
ylim = c(0, maxy),ylab = ylabel, xpd = FALSE, las=1,lwd=2,xaxt=”s”,
cex.axis = 2.0, cex.names = 2.0, cex.lab=2.0, mgp=c(3.5,1,0))
arrows(mp, lowerci ,mp, upperci, code=3, angle=90, length=.1, lwd=2)
abline(0,0,lwd=3)

# Save plot
savePlot(filename=”Occurence_plot_site_type.wmf”,type=”wmf”) #rename as you prefer
### Plot for occupancy: Native wood area

newdat<-expand.grid(
tot_native_wood_area_meters = seq(min(birds[,grep(“wood_area”, colnames(birds))], na.rm=TRUE), max(birds[,grep(“wood_area”, colnames(birds))]) ,length.out=100),
site_type= c(“natural”, “park”, “restoration”),
tot_road_lengths_meters=mean(birds$ tot_road_lengths_meters))
det.predict<-predict(fmList, type=”state”, newdata=newdat, appendData=TRUE)
upperci<-det.predict$upper
lowerci<- det.predict$lower
meanvalues<-det.predict$Predicted
ymax<-1
ymin<-0
xmax<-ceiling(max(birds[,grep(“wood_area”, colnames(birds))], na.rm=TRUE)) #set maximum x value. change “time” to be the text of the variable you are plotting
xmin<-0 #min road was 1
yinterval<-0.2
xinterval<-20
plot(Predicted~ tot_native_wood_area_meters,data=det.predict[1:100,],
type=”n”,ann=F,axes=F,ylim=c(0,ymax)) #sets up the space to do the plot, without drawing the plot
axis(1,at=seq(xmin, xmax,xinterval) ,cex.axis=1.3) #add axis
axis(2,at=seq(ymin,ymax,yinterval),las=2,cex.axis=1.3) #add axis
box(bty=’l’) #what kind of box to draw around plot
title(xlab=”Total Native Woody Veg (ha)”,ylab=”Probability of detection”,cex.lab=1.5) #use title() function to add x and y axis labels. Change the text here to be what
you want you x label to be
points(Predicted~tot_native_wood_area_meters, data=det.predict[1:100,], type=’l’,lty=1, lwd=2) #draw on the predicted line
xs<- c(det.predict$tot_native_wood_area_meters [1:100], sort(det.predict$tot_native_wood_area_meters [1:100],decreasing=TRUE)) #it assumes that it joins first to last
point
ys<- c(upperci[1:100],lowerci[1:100][order(det.predict$tot_native_wood_area_meters [1:100],decreasing=TRUE)])
polygon(xs,ys, col=rgb(0,0,0,0.1), border=NA) #draw on the confidence bands as a polygon

# Save plot
savePlot(filename=”Occurrence_plot_wood.wmf”,type=”wmf”) #update graph name

find the cost of your paper