Sparta provides a range of tools for analysing trends in species occurrence data and is based on the work presented in Isaac et al (2014). The data that is used in these method is 'what where and when'. The 'what' is typically a species name. 'Where' is the location of the observation, sometimes referred to as the site. This is typically a 1km, 2km or 10km grid square but could also be a none regular location such as field sites or counties. 'When' is the time when an observation is made, and the requirements differ between methods. Some methods require a date while others require you to aggregate dates into time periods for comparison.
All of the methods described here require multi species data. This is because they use information across all species to assess biases.
In this vignette we will run through the methods and show how they can be used in reproducible examples.
Installing the package is easy and can be done from CRAN. Alternatively the development version can be installed from GitHub.
NOTE: JAGS must be installed before the R package installation will work. JAGS can be found here - http://sourceforge.net/projects/mcmc-jags/files/JAGS/
# Install the package from CRAN # THIS WILL WORK ONLY AFTER THE PACKAGE IS PUBLISHED # install.packages('sparta') # Or install the development version from GitHub library(devtools) # install_github('biologicalrecordscentre/sparta')
# Once installed, load the package library(sparta)
The functions in sparta cover a range of tasks. Primarily they are focused on analysing trends in species occurrence data while accounting for biases (see Isaac et al, 2014). In this vignette we step through these functions and others so that you can understand how the package works. If you have any questions you can find the package maintainers email address using maintainer('sparta')
, and if you have issues or bugs you can report them here
\pagebreak
Clearly when you are using sparta you will want to use your own data, however perhaps you are only at the planning stage of your project? This code shows you how to create some example data so that you can try out sparta's functionality.
# Create data n <- 8000 # size of dataset nyr <- 50 # number of years in data nSamples <- 200 # set number of dates nSites <- 100 # set number of sites set.seed(125) # set a random seed # Create somes dates first <- as.Date(strptime("1950/01/01", "%Y/%m/%d")) last <- as.Date(strptime(paste(1950+(nyr-1),"/12/31", sep=''), "%Y/%m/%d")) dt <- last-first rDates <- first + (runif(nSamples)*dt) # taxa are set semi-randomly taxa_probabilities <- seq(from = 0.1, to = 0.7, length.out = 26) taxa <- sample(letters, size = n, TRUE, prob = taxa_probabilities) # sites are visited semi-randomly site_probabilities <- seq(from = 0.1, to = 0.7, length.out = nSites) site <- sample(paste('A', 1:nSites, sep=''), size = n, TRUE, prob = site_probabilities) # the date of visit is selected semi-randomly from those created earlier time_probabilities <- seq(from = 0.1, to = 0.7, length.out = nSamples) time_period <- sample(rDates, size = n, TRUE, prob = time_probabilities) myData <- data.frame(taxa, site, time_period) # Let's have a look at the my example data head(myData)
In general this is the format of data you will need for all of the functions in sparta. The taxa and site columns should be characters and the time_period column should ideally be a date but can in some cases be a numeric.
There are many sources of wildlife observation data including GBIF (Global Biodiversity Information Facility) and the NBN gateway (National Biodiversity Network). Both of these repositories have R packages that will allow you to download this type of data straight into your R session (see rgbif and rnbn for details)
It can be useful to have a look at your data before you do any analyses. For example it is important to understand the biases in your data. The function dataDiagnostics
is designed to help with this.
# Run some data diagnostics on our data results <- dataDiagnostics(taxa = myData$taxa, site = myData$site, time_period = myData$time_period, progress_bar = FALSE)
The plot produced shows the number of records for each year in the top plot and the average list length in a box plot at the bottom. List length is the number of taxa observed on a visit to a site, where a visit is taken to be a unique combination of 'where' and 'when'. A trend in the number of observations across time is not uncommon and a formal test for such a trend is performed in the form of a linear model. Trends in the number of records over time are handled by all of the methods presented in sparta in a variety of different ways. Trends in list length are tested in the same manner, and both are returned to the console. A in list length can cause some methods such as the reporting rate methods to fail (see 'LessEffortPerVisit' scenario in Isaac et al (2014)) Unsurprisingly, since this is a random dataset, we have no trend in either the number of records or list length over time. This function also works if we have a numeric for time period such as the year
# Run some data diagnostics on our data, now time_period # is set to be a year results <- dataDiagnostics(taxa = myData$taxa, site = myData$site, time_period = as.numeric(format(myData$time_period, '%Y')), progress_bar = FALSE)
If we want to view these results in more detail we can interrogate the object results
# See what is in results.. names(results) # Let's have a look at the details head(results$RecordsPerYear) head(results$VisitListLength) summary(results$modelRecs) summary(results$modelList)
Telfer's change index is designed to assess the relative change in range size of species between two time periods (Telfer et al, 2002). This is a simple method that is robust but has low power to detect trends where they exist. While this method is designed to compare two time period sparta can take many time periods and will complete all pairwise comparisons.
Our data is not quite in the correct format for Telfer since it is used to compare time periods but our time_period
column is a date. We can fix this by using the date2timeperiod
function.
## Create a new column for the time period # First define my time periods time_periods <- data.frame(start = c(1950, 1960, 1970, 1980, 1990), end = c(1959, 1969, 1979, 1989, 1999)) time_periods # Now use these to assign my dates to time periods myData$tp <- date2timeperiod(myData$time_period, time_periods) head(myData)
As you can see our new column indicates which time period each date falls into with 1 being the earliest time period, 2 being the second and so on. This function will also work if instead of a single date for each record you have a date range
## Create a dataset where we have date ranges Date_range <- data.frame(startdate = myData$time_period, enddate = (myData$time_period + 600)) head(Date_range) # Now assign my date ranges to time periods Date_range$time_period <- date2timeperiod(Date_range, time_periods) head(Date_range)
As you can see in this example when a date range spans the boundaries of your time periods NA is returned.
Now we have our data in the right format we can use the telfer
function to analyse the data. The Telfer index for each species is the standardized residual from a linear regression across all species and is a measure of relative change only as the average real trend across species is obscured (Isaac et al (2014); Telfer et al, 2002).Telfer is used for comparing two time periods and if you have more than this the telfer
function will all pair-wise comparisons.
# Here is our data head(myData) telfer_results <- telfer(taxa = myData$taxa, site = myData$site, time_period = myData$tp, minSite = 2)
We get a warning message indicating that a large number of rows are being removed as duplicates. This occurs since we are now aggregating records into time periods and therefore creating a large number of duplicates.
The results give the change index for each species (rows) in each of the pairwise comparisons of time periods (columns).
head(telfer_results)
The reporting rates models in sparta are all either GLMs or GLMMs with year as a continuous covariate but are flexible, giving the user a number of options for their analyses. These options include the addition of covariates to account for biases in the data including a random site effect and fixed effect of list length.
In Isaac et al (2014) it was shown that reporting rate models can be susceptible to type 1 errors under certain scenarios and that with site and list length covariates the models performed better when the data were biased. These methods were found to out perform simple methods like Telfer.
The common feature among these models is that the quantity under consideration is the 'probability of being recorded'. When binomial models are used (as is the default), it's the 'probability for an average visit' for the Bernoulli version it is the probability of being recorded per time period.
Before undertaking modelling the data can be subset in an effort to remove data that may introduce bias. Model sub-setting was found to reduce power in Isaac et al (2014) but can partially deal with uneven sampling of site. This process can also be used with other methods and is not solely applicable to the reporting rate models.
The first function allows you to subset your data by list length. This works out, for each combination of 'where' and 'when' (a visit), the number of species observed (list length). Any records that to not come from a list that meets your list length criteria are then dropped.
# Select only records which occur on lists of length 2 or more myDataL <- siteSelectionMinL(taxa = myData$taxa, site = myData$site, time_period = myData$time_period, minL = 2) head(myDataL) # We now have a much smaller dataset after subsetting nrow(myData) nrow(myDataL)
We are also able to subset by the number of times a site is sampled. The function siteSelectionMinTP
does this. When time_period is a date, as in this case, minTP is minimum number of years a site must be sampled in for it be included in the subset.
# Select only data from sites sampled in at least 10 years myDataTP <- siteSelectionMinTP(taxa = myData$taxa, site = myData$site, time_period = myData$time_period, minTP = 10) head(myDataTP) # Here we have only lost a small number rows, this is because # many sites in our data are visited in a lot of years. Those # rows that have been removed are duplicates nrow(myData) nrow(myDataTP)
As you can see in the above example minTP specifies the number of years a site must be sampled in order to be included. However, our dataset is very well sampled so we might be interested in another measure of time. For example, you might want only sites that have been observed in at least 60 months. Let's see how this could be done.
# We need to create a new column to represent unique months # this could also be any unit of time you wanted (week, decade, etc.) # This line returns a unique character for each month unique_Months <- format(myData$time_period, "%B_%Y") head(unique_Months) # Week could be done like this, see ?strptime for more details unique_Weeks <- format(myData$time_period, "%U_%Y") head(unique_Weeks) # Now lets subset to records found on 60 months or more myData60Months <- siteSelectionMinTP(taxa = myData$taxa, site = myData$site, time_period = unique_Months, minTP = 60) head(myData60Months) # We could merge this back with our original data if # we need to retain the full dates myData60Months <- merge(myData60Months, myData$time_period, all.x = TRUE, all.y = FALSE, by = "row.names") head(myData60Months) nrow(myData) nrow(myData60Months)
Following the method in Roy et al (2012) we can combine these two functions to subset both by the length of lists and by the number of years that sites are sampled. This has been wrapped up in to the function siteSelection
which takes all the arguments of the previous two functions plus the argument LFirst
which indicates whether the data should be subset by list length first (TRUE
) or second (FALSE
).
# Subset our data as above but in one go myDataSubset <- siteSelection(taxa = myData$taxa, site = myData$site, time_period = myData$time_period, minL = 2, minTP = 10, LFirst = TRUE) head(myDataSubset) nrow(myDataSubset)
Once you have subset your data using the above functions (or perhaps not at all) the reporting rate models can be applied using the function reportingRateModel
. This function offers flexibility in the model you wish to fit, allowing the user to specify whether list length and site should be used as covariates, whether over-dispersion should be used, and whether the family should be binomial or Bernoulli. A number of these variants are presented in Isaac et al (2014). While multi-species data is required it is not nessecary to model all species. In fact you can save a significant amount of time by only modelling the species you are interested in.
# Run the reporting rate model using list length as a fixed effect and # site as a random effect. Here we only model a few species. system.time({ RR_out <- reportingRateModel(taxa = myData$taxa, site = myData$site, time_period = myData$time_period, list_length = TRUE, site_effect = TRUE, species_to_include = c('e','u','r','o','t','a','s'), overdispersion = FALSE, family = 'Bernoulli', print_progress = TRUE) }) # Let's have a look at the data that is returned str(RR_out) # We could plot these to see the species trends with(RR_out, # Plot graph {plot(x = 1:7, y = year.estimate, ylim = range(c(year.estimate - year.stderror, year.estimate + year.stderror)), ylab = 'Year effect (+/- Std Dev)', xlab = 'Species', xaxt = "n") # Add x-axis with species names axis(1, at = 1:7, labels = species_name) # Add the error bars arrows(1:7, year.estimate - year.stderror, 1:7, year.estimate + year.stderror, length = 0.05, angle = 90, code = 3)} )
The returned object is a data frame with one row per species. Each column gives information on an element of the model output including covariate estimates, standard errors and p-values. This object also has some attributes giving the year that was chosen as the intercept, the number of visits in the dataset and the model formula used.
These models can take a long time to run when your data set is large or you have a large number of species to model. To make this faster it is possible to parallelise this process across species which can significantly improve your run times. Here is an example of how we would parallelise the above example using hte R package snowfall.
# Load in snowfall library(snowfall) # I have 4 cpus on my PC so I set cpus to 4 # when I initialise the cluster sfInit(parallel = TRUE, cpus = 4) # Export my data to the cluster sfExport('myData') # I create a function that takes a species name and runs my models RR_mod_function <- function(taxa_name){ library(sparta) RR_out <- reportingRateModel(species_to_include = taxa_name, taxa = myData$taxa, site = myData$site, time_period = myData$time_period, list_length = TRUE, site_effect = TRUE, overdispersion = FALSE, family = 'Bernoulli', print_progress = FALSE) } # I then run this in parallel system.time({ para_out <- sfClusterApplyLB(c('e','u','r','o','t','a','s'), RR_mod_function) }) # Name each element of this output by the species RR_out_combined <- do.call(rbind, para_out) # Stop the cluster sfStop() # You'll see the output is the same as when we did it serially but the # time taken is shorter. Using a cluster computer with many more than # 4 cores can greatly reduce run time. str(RR_out_combined)
Using these functions it is possible to recreate the 'Well-sampled sites' method that is presented in Roy et al (2012) and Thomas et al (2015). This is made available in the function WSS
which is a simple wrapper around siteSelection
and reportingratemodel
. In this variant the data is subset by list length and the number of years each site was sampled before being run in a GLMM with site as a random effect.
# Run our data through the well-sampled sites function # This time we run all species WSS_out <- WSS(taxa = myData$taxa, site = myData$site, time_period = myData$time_period, minL = 2, minTP = 10, print_progress = FALSE) # The data is returned in the same format as from reportingRateModel str(WSS_out) # We can plot these and see that we get different results to our # previous analysis since this time the method includes subsetting with(WSS_out[1:10,], # Plot graph {plot(x = 1:10, y = year.estimate, ylim = range(c(year.estimate - year.stderror, year.estimate + year.stderror)), ylab = 'Year effect (+/- Std Dev)', xlab = 'Species', xaxt="n") # Add x-axis with species names axis(1, at=1:10, labels = species_name[1:10]) # Add the error bars arrows(1:10, year.estimate - year.stderror, 1:10, year.estimate + year.stderror, length=0.05, angle=90, code=3)} )
Occupancy models were found by Isaac et al (2014) to be one of the best tools for analysing species occurrence data typical of citizen science projects, being both robust and powerful. This method models the occupancy process separately from the detection process, but we will not go in to the details of the model here since there is a growing literature about occupancy models, how and when they should be used. Here we focus on how the occupancy model discussed in Isaac et al 2014 is implemented in sparta
.
This function works in a very similar fashion to that of the previous functions we have discussed. The data it takes is 'What, where, when' as in other functions, however here we have the option to specify which species we wish to model. This feature has been added as occupancy models are computationally intensive. The parameters of the function allow you control over the number of iterations, burnin, thinning, the number of chains, the seed and for advanced users there is also the possibility to pass in your own BUGS script.
# Here is our data str(myData) # Run an occupancy model for three species # Here we use very small number of iterations # to avoid a long run time system.time({ occ_out <- occDetModel(taxa = myData$taxa, site = myData$site, survey = myData$time_period, species_list = c('a','b','c','d'), write_results = FALSE, n_iterations = 200, burnin = 15, n_chains = 3, thinning = 3, seed = 123) }) # Lets look at the results ## The object returned is a list with one element for each species names(occ_out) # Each of these is an object of class 'occDet' class(occ_out$a) # Inside these elements is the information of interest names(occ_out$a) # Of particular interest to many users will be the summary # data in the BUGSoutput head(occ_out$a$BUGSoutput$summary) # We have included a plotting feature for objects of class # occDet which provides a useful visualisation of the trend # in occupancy over time plot(occ_out$a)
Here we have run a small example but in reality these models are usually run for many thousands of iterations, making the analysis of more than a handful of species impractical. For those with access to the necessary facilities it is possible to parallelise across species. To do this we use a pair of functions that are used internally by occDetModel
. These are formatOccData
which is used to format our occurrence data into the format needed by JAGS, and occDetFunc
, the function which undertakes the modelling.
# First format our data formattedOccData <- formatOccData(taxa = myData$taxa, site = myData$site, survey = myData$time_period) # This is a list of two elements names(formattedOccData)
formatOccData
returns a list of length 2; the first element 'spp_vis' is a data.frame with visit (unique combination of site and time period) in the first column and taxa for all the following columns. Values in taxa columns are either TRUE or FALSE depending on whether they were observed on that visit.
# Lets have a look at spp_vis head(formattedOccData$spp_vis[,1:5])
The second element ('occDetData') is a data frame giving the site, list length (the number of species observed on a visit) and year for each visit.
# Lets have a look at occDetData head(formattedOccData$occDetdata)
With our data in the correct format this can now go into the modelling function
# Use the occupancy modelling function to parrellise the process # Here we are going to use the package snowfall library(snowfall) # I have 4 cpus on my PC so I set cpus to 4 # when I initialise the cluster sfInit(parallel = TRUE, cpus = 4) # Export my data to the cluster sfExport('formattedOccData') # I create a function that takes a species name and runs my model occ_mod_function <- function(taxa_name){ library(sparta) occ_out <- occDetFunc(taxa_name = taxa_name, n_iterations = 200, burnin = 15, occDetdata = formattedOccData$occDetdata, spp_vis = formattedOccData$spp_vis, write_results = FALSE, seed = 123) } # I then run this in parallel system.time({ para_out <- sfClusterApplyLB(c('a','b','c','d'), occ_mod_function) }) # Name each element of this output by the species for(i in 1:length(para_out)) names(para_out)[i] <- para_out[[i]]$SPP_NAM # Stop the cluster sfStop() # This takes about half the time of the # serial version we ran earlier, and the resulting object # is the same (since we set the random seed to be the same # in each) head(para_out$a$BUGSoutput$summary) plot(para_out$a)
This same approach can be used on cluster computers, which can have hundreds of processors, to dramatically reduce run times.
The frescalo method is outlined in Hill (2012) and is a means to account for both spatial and temporal bias. This method was shown by Isaac et al (2014) to be a good method for data that is aggregated into time periods such as when comparing atlases. The frescalo method is run using a .exe, you will need to download this file by visiting this link - https://github.com/BiologicalRecordsCentre/frescalo. Once you have downloaded the .exe make a note of the directory you have placed it in, we will need that in a moment.
Again we will assume that your data is in a 'what, where, when' format similar to that we used in the previous method:
head(myData)
Frescalo's requirements in terms of data structure and types is a little different to that we have seen in other functions. Firstly the entire data.frame is passed in as an argument called Data
, and the column names of your various elements (taxa, site, etc) are given as other arguments. Secondly frescalo requires that the 'when' component is either a column of year or two columns, one of 'start date' and one of 'end date'. Our data as presented above does not fit into this format so first we must reformat it. In our situation the simplest thing to do is to add a column giving the year. Since frescalo aggregates across time periods (often decades or greater) this loss of temporal resolution is not an issue.
# Add a year column myData$year <- as.numeric(format(myData$time_period, '%Y')) head(myData)
Now we have our data in the correct format for frescalo there is one other major component we need, a weights file. You can find out more about the weights file and what it is used for in the original paper (Hill, 2012). In short the weights file outlines the similarity between sites in your dataset. This information is used to weight the analysis of each site accordingly. If you are undertaking this analysis in the UK at 10km square resolution there are some built in weights files you can use. Some of these weights files use the UK landcover map instead of floristic similarity (as used in Hill (2012)). You can find out more about these in the frescalo help file.
For the sake of demonstration let us assume that you do not have a weights file for your analysis, or that you want to create your own. To create a weights file you need two things, a measure of physical distance between your sites and a measure of similarity. In the original paper this similarity measure was floristic similarity, but it could also be habitat similarity or whatever is relevant for the taxa you are studying. In this example I have a table of distances and of land cover proportions at each site
# I'm going to create some made up data mySites <- unique(myData$site) # Build a table of distances myDistances <- merge(mySites, mySites) # add random distances myDistances$dist <- runif(n = nrow(myDistances), min = 10, max = 10000) # to be realistic the distance from a site to itself should be 0 myDistances$dist[myDistances$x == myDistances$y] <- 0 # Build a table of attributes # This can be done in numerous ways, here is one example # Lets say I have habitat data for all my sites myHabitatData <- data.frame(site = mySites, grassland = runif(length(mySites), 0, 1), woodland = runif(length(mySites), 0, 1), heathland = runif(length(mySites), 0, 1), urban = runif(length(mySites), 0, 1), freshwater = runif(length(mySites), 0, 1)) # This pretend data is supposed to be proportional cover so lets # make sure each row sums to 1 multiples <- apply(myHabitatData[,2:6], 1, sum) for(i in 1:length(mySites)){ myHabitatData[i,2:6] <- myHabitatData[i,2:6]/multiples[i] }
# Here is the distance table head(myDistances) # Here is our habitat data head(myHabitatData) # With our distance and habitat tables in hand we can # use the createWeights function to build our weights file # I have changed the defualts of dist_sub and sim_sub since # we have a very small example dataset of only 50 sites myWeights <- createWeights(distances = myDistances, attributes = myHabitatData, dist_sub = 20, sim_sub = 10) head(myWeights)
The createWeights
function follows the procedure outlined in Hill (2012) for creating weights and more information can be found in the help file of the function. With our data and weights file we are now ready to proceed with frescalo. As with other functions frescalo can take a range of additional arguments which you can see by entering ?frescalo
at the console, here we will do a minimal example.
# First we need to enter the location where we placed the .exe # In my case I saved it to a folder within the source package directory. Note that this installer has been filtered out in the .gitignore to reduce size, and therefore you will need to install it separately myFrescaloPath <- file.path(getwd(), "vignette_frescalo.exe") # I then want to set up the time periods I want to analyse # Here I say I want to compare 1980-89 to 1990-99 myTimePeriods <- data.frame(start = c(1980, 1990), end = c(1989, 1999)) head(myTimePeriods) # I also need to specify where I want my results to be saved # I'm going to save it in a folder in my working directory myFolder <- 'myFolder' # Simple run of frescalo frescalo_results <- frescalo(Data = myData, frespath = myFrescaloPath, time_periods = myTimePeriods, site_col = 'site', sp_col = 'taxa', year = 'year', Fres_weights = myWeights, sinkdir = myFolder)
We get a warning from this analysis that our value of phi is too low. In this case this is because our simulated data suggests every species is found on every site in our time periods. This is a little unrealistic but should you get a similar warning with your data you might want to consult Hill (2012) and change your input value of phi.
The object that is returned (frescalo_results
in my case) is an object of class frescalo
. this means there are a couple of special methods we can use with it.
# Using 'summary' gives a quick overview of our data # This can be useful to double check that your data was read in correctly summary(frescalo_results) # Using 'print' we get a preview of the results print(frescalo_results) # There is a lot of information here and you can read more about # what these data mean by looking at the frescalo help file # The files detailed in paths are also in the object returned frescalo_results$paths names(frescalo_results) # However we additionally get some model results in our returned object # under '$lm_stats'
The results from frescalo may seem complex at first and I suggest reading the Value section of the frescalo help file for details. In brief: frescalo_results$paths
lists the file paths of the raw data files for $log
, $stat
, $freq
and $trend
, in that order. frescalo_results$trend
is a data.frame providing the list of time factors (a measure of probability of occurrence relative to benchmark species) for each species-timeperiod. frescalo_results$stat
is a data.frame giving details about sites such as estimated species richness. frescalo_results$freq
is a data.frame of the species frequencies, that is the probabilities that a species was present at a certain location. frescalo_results$log
, a simple report of the console output from the .exe. frescalo_results$lm_stats
is a data.frame giving the results of a linear regression of Tfactors for each species when more than two time periods are used. If only 2 time periods are used (as in our example) the linear modeling section of this data.frame is filled with NAs and a z-test is performed instead (results are given in the last columns).
# Lets look at some results fo the first three species frescalo_results$lm_stats[1:3, c('NAME','Z_VAL','SIG_95')] # None of these have a significant change using a z-test # Lets look at the raw data frescalo_results$trend[frescalo_results$trend$Species %in% c('a', 'b', 'c'), c('Species', 'Time', 'TFactor', 'StDev')] # We can see from these results that the big standard deviations on # the tfactor values means there is no real difference between the # two time periods
If your data are from the UK and sites are given as grid referenes that there is functionality to plot a simple output of your results
# This only works with UK grid references # We can load an example dataset from the UK data(unicorns) head(unicorns) # Create a new time period range myTimePeriods <- data.frame(start= c(1968, 2001), end = c(2000, 2023)) # Now run frescalo using hte built in weights file unicorn_results <- frescalo(Data = unicorns, frespath = myFrescaloPath, time_periods = myTimePeriods, site_col = "site", sp_col = "species", start_col = "start_date", end_col = 'end_date', sinkdir = myFolder)
It is worth noting the console output here. We get a warning telling us that I have some data from a site that is not in my weights file, so I might want to investigate that and add the site to my weights file. We will ignore it for now. The second warning tells us that the sinkdir
that we gave already has frescalo output in it. The function has got around this by renaming the output. We finally got a long list of all the species as their data were compiled internally.
Now for the plotting.
plot(unicorn_results)
Each panel of the plot gives different information for your results. The top left plot shows the observed number of species at each site (given in unicorn_results$stat$No_spp
), this can be contrasted with the top right plot which gives the estimated number of species after accounting for recording effort (given in unicorn_results$stat$Spnum_out
). Recording effort is presented in the bottom left panel - low values of alpha (white) show areas of high recording effort (given in unicorn_results$stat$Alpha
), and a summary of the species trends are given in the bottom right (given in unicorn_results$lm_stats
). In this case there is a skew towards species increasing, however some of these may be non-significant, this could be explored in more detail be referring to unicorn_results$lm_stats
.
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.