Introduction-to-doremi"

 library(ggplot2)
 library(data.table)
 library(doremi)
set.seed(1)
knitr::opts_chunk$set(
  collapse = TRUE,
  warning = FALSE,
  fig.align="center",
  fig.width = 5,
  fig.height = 5,
  comment = "#>"
)

The main purpose of the Dynamics Of Return to Equilibrium during Multiple Inputs (doremi) package is to provide methods to estimate the dynamics parameters of self-regulated homeostatic systems experiencing multiple excitations. To do so, doremi provides functions to generate solutions of a first or second order differential equation, functions to estimate the derivatives of a given variable, and functions to estimate the coefficients of a first or second order differential equations with constant coefficients using a two-step estimation method.

In this introduction vignette, you will find in the examples section two examples of analysis that can be performed with doremi, a short presentation of what a first or second order differential equation is in the first or second order differential equation sections, and a list of the functions included in the package in the last section.

The detail of the analysis when supposing that the studied variable follows a first or second order differential equation is detailed in the associated vignettes for the first or second. The functions used to calculate the derivatives are detailed in the derivatives vignette.

EXAMPLES {#example}

Two examples of analysis are shown below. Detailed use of the functions can be found in the vignettes for the analysis based on a first or second order differential equation. Let us consider a variable following a first order differential equation and set out of equilibrium by a given excitation function. It can be generated as follows:

# create a time vector
time <- 0:90
# create an excitation mechanism
exc <- rep(c(0,1,0),c(11,30,50))

# generate simulated data
set.seed(123)
variable <- generate.panel.1order(time = time,
                           excitation = exc,
                           tau = 5,
                           nind = 1,
                           intranoise = 1)

plot(variable)

The green dots are the noisy variable following a first order differential equation model that we want to analyze. The blue curve is the same without noise (i.e. the solution of the differential equation for the given parameters). The analysis can be simply performed by calling a single function:

est_result <- analyze.1order(data = variable,
                             input = "excitation",
                             time = "time",
                             signal = "signal",
                             dermethod = "gold",
                             derparam = 13)

The results can be plotted easily:

plot(est_result)+
  ggplot2::geom_line(data = variable,aes(time,signalraw,color = "underlying model"))+
  ggplot2::labs(y = "",
       title = "estimation of first order\n differential equation parameters")

The blue curve is the curve reconstructed from the coefficients estimated by the analysis, to be compared to the magenta curve, representing the model underlying the noisy data analyzed. The est_result object contains all estimated coefficients.

The same can be performed for a noisy signal following a second order differential equation driven by the same excitation. Let us generate noisy data driven by the excitation and governed by a second order differential equation:

set.seed(123)
variable2 <- generate.panel.2order(time = time,
                                  excitation = exc,
                                  period = 30,
                                  y0 = 0,
                                  xi = 0.1,
                                  nind = 1,
                                  intranoise = 0.8)

plot(variable2)

We can perform the analysis based on the second order differential equation:

est_result2 <- analyze.2order(variable2,
                              input = "excitation",
                             time = "time",
                             signal = "signal",
                             dermethod = "glla",
                             derparam = 14)

plot(est_result2)+
  ggplot2::geom_line(data = variable2,aes(time,signalraw,color = "underlying model"))+
  ggplot2::labs(color = "",
       y = "",
       title = "estimation of second order \ndifferential equation parameters")

The first step of the analysis consists in estimating the time derivatives of the variable. Several methods are proposed. They all depend on the number of point they consider (the embedding number) or on a smoothing parameter. More detail can be found in the section detailing the package functions and in the vignette associated. Once the derivatives are estimated, the constant coefficients of the first or second order differential equation can be estimated by a simple linear mixed-effects (multilevel) regression. As the embedding number or smoothing parameter affects the estimation of the derivatives, it will also impact the quality of the estimation of the differential equation coefficients. The package doremi provides a function optimum_param which estimates the optimum embedding number (gold/glla) or smoothing parameter (fda) by varying the latter in a range provided as input and keeping the parameter that produces the optimal estimation.

FIRST ORDER DIFFERENTIAL EQUATION MODEL {#first}

We will shortly describe the first order differential equation and its coefficients here. A more detailed explanation of the possible analysis for a variable following a first order differential equation can be found in this vignette.

A first order differential equation with constant coefficient for the variable $y$ function of the variable $t$ (here the time) is:

\begin{equation} \frac{1}{\gamma} \dot{y}(t) + y(t) = kU(t) + y_{eq} \label{eq1} \end{equation}

Where:

The coefficients are:

time <- seq(0,100,0.1)
exc <- as.numeric(time >= 10 & time <= 40)
variable <- generate.1order(time = time,
                         excitation = exc,
                         y0 = 0,
                         tau = 10,
                         k = 1)

ggplot2::ggplot(data = as.data.table(variable)) +
ggplot2::ggtitle( "First order differential equation solution")+
  ggplot2::geom_line(ggplot2::aes(t,y, colour = "variable"))+
  ggplot2::geom_line(ggplot2::aes(t,exc, colour = "Excitation"),size = 1.5)+
  ggplot2::geom_hline(yintercept=0.63*max(variable$y), linetype="dashed", colour = "gray")+
  ggplot2::geom_hline(yintercept=0.37*max(variable$y), linetype="dashed", colour = "gray")+
  ggplot2::geom_vline(xintercept=variable$t[variable$y==max(variable$y)], colour = "gray")+
  ggplot2::geom_vline(xintercept=50, colour = "gray")+
  ggplot2::geom_vline(xintercept=19, colour = "gray")+
  ggplot2::geom_vline(xintercept=10, colour = "gray")+
  ggplot2::annotate("segment", x = 10, xend = 19, y = -0.1, yend = -0.1, colour = "dark green", size = 1)+
  ggplot2::annotate("text", x = 15, y = -0.2, label = "tau", parse = TRUE, colour = "dark green")+

  ggplot2::annotate("segment", x = 40, xend = 50, y = -0.1, yend = -0.1, colour = "dark green", size = 1)+
  ggplot2::annotate("text", x = 45, y = -0.2, label = "tau", parse = TRUE, colour = "dark green")+

  ggplot2::annotate("text", x = 75, y = 0.7, label = "63% diff. max and eq. value", colour = "gray")+
  ggplot2::annotate("text", x = 75, y = 0.3*max(variable$y), label = "37% diff. max and eq. value", colour = "gray")+
  ggplot2::labs(x = "Time (arb. unit)",
           y = "variable (arb. unit)",
           colour = "Legend")+
  ggplot2::theme_bw() + 
  ggplot2::theme(legend.position = "top", plot.title = ggplot2::element_text(hjust = 0.5))

Once the derivative estimated, doremi performs a linear mixed-effect regression to estimate these parameters. The estimated variable can be reconstructed for each individual using the numerical estimation provided by the function ode from the package deSolve.

SECOND ORDER DIFFERENTIAL EQUATION MODEL {#second}

We will in this section shortly describe the second order differential equation and its coefficients. A more detailed explanation of the possible analysis for a variable following a second order differential equation can be found in this vignette. The differential equation considered in this case is the following:

\begin{equation} \frac{d^2y}{dt} + 2\xi\omega_{n}\frac{dy}{dt} + \omega_{n}^2 y = \omega_{n}^2(y_{eq} + k U(t)) \label{eq3} \end{equation}

Where:

And regarding the coefficients:

Once the derivative estimated, doremi performs a linear mixed-effect regression to estimate these parameters. The estimated variable can be reconstructed for each individual using the numerical estimation provided by the function ode from the package deSolve.

We can use generate.2order to generate a solution of this differential equation, for a given $U(t)$, a period of 15, and a damping factor of 0.2:

time <- 0:130
excitation <- c(rep(0,30),rep(1,50),rep(0,51))
variable <- generate.2order(time = time,
                            excitation = excitation,
                            xi = 0.2,
                            period = 15,
                            k = 1)

ggplot2::ggplot(variable)+
  ggplot2::ggtitle( "Second order differential equation solution")+
  ggplot2::geom_line(aes(t,y,color = "variable"))+
  ggplot2::geom_line(aes(time,excitation,color = "excitation"))+
  ggplot2::geom_vline(xintercept=80, colour = "gray")+
  ggplot2::annotate("text", x = 100, y = 1, label = "DLO", parse = TRUE, colour = "gray")+
  ggplot2::labs(x = "Time (arb. unit)",
           y = "variable (arb. unit)",
           colour = "Legend")+
  ggplot2::theme_bw() + 
  ggplot2::theme(legend.position = "top", plot.title = ggplot2::element_text(hjust = 0.5))

FUNCTIONS INCLUDED ON THE PACKAGE {#functions}

The doremi package contains three types of functions:

Simulation functions

These are functions that allow the user to generate data following the first/second order differential equation models (i.e. solution of equations (1)/(3)). More specifically:

Analysis functions

These functions allow to analyze a set of data and verify how close it is to being a solution of a first/second order differential equation with constant coefficients by following a two-step estimation method (derivative estimation first and then estimation of the coefficients through a multilevel regression).

*GOLD- Generalized Orthogonal Local Derivative, method described in \href{https://doi.org/10.1080/00273171.2010.498294}{Deboeck (2010)}. The code available on this paper was extracted and adapted for non-constant time steps. This method allows calculating over a number of measurement points (called the embedding number) the first and second derivatives (or higher, depending on the order set as input parameter) with errors uncorrelated with the variable.

*GLLA- Generalized Local Linear Approximation, method described \href{https://doi.org/10.4324/9780203864746}{Boker et al.(2010)}. This method allows to estimate the derivatives over a number of measurement points called the embedding number assuming an equally spaced time series.

*FDA- Functional Data Analysis. This method creates a function that fits the data and then derives it to evaluate the derivatives in those same points.

The analyze.1order/analyze.2order functions use one of these three methods according to what is set as input and then, once the derivatives are estimated, they are used as a known term in the multilevel regression so that the coefficients of the equation are the only unknowns and can be estimated through the fit. Once the coefficients are estimated, the estimated 1st order/2nd order variable is generated from these and the R2 is calculated and provided as a result. The estimation of the derivatives, summary of the regression, coefficients found for each individual (random coefficients of the regression) and for the group (fixed coefficients) are also provided.

derivative functions



Try the doremi package in your browser

Any scripts or data that you put into this service are public.

doremi documentation built on Jan. 29, 2021, 5:06 p.m.