spatial-ews / spatialwarnings Goto Github PK

fitpsd.R et patchsizedistr_ews.R

patchsizedistr_ews.R seems to be obsolete.

• Add the documentation.

May be we have also to think to the structure of this function ? Cumulative patch size distribution could be separated from power law fitting function ?

Add a workflow functions (spectral_spews + print/summary/plot methods) that informs the use about :

• changes in shape of powerspectrum
• summarise those changes in SDR

In more details :

• Add function that computes the result
• Summary method for the result
  • single run
  • list run
• Plot method for the result displaying a trend
  • list run
• Print method for the result
  • single run
  • list run

Make sure the spectral ews workflow is documented the following functions should have their respective roxygen comments:

• spectral_spews (main function)
  • plot methods
  • plot_spectrum methods
• indictest.spectral_spews_* (assess significance)
  • plot methods
  • plot_spectrum methods

There are a couple of data files in the original repo, but no documentation on how to interpret them.

Could we add a couple of valid binary matrices as data objects that we can use for writing examples?

As put in http://r-pkgs.had.co.nz/vignettes.html :

A vignette is like a book chapter or an academic paper: it can describe the problem that your package is designed to solve, and then show the reader how to solve it. A vignette should divide functions into useful categories, and demonstrate how to coordinate multiple functions to solve problems. Vignettes are also useful if you want to explain the details of your package. For example, if you have implemented a complex statistical algorithm, you might want to describe all the details in a vignette so that users of your package can understand what’s going on under the hood, and be confident that you’ve implemented the algorithm correctly.

Following this idea, I suggest dividing our vignette into several parts :

• Intro: General context and background. What are indicators ? In which context should somebody have a look at them ?
• Which indicators are implemented in the package ? (in particular: we do not address the irregular patterns cases).
• Generic EWS
  • Intro: The theory behind it + little description of each indicator
  • A workflow example
• Patch-based EWS (psd, psdfit)
  • Intro
  • A workflow example ( + some ideas to converting continuous data to discrete)
• Spectral EWS (spectral density ratio, powerspectrum)
  • Intro, the rationale behind it
  • A workflow example

Some R script contain source() command. We don't need those commands since when the package is loaded, all the scripts of the R folder are sourced.

A grep command in the repository gives that 20 source() command are present, we need to delete them.

To work locally, you can install the package by doing

install_github("fdschneider/spatial_warnings")
library(spatialwarnings)

Normally, it should work and you don't need to source manually individual scripts.

If you modify a function locally and you want to test it, you can run devtools::load_all() which load all the scripts which are in the R/ folder.

1. Check 2d phase diagram axes as per example.R code. Are the axes swapped?
2. Not all replicates for d values very close to transition in the continuous phase have non-zero densities. Patch Size distributions to be calculated only on the basis of finite density snapshots.

Create a function that will compute all the indicators and produce a summarizable/plottable object.

# We need a catchy short name for that
all_indicators <- function(data, # `binary_matrix` object
                           indicators = list(variance = indicator_variance, psd = indicator_psd, ...), # list of functions
                           ?)

It should return an object of class indicators with the following structure:

list(variance = list(value = xx, ...), # OR variance = xx
     skewness = list(value = xx, ...), 
     psd = `psd object`, 
     ...) # Component names are taken from the list of functions above.
(class: c('spatial_indicators', 'list'))

Numerical indicators (e.g. variance, skewness, moran's I) should have a value= in their output list so they get recognized by the summarize function.

We need to wrap this into a valid R package.

The globally accepted input should be a binary matrix, A.

• Any factorial matrix F with the focal state "+" can be easily converted into such by A <- F == "+".
• Same works for integer matrices, I. A <- I == 1

Many of the present functions can be run on that directly. Others need patch counts. Thus a core function label()should provide a fast labelling of patches and return a matrix with integer patch IDs.

Then a secondary function could take in the binary matrix as input and return a vector of patch sizes and/or a table of cumulative patch sizes.

Make sure the spectral ews workflow works: the following functions should be implemented and working.

• spectral_spews (main function)
  • print methods
  • summary methods
  • as.data.frame methods
  • plot methods
  • plot_spectrum methods
• indictest.spectral_spews_* (assess significance)
  • print methods
  • summary methods
  • as.data.frame methods
  • plot methods
  • plot_spectrum methods

For all of them we compute a null distribution of indicator values based on shuffled matrices, but the way we test against this distribution might have to differ from indicator to indicator:

• Moran's I should always be higher so we just do a one-tailed test for higher than expected values
• Variance should always be higher in the coarse-grained observed data than in the coarsed-grained null dat (so one-tailed test too). However, the significance can decrease just before the tipping point (e.g. observed in Sabiha's runs for the pq model).
• Skewness can increase or decrease depending on whether the shift is upwards or downwards. In all cases, the absolute value of skewness should always increase before a critical point. This absolute value skewness is expected to be lower on a null coarsed-grain matrix (??) than observed.

These functions will be available to users as they can be useful, and need documentation:

• rspectrum (compute the r-spectrum)
• coarse_grain
• patchsizes (get patch size distribution of a binary matrix)
• as.binary_matrix (convert to binary matrices)

Currently there is no standard naming for the output of the indicator_* functions: we need that so the summary function can produce the right output for each of them (e.g. a table with values for numerical indicators, see #19).

Numerical indicators (e.g. variance, skewness, etc.)

In particular, all indicators that return "a single number" and not a complex object (e.g. a curve fit) should return EITHER :

• a single numeric value (numeric vector of length 1) OR
• a list with (but not necessarily limited to) a component named value

Here is a list of the numerical functions that need to be modified :

• indicator_moran
• indicator_skewness (/!\ [Alex] I'm currently cleaning so please wait a bit)
• indicator_variance (/!\ [Alex] I'm currently cleaning so please wait a bit)
• indicator_largestpatch

Non numerical indicators (e.g. patch-size distr, etc)

These functions should explicitely add a specific class to their output. For example, indicator_fitpsd sets the class of its output to psdfit.

There should be a summary method available for this class that will be called within the global summary of all the indicators. (e.g. summary.psdfit).

These ones are :

• indicator_psdfit (mostly done, a good example)
• indicator_cumpsd
• indicator_corrfunc
• indicator_powerspectrum
• indicator_fractalgeometry

I might have forgotten some indicators so please add them below :o)

After running all the indicators the user is likely to want to export the results to a more standard object like a data.frame (so it can do beautiful plots in ggplot for example).

To do this we need functions that take the output object:

list(variance = list(value = xx, ...), # OR variance = xx
     skewness = list(value = xx, ...), 
     psd = `psd object`, 
     ...) 
(class: c('spatial_indicators', 'list'))

and convert it to something like this

> as.data.frame(indic_object) 

data.frame of x rows and x columns
         replicate value
variance 1            xx
variance 2            yy
variance 3            zz
skewness 1            aa
skewness 2            bb
(etc.)

class: c('data.frame', 'list')

Non numerical indicators cannot fit into this data structure so I suggest dicarding them.

Make sure the patch-based ews workflow is documented: the following functions should have their respective roxygen comments:

• patchdistr_spews (main function)
  • plot methods
  • predict methods

This function should print the ic_value ~ coarse_graining size plot so the user can decide on an optimal coarse-graining size.

When the input is several matrices along a gradient, then we don't know on which snapshot this curve should be displayed. Moreover they are not comparable across values of the gradient as density changes.

A rule of thumb could be using (1) something not too close to the tipping point and (2) with a density that is not too high.

If we plan to submit our package on the CRAN, we need to pass all the required automatic tests. CRAN reviewers are very strict about this so all NOTEs/WARNINGs/ERRORs must be squashed.

The package can be checked within R using the devtools::check() function. Pass the argument cran=TRUE to run a CRAN-like check.

Note that this is not to be confused with the unit-tests run using test(). What check() does is checking whether the package will be runnable on other computers and whether it follows the CRAN guidelines. (check() run the unit tests though and stops if those fail).

Example output:

* using log directory ‘<snip>/spatial_warnings.Rcheck’
* using R version 3.2.1 (2015-06-18)
<snip>
Status: 1 ERROR, 10 WARNINGs, 6 NOTEs
See

Error:
unexpected end of input in spatial_indicators_main.R

e.g. This code produces an error

ex <- summary(generic_spews(forestdat$matrices), null_replicates = 0)
print(ex)

Make sure the generic ews workflow works: the following functions should be implemented and working.

• generic_spews (main function)
  • print methods
  • summary methods
  • as.data.frame methods
  • plot methods
• indictest.generic_spews_* (assess significance)
  • print methods
  • summary methods
  • as.data.frame methods
  • plot methods

When using detrending, means can be negative due to floating point errors, so a minus sign is added in front of the number => an extra column is added and the layourt breaks. Similarly, plots show a bogus trend.

Solution : set means to zero when detrending has been used.

Reproducible example:

a <- generic_spews(forestdat$matrices, detrend = TRUE)
a.summary <- summary(a)

a.summary # broken layout
plot(a.summary) # broken mean trend

input is a binary matrix or a list of binary matrices (using S3 methods)
See discussion at fdschneider/caspr#4

Indicator functions should return a warning when a matrix looks like a binary matrix while not being one as it has big consequences on the way indicators are computed.

Maybe the right test could be length(unique(states)) <= 1

• Change the structure of the fonction in the standard way as it was proposed by Alex.
• Add the documentation

• Add documentation

This function is called in the case of spectral analysis, right ?

• Find reference of the index.
• Change the structure of the fonction in the standard way as it was proposed by Alex.
• Add the documentation

The package needs to do its best try at converting data to binary matrices. The best way to do this I think is to create a binary_matrix class and as.binary_matrix.* S3 methods :

• as.binary_matrix.data.frame
• as.binary_matrix.matrix (for non-binary matrices)
• as.binary_matrix.list (for lists of stuff, determines what kind of objects the list is made of and apply the appropriate as.* function to them)
• as.binary_matrix.ca_run (not to be defined in this package but in caspr)
• ?

We should use the following ranges:

• remove the zero frequency
• low_range: lower 20%
• high_range: higher 20%

Having all the single-indicator functions, we need a user-friendly function now that computes all the indicators and reports about their values. I suggest following the usual pattern of creating a new object class "shift_indicator" and add usual methods.

A user workflow could look this way (a bit similar to a glm workflow basically) :

1. convert original data to binary matrix(/ces) (use as.* function family)
2. run the indicators function to create the result object containing all the appropriate indicator values (naming is to be discussed)
3. using the summary function on the indicators object, see what indicators show an interesting value
4. plot diagnostics

What I have in mind for a summary function is something with an output like :

 > summary(glm(Sepal.Length ~ Sepal.Width + Species, data = iris))

Call:
glm(formula = Sepal.Length ~ Sepal.Width + Species, data = iris)

Deviance Residuals: 
     Min        1Q    Median        3Q       Max  
-1.30711  -0.25713  -0.05325   0.19542   1.41253  

Coefficients:
                  Estimate Std. Error t value Pr(>|t|)    
(Intercept)         2.2514     0.3698   6.089 9.57e-09 ***
Sepal.Width         0.8036     0.1063   7.557 4.19e-12 ***
Speciesversicolor   1.4587     0.1121  13.012  < 2e-16 ***
Speciesvirginica    1.9468     0.1000  19.465  < 2e-16 ***

---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Thoughts ? Ideas for the plot function ?

Use one of the templates provided in [1] to add a function that returns the largest patch size (with user documentation).

[1] https://github.com/fdschneider/spatial_warnings/blob/master/R/indicator_example.R

input should be a vector of patch sizes.
returns a list of fitting results for all candidate distributions and a general summary of the best fitting function. see issue fdschneider/caspr#2

Some indicators will not run on the example datasets :

• indicator_powerspectrum (cannot find function rspec_ews)
• indicator_corrfunc (missing value where TRUE/FALSE needed)

Here are the bifurcation diagrams with two cross-sections (Figure in the same style as Flo's Fig. 2).

In top row plots, areas shaded in gray indicate that density of mussels/veg/trees is zero. Contour lines indicate the density.

I also put this figure in the dropbox with the corresponding code (in ./figures/bifurcation_diagrams)

• Find reference of the index.
• Change the structure of the fonction in the standard way as it was proposed by Alex.
• Add the documentation

• Change the structure of the fonction in the standard way as it was proposed by Alex.
• Add the documentation

Most of the code is shared between those three functions, most of the code can be factorized. (not a high priority though)

Right now skewness returns NA if the vector of values has a std. deviation of zero. It should return zero instead.

Add a check in coarse_grain() to catch this ? (or maybe do it at higher level to avoid exceptions in highly-repeated code).

Make sure the generic ews workflow is documented the following functions should have their respective roxygen comments:

• generic_spews (main function)
  • plot methods
• indictest.generic_spews_* (assess significance)
  • plot methods

• Find publication of the index. I have already the paper of Moran (1950) but it's quite hard to read. There are others publications on this index.
• Change the structure of the fonction in the standard way as it was proposed by Alex.
• Add the documentation

The model results are starting to arrive !

• Grazing model (~220 MB, ~3GB in memory)
• Musselbed model
• Forestgap model (~88MB)

As usual these files contain two variables: upper_branch and lower_branch. Both are lists of three components :

1. DatBif
2. snaps
3. time_series (new!) containing the full time series of global and local cover

access them with upper_branch[[1]] (or 2 or 3)

It would be great to have the list of all the indicator functions that we need for the package in this respository.

• Change the structure of the fonction in the standard way as it was proposed by Alex.
• Add the documentation

It seems a function to do coarsed grain

• Change the structure of the fonction in the standard way as it was proposed by Alex.
• Add the documentation

Make sure the generic ews workflow works: the following functions should be implemented and working.

• patchdistr_spews (main function)
  • print methods
  • summary methods
  • as.data.frame methods
  • plot methods
  • plot_distr methods

It is unclear what the indictest method should be (if it is needed).