Land-cover use

table(hr$EOSD.class, hr$Used)

##                      
##                       FALSE TRUE
##   Broadleaf-dense        10    3
##   Broadleaf-open          5    6
##   Bryoids                 2    3
##   Coniferous-dense      142 1373
##   Coniferous-open       235 3174
##   Coniferous-sparse      58  460
##   Developed               2    0
##   Exposed/Barren Land     1    0
##   Herbs                   0    0
##   Mixedwood-dense        48  164
##   Mixedwood-open         12   72
##   Recent Burn            64  650
##   Rock/Rubble             0    0
##   Shadow                  0    0
##   Shrub Low               8   32
##   Shrub Tall              5    0
##   Water                  91   91
##   Wetland-herb           33  682
##   Wetland-shrub          76 1404
##   Wetland-treed          53  724

par(mar = c(3,10,2,2), cex.axis = 0.8)
barplot(t(table(hr$EOSD.class, hr$Used)), 
        beside = TRUE, horiz = TRUE, las = 1)

counts.table <- table(hr$EOSD.class, hr$Used)
counts.df <- as.data.frame.matrix(counts.table, 
                                  row.names = row.names(counts.table)) 
names(counts.df) <- c("Available", "Used")
counts.df$Total <- rowSums(counts.df)
counts.df$EOSD <- row.names(counts.df)

counts.df <- mutate(counts.df, 
                    Used = Used/sum(Used),
                                        Available = Available/sum(Available))   %>% 
                                        plyr::arrange(Available) 

par(mar = c(4,15,2,2), cex = 0.8)
cols <- c("grey", "darkgreen")
barplot(t(counts.df[,c("Available", "Used")]), 
                    beside = TRUE, col = cols, 
          horiz = TRUE, las = 1, bor = NA,  
                    names.arg = counts.df$EOSD)
legend("bottomright", fill = cols, legend = c("available", "used"))

Relative use plots

Ratios of use versus available counts can also be useful, e.g.:

counts.df <- mutate(counts.df, rel.preference = Used / Available) %>% subset(!is.na(rel.preference) & Used > 0) %>% arrange(rel.preference)
par(mar = c(4,10,2,2), cex = 0.8, las = 1)
with(counts.df, barplot(rel.preference, width = Available, horiz = TRUE, names.arg = EOSD))
abline(v = 1, col = 2, lwd = 2, lty = 3)

The widths of the bars are proportional to the availability in the landscape, those that are greater than 1 indicate greater use than expected (i.e. possible preference), those that are less than one indicate - possible - avoidance. The RSF will help us quantify these results.

Note, with a little but of fussing, we can actually match these bars to the colors of the original EOSD map:

require(raster); eosd <- raster("_data/rasters/EOSD_withfire")
counts.df <- merge(counts.df, levels(eosd)[[1]][,c("Description","Color")], by.x = "EOSD", by.y = "Description")
par(mar = c(4,10,2,2), cex = 0.8, las = 1)
with(counts.df, 
     barplot(rel.preference, width = Available, 
             horiz = TRUE, names.arg = EOSD,
             col = Color))
abline(v = 1, col = 2, lwd = 2, lty = 3)

Continuous covariates

Boxplot are the most obvious way to look for differences:

boxplot(DEM ~ Used, data = hr, main = "Elevation (m)")
boxplot(Road.Density ~ Used, data = hr, "Road Density (m/m^2)")

But it’s hard to see any really compelling differences.

Interesting - always - to look at differences against individuals:

boxplot(DEM ~ nick_name, data = subset(hr, Used))
boxplot(Road.Density ~ nick_name, data = subset(hr, Used))

I often like looking at density plots or violinplots to compare continuous covariates across use.

par(mfrow = c(1,2), bty = "l")
with(hr, {
    plot(density(DEM[!Used]), xlim = c(0,800), main = "elevation distrubution", 
         xlab = "DEM", col = "grey", lwd = 2, ylim = c(0,0.012))
    lines(density(DEM[Used], na.rm = TRUE), col = 2, lwd = 2)
    abline(v = tapply(DEM, Used, median), col = c("grey", "red"), lty = 3)

    plot(density(sqrt(Road.Density[Used]), na.rm = TRUE), main = "road densities", 
         xlab = "sqrt(Road density)", col = "grey", lwd = 2)
    lines(density(sqrt(Road.Density[!Used]), na.rm = TRUE), col = 2, lwd = 2)
    abline(v = tapply(Road.Density, Used, median), col = c("grey", "red"), lty = 3)

    
    legend("topright", col = c("red", "grey"), legend = c("used", "available"), lty = 1, lwd = 2)
    })

Note, road densities are extremely skewed, but also contain a lot of zeros, so we take the square root (rather than log) just to fit them on the figure.

There are few very obvious difference (e.g. in medians - the dashed lines above) between the two distributions. Here are some quick summary statistics:

ddply(hr, "Used", summarize, DEM = t(quantile(DEM))) 

##    Used   DEM.0%  DEM.25%  DEM.50%  DEM.75% DEM.100%
## 1 FALSE 146.0664 239.6714 291.6505 323.6667 886.6441
## 2  TRUE 155.2552 228.1318 280.1572 313.0000 813.1412

ddply(hr, "Used", summarize, DRoad = t(quantile(Road.Density)))

##    Used   DRoad.0%  DRoad.25%  DRoad.50%  DRoad.75% DRoad.100%
## 1 FALSE 0.00000000 0.00000000 0.00000000 0.17517307 3.35663557
## 2  TRUE 0.00000000 0.00000000 0.02342629 0.17865342 2.04081440

It is important also to make sure that our covariates aren’t too correlated:

plot(hr$DEM, jitter(sqrt(hr$Road.Density), 200), 
     col = rgb(hr$Used,0,0,.2), pch = 19, cex = 0.5)

Red points, in the above, are “used” locations. I added the jitter to see all the 0 density points. But in any case, correlation shouldn’t be too much of a problem

Seasonal / individual breakdowns

Here are some fun, quick visualizations of the distribution of caribou use vs. availability across seasons for our two continuous covariates.

Elevation:

ggplot(hr, aes(Used, DEM, fill = factor(month))) + geom_violin(draw_quantiles = c(.25,.5,.75)) 

There appear to be some seasonal differences, but also a sharp split in the summer months between high and low elevations. This almost certainly reflects individual differences, which we can quickly see here:

require(ggthemes)
ggplot(hr, 
  aes(timestamp, DEM, factor = nick_name, color = nick_name)) + geom_path() + theme_few()

This is actually a super-interesting plot, since it indicates that quite a few animals have a strong seasonal range shift with some altitudinal site fidelity in the winter, but less so in the summer. It also shows exactly where one rather extreme altitudinal movement occurs near the beginning for HR110.

Here’s an example of an (arguably - extreme) piped workflow that (a) subsets the portion of data that we’re interested in, (b) extracts the columns we’re interested in, (c) “melts” the data frame, i.e. repeats it three times, once each for longitude, latitude and elevation, and finally (d) plots it along the three “facets” of interest:

require(reshape2)
hr %>% subset(Used & DEM < 800) %>% 
  magrittr::extract(,c("longitude","latitude","DEM","nick_name","timestamp", "Used")) %>% 
  melt(c("nick_name","timestamp","Used")) %>%  
  ggplot(aes(timestamp, value, color = nick_name)) +
    facet_wrap(.~variable, scales = "free_y", nrow = 3) + 
    geom_path() + theme_few()

The altitudinal migration really is clearest here - and something eventually to think about when modeling elevation.

Here’s are some equivalent road density violinplots against season:

ggplot(hr %>% mutate(Road.Density = Road.Density + 1e-6), 
       aes(Used, Road.Density, fill = factor(month))) + geom_violin(draw_quantiles = c(.25,.5,.75)) + scale_y_log10() 

It is harder to pick out any strong patterns here, though there might be greater avoidance of roads in the summer than in the late winter. Note, the bottom of the barbells are artificially set to 1e-6.