In BiologicalRecordsCentre/sparta: Trend Analysis for Unstructured Data

This is a visualiation of an occupancy model produced using the r-package sparta. For more information of sparta visit https://github.com/biologicalrecordscentre/sparta

knitr::opts_chunk$set(echo = TRUE, fig.align = 'center')
results <- readRDS(params$dataFile)
require(sparta)
require(R2jags)

Table of contents

1. Basic Information
2. Species trend
3. Traceplots - Annnual occupancy estimates
4. Traceplots - Annnual detectability estimates
5. Traceplots - Other parameters
6. Rhat values

Basic Information {#Basic_Information}

cat(paste('Species:', results$SPP_NAME, '\n'))
cat(paste('Year range:', results$min_year, '-', results$max_year, '\n'))
cat(paste('Iterations:', results$n.iter, '\n'))
cat(paste('Chains:', results$BUGSoutput$n.chains, '\n'))
cat(paste('Burn in:', results$BUGSoutput$n.burnin, '\n'))
cat(paste('Thinning:', results$BUGSoutput$n.thin, '\n'))
cat(paste('Number of sites:', results$nsites, '\n'))
cat(paste('Number of visits:', results$nvisit, '\n'))
cat(paste('Number of sites with records of', paste0(results$SPP_NAME,  ':'), results$species_sites), '\n')
cat(paste('Number of observations of', paste0(results$SPP_NAME,  ':'), sum(results$species_observations), '\n'))

Species trend {#Species_trend}

plot(results)

if('regions' %in% names(results)){
  for(i in results$regions){
    print(plot(results, reg_agg = i, main = i))
  }
}
if('region_aggs' %in% names(results)){
  for(i in names(results$region_aggs)){
    print(plot(results, reg_agg = i, main = i))
  }
}

Traceplots - Annnual occupancy estimates {#Traceplots_Annnual_occupancy_estimates}

array_sim <- results$BUGSoutput$sims.array
comb.samples <- mcmc.list(
  lapply(1:results$BUGSoutput$n.chains,
         FUN = function(x, array_sim){
           year_ests <- colnames(array_sim[,x,])[grepl('^psi.fs\\[', colnames(array_sim[,x,]))]
           ar_temp <- array_sim[ , x, year_ests]
           colnames(ar_temp) <- paste('Occupancy - Year', 1:ncol(ar_temp))
           as.mcmc(ar_temp)
           },
         array_sim = array_sim)
)
plot(comb.samples)

Traceplots - Annnual detectability estimates {#Traceplots_Annnual_detectability_estimates}

array_sim <- results$BUGSoutput$sims.array
comb.samples <- mcmc.list(
  lapply(1:results$BUGSoutput$n.chains,
         FUN = function(x, array_sim){
           year_ests <- colnames(array_sim[,x,])[grepl('^alpha.p\\[', colnames(array_sim[,x,]))]
           ar_temp <- array_sim[ , x, year_ests]
           colnames(ar_temp) <- paste('Detectability - Year', 1:ncol(ar_temp))
           as.mcmc(ar_temp)
           },
         array_sim = array_sim)
)
plot(comb.samples)

Traceplots - Other parameters {#Traceplots_Other_parameters}

array_sim <- results$BUGSoutput$sims.array
comb.samples <- mcmc.list(
  lapply(1:results$BUGSoutput$n.chains,
         FUN = function(x, array_sim){
           params_other <- colnames(array_sim[,x,])[!(grepl('^alpha.p\\[', colnames(array_sim[,x,])) | grepl('^psi.fs\\[', colnames(array_sim[,x,])))]
           ar_temp <- array_sim[ , x, params_other]
           # colnames(ar_temp) <- paste('Detectability - year', 1:ncol(ar_temp))
           as.mcmc(ar_temp)
           },
         array_sim = array_sim)
)
plot(comb.samples,
     density = FALSE,
     smooth = FALSE,
     omi = c(1,1,1,1),
     # yaxs = 'i',
     lty = 1)

Rhat values {#Rhat_values}

cols <- vals <- rev(results$BUGSoutput$summary[,'Rhat'])
cols[vals <= 1.01] <- 'green'
cols[vals <= 1.1 & vals > 1.01] <- 'yellow'
cols[vals > 1.1] <- 'red'
label_max <- max(4.1, max(nchar(names(vals)))/1.8)
par(mar = c(5, label_max, 4, 2))
barplot(vals,
        horiz = TRUE,
        col = cols,
        las = 1,
        xlim = c(0.9, max(vals) + 0.2),
        xlab = 'Rhat',
        main = 'Rhat values\nThresholds are at 1.01 and 1.1',
        width = 1,
        offset = 0,
        yaxs = 'i',
        xpd = FALSE)
abline(v = 1.1, col = 'black', lty = 5)
# text(labels = '1.1',
#      pos = 4,
#      x = 1.1,
#      y = length(vals) + 1)
abline(v = 1.01, col = 'black', lty = 5)
# text(labels = '1.01',
#      pos = 2,
#      x = 1,
#      y = (length(vals) + 1) * 1.25)

BiologicalRecordsCentre/sparta documentation built on April 22, 2024, 2:34 p.m.