Multi-layered linear stochastic differential equation inference
This page presents the supplementary material and software used for
multi-layered analysis using vectorial linear
stochastic differential equations. The reason for using stochastic
differential equation is to facilitate studying irregularly spaced
observations like those we find in our study dataset of micro-fossils.
This work was performed in
conjuction with doing the research for and writing the manuscript
"Phenotypic Evolution studied by Layered Stochastic Differential Equations" by
T. Reitan, T. Schweder and J. Henderiks.
Originally, this was used for modelling phenotypic
evolution towards a optima that itself could be a stochastic process
depending responding to layers below it (climate or primary optimum).
However, it may well be that the method can be used for other
applications as well, with a little tweaking.
Current versions of the analysis code:
- An R package,
for the analyzis has been made, which uses the same
underlying C++ code (layeranalyzer.cpp).
It can be installed with the R code:
and used after invoking library(layeranalyzer).
A more official (less frequently updated but thus more stable)
version lies here:
- layeranalyzer.cpp. This is the
current source version of the program (replacing "layer_analyzer_complex"
see below). This version does not depend on any library but the
Lapack library (for linear algebra operations), which typically can be
found in an R installation. The program can then typically be
compiled on a Linux machine like this:
"g++ -DMAIN -DENGLISH_LANGUAGE -I/usr/include/R -I/usr/include/R/R_ext -o layeranalyzer layeranalyzer.cpp -lm -llapack". If the Lapack library files or header files
are somewhere else on your computer, you should change the "-I" (which
tells the compiler where to look for header files) and "-L" (which tells
the compiler where to search for library files other than in the
This has already been tested on multiple platforms. See "help" texts for more
on usage. See below in this text for some examples.
- A short powerpoint presentation can be found
here, with a PDF version
here. Made for the Seventh Conference on
Bayesian Inference in Stochastic Processes in Madrid, Spain,
- A simpler presentation made for CEES
of the basic modelling framework and how it connects to the application
can be found here.
- A how-to concerning the program usage (an older version called
'layer_analyzer_complex' is used, but it has the same input structure as the
current 'layeranalyzer' program), with a short intro to
continuous time modelling is found
Papers on theory and application:
Reitan, T., Schweder, T., Henderiks, J. (2012)
Phenotypic Evolution studied by Layered Stochastic
Annals of Applied Statistics, Volume 6 (4): 1531-1551.
- supplementary.pdf: A
supplementary to the article "Phenotypic Evolution studied by
Layered Stochastic Differential Equations", which
goes more into detail concerning the algorithmical aspects
of the inference.
- Liow, L.H., Reitan, T., Harnik, P.G. (2015)
Ecological interactions on macroevolutionary time scales;
clams and brachiopods are more than ships that pass in the night
Ecology Letters, Volume 18(10): 1030-1039.
- Reitan, T., Liow, L. H. (2017)
An unknown Phanerozoic driver of brachiopod extinction rates unveiled
by multivariate linear stochastic differential equations.
Paleobiology, Published online.
R example code:
PS: All example datasets are provided in the package, so that the
file reading is not necessary.
- analyze_malta_data.R: Analyzis of reed
warbler body size, provided by
Camilla Lo Cascio Sætre, published in
and normalize_harelynx.R: Analyzis
for hare/lynx data from the Hudson Bay Company. The data is first
normalized (using normalize_harelynx.R) before analysis is performed, in
normalize_harelynx.R. This script is
not strictly necessary to run, for the
script, since the normalized datasets ('hare.norm' and 'lynx.norm')
are provided in the package. But the transformation
function is defined in the first sccript and is used for plotting
in the second script.
The code does a brief foray into the process structure of each
time series separately, before doing a "simple" connection analysis
with the simplest process structure available. Some residual analysis
is then performed, which reveals that changes in the process
structure may be warranted. The process structure for lynx is
then expanded and a new connection analysis is performed.
Please not that running the code takes some time!
. Analysis of bivalve+brachiopod diversifications rates.
This analysis seeks to check the model deemed best in
Reitan & Liow (2017), using the R package. Some extra residual analysis
is also done.
Example datasets (see in addition the example code):
Simulations can be used for testing the behaviour of the analysis.
Updates and older versions
- 05/01-2019: R package added:
- 23/11-2017: New version of the analysis program simply
called "layeranalyzer" added. This version does not depend on
anything but the Lapack library. All parts of the Hydrasub library
that were needed, has been included in the program code itself.
This should make compilation easier and allow for compilation on
more platforms. The change was made in order to prepare the groun for
an R package version.
- 25/2-2015: Debug on layer_analyzer_complex has been performed
on the option for linear time trends and the option has been tested on
simulated files. The input for multiple files has been shortened
down to avoid using having dual ways of specifying the same options.
Run layer_analyzer_complex without input arguments to see the changes.
- 21/3-2013: The program (layer_analyzer_complex) can now
subtract a periodic
pattern (sine/cosine) from the observations, thus allowing seasonality or
other strict periodic phenomena to be removed. The length of
the period is specified. Multiple periodicities can be specified,
making ti possible to handle arbitrarily complex period patterns.
The coefficients in fron of the sine/cosine terms enter in as parameters.
In addition, some extra debug info was added as options. Also
the parallel tempering was given an extra option, where the
multiplicative steps of the 'temperature' now can be specified.
Also, starting parameters for the MCMC chain can now be set.
A short intro to using the program was made,
here. The into does not show
all the options. Call the program without input in order to see that.
- 21/5-2012: The program layer_analyzer_complex is now
able to deal with multiple correlated series, using an
extra file specifying the correlation in the observations. Thus
allometry in biology is dealt with in the model. Correlations
as between series as well as causal pathways are possible. This
machinery is currently being tested on microfossil data. The
results so far look encouraging.
- 23/3-2012: A large expansion that will later replace
the current layer_analyzer is now made available as a
separate program named layer_analyzer_complex. The program deals
with separately specified feedback loops and also with the resulting
problem of complex eigenvalues.
- 7/2-2012: New code for dealing with pair-wise inter-regional
correlations in a specified layer. (Note that there should be no
more than 6 or at worst 7 sites for the methodology to be
efficient. If there's more, the resources to the initial Monte Carlo
estimation of the probability of drawing positive covariance
structures needs to be adjusted in the code. Either that or a
completely different prior for the correlation structure, like the
inverse-Wishart must be implemented.) Also, an indicator on/off
option for allowing the analysis to switch off inter-regional
correlations for some sites.
- 10/2-2012: Bug fix for 'layer_analyzer' concerning the prior
identification strategies in combination with itnegrated layers.
Also, test code for importance sampling added.
- 17/1-2012: Bug fix on 'layer_analyzer' on
the priors for hierarchical models with
identification restriction. Possibility of switching off
identification restrictions. Realizations from smoothing estimates
as well as simulations based on given parameters also implemented
as options. Possibility of plugging in secondary data sources that
will also be modelled using linear SDEs. The point is to make
'layer_analyzer' a single tool for all practical linear SDE problems.
The latest version of the analysis program and R package uses only
the Lapack library.
Older versions use the Hydrasub library for
support. This library in turn uses
(Gnu Scientific Library). Plotting programs, such as 'vvgraph' and
'histogramme' are used for displaying the results. PS: If plotting in
the program (rather than the R package) is to be used, then even
the latest version requires the 'vvgraph' and 'histogramme' programs
from the Hydrasub library. However, there are options for not doing
plots and for sending plots to file instead.
layer_analyzer but is able to deal with feedback loops and complex
eigenvalues (thus making this a tractable way of continuous time
ecological modelling). May run slower on such problems, but should run just
as fast as usual on models without feedback loops. Also able to
deal with non-stationary (exploding) time series, as long as
an initial value treatment is specified. Also, multiple measurements
at the same time possble (as long as they don't belong to the same
site and the same series). Thus phenotypical traits with
allometry can be handled. Replaced the old "layer_analyzer"
as the default analysis program, but has now been replaced by the
program "layeranalyzer" (see above, and sorry for the lack of imagination),
which doesn not use the Hydrasub library and has replaced the use of the
GSL library for the more universally available "lapack" library.
layer_analyzer.H: A generalized
layered analyzer, though has been replaced by an even
more generalized program (now called layer_analyzer_complex),
which allows for feedback loops (something this version does not).
It does Bayesian analysis (and also
a subsequent ML analysis, if that is wanted) on a
layered stochastic differential equation set. More than one
time series can be analyzed at a time, with options for letting
a layer belonging to one series affect a layer belonging to another.
A Perl-script for traversing a
host of models (710 in number)
in order to find the posterior model probabilities and find the most
probable model. Uses 'layer_analyzer'.
A bundle of the supplementary text,
the dataset and
the main source codes (layer_analyzer).
(as well as the previously used underlying program library
Hydrasub) is licenced as
meaning that anyone can modify and use it
for other purposes than the original intent.
Up to my home page.
Last modified: Fri Aug 2 15:26:00 CET 2019