analyze.2order {doremi}

DOREMI second order analysis function

Description

analyze.2order estimates the coefficients of a second order differential equation of the form:

\frac{d^2y}{dt}(t) + 2\xi\omega_{n}\frac{dy}{dt}(t) + \omega_{n}^2 (y - y_{eq}) = \omega_{n}^2 k*u(t)

Where y(t) is the individual's signal, \frac{dy}{dt}(t) is its first derivative,\frac{d^2y}{dt}(t) its second derivative and u(t) is the excitation. The function estimates the coefficients 2\xi\omega_{n}, \omega_{n}^2 k and \omega_{n}^2 y_{eq}, from which the oscillation period T, the damping ratio \xi, the equilibrium y_{eq} value and the gain k can be derived. Th estimation is based on a two step procedure: the first step consists in estimating the derivatives to then estimate in a second step the differential equation coefficients through a linear mixed-effect model. Three different method to estimate the derivative during the first step are proposed.

Usage

analyze.2order(
  data,
  id = NULL,
  input = NULL,
  time = NULL,
  signal,
  dermethod = "gold",
  derparam = 3,
  order = 2,
  verbose = FALSE
)

Arguments

data

Is a data frame containing at least one column, that is the signal to be analyzed.

id

Is a CHARACTER containing the NAME of the column of data containing the identifier of the individual. If this parameter is not entered when calling the function, a single individual is assumed and a linear regression is done instead of the linear mixed-effects regression.

input

Is a CHARACTER or a VECTOR OF CHARACTERS containing the NAME(s) of data column(s) containing the EXCITATION vector(s). If this parameter is not entered when calling the function, the excitation is assumed to be unknown. In this case, the linear mixed-effect regression is still carried out but no coefficient is calculated for the excitation term. If no excitation term is supplied, one of the initial conditions is different from 0 (signal or derivative) and xi<1 the function will estimate a damped linear oscillator (DLO)

time

Is a CHARACTER containing the NAME of the column of data containing the time vector. If this parameter is not entered when calling the function, it is assumed that time steps are of 1 unit and the time vector is generated internally in the function.

signal

Is a CHARACTER containing the NAME of the column of the data frame containing the SIGNAL to be studied.

dermethod

is the derivative estimation method. The following methods are available: "gold","glla" and "fda"

derparam

If dermethod "glla" or "gold" are chosen, it is the embedding number, a positive integer containing the number of points to be used for the calculation of the derivatives. Its value by default is 3 as at least three points are needed for the calculation of the second derivative. If dermethod "fda" is chosen, this parameter is spar, the parameter related to the smoothing parameter lambda used in the penalization function of the function smooth.spline to estimate the derivative via splines (Functional Data Analysis)

order

is the maximum order of the derivative to estimate. Using an order higher than that of the maximum derivative to estimate (1 in first order differential equations and 2 in second order differential equations), for instance, order=4 might enhance derivative estimation (see doi: 10.1080/00273171.2015.1123138Chow et al.(2016))

verbose

Is a boolean that displays status messages of the function when set to 1.

Details

The analysis performs the following linear mixed-effects regression:

\ddot{y_{ij}} \sim (b_{0} +b_{0j}) + (b_{1}+b_{1j}) \dot{y_{ij}} + (b_{2}+b_{2j}) y_{ij} + (b_{3}+b_{3j}) u_{ij} + e_{ij}

with i accounting for the time and j for the different individuals. e_{ij} are the residuals, \dot{y_{ij}} is the first derivative and \ddot{y_{ij}} the second derivative calculated on embedding points, and y_{ij} and u_{ij} are the signal and the excitation averaged on embedding points. The fixed effect coefficients estimated from the regression are:

• xi2omega: b_1 (2\xi\omega^2 from the differential equation)

• omega2: b_2 (\omega^2 from the differential equation)

• komega2: b_3 (k\omega^2 from the differential equation)

• yeqomega2: b_0 (\omega^2 y_{eq} from the differential equation)

The coefficients derived to characterize the signal are calculated as follows:

• the oscillation period T = \sqrt{ \frac{2\pi}{ b_2 }}

• the damping factor xi: xi = \frac{b_1}{2*\sqrt{b_2}}

• the gain, k: k = \frac{b_3}{b_1}. It is the proportionality between the stationary increase of the signal and the excitation increase that caused it.

• the equilibrium value, yeq: yeq = \frac{b_0}{b_2}. It is the stationary value reached in the absence of excitation.

The estimation is performed using the function lmer if there are several individuals or lm if there is only one. With the above estimated parameters, the estimated signal is reconstructed for each individual by using the generate.2order function of this package (based on deSolve's ode function). The function returns five objects:

1. data- A data.frame including the input data, the intermediate calculations used to prepare the variables for the fit and the estimated trajectories for each individual.

   1. signal_rollmean - calculation of the o order derivative of the signal over embedding points.

   2. signal_derivate1 - calculation of the first derivative of the signal with the chosen method in embedding points/with smoothing parameter spar

   3. signal_derivate2 - calculation of the second derivative of the signal with the chosen method in embedding points/with smoothing parameter spar

   4. time_derivate - calculation of the moving average of the time vector over embedding points.

   5. input_rollmean - calculation of the moving average of the excitation vector over embedding points.

   6. estimated - the estimated signal calculated using deSolve's ode function with a second order model, the excitation provided as input and the coefficients obtained from the fit.

2. resultid- A data.frame including for each individual, listed by id number, the coefficients calculated (thus fixed + random component)

3. resultmean- A data.frame including the fixed effects of the coefficients mentioned above with their standard error, the coefficients characterizing the signal shape (i.e. the period, the damping factor, the gain and the equilibrium value), and the R2 resulting from the estimation

4. regression- A list containing the summary of the linear mixed-effects regression.

   As seen in the Description section, the print method by default prints only the resultmean element. Each one of the other objects can be accessed by indicating $ and their name after the result, for instance, for a DOREMI object called "result", it is possible to access the regression summary by typing result$regression.

5. derparam - contains the embedding number used to generate the results (if the derivative estimation method chosen is "glla" or "gold") or the smoothing parameter spar if the chosen method is fda

Value

Returns a summary of the fixed effect coefficients (see details)

See Also

calculate.gold, calculate.glla, calculate.fda to compute the derivatives, for details on embedding/spar See generate.2order for the generation of the solution of the second order differential equation.

Examples

# generate a panel of oscillating signals
test   <- generate.panel.2order(time = 0:100,
                              excitation = as.numeric(0:100>50),
                              period = 15,
                              xi = 0.3,
                              k = 2,
                              internoise = 0.2,
                              intranoise = 0.3,
                              nind = 10)

# plot the signal to analyze
plot(test)


# analyze them
res <- analyze.2order(data = test,
                      id = "id",
                      input = "excitation",
                      time =  "time",
                      signal = "signal",
                      derparam = 13)
res
plot(res)