Riverbird Habitat Use Data Analysis

John McLaughlin

1 March 2020



Commands and code for data analysis using the R software environment.

This document contains raw R commands, suitable for copying and pasting directly into R.

Analysis below is almost entirely automated.

Users need enter their own numbers in only three places:

(1) Number of riverbird detections in each habitat type:

     e.g., if 12 birds detected at forested riverbanks, line below would be:  fu <- 12

(2) Number of each habitat type detected during 5-minute interval habitat survey.

     e.g., if 35 points had forested banks, line below would be: fa <- 35

(3) Standard errors for each habitat type coefficient. These values are given

     in each model summary, generated by commands below. 

     If standard errors are not entered, commands below use zero as default values.

Commands below generate numerical results, and plot resource selection function values

  for each habitat type -- including standard errors.

This analysis should be repeated for each riverbird species.

______________________________________________________________________________



# This is a comment. R ignores anything that follows the pound sign (#).

# To learn more about R, including how to download your own copy of R, go to

#    http://cran.R-project.org



# Below are tallies of riverbird detections in each riverbank habitat type,

#    followed by tallies of habitat types recorded at 5-minute intervals along river. 

# Repeat data entry and analysis for each riverbird species.



# Enter number of riverbird detections in each habitat 

#    (First letter denotes habitat type; second letter ("u") stands for use.)

bu <-     # enter no. riverbird detections on bare substrate here

hu <-     # enter no. riverbird detections in herbaceous vegetation here

su <-     # enter no. detections at/in shrubs here

fu <-     # enter no. detections at forest here

lu <-     # enter no. detections at Large woody debris (LWD) here



# Enter number of each habitat type recorded at 5-minute intervals = habitat "available"

#    (First letter denotes habitat type; second letter ("a") stands for "available.")

ba <-     # enter no. points w/ bare substrate here

ha <-     # enter no. points w/ herbaceous vegetation here

sa <-     # enter no. points w/ shrubby banks here

fa <-     # enter no. points w/ forested banks here

la <-     # enter no. points w/ LWD here



# Now construct response (riverbird use) and predictor (habitat) variables 

#     for all sampled locations -- both riverbird and habitat sampling.

# For riverbird use, a value of 1 indicates riverbird detected there; 0 indicates non-detection.



use.spp <- c(rep(1,bu), rep(0,ba), rep(1,hu), rep(0,ha), rep(1,su), rep(0,sa), rep(1,fu), rep(0,fa), rep(1,lu), rep(0,la))



# For habitat variables, 1 indicates location was that habitat type; 0 indicates a different habitat.

# Each of those bare habitat locations are assigned value "1" for the "bare" category.

# Assign 1 and 0 to remaining habitat categories in similar way. 

# Variable assignments below include 10 riverbird-habitat combinations, to facilitate error-checking.

# Assignments below use R's repeat function, rep(), to simplify data entry.



bare <- c(rep(1,bu), rep(1,ba), rep(0,hu), rep(0,ha), rep(0,su), rep(0,sa), rep(0,fu), rep(0,fa), rep(0,lu), rep(0,la))

herb <- c(rep(0,bu), rep(0,ba), rep(1,hu), rep(1,ha), rep(0,su), rep(0,sa), rep(0,fu), rep(0,fa), rep(0,lu), rep(0,la))

shrub <- c(rep(0,bu), rep(0,ba), rep(0,hu), rep(0,ha), rep(1,su), rep(1,sa), rep(0,fu), rep(0,fa), rep(0,lu), rep(0,la))

forest <- c(rep(0,bu), rep(0,ba), rep(0,hu), rep(0,ha), rep(0,su), rep(0,sa), rep(1,fu), rep(1,fa), rep(0,lu), rep(0,la))

lwd <- c(rep(0,bu), rep(0,ba), rep(0,hu), rep(0,ha), rep(0,su), rep(0,sa), rep(0,fu), rep(0,fa), rep(1,lu), rep(1,la))



# Check that Riverbird variable and all habitat variables contain the same number of values.

length(use.spp)

length(bare)

length(herb)

length(shrub)

length(forest)

length(lwd)



# ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~

# Fit Generalized Linear Models, glm(), for each categorical habitat variable



spp.bare <- glm(use.spp ~ bare, family=binomial)

summary(spp.bare)



spp.herb <- glm(use.spp ~ herb, family=binomial)

summary(spp.herb)



spp.shrub <- glm(use.spp ~ shrub, family=binomial)

summary(spp.shrub)



spp.for <- glm(use.spp ~ forest, family=binomial)

summary(spp.for)



spp.lwd <- glm(use.spp ~ lwd, family=binomial)

summary(spp.lwd)



# ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~



# Calculate resource selection function (RSF) values for each habitat 

# RSF values for each habitat are glm coefficients, exponentiated (= e^coef)



spp.coef <- c(spp.bare$coef[2], spp.herb$coef[2], spp.shrub$coef[2], 

    spp.for$coef[2], spp.lwd$coef[2])



# RSF std errors

spp.err <- c(0, 0, 0, 0, 0)    # Replace zeros with coefficient std.error values

#    Standard error values for each habitat type are given in model summaries above



# Function to plot RSF values and standard errors for each habitat type.

# Function adapted from Crawley MJ (2007) The R Book. Wiley. p.56.



error.bar5 <- function(yv, z, nn, yx){

# yv=bar heights, z=error bar heights, nn=labels for each bar, yx=y axis label

xv <- 

barplot(exp(yv), ylim=c(0,(max(exp(yv+z)))),names=nn, ylab=yx)

g=(max(xv)-min(xv))/50

for (i in 1:length(xv)) {

  lines(c(xv[i], xv[i]), c(exp(yv[i]+z[i]),exp(yv[i]-z[i]))) 

  lines(c(xv[i]-g, xv[i]+g), c(exp(yv[i]+z[i]),exp(yv[i]+z[i]))) 

  lines(c(xv[i]-g, xv[i]+g), c(exp(yv[i]-z[i]),exp(yv[i]-z[i])))  

  }

}



spp.habitat <- c("Bare", "Herbaceous", "Shrub", "Forest", "LWD") # x-axis labels

spp.label <- "Riverbird RSF"    # y-axis label: replace "Riverbird" with species name



# Plot RSF values and SE



error.bar5(spp.coef, spp.err, spp.habitat, spp.label)

# Add dashed horizontal line at y=1, or no selection; use in same proportion as availability

lines(c(0,18), c(1,1), lty=2)