Package {layeranalyzer}

Objects of class 'layered', as returned from the 'layer.analyzer' method.

A set of 'layered' objects returned from the layer.analyzer function. This set represents the model analyses performed and should contain information used for model comparison (marginal likelihoods for Bayesian analyses, AIC/BIC/AICc for classic analyses).

A numeric array describing the prior probabilities of each models sent to the function. The size of the array needs to match the number of models sent to the function. If prior probability array, 'p0', is not given each model is default given the same prior probabiliuty, but see item 'first.is.nullhypothesis'.

If prior probabilities (item 'p0') is not given, the prior probabilities are default equal for all models. However, if first.is.nullhypothesis is set to 'TRUE', the first model is given 50% of the prior probability, while the rest of the models shares the remaining 50% equally.

If classic ML analysis is performed instead of Bayesian analysis, "AIC","BIC" or "AICc" will be used for model comparison instead of the Bayesian marginal likelihood. The item 'ML.IC' determines which of these three model comparison criteria are to be used.

A set of 'layer.series.structure' objects that represents the set of time series that are to be analyzed and how the underlying processes of each such series is to be modelles. (See 'layer.series.structure'). At least one such object must be given.

num.MCMC is an integer specifying the number of Markov chain Monte Carlo (MCMC) samples that the analysis will rest on. Default 1000. More MCMC samples are run internally, depending on the spacing (default 10) and burnin (default 1000).

Specifies the number of MCMC iteration between each sample that is used (default 10). A higher number means the samples will be less correlated but requires proptionally more computer time.

Markov chain Monte Carlo sampling takes some iterations to converge to the posterior parameter distribution. The burn-in phase is the set of iterations performed where the samples are then discarded, so as to avoid including samples before convergence. Convergence can be check by for instance 'gelman.diag' in the 'coda' package, and the MCMC samples themselves can be gotten by using the option 'mcmc'. Still, stability of the results can be quite a good enough measure for most purposes. Default is 10000 iterations. PS: The burn-in algorithm here includes two steps where the random walk MCMC variance is adjusted in order to optimize the efficiency of the algorithm.

The number of tempering chains (default 1, meaning it is just the MCMC chain itself that is beeing run). Parallel tempering (Geyer 1991) is a method for overcoming local optima and can thus increase the stability of the results. The way it is done is that in addition to the MCMC chain for the posterior distribution, parallel chains sampling from the distribution proportional to exp(-log(likelihood*prior)/T) where T is a "temperature" larger than one (this is in practise a "smoothed" version of the posterior), is also sampled from. Swaps between neighbouring chains are then allowed, making it possible for the MCMC chain for the posterior distribution to find a new optima.

As the code is not parallel, the computer time will increase proportional to the number of tempering chains. The set of "temperatures" is set to (1, T.ground, T.ground^2, ..., T.ground^(numtemp-1)), where 'T.ground' is default set to 1.5. If no swapping is done (this can be seen if 'silent.mode' is set to 'FALSE'), this "ground temperature" must be lowered, in order to facilitate the swaps.

Logical variable. Only makes sense for Bayesian analysis (which is default). If set to 'TRUE' (default), the Bayesian model (marginal) likelihood is estimated using an importance sampler. If set to 'FALSE', this will not happen, which saves time but means that Bayesian model comparison will not be possible.

Logical variable. If set to 'FALSE' (default), Bayesian analysis is performed. If set to 'TRUE', a classic maximum likelihood (ML) is performed, where the starting points of a set of optimizations is drawns from the MCMC samples. (This means that Bayesian MCMC sampling is performed anyway, but that the results from this is not focused on in the summary). The number of optimizations is determined by option 'maximum.likelihood.numstart', which is default set to 10. Estimation of the Bayesian model (marginal) likelihood is switched off when this option is used.

The number of ML optimizations performed, it applicable. (If 'do.maximum.likelihood' is set to 'TRUE'.) Default:10.

In layeranalyzer, ML optimizations are per default (maximum.likelihood.strategy="MCMC-from-model") started from MCMC samples,a s it was found that completely randomly drawn starting points worked poorly for large models. However, MCMC convergence is no guarantee either, especially for models with multiple datasets and connections between them. 'maximum.likelihood.strategy' allows for using alternative startegies for optimization in connection models by MCMC sampling from a model that is the same in structure but without connections and using these samples as starting points for the optimization instead, which is specified by setting maximum.likelihood.strategy="MCMC-from-nonconnected". The third possibility is to set maximum.likelihood.strategy="MCMC-from-hybrid", which alternates betwen using MCMC samples from the connection model and from MCMC samples for the non-connected model.

If set to 'FALSE', the routine will show a lot debug information while running. This can be useful for tempering purposes, since the number of tempering swaps are shown in this debug information.

If set to 'TRUE', this triggers the printing of more specific debug information concerning the burn-in phase.

If set to 'TRUE', this triggers the printing of more specific debug information concerning the likelihood calculation phase. (PS: This means a lot of printing.)

Determines how to handle the identifiability problem in Bayesian multi-layer analysis. As mentioned in the supplementary of Reitan et al. (2012), one can switch the characteristic times of neighbouring layers and reorganize the stochastic contributions in such a way that the same top layer process takes place. This can be solved by requiring that the lower the layer, the larger the characteristic time. However, simply requiring this means the prior distribution of the characteristic times get shortened for multi-layered models, as compared to one-layered model. As long as the data is within the range of these shortened prior distributions, this puts the one-layer model at a disadvantage when it comes to model comparison. There are various solutions to this problem, and which solution is used is determined by the variable 'id.strategy'.

The variable 'id.strategy' should be an integer between 0 and 4, default 2. 0 - No identification treatment. (Default) PS: Can and even should yield multimodal characteristic times 1 - Keep upper characteristic time prior. Add lognormally to beneath-lying characteristic times. 2 - Keep lower characteristic time prior. Substract lognormally to above-lying characteristic times. 3 - Keep lower characteristic time prior. Cut depending on that on the above-lying characteristic times. 4 - Keep upper characteristic time prior. Cut depending on that on the below-lying characteristic times.

If set to 'TRUE', instead of reporting the stochastic contribution size in the stochastic differential equations, the stationary standard deviation is reported. The stationary standard deviation is the standard deviation of the process state at any given time after convergence, unconditional on previous states. In a multi-layered model, this is to be interpreted as the stationary standard deviation if this layer when let alone (with no other process affecting it). This is calculated as s*sqrt(characteristic time/2), where s is the stochastic contribution size. PS: This option does not make any sense if the bottom layer is a Wiener process, as that process is not stationary.

This determines how the parallel tempering chains are to be defined. Parallel tempering (Geyer 1991) is a method for overcoming local optima and can thus increase the stability of the results. The way it is done is that in addition to the MCMC chain for the posterior distribution, parallel chains sampling from the distribution proportional to exp(-log(likelihood*prior)/T) where T is a "temperature" larger than one (this is in practise a "smoothed" version of the posterior), is also sampled from. Swaps between neighbouring chains are then allowed, making it possible for the MCMC chain for the posterior distribution to find a new optima.

As the code is not parallel, the computer time will increase proportional to the number of tempering chains. The set of "temperatures" is set to (1, T.ground, T.ground^2, ..., T.ground^(numtemp-1)), where 'T.ground' is default set to 1.5. If no swapping is done (this can be seen if 'silent.mode' is set to 'FALSE'), this "ground temperature" must be lowered, in order to facilitate the swaps.

If set to 'TRUE', reports half-lives rather than characteristic times. Half-lives are the time the auto-correlation of an OU process drops to 1/2, or the time it takes for a perturbation from the expected value to drop to 1/2 the original size in expected value. (As opposed to dropping to exp(-1) for characteristic times.) Half-life=log(2)*characteristic time. The auto-correlation for an OU process (and thus for the inner dynamics of any linear layer) is exp(-diff.time/characteristic.time)=exp(-log(2)*diff.time/half.life)= (1/2)^(diff.time/half.life), where diff.time is the time difference between two process states.

If set to 'TRUE', returns the MCMC samples as an 'mcmc' object within the return object. The name of the list object will also be 'mcmc'.

If multiple time series are given, one can specify causal connections from underlying processes belonging to such time series to the underlying processes belonging to some other time series. This specification should come in the form of an 4 x n matrix, where n is the number of causal connections. For each column (i.e. for each connection), row 1 specifies the number of the cause time series (same as the ordering of the input 'layer.series.structure' input to the routine), row 2 specifies the layer of the cause process, row 3 specifies the number of the effect time series and row 4 specifies the layer of the effect process. So if one has two time series, X and Y, both being modelled as having 3 layers, one can specify a causal connection from layer 3 in X to layer 1 in Y by the matrix "causal<-matrix(c(1,3,2,1),ncol=1)".

Note that it is possible to specify a causal connection even with only one time series, if a multi-layered process is considered. So if X is two-layered, "causal<-matrix(c(1,1,1,2),ncol=1)" will work. As layers are defined so that there is always a causal connection from layer 'k' to layer 'k-1' (with strength=1 for identifiability purposes), the reverse specification ("causal<-matrix(c(1,2,1,1),ncol=1)") will not make sense.

Specifies a set of causal connections going both ways with equal strength. For more on the usage, see option 'causal', but note that symmetric causal connection between one layer and the next does not make sense, as the connection strength from layer 'k' to layer 'k-1' is always set to "1" for identifiability purposes. Note that which series is first and which is second does not matter, so the matrix "causal.symmetric<-matrix(c(2,1,1,3),ncol=1)" would specify the same model as "causal.symmetric<-matrix(c(1,3,2,1),ncol=1)".

A 4 x n matrix which specifies the correlative connections between the processes, where n is the number of such connections. For each column (i.e. for each connection), row 1 specifies the number of the first time series (same as the ordering of the input 'layer.series.structure' input to the routine), row 2 specifies the layer of the first process, row 3 specifies the number of the second time series and row 4 specifies the layer of the second process. So if one has two time series, X and Y, both being modelled as having 3 layers, one can specify a correlative connection between layer 3 in X and layer 1 in Y by the matrix "corr<-matrix(c(1,3,2,1),ncol=1)". Note that which series is first and which is second does not matter, so the matrix "corr<-matrix(c(2,1,1,3),ncol=1)" would specify the exact same model.

A list that specifies whether smoothing samples are to be returned and if so how. Smoothing samples are a posteriori samples of the process states at measurement times and possibly also other times. Smoothing samples can be turned on by setting smoothing.specs$do.smoothing=TRUE. If so, the smoothing.specs$smoothing.time.diff and smoothing.specs$num.smooth.per.mcmc must also be given. If smoothing.specs$smoothing.time.diff is set to "0", only the states at the measurement times will be inferred. If smoothing.specs$smoothing.time.diff is set to a positive real value, the state at the regularly spaced time points (smoothing.specs$smoothing.start, smoothing.specs$smoothing.start+smoothing.specs$smoothing.time.diff, smoothing.specs$smoothing.start+2*smoothing.specs$smoothing.time.diff, ... smoothing.specs$smoothing.end) will also be reported. Note that this means that smoothing.specs$smoothing.start and smoothing.specs$smoothing.end must also be given. Also, if smoothing samples are to be given, smoothing.specs$num.smooth.per.mcmc must be given. This is the number of smoothing samples fetched per MCMC sample. Lastly, smoothing.specs$do.return.smoothing.sample determines if the smoothing samples are to be returned or not. If smoothing.specs$do.smoothing=TRUE, it rarely makes sense to set smoothing.specs$do.return.smoothing.sample<-FALSE, so if this is not given, the algorithm sets smoothing.specs$do.return.smoothing.sample<- smoothing.specs$do.smoothing.

A list that specifies whether realization samples are to be returned and if so how. Realization samples are a posteriori samples of the process itself (including the correlations between time points) at measurement times and possibly also other times. Realization samples can be turned on by setting realization.specs$do.realizations=TRUE. If so, the realization.specs$realization.time.diff and realization.specs$num.realizations must also be given. If realization.specs$realization.time.diff is set to "0", only the states at the measurement times will be inferred. If realization.specs$realization.time.diff is set to a positive real value, the state at the regularly spaced time points (realization.specs$realization.start, realization.specs$realization.start+realization.specs$realization.time.diff, realization.specs$realization.start+2*realization.specs$realization.time.diff, ... realization.specs$realization.end) will also be reported. Note that this means that realization.specs$realization.start and realization.specs$realization.end must also be given. Also, if realization samples are to be given, realization.specs$num.realizations must be given. This is the number of realization samples fetched in total. (This cannot be more than the number of MCMC samples. the reason it is not equal to that, is that when censoring is in use, not all MCMC samples will result in valid realizations.)

Lastly, realization.specs$strategy determines whether censoring of realizations is to take place. Censoring has been used when the process in question represented sedimentation, in order to avoid realizations such that the process would have erased one or more of the measurements. Thus if X is sedimentation depth and X_measurement < X_later_time_point for any combination of measurement or later time point, a censoring would take place. realization.specs$strategy="N" means no censoring. realization.specs$strategy="A" means ascending censoring (censors realizations that anywhere ascends beyond the value of the previous measurement point and realization.specs$strategy="D" means descending censoring (censors realizations that anywhere descends beyond the value of the previous measurement point).

If set to 'TRUE', the standardized residuals for all observed series as a matrix called 'standardized.residuals' (one column for each measured series) and the returns the time points for the residuals ('residuals.time'). Also returns the prior expected values from the Kalman filter ('prior.expected.values'), which should be independent of the residuals if the model assumptions are correct. Note that the residuals are calculated for the Bayesian median parameters if Bayesian analysis is used and maximum likelihood estimated parameters if classic (ML) analysis is used. Also note that the residuals for each observed series is the measurement minus the expected value divided by the uncertainty from the Kalman *filter* not the Kalman smoother (which is used for process inference), as it is the deviations in the Kalman filter that are assumed to be independent from measurement to measurement.

There are sanity checks for each measurement point handled during the log-likelihood calculation, in order to avoid errors in the calculation due to numerical errors. However, some numerical inaccuracies will always exist, so some laxness in these checks must be allowed. Varianles that should be equal might not be exactly equal. And covariance matrices, which should be symmetric, might not be exactly symmetric or be so nearly singular that numerical errors destroys positive definiteness (only checked for levels "low" and "moderate"). Too much laxness and erroneous log-likelihoods may be used. Too much rigidity, i.e. too little laxness and parameter values may needlessly be excluded and the optimal estimation is not achieved. Default laxness is "high", but options are also "low" (very rigid), "moderate" and "complete" (sanity checks are not performed).

PS: The test covariance matrix test for positive definiteness (positive eigenvalues) do take quite a lot of computational resources. If you run with loglik.laxness="complete" or "high" (default), you might run things as much as 3 times faster! But unless you run some fairly simple models, be sure to check both parameter estimates and the stability of the log-likelihood. (Perhaps compare it to a simpler model). loglik.laxness="high" (default), doesn't cost much more than loglik.laxness="complete" (nothing is checked), but doesn't check for positive definiteness.

External series are time series that are treated as explanatory variables for the time series that are analyzed. They are not themselves analyzed statistically. Instead they are taken as is, with the assumption that what happenes between time points can be exppromately attained through linear interpolation. Since the stochastic nature of these series are not analyzed, only causal links are possible. Note that if external time series are given, one must also specify which series and layers they go to (see the next two entries.) Note that external data series are assumed global, and this connect equally to all sites (if that is a thing in your datasets).

An vector of integers that should specify which external time series are causally connected to which observed series, numbered from 1 to the number of series. The length of this vector should match the number of external time series.

An vector of integers that should specify which external time series are causally connected to which layers of the observed series, numbered from 1 to the number of layers in the series specified in 'external.series.conneciton'. The length of this vector should match the number of external time series.

A matrix that specifies the distance between each pair of sites. Only relevant when on series uses the 'distance.correlated.sigma' option to specify between-site correlation based on distance and a distance correlation falloff parameter. But in such a case, the site distance matrix is very necessary! (Keep in mind, sites are numbered from 0 to number of sites-1, so element i,j in the matrix should tell the distance between site i-1 and site j-1. Also keep in mind that the size of the matrix must match the number of sites in the dataset. Also, diagonal elements are the distance from a specific site to itself and so shoulød be 0. Also, the matrix must be square and symmetric, since the distance from site i-1 to site j-1 is the distance from site j-1 to site i-1.)

Alternative connection specification through an array of integers between 0 and 4, one for each possible pairing of two processes (series+layer) belonging to different series (see 'layer.pairs'). For each pair, 0 means no connection between the two, 1 means correlative connection between the two, 2 means causal connection from the first process in the pair to the second, 3 means causal connection from the second process in the pair to the first and 4 means causal connection going both ways between the pair of processes.

Used instead of 'causal' and 'corr' if one wants a compact (but less user-friendly) way to specify connections. This alternative way of specifying a connection was introduced in order to make it easy to produce different connection models by simply traversing the possible connection between each pair of processes. If one use a script to traverse all possible connections models, reads a script and corrects the connection specification, one has a tool for outputing one script for each connection models. These scripts can then possibly be run in parallel on a computer with many cores or a computing cluster.

See the function 'layer.pairs' and 'layer.connections' for ways to make sure that you use this specification correctly.

A list containing 'layer.series.structure' objects that represents the set of time series that are to be analyzed and how the underlying processes of each such series is to be modelles. (See 'layer.series.structure'). At least one such object must be given.

num.MCMC is an integer specifying the number of Markov chain Monte Carlo (MCMC) samples that the analysis will rest on. Default 1000. More MCMC samples are run internally, depending on the spacing (default 10) and burnin (default 1000).

Specifies the number of MCMC iteration between each sample that is used (default 10). A higher number means the samples will be less correlated but requires proptionally more computer time.

Markov chain Monte Carlo sampling takes some iterations to converge to the posterior parameter distribution. The burn-in phase is the set of iterations performed where the samples are then discarded, so as to avoid including samples before convergence. Convergence can be check by for instance 'gelman.diag' in the 'coda' package, and the MCMC samples themselves can be gotten by using the option 'mcmc'. Still, stability of the results can be quite a good enough measure for most purposes. Default is 10000 iterations. PS: The burn-in algorithm here includes two steps where the random walk MCMC variance is adjusted in order to optimize the efficiency of the algorithm.

The number of tempering chains (default 1, meaning it is just the MCMC chain itself that is beeing run). Parallel tempering (Geyer 1991) is a method for overcoming local optima and can thus increase the stability of the results. The way it is done is that in addition to the MCMC chain for the posterior distribution, parallel chains sampling from the distribution proportional to exp(-log(likelihood*prior)/T) where T is a "temperature" larger than one (this is in practise a "smoothed" version of the posterior), is also sampled from. Swaps between neighbouring chains are then allowed, making it possible for the MCMC chain for the posterior distribution to find a new optima.

As the code is not parallel, the computer time will increase proportional to the number of tempering chains. The set of "temperatures" is set to (1, T.ground, T.ground^2, ..., T.ground^(numtemp-1)), where 'T.ground' is default set to 1.5. If no swapping is done (this can be seen if 'silent.mode' is set to 'FALSE'), this "ground temperature" must be lowered, in order to facilitate the swaps.

Logical variable. Only makes sense for Bayesian analysis (which is default). If set to 'TRUE' (default), the Bayesian model (marginal) likelihood is estimated using an importance sampler. If set to 'FALSE', this will not happen, which saves time but means that Bayesian model comparison will not be possible.

Logical variable. If set to 'FALSE' (default), Bayesian analysis is performed. If set to 'TRUE', a classic maximum likelihood (ML) is performed, where the starting points of a set of optimizations is drawns from the MCMC samples. (This means that Bayesian MCMC sampling is performed anyway, but that the results from this is not focused on in the summary). The number of optimizations is determined by option 'maximum.likelihood.numstart', which is default set to 10. Estimation of the Bayesian model (marginal) likelihood is switched off when this option is used.

The number of ML optimizations performed, it applicable. (If 'do.maximum.likelihood' is set to 'TRUE'.) Default:10.

If set to 'FALSE', the routine will show a lot debug information while running. This can be useful for tempering purposes, since the number of tempering swaps are shown in this debug information.

If set to 'TRUE', this triggers the printing of more specific debug information concerning the burn-in phase.

If set to 'TRUE', this triggers the printing of more specific debug information concerning the likelihood calculation phase. (PS: This means a lot of printing.)

Determines how to handle the identifiability problem in Bayesian multi-layer analysis. As mentioned in the supplementary of Reitan et al. (2012), one can switch the characteristic times of neighbouring layers and reorganize the stochastic contributions in such a way that the same top layer process takes place. This can be solved by requiring that the lower the layer, the larger the characteristic time. However, simply requiring this means the prior distribution of the characteristic times get shortened for multi-layered models, as compared to one-layered model. As long as the data is within the range of these shortened prior distributions, this puts the one-layer model at a disadvantage when it comes to model comparison. There are various solutions to this problem, and which solution is used is determined by the variable 'id.strategy'.

The variable 'id.strategy' should be an integer between 0 and 4, default 2. 0 - No identification treatment. (Default) PS: Can and even should yield multimodal characteristic times 1 - Keep upper characteristic time prior. Add lognormally to beneath-lying characteristic times. 2 - Keep lower characteristic time prior. Substract lognormally to above-lying characteristic times. 3 - Keep lower characteristic time prior. Cut depending on that on the above-lying characteristic times. 4 - Keep upper characteristic time prior. Cut depending on that on the below-lying characteristic times.

If set to 'TRUE', instead of reporting the stochastic contribution size in the stochastic differential equations, the stationary standard deviation is reported. The stationary standard deviation is the standard deviation of the process state at any given time after convergence, unconditional on previous states. In a multi-layered model, this is to be interpreted as the stationary standard deviation if this layer when let alone (with no other process affecting it). This is calculated as s*sqrt(characteristic time/2), where s is the stochastic contribution size. PS: This option does not make any sense if the bottom layer is a Wiener process, as that process is not stationary.

This determines how the parallel tempering chains are to be defined. Parallel tempering (Geyer 1991) is a method for overcoming local optima and can thus increase the stability of the results. The way it is done is that in addition to the MCMC chain for the posterior distribution, parallel chains sampling from the distribution proportional to exp(-log(likelihood*prior)/T) where T is a "temperature" larger than one (this is in practise a "smoothed" version of the posterior), is also sampled from. Swaps between neighbouring chains are then allowed, making it possible for the MCMC chain for the posterior distribution to find a new optima.

As the code is not parallel, the computer time will increase proportional to the number of tempering chains. The set of "temperatures" is set to (1, T.ground, T.ground^2, ..., T.ground^(numtemp-1)), where 'T.ground' is default set to 1.5. If no swapping is done (this can be seen if 'silent.mode' is set to 'FALSE'), this "ground temperature" must be lowered, in order to facilitate the swaps.

If set to 'TRUE', reports half-lives rather than characteristic times. Half-lives are the time the auto-correlation of an OU process drops to 1/2, or the time it takes for a perturbation from the expected value to drop to 1/2 the original size in expected value. (As opposed to dropping to exp(-1) for characteristic times.) Half-life=log(2)*characteristic time. The auto-correlation for an OU process (and thus for the inner dynamics of any linear layer) is exp(-diff.time/characteristic.time)=exp(-log(2)*diff.time/half.life)= (1/2)^(diff.time/half.life), where diff.time is the time difference between two process states.

If set to 'TRUE', returns the MCMC samples as an 'mcmc' object within the return object. The name of the list object will also be 'mcmc'.

If multiple time series are given, one can specify causal connections from underlying processes belonging to such time series to the underlying processes belonging to some other time series. This specification should come in the form of an 4 x n matrix, where n is the number of causal connections. For each column (i.e. for each connection), row 1 specifies the number of the cause time series (same as the ordering of the input 'layer.series.structure' input to the routine), row 2 specifies the layer of the cause process, row 3 specifies the number of the effect time series and row 4 specifies the layer of the effect process. So if one has two time series, X and Y, both being modelled as having 3 layers, one can specify a causal connection from layer 3 in X to layer 1 in Y by the matrix "causal<-matrix(c(1,3,2,1),ncol=1)".

Note that it is possible to specify a causal connection even with only one time series, if a multi-layered process is considered. So if X is two-layered, "causal<-matrix(c(1,1,1,2),ncol=1)" will work. As layers are defined so that there is always a causal connection from layer 'k' to layer 'k-1' (with strength=1 for identifiability purposes), the reverse specification ("causal<-matrix(c(1,2,1,1),ncol=1)") will not make sense.

Specifies a set of causal connections going both ways with equal strength. For more on the usage, see option 'causal', but note that symmetric causal connection between one layer and the next does not make sense, as the connection strength from layer 'k' to layer 'k-1' is always set to "1" for identifiability purposes. Note that which series is first and which is second does not matter, so the matrix "causal.symmetric<-matrix(c(2,1,1,3),ncol=1)" would specify the same model as "causal.symmetric<-matrix(c(1,3,2,1),ncol=1)".

A 4 x n matrix which specifies the correlative connections between the processes, where n is the number of such connections. For each column (i.e. for each connection), row 1 specifies the number of the first time series (same as the ordering of the input 'layer.series.structure' input to the routine), row 2 specifies the layer of the first process, row 3 specifies the number of the second time series and row 4 specifies the layer of the second process. So if one has two time series, X and Y, both being modelled as having 3 layers, one can specify a correlative connection between layer 3 in X and layer 1 in Y by the matrix "corr<-matrix(c(1,3,2,1),ncol=1)". Note that which series is first and which is second does not matter, so the matrix "corr<-matrix(c(2,1,1,3),ncol=1)" would specify the exact same model.

A list that specifies whether smoothing samples are to be returned and if so how. Smoothing samples are a posteriori samples of the process states at measurement times and possibly also other times. Smoothing samples can be turned on by setting smoothing.specs$do.smoothing=TRUE. If so, the smoothing.specs$smoothing.time.diff and smoothing.specs$num.smooth.per.mcmc must also be given. If smoothing.specs$smoothing.time.diff is set to "0", only the states at the measurement times will be inferred. If smoothing.specs$smoothing.time.diff is set to a positive real value, the state at the regularly spaced time points (smoothing.specs$smoothing.start, smoothing.specs$smoothing.start+smoothing.specs$smoothing.time.diff, smoothing.specs$smoothing.start+2*smoothing.specs$smoothing.time.diff, ... smoothing.specs$smoothing.end) will also be reported. Note that this means that smoothing.specs$smoothing.start and smoothing.specs$smoothing.end must also be given. Also, if smoothing samples are to be given, smoothing.specs$num.smooth.per.mcmc must be given. This is the number of smoothing samples fetched per MCMC sample. Lastly, smoothing.specs$do.return.smoothing.sample determines if the smoothing samples are to be returned or not. If smoothing.specs$do.smoothing=TRUE, it rarely makes sense to set smoothing.specs$do.return.smoothing.sample<-FALSE, so if this is not given, the algorithm sets smoothing.specs$do.return.smoothing.sample<- smoothing.specs$do.smoothing.

A list that specifies whether realization samples are to be returned and if so how. Realization samples are a posteriori samples of the process itself (including the correlations between time points) at measurement times and possibly also other times. Realization samples can be turned on by setting realization.specs$do.realizations=TRUE. If so, the realization.specs$realization.time.diff and realization.specs$num.realizations must also be given. If realization.specs$realization.time.diff is set to "0", only the states at the measurement times will be inferred. If realization.specs$realization.time.diff is set to a positive real value, the state at the regularly spaced time points (realization.specs$realization.start, realization.specs$realization.start+realization.specs$realization.time.diff, realization.specs$realization.start+2*realization.specs$realization.time.diff, ... realization.specs$realization.end) will also be reported. Note that this means that realization.specs$realization.start and realization.specs$realization.end must also be given. Also, if realization samples are to be given, realization.specs$num.realizations must be given. This is the number of realization samples fetched in total. (This cannot be more than the number of MCMC samples. the reason it is not equal to that, is that when censoring is in use, not all MCMC samples will result in valid realizations.)

Lastly, realization.specs$strategy determines whether censoring of realizations is to take place. Censoring has been used when the process in question represented sedimentation, in order to avoid realizations such that the process would have erased one or more of the measurements. Thus if X is sedimentation depth and X_measurement < X_later_time_point for any combination of measurement or later time point, a censoring would take place. realization.specs$strategy="N" means no censoring. realization.specs$strategy="A" means ascending censoring (censors realizations that anywhere ascends beyond the value of the previous measurement point and realization.specs$strategy="D" means descending censoring (censors realizations that anywhere descends beyond the value of the previous measurement point).

If set to 'TRUE', the standardized residuals for all observed series as a matrix called 'standardized.residuals' (one column for each measured series) and the returns the time points for the residuals ('residuals.time'). Also returns the prior expected values from the Kalman filter ('prior.expected.values'), which should be independent of the residuals if the model assumptions are correct. Note that the residuals are calculated for the Bayesian median parameters if Bayesian analysis is used and maximum likelihood estimated parameters if classic (ML) analysis is used. Also note that the residuals for each observed series is the measurement minus the expected value divided by the uncertainty from the Kalman *filter* not the Kalman smoother (which is used for process inference), as it is the deviations in the Kalman filter that are assumed to be independent from measurement to measurement.

If set true, performs smoothing (process inference) on a previously performed MCMC run, defined in the subsecquent item, 'previous.run'. Does not perform an new MCMC analysis. Note that the new data can still be given, so that the inference is for the process parameters as revealed by the previously run MCMC sampling, but with data equal to the new input. Note that the 'smoothing.specs' option also has to be set, in order for this to make sense.

A 'layered' object representing a previous run analysis, where MCMC samples have been saved (using the 'mcmc=True' option). Only makes sense if the 'smooth.previous.run' has been set to TRUE, otherwise the object is irrelevant.

Mode literally set the mode for the call to the C++ code. Options are "Bayes","ML-from-MCMC","ML-from-input","ML-from-hybrid", "Loglik-from-input","Smooth-from-input" and "Numpar". "Bayes" is Bayesian analysis, "ML-from-MCMC" runs ML optimizations starting from an MCMC run. "ML-from-input" instead starts the ML optimizations from the MCMC samples in a previous run (this also allows for explicitely stating the starting points, as one can set the "MCMC" samples from the previous run), while "ML-from-hybrid" starts alternating between MCMC samples fromt he current analysis and the previous run. "Loglik-from-input" is used by the layer.param.loglik(s) methods. "Smooth-from-input" is used by the layer.predict.mcmc methods. "Numpar" is used for extracting just the number of parameters (and the parameter names), which is useful if you need to know the number of parameters before running an anlysis (used in 'layer.analyzer' in case the user wants to start ML optimizations from a previous run without connections). Note that mode needs to match input. So for instance, "Bayes" and "ML-from-MCMC" have no need or want for a previous run, while "ML-from-input","ML-from-hybrid","Loglik-from-input" and "Smooth-from-input" do.

There are sanity checks for each measurement point handled during the log-likelihood calculation, in order to avoid errors in the calculation due to numerical errors. However, some numerical inaccuracies will always exist, so some laxness in these checks must be allowed. Varianles that should be equal might not be exactly equal. And covariance matrices, which should be symmetric, might not be exactly symmetric or be so nearly singular that numerical errors destroys positive definiteness (only checked for levels "low" and "moderate"). Too much laxness and erroneous log-likelihoods may be used. Too much rigidity, i.e. too little laxness and parameter values may needlessly be excluded and the optimal estimation is not achieved. Default laxness is "high", but options are also "low" (very rigid), "moderate" and "complete" (sanity checks are not performed).

PS: The test covariance matrix test for positive definiteness (positive eigenvalues) do take quite a lot of computational resources. If you run with loglik.laxness="complete" or "high" (default), you might run things as much as 3 times faster! But unless you run some fairly simple models, be sure to check both parameter estimates and the stability of the log-likelihood. (Perhaps compare it to a simpler model). loglik.laxness="high" (default), doesn't cost much more than loglik.laxness="complete" (nothing is checked), but doesn't check for positive definiteness.

External series are time series that are treated as explanatory variables for the time series that are analyzed. They are not themselves analyzed statistically. Instead they are taken as is, with the assumption that what happenes between time points can be exppromately attained through linear interpolation. Since the stochastic nature of these series are not analyzed, only causal links are possible. Note that if external time series are given, one must also specify which series and layers they go to (see the next two entries.) Note that external data series are assumed global, and this connect equally to all sites (if that is a thing in your datasets).

An vector of integers that should specify which external time series are causally connected to which observed series, numbered from 1 to the number of series. The length of this vector should match the number of external time series.

An vector of integers that should specify which external time series are causally connected to which layers of the observed series, numbered from 1 to the number of layers in the series specified in 'external.series.conneciton'. The length of this vector should match the number of external time series.

A matrix that specifies the distance between each pair of sites. Only relevant when on series uses the 'distance.correlated.sigma' option to specify between-site correlation based on distance and a distance correlation falloff parameter. But in such a case, the site distance matrix is very necessary! (Keep in mind, sites are numbered from 0 to number of sites-1, so element i,j in the matrix should tell the distance between site i-1 and site j-1. Also keep in mind that the size of the matrix must match the number of sites in the dataset. Also, diagonal elements are the distance from a specific site to itself and so shoulød be 0. Also, the matrix must be square and symmetric, since the distance from site i-1 to site j-1 is the distance from site j-1 to site i-1.)

A set of 'layer.series.structure' objects that represents the set of time series that are to be analyzed and how the underlying processes of each such series is to be modelles. (See 'layer.series.structure'). At least two such objects must be given.

An array of integers between 0 and 4, one for each possible pairing of two processes (series+layer) belonging to different series (see 'layer.pairs'). For each pair, 0 means no connection between the two, 1 means correlative connection between the two, 2 means causal connection from the first process in the pair to the second, 3 means causal connection from the second process in the pair to the first and 4 means causal connection going both ways between the pair of processes.

A list containing 'layer.series.structure' objects that represents the set of time series that are to be analyzed and how the underlying processes of each such series is to be modelles. (See 'layer.series.structure'). At least two such objects must be given.

An array of integers between 0 and 4, one for each possible pairing of two processes (series+layer) belonging to different series (see 'layer.pairs'). For each pair, 0 means no connection between the two, 1 means correlative connection between the two, 2 means causal connection from the first process in the pair to the second, 3 means causal connection from the second process in the pair to the first and 4 means causal connection going both ways between the pair of processes.

The time points of the time series in chronological (ascending) order. Does not need

The corresponding measurements of the time series (or the average measurement if more than one goes into the value, if 'num.meas.per.value' is given and larger than 1). Must necessarily be the same size as 'time.points'.

The name of the time series. Will be used in assigning parameter names for parameters associated with this series.

The individual measurement standard deviations of each measurement (if applicable). If this does not exist, the model will estimate an overall measurement standard deviation belonging to all measurements. If 'num.meas.per.value', these are the standard deviations of of the individual measurements going into each value. Must necessarily be the same size as 'time.points', if given.

The number of each individual measurement that goes into each value in the time series. Must necessarily be the same size as 'time.points', if given.

The site number of each measurement, if applicable. Should be labelled from 0 to <number of sites>-1. Must necessarily be the same size as 'time.points', if given.

The name of the semi-colon-separated csv file that contains the prior specification.

A set of 'layer.series.structure' objects that represents the set of time series that are to be analyzed and how the underlying processes of each such series is to be modelles. (See 'layer.series.structure'). At least two such objects must be given.

A list containing 'layer.series.structure' objects that represents the set of time series that are to be analyzed and how the underlying processes of each such series is to be modelles. (See 'layer.series.structure'). At least two such object must be given.

The 'analysis' variable must be a 'layered' object. I.e. a 'layer.analyzer' analysis needs to be performed first. This variable is not an option. If it is NULL, the function will complain. it is from here the function gets the model specification (including conneciton elements in the case of multiple time series). Thus a small analysis should have been performed alread. Alternatively, one can try to fill in a 'layered' object by hand, but that takes effort.

Specifies the parameter values for the log-likelihood calculation. If this is not specified (NULL), then the parameter estimates from the analysis is used instead. For the order of the parameters, see for instance analysis$parameter.names. Note that the scale is the original scale, not the reparametrized scale used internally (where for instance stochastic contributions, characteristic times and observational noise level is log-trnasformed, and correlation is logit-transformed.) It is possible to set some of the parameter values to NA, in which case these will be optimized over.

If set to FALSE, shows debug info for each iteration in the Kalman filter.

If some parameter values are set to NA, optimization is performed on those parameters. 'num.optim' specifies the number of optimizations performed on each such parameter (since optimization is performed by hill-climbing and can thus converge to different optima from different starting positions).

Specifies that an MCMC analysis should take palce before doing the log-likleihood calculation(s). Only a point if some of the parameters are unspecified. If so, and do.preanalysis.mcmc is set to TRUE, then the first optimization will use the parameter values from closest MCMC samples closest to the specified parts of the parameter set.

num.MCMC is an integer specifying the number of Markov chain Monte Carlo (MCMC) samples that the pre-analysis (see do.preanalysis.mcmc) will rest on. Default 1000. More MCMC samples are run internally, depending on the spacing (default 10) and burnin (default 1000).

Specifies the number of MCMC iteration (see do.preanalysis.mcmc) between each sample that is used (default 10). A higher number means the samples will be less correlated but requires proptionally more computer time.

Markov chain Monte Carlo sampling (see do.preanalysis.mcmc) takes some iterations to converge to the posterior parameter distribution.

The number of tempering chains (default 1, meaning it is just the MCMC chain (see do.preanalysis.mcmc)itself that is beeing run). Parallel tempering (Geyer 1991) is a method for overcoming local optima and can thus increase the stability of the results.

The 'analysis' variable must be a 'layered' object. I.e. a 'layer.analyzer' analysis needs to be performed first. This variable is not an option. If it is NULL, the function will complain. it is from here the function gets the model specification (including conneciton elements in the case of multiple time series). Thus a small analysis should have been performed alread. Alternatively, one can try to fill in a 'layered' object by hand, but that takes effort.

Specifies the set of parameter value sets for the log-likelihood calculations. For the order of the parameters, see for instance analysis$parameter.names. Note that the scale is the original scale, not the reparametrized scale used internally (where for instance stochastic contributions, characteristic times and observational noise level is log-trnasformed, and correlation is logit-transformed.) It is possible to set some of the parameter values to NA, in which case these will be optimized over.

If set to FALSE, shows some debug info.

If some parameter values are set to NA, optimization is performed on those parameters. 'num.optim' specifies the number of optimizations performed on each such parameter (since optimization is performed by hill-climbing and can thus converge to different optima from different starting positions).

Specifies that an MCMC analysis should take palce before doing the log-likleihood calculation(s). Only a point if some of the parameters are unspecified. If so, and do.preanalysis.mcmc is set to TRUE, then the first optimization will use the parameter values from closest MCMC samples closest to the specified parts of the parameter set.

num.MCMC is an integer specifying the number of Markov chain Monte Carlo (MCMC) samples that the pre-analysis (see do.preanalysis.mcmc) will rest on. Default 1000. More MCMC samples are run internally, depending on the spacing (default 10) and burnin (default 1000).

Specifies the number of MCMC iteration (see do.preanalysis.mcmc) between each sample that is used (default 10). A higher number means the samples will be less correlated but requires proptionally more computer time.

Markov chain Monte Carlo sampling (see do.preanalysis.mcmc) takes some iterations to converge to the posterior parameter distribution.

The number of tempering chains (default 1, meaning it is just the MCMC chain (see do.preanalysis.mcmc)itself that is beeing run). Parallel tempering (Geyer 1991) is a method for overcoming local optima and can thus increase the stability of the results.

A set of 'layer.data.series' objects that represents the set of time series that are to be predicted on. Note that these must match the ones used for running 'layer.analyzer' when it comes to number of time series and order. (See 'layer.data.series' for how they are defined).

The 'analysis' variable must be a 'layer.analyzer' object. I.e. a 'layer.analyzer' analysis needs to be performed first. This variable is not an option. If it is NULL, the function will complain.

If 'smoothing.time.diff' is set to "0", only the states at the measurement times will be inferred. (PS: you can get the function to perform process state inference for unmeasured time points by injecting them into the time series with' missing values). If 'smoothing.time.diff' is set to a positive real value, the state at the regularly spaced time points (smoothing.start, smoothing.start+smoothing.time.diff, smoothing.start+smoothing.time.diff, ... smoothing.end) will also be reported. Note that this means that 'smoothing.start' and 'smoothing.end' must also be given.

Start time for extra smoothing points.

End time for extra smoothing points.

The difference between predicted and actual measurements are returned. If set to 'TRUE', the standardized residuals for all observed series as a matrix called 'standardized.residuals' (one column for each measured series) and the returns the time points for the residuals ('residuals.time'). Also returns the prior expected values from the Kalman filter ('prior.expected.values'), which should be independent of the residuals if the model assumptions are correct. Note that the residuals are calculated for the Bayesian median parameters if Bayesian analysis is used and maximum likelihood estimated parameters if classic (ML) analysis is used. Also note that the residuals for each observed series is the measurement minus the expected value divided by the uncertainty from the Kalman *filter* not the Kalman smoother (which is used for process inference), as it is the deviations in the Kalman filter that are assumed to be independent from measurement to measurement.

A set of 'layer.data.series' objects that represents the set of time series that are to be predicted on. Note that these must match the ones used for running 'layer.analyzer' when it comes to number of time series and order. (See 'layer.data.series' for how they are defined).

The 'analysis' variable must be a 'layer.analyzer' object. I.e. a 'layer.analyzer' analysis needs to be performed first. This variable is not an option. If it is NULL, the function will complain.

If 'smoothing.time.diff' is set to "0", only the states at the measurement times will be inferred. (PS: you can get the function to perform process state inference for unmeasured time points by injecting them into the time series with' missing values). If 'smoothing.time.diff' is set to a positive real value, the state at the regularly spaced time points (smoothing.start, smoothing.start+smoothing.time.diff, smoothing.start+smoothing.time.diff, ... smoothing.end) will also be reported. Note that this means that 'smoothing.start' and 'smoothing.end' must also be given.

Start time for extra smoothing points.

End time for extra smoothing points.

This is the number of smoothing process state samples fetched per MCMC sample.

'do.return.smoothing.sample' determines if the raw smoothing samples are to be returned or not, or only the mean and lower/upper 95% credibility band of the processes. These represents samples of the individual process states given the new data. Note that these are not process realizations, as the autocorrelation between one time point and the next is not taken into account when sampling. Process realizations may perhaps be added later.

A list item of 'layer.data.series' objects that represents the set of time series that are to be predicted on. Note that these must match the ones used for running 'layer.analyzer' when it comes to number of time series and order. (See 'layer.data.series' for how they are defined).

The 'analysis' variable must be a 'layer.analyzer' object. I.e. a 'layer.analyzer' analysis needs to be performed first. This variable is not an option. If it is NULL, the function will complain.

If 'smoothing.time.diff' is set to "0", only the states at the measurement times will be inferred. (PS: you can get the function to perform process state inference for unmeasured time points by injecting them into the time series with' missing values). If 'smoothing.time.diff' is set to a positive real value, the state at the regularly spaced time points (smoothing.start, smoothing.start+smoothing.time.diff, smoothing.start+smoothing.time.diff, ... smoothing.end) will also be reported. Note that this means that 'smoothing.start' and 'smoothing.end' must also be given.

Start time for extra smoothing points.

End time for extra smoothing points.

This is the number of smoothing process state samples fetched per MCMC sample.

'do.return.smoothing.sample' determines if the raw smoothing samples are to be returned or not, or only the mean and lower/upper 95% credibility band of the processes. These represents samples of the individual process states given the new data. Note that these are not process realizations, as the autocorrelation between one time point and the next is not taken into account when sampling. Process realizations may perhaps be added later.

The difference between predicted and actual measurements are returned. If set to 'TRUE', the standardized residuals for all observed series as a matrix called 'standardized.residuals' (one column for each measured series) and the returns the time points for the residuals ('residuals.time'). Also returns the prior expected values from the Kalman filter ('prior.expected.values'), which should be independent of the residuals if the model assumptions are correct. Note that the residuals are calculated for the Bayesian median parameters if Bayesian analysis is used and maximum likelihood estimated parameters if classic (ML) analysis is used. Also note that the residuals for each observed series is the measurement minus the expected value divided by the uncertainty from the Kalman *filter* not the Kalman smoother (which is used for process inference), as it is the deviations in the Kalman filter that are assumed to be independent from measurement to measurement.

An array with 2 elements that specifies the low and high end of the 95% credibility band of mu, the expected value for the process.

An array with 2 elements that specifies the low and high end of the 95% credibility band of dt, the characteristic time (or the half-life if that is specified).

An array with 2 elements that specifies the low and high end of the 95% credibility band of sigma, the stochastic contribution.

An array with 2 elements that specifies the low and high end of the 95% credibility band of the stationary standard deviation. If not specified, the specifications for the prior of sigma is copied.

An array with 2 elements that specifies the low and high end of the 95% credibility band of init, the initical value of the process.

An array with 2 elements that specifies the low and high end of the 95% credibility band of lin, the linear trend (if applicable).

An array with 2 elements that specifies the low and high end of the 95% credibility band of beta, the strength of causal links (if applicable).

An array with 2 elements that specifies the low and high end of the 95% credibility band of obs, the standard deviation of the observational errors.

If distance-based spatial correlation on sites is wanted, there will be a parameter describing how the correlation decreases (exponentially) with distance, the distance correlation falloff parameter. The prior distribution for this is log-normal with 95 if not specified is set extremely wide, going from 1e-6 to 1e+6.

If set to 1, the analysis will log-transform the data before analysis. If set to 2, the analysis will assume the data has already been log-transformed.

A 'layer.data.series' object representing a measurement time series.

The number of layers in the model (hidden layers + measured layer). Setting numlayers=1 and nothing more means modellings the process behind the time series as an OU process. Note that increasing the number of layers can dramatically increase the computer resources needed (by as much as the cube of the number of layers), since matrix operations are needed in order to update knowledge of all layered processes.

If set to 'TRUE', this indicates that a linear time trend is to be included. The bottom layer must be an OU process if so (i.e. 'no.pull' must be set to "False").

An array of integers specifying which layers are the time integral of the layer below it. In differential equation terms, that means dX_n=X_(n+1) where 'n' is the layer number in question. This option was introduced in order to model sedimention height by letting it be the time integral of sedimentation rate and letting sedimentation rate be an OU process. Usually, only one time integral layer is needed, but the analysis allows for multiple such layers. Note that the bottom layer cannot be a time integral layer, as it is not tracking anything, so 'time.integral' must be less that 'numlayers'.

If set to 'TRUE', it specifies that the bottom layers has no pull term, i.e. that it is a Wiener process (random walk) rather than an OU process. When such models are included, it pays to use the initial value options, see the 'init' options in this routine. If not, a very wide generic prior for the initial value is used, giving the Wiener process model an unfair disadvantage.

An array of integers indicating which layers are modelled to have no stochastic contributions. That means each layer specified will be a deterministic tracking of the layer below and will in practise be a smoothed version of the underlying process. because of this, it does not make sense setting 'no.sigma' to 'numlayers', as the bottom layers needs to have stoachasticity in order to have any dynamics at all.

Should be a 'layer.prior' object that specifies the prior distribution of the parameters, see the 'layer.prior' documentation. If not given, a standard prior is used. (Warning: This will not be appropriate for all circumstances).

If given, should specify an array of period lengths, in order to include in the model periodic sine/cosine fluctuations (such as seasonality). Regression parameters for the sine/cosine terms are estimated.

If set to 'TRUE', the model includes a parameter for the initial values of each process (one for each combination of site and layer). The initial time is that of the first measurement (unless 'init.time' is specified).

If set, the model includes a parameter for the initial values of each process (one for each combination of site and layer). The initial time is that of 'init.time'.

Only makes sense for time series with multiple sites. If set, an initial value treatment will assume the same initial value for all sites (if applicable).

If set, an initial value treatment will assume the same initial value for all layers (if applicable, i.e. if 'numlayers'>1).

If given, should be an array of two elements, initial time and initial value.

Only makes sense for time series with multiple sites. If set to 'TRUE', specifies that the process of each site should have a separate expected value (mu).

Only makes sense for time series with multiple sites. If set to 'TRUE', specifies that the process of each site should have a separate linear time trend.

Only makes sense for time series with multiple sites. Should be an integer array specifying which layers have separate pulls (and thus characteristic times/halflives and stationary standard deviations) for each site.

Only makes sense for time series with multiple sites. Should be an integer array specifying which layers have separate stochastic noise size (and thus stationary standard deviations) for each site.

Only makes sense for time series with multiple sites. Should be an array of integers specifying which layers have correlated stochastic contributions between the sites. The same correlation is assumed between any two sites.

Only makes sense for time series with multiple sites. Should be an array of integers specifying which layers have correlated stochastic contributions between the sites. Correlation between a pair of sites depend on the distance between them, through a distance falloff parameter. A distance matrix needs to be supplied!

Only makes sense for time series with multiple sites. Should be an array of integers specifying which layers have correlated stochastic contributions between the sites. Different correlations are allowed between different pairs of sites.

Only makes sense for time series with multiple sites. Should be an array of integers specifying which layers have perfect correlated stochastic contributions (meaning correlation=1, i.e. the same stochastic contribution for all sites) between the sites. Can introduce numerical difficulties.

Only makes sense for time series with multiple sites. Should be an array of integers specifying which layers have correlated stochastic contributions between the sites. However, the sites will be put into two groups, with no correlation between sites beloning to one group and the other. The grouping will be inferred. Not well tested.

Only makes sense for time series with multiple sites. Should be an array of integers specifying which layers have a grouping between sites with and without correlation to the stochastic contributions of other sites. Not well tested.

Only makes sense for time series with multiple sites. Attempts to find two groups of sites for the stochastic contribution sizes for the given set of layers, as specified by an array of integers.

Only makes sense for time series with multiple sites. Attempts to find two groups of sites for the pull (characteristic times/halftimes) for the given set of layers, as specified by an array of integers.

Only makes sense for time series with multiple sites. If set to 'True', attempts to find two groups of sites for the expected value (mu) of the processes.

Only makes sense for time series with multiple sites. If set to 'True', attempts to find two groups of sites for the linear time trend.

If set to 'TRUE', allows for positive pull parameters in the linear stochastic equations. Positive pulls means the process is prone to "explode", but for a short time (typically not much more than the characteristic time), such a process could none the less make sense.

A 'layer.series.structure' object that should be described.

Objects of class 'layered', as returned from the 'layer.analyzer' function.

The number of objects to be combined.

The span of the possible object values, going from 0,1,...,base.number. Thus, base.number=2 (default) only allows "0" and "1".

Objects of class 'layered', as returned from the 'layer.analyzer' function.

A 'layer.series.structure' object for which the description should be printed.

Not used, only there to keep complicance with print statements.

An 'summary.layered' object, as returned from the 'summary.layered' method.

Not used. Only included for compatibility with other other 'print.summary' methods.

Name of the input file.

The name given to the series. Will be used in assigning parameter names for parameters associated with this series.

This should specify the content of each column in the file. The file must contain a "time" and a "value" column, as that is the bare minimum in order to specify a time series. In addition, the a measurement standard deviation column can be specified by "std.dev". If in addition a "num.meas.per.value" column exists, the standard error of each mueasurement will be taken to be the content of the "std.dev" column divied by the square root of the content of the "num.meas.per.value" column. A "site" column can also be added, if the time series contains the same typ of measurements from different sites (note that regional options then becomes a possibility). The site index must run from 0 to the number of sites minus one.

This and the rest of the options are just sent to "read.table", so the help text has just been lifted from the manual of that function. You should check that manual for an updated version of this help text.

A logical value indicating whether the file contains the names of the variables as its first line. If missing, the value is determined from the file format: ‘header’ is set to ‘TRUE’ if and only if the first row contains one fewer field than the number of columns.

the character used in the file for decimal points.

the field separator character. Values on each line of the file are separated by this character. If ‘sep = ""’ (the default for ‘read.table’) the separator is ‘white space’, that is one or more spaces, tabs, newlines or carriage returns.

the set of quoting characters. To disable quoting altogether, use ‘quote = ""’. See ‘scan’ for the behaviour on quotes embedded in quotes. Quoting is only considered for columns read as character, which is all of them unless ‘colClasses’ is specified.

string indicating how to convert numbers whose conversion to double precision would lose accuracy, see ‘type.convert’. Can be abbreviated. (Applies also to complex-number inputs.)

a vector of row names. This can be a vector giving the actual row names, or a single number giving the column of the table which contains the row names, or character string giving the name of the table column containing the row names.

If there is a header and the first row contains one fewer field than the number of columns, the first column in the input is used for the row names. Otherwise if ‘row.names’ is missing, the rows are numbered.

Using ‘row.names = NULL’ forces row numbering. Missing or ‘NULL’ ‘row.names’ generate row names that are considered to be ‘automatic’ (and not preserved by ‘as.matrix’).

a vector of optional names for the variables. The default is to use ‘"V"’ followed by the column number.

a character vector of strings which are to be interpreted as ‘NA’ values. Blank fields are also considered to be missing values in logical, integer, numeric and complex fields. Note that the test happens _after_ white space is stripped from the input, so ‘na.strings’ values may need their own white space stripped in advance.

character. A vector of classes to be assumed for the columns. If unnamed, recycled as necessary. If named, names are matched with unspecified values being taken to be ‘NA’.

Possible values are ‘NA’ (the default, when ‘type.convert’ is used), ‘"NULL"’ (when the column is skipped), one of the atomic vector classes (logical, integer, numeric, complex, character, raw), or ‘"factor"’, ‘"Date"’ or ‘"POSIXct"’. Otherwise there needs to be an ‘as’ method (from package ‘methods’) for conversion from ‘"character"’ to the specified formal class.

Note that ‘colClasses’ is specified per column (not per variable) and so includes the column of row names (if any).

integer: the maximum number of rows to read in. Negative and other invalid values are ignored.

integer: the number of lines of the data file to skip before beginning to read data.

logical. If ‘TRUE’ then the names of the variables in the data frame are checked to ensure that they are syntactically valid variable names. If necessary they are adjusted (by ‘make.names’) so that they are, and also to ensure that there are no duplicates.

logical. If ‘TRUE’ then in case the rows have unequal length, blank fields are implicitly added. See ‘Details’.

logical. Used only when ‘sep’ has been specified, and allows the stripping of leading and trailing white space from unquoted ‘character’ fields (‘numeric’ fields are always stripped). See ‘scan’ for further details (including the exact meaning of ‘white space’), remembering that the columns may include the row names.

logical: if ‘TRUE’ blank lines in the input are ignored.

character: a character vector of length one containing a single character or an empty string. Use "" to turn off the interpretation of comments altogether.

logical. Should C-style escapes such as "\n" be processed or read verbatim (the default)? Note that if not within quotes these could be interpreted as a delimiter (but not as a comment character). For more details see ‘scan’.

logical: if "TRUE", "scan" will flush to the end of the line after reading the last of the fields requested. This allows putting comments after the last field.

logical: should character vectors be converted to factors? Note that this is overridden by ‘as.is’ and ‘colClasses’, both of which allow finer control.

the default behavior of ‘read.table’ is to convert character variables (which are not converted to logical, numeric or complex) to factors. The variable ‘as.is’ controls the conversion of columns not otherwise specified by ‘colClasses’. Its value is either a vector of logicals (values are recycled if necessary), or a vector of numeric or character indices which specify which columns should not be converted to factors.

Note: to suppress all conversions including those of numeric columns, set ‘colClasses = "character"’.

Note that ‘as.is’ is specified per column (not per variable) and so includes the column of row names (if any) and any columns to be skipped.

character string: if non-empty declares the encoding used on a file (not a connection) so the character data can be re-encoded. See the ‘Encoding’ section of the help for ‘file’, the ‘R Data Import/Export Manual’ and ‘Note’.

encoding to be assumed for input strings. It is used to mark character strings as known to be in Latin-1 or UTF-8 (see ‘Encoding’): it is not used to re-encode the input, but allows R to handle encoded strings in their native encoding (if one of those two). See ‘Value’ and ‘Note’.

character string: if ‘file’ is not supplied and this is, then data are read from the value of ‘text’ via a text connection. Notice that a literal string can be used to include (small) data sets within R code.

logical: should nuls be skipped?

A set of 'layer.series.structure' objects that represents the set of time series that are to be analyzed and how the underlying processes of each such series is to be modelles. (See 'layer.series.structure'). At least one such object must be given. 'traverse.standalone.layered' can be used for settling the structure of each time series.

num.MCMC is an integer specifying the number of Markov chain Monte Carlo (MCMC) samples that the analysis will rest on. Default 1000. More MCMC samples are run internally, depending on the spacing (default 10) and burnin (default 1000).

Specifies the number of MCMC iteration between each sample that is used (default 10). A higher number means the samples will be less correlated but requires proptionally more computer time.

Markov chain Monte Carlo sampling takes some iterations to converge to the posterior parameter distribution. The burn-in phase is the set of iterations performed where the samples are then discarded, so as to avoid including samples before convergence. Convergence can be check by for instance 'gelman.diag' in the 'coda' package, and the MCMC samples themselves can be gotten by using the option 'mcmc'. Still, stability of the results can be quite a good enough measure for most purposes. Default is 10000 iterations. PS: The burn-in algorithm here includes two steps where the random walk MCMC variance is adjusted in order to optimize the efficiency of the algorithm.

The number of tempering chains (default 1, meaning it is just the MCMC chain itself that is beeing run). Parallel tempering (Geyer 1991) is a method for overcoming local optima and can thus increase the stability of the results. The way it is done is that in addition to the MCMC chain for the posterior distribution, parallel chains sampling from the distribution proportional to exp(-log(likelihood*prior)/T) where T is a "temperature" larger than one (this is in practise a "smoothed" version of the posterior), is also sampled from. Swaps between neighbouring chains are then allowed, making it possible for the MCMC chain for the posterior distribution to find a new optima.

As the code is not parallel, the computer time will increase proportional to the number of tempering chains. The set of "temperatures" is set to (1, T.ground, T.ground^2, ..., T.ground^(numtemp-1)), where 'T.ground' is default set to 1.5. If no swapping is done (this can be seen if 'silent.mode' is set to 'FALSE'), this "ground temperature" must be lowered, in order to facilitate the swaps.

Logical variable. If set to 'FALSE' (default), Bayesian analysis is performed. If set to 'TRUE', a classic maximum likelihood (ML) is performed, where the starting points of a set of optimizations is drawns from the MCMC samples. (This means that Bayesian MCMC sampling is performed anyway, but that the results from this is not focused on in the summary). The number of optimizations is determined by option 'maximum.likelihood.numstart', which is default set to 10. Estimation of the Bayesian model (marginal) likelihood is switched off when this option is used.

The number of ML optimizations performed, it applicable. (If 'do.maximum.likelihood' is set to 'TRUE'.) Default:10.

If set to 'FALSE', the routine will show a lot debug information while running. This can be useful for tempering purposes, since the number of tempering swaps are shown in this debug information.

If set to 'TRUE', this triggers the printing of more specific debug information concerning the burn-in phase.

If set to 'TRUE', this triggers the printing of more specific debug information concerning the likelihood calculation phase. (PS: This means a lot of printing.)

If set to 'TRUE'(default), the connection pairs and the causal and correlative connections of each model.

If set to 'TRUE', the analyses themselves are not performed. Only the connection pairs and the causal and correlative connections of each model are shown. This can be useful if one wants to know how many models are are to be analyzed.

Determines how to handle the identifiability problem in Bayesian multi-layer analysis. As mentioned in the supplementary of Reitan et al. (2012), one can switch the characteristic times of neighbouring layers and reorganize the stochastic contributions in such a way that the same top layer process takes place. This can be solved by requiring that the lower the layer, the larger the characteristic time. However, simply requiring this means the prior distribution of the characteristic times get shortened for multi-layered models, as compared to one-layered model. As long as the data is within the range of these shortened prior distributions, this puts the one-layer model at a disadvantage when it comes to model comparison. There are various solutions to this problem, and which solution is used is determined by the variable 'id.strategy'.

The variable 'id.strategy' should be an integer between 0 and 4, default 2. 0 - No identification treatment. (Default) PS: Can and even should yield multimodal characteristic times 1 - Keep upper characteristic time prior. Add lognormally to beneath-lying characteristic times. 2 - Keep lower characteristic time prior. Substract lognormally to above-lying characteristic times. 3 - Keep lower characteristic time prior. Cut depending on that on the above-lying characteristic times. 4 - Keep upper characteristic time prior. Cut depending on that on the below-lying characteristic times.

If set to 'TRUE', instead of reporting the stochastic contribution size in the stochastic differential equations, the stationary standard deviation is reported. The stationary standard deviation is the standard deviation of the process state at any given time after convergence, unconditional on previous states. In a multi-layered model, this is to be interpreted as the stationary standard deviation if this layer when let alone (with no other process affecting it). This is calculated as s*sqrt(characteristic time/2), where s is the stochastic contribution size. PS: This option does not make any sense if the bottom layer is a Wiener process, as that process is not stationary.

This determines how the parallel tempering chains are to be defined. Parallel tempering (Geyer 1991) is a method for overcoming local optima and can thus increase the stability of the results. The way it is done is that in addition to the MCMC chain for the posterior distribution, parallel chains sampling from the distribution proportional to exp(-log(likelihood*prior)/T) where T is a "temperature" larger than one (this is in practise a "smoothed" version of the posterior), is also sampled from. Swaps between neighbouring chains are then allowed, making it possible for the MCMC chain for the posterior distribution to find a new optima.

As the code is not parallel, the computer time will increase proportional to the number of tempering chains. The set of "temperatures" is set to (1, T.ground, T.ground^2, ..., T.ground^(numtemp-1)), where 'T.ground' is default set to 1.5. If no swapping is done (this can be seen if 'silent.mode' is set to 'FALSE'), this "ground temperature" must be lowered, in order to facilitate the swaps.

If set to 'TRUE', reports half-lives rather than characteristic times. Half-lives are the time the auto-correlation of an OU process drops to 1/2, or the time it takes for a perturbation from the expected value to drop to 1/2 the original size in expected value. (As opposed to dropping to exp(-1) for characteristic times.) Half-life=log(2)*characteristic time. The auto-correlation for an OU process (and thus for the inner dynamics of any linear layer) is exp(-diff.time/characteristic.time)=exp(-log(2)*diff.time/half.life)= (1/2)^(diff.time/half.life), where diff.time is the time difference between two process states.

If set to 'TRUE', returns the MCMC samples as an 'mcmc' object within the return object. The name of the list object will also be 'mcmc'.

If set to 'FALSE', no causal connection models will be analyzed, only correlative ones.

If set to 'FALSE', no correlative conneciton models will be analyzed, only causal ones. (Incompatible with 'allow.causal=FALSE').

If set to 'FALSE', no direct causal loops will be examined, i.e. causal links from series 1, layer 1 to series 2, layer 2 and another liunk back again. (Does not rule out more elaborate causal loops.)

If prior probabilities (item 'p0') is not given, the prior probabilities are default equal for all models. However, if first.is.nullhypothesis is set to 'TRUE', the first model is given 50% of the prior probability, while the rest of the models shares the remaining 50% equally. This variable is sent to 'compared.layered' for model comparison.

If classic ML analysis is performed instead of Bayesian analysis, "AIC","BIC" or "AICc" will be used for model comparison instead of the Bayesian marginal likelihood. The item 'ML.IC' determines which of these three model comparison criteria are to be used. This variable is sent to 'compared.layered' for model comparison.

An 'layered' object, as returned from 'layer.analyzer'.

Not used. Only included for compatibility with other other 'summary' methods.

A set of 'layer.series.structure' objects that represents the set of time series that are to be analyzed and how the underlying processes of each such series is to be modelles. (See 'layer.series.structure'). At least one such object must be given. 'traverse.standalone.layered' can be used for settling the structure of each time series.

num.MCMC is an integer specifying the number of Markov chain Monte Carlo (MCMC) samples that the analysis will rest on. Default 1000. More MCMC samples are run internally, depending on the spacing (default 10) and burnin (default 1000).

Specifies the number of MCMC iteration between each sample that is used (default 10). A higher number means the samples will be less correlated but requires proptionally more computer time.

Markov chain Monte Carlo sampling takes some iterations to converge to the posterior parameter distribution. The burn-in phase is the set of iterations performed where the samples are then discarded, so as to avoid including samples before convergence. Convergence can be check by for instance 'gelman.diag' in the 'coda' package, and the MCMC samples themselves can be gotten by using the option 'mcmc'. Still, stability of the results can be quite a good enough measure for most purposes. Default is 10000 iterations. PS: The burn-in algorithm here includes two steps where the random walk MCMC variance is adjusted in order to optimize the efficiency of the algorithm.

The number of tempering chains (default 1, meaning it is just the MCMC chain itself that is beeing run). Parallel tempering (Geyer 1991) is a method for overcoming local optima and can thus increase the stability of the results. The way it is done is that in addition to the MCMC chain for the posterior distribution, parallel chains sampling from the distribution proportional to exp(-log(likelihood*prior)/T) where T is a "temperature" larger than one (this is in practise a "smoothed" version of the posterior), is also sampled from. Swaps between neighbouring chains are then allowed, making it possible for the MCMC chain for the posterior distribution to find a new optima.

As the code is not parallel, the computer time will increase proportional to the number of tempering chains. The set of "temperatures" is set to (1, T.ground, T.ground^2, ..., T.ground^(numtemp-1)), where 'T.ground' is default set to 1.5. If no swapping is done (this can be seen if 'silent.mode' is set to 'FALSE'), this "ground temperature" must be lowered, in order to facilitate the swaps.

Logical variable. If set to 'FALSE' (default), Bayesian analysis is performed. If set to 'TRUE', a classic maximum likelihood (ML) is performed, where the starting points of a set of optimizations is drawns from the MCMC samples. (This means that Bayesian MCMC sampling is performed anyway, but that the results from this is not focused on in the summary). The number of optimizations is determined by option 'maximum.likelihood.numstart', which is default set to 10. Estimation of the Bayesian model (marginal) likelihood is switched off when this option is used.

The number of ML optimizations performed, it applicable. (If 'do.maximum.likelihood' is set to 'TRUE'.) Default:10.

If set to 'FALSE', the routine will show a lot debug information while running. This can be useful for tempering purposes, since the number of tempering swaps are shown in this debug information.

If set to 'TRUE', this triggers the printing of more specific debug information concerning the burn-in phase.

If set to 'TRUE', this triggers the printing of more specific debug information concerning the likelihood calculation phase. (PS: This means a lot of printing.)

If set to 'TRUE'(default), the connection pairs and the causal and correlative connections of each model.

If set to 'TRUE', the analyses themselves are not performed. Only the connection pairs and the causal and correlative connections of each model are shown. This can be useful if one wants to know how many models are are to be analyzed.

Determines how to handle the identifiability problem in Bayesian multi-layer analysis. As mentioned in the supplementary of Reitan et al. (2012), one can switch the characteristic times of neighbouring layers and reorganize the stochastic contributions in such a way that the same top layer process takes place. This can be solved by requiring that the lower the layer, the larger the characteristic time. However, simply requiring this means the prior distribution of the characteristic times get shortened for multi-layered models, as compared to one-layered model. As long as the data is within the range of these shortened prior distributions, this puts the one-layer model at a disadvantage when it comes to model comparison. There are various solutions to this problem, and which solution is used is determined by the variable 'id.strategy'.

The variable 'id.strategy' should be an integer between 0 and 4, default 2. 0 - No identification treatment. (Default) PS: Can and even should yield multimodal characteristic times 1 - Keep upper characteristic time prior. Add lognormally to beneath-lying characteristic times. 2 - Keep lower characteristic time prior. Substract lognormally to above-lying characteristic times. 3 - Keep lower characteristic time prior. Cut depending on that on the above-lying characteristic times. 4 - Keep upper characteristic time prior. Cut depending on that on the below-lying characteristic times.

If set to 'TRUE', instead of reporting the stochastic contribution size in the stochastic differential equations, the stationary standard deviation is reported. The stationary standard deviation is the standard deviation of the process state at any given time after convergence, unconditional on previous states. In a multi-layered model, this is to be interpreted as the stationary standard deviation if this layer when let alone (with no other process affecting it). This is calculated as s*sqrt(characteristic time/2), where s is the stochastic contribution size. PS: This option does not make any sense if the bottom layer is a Wiener process, as that process is not stationary.

This determines how the parallel tempering chains are to be defined. Parallel tempering (Geyer 1991) is a method for overcoming local optima and can thus increase the stability of the results. The way it is done is that in addition to the MCMC chain for the posterior distribution, parallel chains sampling from the distribution proportional to exp(-log(likelihood*prior)/T) where T is a "temperature" larger than one (this is in practise a "smoothed" version of the posterior), is also sampled from. Swaps between neighbouring chains are then allowed, making it possible for the MCMC chain for the posterior distribution to find a new optima.

As the code is not parallel, the computer time will increase proportional to the number of tempering chains. The set of "temperatures" is set to (1, T.ground, T.ground^2, ..., T.ground^(numtemp-1)), where 'T.ground' is default set to 1.5. If no swapping is done (this can be seen if 'silent.mode' is set to 'FALSE'), this "ground temperature" must be lowered, in order to facilitate the swaps.

If set to 'TRUE', reports half-lives rather than characteristic times. Half-lives are the time the auto-correlation of an OU process drops to 1/2, or the time it takes for a perturbation from the expected value to drop to 1/2 the original size in expected value. (As opposed to dropping to exp(-1) for characteristic times.) Half-life=log(2)*characteristic time. The auto-correlation for an OU process (and thus for the inner dynamics of any linear layer) is exp(-diff.time/characteristic.time)=exp(-log(2)*diff.time/half.life)= (1/2)^(diff.time/half.life), where diff.time is the time difference between two process states.

If set to 'TRUE', returns the MCMC samples as an 'mcmc' object within the return object. The name of the list object will also be 'mcmc'.

If set to 'FALSE', no causal connection models will be analyzed, only correlative ones.

If set to 'FALSE', no correlative conneciton models will be analyzed, only causal ones. (Incompatible with 'allow.causal=FALSE'.

If set to 'FALSE', no direct causal loops will be examined, i.e. causal links from series 1, layer 1 to series 2, layer 2 and another liunk back again. (Does not rule out more elaborate causal loops.)

External series are time series that are treated as explanatory variables for the time series that are analyzed. They are not themselves analyzed statistically. Instead they are taken as is, with the assumption that what happenes between time points can be exppromately attained through linear interpolation. Since the stochastic nature of these series are not analyzed, only causal links are possible. Note that if external time series are given, one must also specify which series and layers they go to (see the next two entries.) Note that external data series are assumed global, and this connect equally to all sites (if that is a thing in your datasets).

An vector of integers that should specify which external time series are causally connected to which observed series, numbered from 1 to the number of series. The length of this vector should match the number of external time series.

An vector of integers that should specify which external time series are causally connected to which layers of the observed series, numbered from 1 to the number of layers in the series specified in 'external.series.conneciton'. The length of this vector should match the number of external time series.

A matrix that specifies the distance between each pair of sites. Only relevant when on series uses the 'distance.correlated.sigma' option to specify between-site correlation based on distance and a distance correlation falloff parameter. But in such a case, the site distance matrix is very necessary! (Keep in mind, sites are numbered from 0 to number of sites-1, so element i,j in the matrix should tell the distance between site i-1 and site j-1. Also keep in mind that the size of the matrix must match the number of sites in the dataset. Also, diagonal elements are the distance from a specific site to itself and so shoulød be 0. Also, the matrix must be square and symmetric, since the distance from site i-1 to site j-1 is the distance from site j-1 to site i-1.)

A single time series in the form of a 'layer.data.series' object. This is the input time series for the analysis.

The maximal number of layers (hidden and observed) that is allowed in the model traversal. If any of the switches for how each layer is handled (such as 'regional.options=TRUE', 'time.integrals.possible=TRUE', 'allow.one.feedback.loop=TRUE' and 'allow.deterministic.layers=TRUE') are in use, the number of models to traverse will quickly increase with this maximal number of layers.

Gives a quick description of each model traversed and shows the progress.

If set to 'TRUE', allows for one feedback loop from a higher to a lower layer (so 'against' the normal flow of causality in these models). This can provide an explanation for quasi-cyclical behaviour.

If set to 'FALSE' allows for other models than OU (Ornstein-Uhlenbeck) at the lowest layer, namely also the Wiener process (also called Random Walk) and linear trend plus OU. If set to 'TRUE*, only OU is allowed at the bottom layer.

Only in use when 'just.stationary=FALSE'. Switches off looking at the Wiener process (Random Walk) in the lowest layer if set to 'TRUE'.

If set to 'TRUE' switches on traversal of the different regional options. This only makes sense if the number of sites is larger than one. Also, it can massively expland the number of models traversed. Not implemented yet!

If set to 'TRUE' allows for time integral layers, i.e. the existence of one or more layers which are the time integral of the layer below. This of course only makes sense if the number of layers is alrger than one.

If set to 'TRUE', allows for traversing models where there is no stochasitcity (i.e. determinism) on some layers. This can easily be the case as long as the number of layers is larger than one. (The lowest layer cannot be deterministic, as that just means a flat line.)

num.MCMC is an integer specifying the number of Markov chain Monte Carlo (MCMC) samples that the analysis will rest on. Default 1000. More MCMC samples are run internally, depending on the spacing (default 10) and burnin (default 1000).

Specifies the number of MCMC iteration between each sample that is used (default 10). A higher number means the samples will be less correlated but requires proptionally more computer time.

Markov chain Monte Carlo sampling takes some iterations to converge to the posterior parameter distribution. The burn-in phase is the set of iterations performed where the samples are then discarded, so as to avoid including samples before convergence. Convergence can be check by for instance 'gelman.diag' in the 'coda' package, and the MCMC samples themselves can be gotten by using the option 'mcmc'. Still, stability of the results can be quite a good enough measure for most purposes. Default is 10000 iterations. PS: The burn-in algorithm here includes two steps where the random walk MCMC variance is adjusted in order to optimize the efficiency of the algorithm.

The number of tempering chains (default 1, meaning it is just the MCMC chain itself that is beeing run). Parallel tempering (Geyer 1991) is a method for overcoming local optima and can thus increase the stability of the results. The way it is done is that in addition to the MCMC chain for the posterior distribution, parallel chains sampling from the distribution proportional to exp(-log(likelihood*prior)/T) where T is a "temperature" larger than one (this is in practise a "smoothed" version of the posterior), is also sampled from. Swaps between neighbouring chains are then allowed, making it possible for the MCMC chain for the posterior distribution to find a new optima.

As the code is not parallel, the computer time will increase proportional to the number of tempering chains. The set of "temperatures" is set to (1, T.ground, T.ground^2, ..., T.ground^(numtemp-1)), where 'T.ground' is default set to 1.5. If no swapping is done (this can be seen if 'silent.mode' is set to 'FALSE'), this "ground temperature" must be lowered, in order to facilitate the swaps.

Logical variable. If set to 'FALSE' (default), Bayesian analysis is performed. If set to 'TRUE', a classic maximum likelihood (ML) is performed, where the starting points of a set of optimizations is drawns from the MCMC samples. (This means that Bayesian MCMC sampling is performed anyway, but that the results from this is not focused on in the summary). The number of optimizations is determined by option 'maximum.likelihood.numstart', which is default set to 10. Estimation of the Bayesian model (marginal) likelihood is switched off when this option is used.

The number of ML optimizations performed, it applicable. (If 'do.maximum.likelihood' is set to 'TRUE'.) Default:10.

Determines how to handle the identifiability problem in Bayesian multi-layer analysis. As mentioned in the supplementary of Reitan et al. (2012), one can switch the characteristic times of neighbouring layers and reorganize the stochastic contributions in such a way that the same top layer process takes place. This can be solved by requiring that the lower the layer, the larger the characteristic time. However, simply requiring this means the prior distribution of the characteristic times get shortened for multi-layered models, as compared to one-layered model. As long as the data is within the range of these shortened prior distributions, this puts the one-layer model at a disadvantage when it comes to model comparison. There are various solutions to this problem, and which solution is used is determined by the variable 'id.strategy'.

The variable 'id.strategy' should be an integer between 0 and 4, default 2. 0 - No identification treatment. (Default) PS: Can and even should yield multimodal characteristic times 1 - Keep upper characteristic time prior. Add lognormally to beneath-lying characteristic times. 2 - Keep lower characteristic time prior. Substract lognormally to above-lying characteristic times. 3 - Keep lower characteristic time prior. Cut depending on that on the above-lying characteristic times. 4 - Keep upper characteristic time prior. Cut depending on that on the below-lying characteristic times.

If set to 'TRUE', instead of reporting the stochastic contribution size in the stochastic differential equations, the stationary standard deviation is reported. The stationary standard deviation is the standard deviation of the process state at any given time after convergence, unconditional on previous states. In a multi-layered model, this is to be interpreted as the stationary standard deviation if this layer when let alone (with no other process affecting it). This is calculated as s*sqrt(characteristic time/2), where s is the stochastic contribution size. PS: This option does not make any sense if the bottom layer is a Wiener process, as that process is not stationary.

This determines how the parallel tempering chains are to be defined. Parallel tempering (Geyer 1991) is a method for overcoming local optima and can thus increase the stability of the results. The way it is done is that in addition to the MCMC chain for the posterior distribution, parallel chains sampling from the distribution proportional to exp(-log(likelihood*prior)/T) where T is a "temperature" larger than one (this is in practise a "smoothed" version of the posterior), is also sampled from. Swaps between neighbouring chains are then allowed, making it possible for the MCMC chain for the posterior distribution to find a new optima.

As the code is not parallel, the computer time will increase proportional to the number of tempering chains. The set of "temperatures" is set to (1, T.ground, T.ground^2, ..., T.ground^(numtemp-1)), where 'T.ground' is default set to 1.5. If no swapping is done (this can be seen if 'silent.mode' is set to 'FALSE'), this "ground temperature" must be lowered, in order to facilitate the swaps.

If set to 'TRUE', reports half-lives rather than characteristic times. Half-lives are the time the auto-correlation of an OU process drops to 1/2, or the time it takes for a perturbation from the expected value to drop to 1/2 the original size in expected value. (As opposed to dropping to exp(-1) for characteristic times.) Half-life=log(2)*characteristic time. The auto-correlation for an OU process (and thus for the inner dynamics of any linear layer) is exp(-diff.time/characteristic.time)=exp(-log(2)*diff.time/half.life)= (1/2)^(diff.time/half.life), where diff.time is the time difference between two process states.

If set to 'TRUE', returns the MCMC samples as an 'mcmc' object within the return object. The name of the list object will also be 'mcmc'.

If set to 'TRUE', the model includes a parameter for the initial values of each process (one for each combination of site and layer). The initial time is that of the first measurement (if not 'init.time' is specified).

If set, the model includes a parameter for the initial values of each process (one for each combination of site and layer). The initial time is that of 'init.time'.

Only makes sense for time series with multiple sites. If set, an initial value treatment will assume the same initial value for all sites (if applicable).

If set, an initial value treatment will assume the same initial value for all layers (if applicable, i.e. if 'numlayers'>1).

If given, should be an array of two elements, initial time and initial value.

If set to 'TRUE', allows for positive pull parameters in the linear stochastic equations. Positive pulls means the process is prone to "explode", but for a short time (typically not much more than the characteristic time), such a process could none the less make sense.

If given, should specify an array of period lengths, in order to include in the model periodic sine/cosine fluctuations (such as seasonality). Regression parameters for the sine/cosine terms are estimated.

An object specifying the prior distribution, in the form of a 'layer.prior' object.

A matrix that specifies the distance between each pair of sites. Only relevant when on series uses the 'distance.correlated.sigma' option to specify between-site correlation based on distance and a distance correlation falloff parameter. But in such a case, the site distance matrix is very necessary! (Keep in mind, sites are numbered from 0 to number of sites-1, so element i,j in the matrix should tell the distance between site i-1 and site j-1. Also keep in mind that the size of the matrix must match the number of sites in the dataset. Also, diagonal elements are the distance from a specific site to itself and so shoulød be 0. Also, the matrix must be square and symmetric, since the distance from site i-1 to site j-1 is the distance from site j-1 to site i-1.)