Package 'dynr'

Description

Intensive longitudinal data have become increasingly prevalent in various scientific disciplines. Many such data sets are noisy, multivariate, and multi-subject in nature. The change functions may also be continuous, or continuous but interspersed with periods of discontinuities (i.e., showing regime switches). The package 'dynr' (Dynamic Modeling in R) is an R package that implements a set of computationally efficient algorithms for handling a broad class of linear and nonlinear discrete- and continuous-time models with regime-switching properties under the constraint of linear Gaussian measurement functions. The discrete-time models can generally take on the form of a state-space or difference equation model. The continuous-time models are generally expressed as a set of ordinary or stochastic differential equations. All estimation and computations are performed in C, but users are provided with the option to specify the model of interest via a set of simple and easy-to-learn model specification functions in R. Model fitting can be performed using single-subject time series data or multiple-subject longitudinal data. Ou, Hunter, & Chow (2019) <doi:10.32614%2FRJ-2019-012> provided a detailed introduction to the interface and more information on the algorithms.

Details

The DESCRIPTION file:

Index of help topics:

EMG                     Single-subject time series of facial
                        electromyography data
EMGsim                  Simulated single-subject time series to capture
                        features of facial electromyography data
ExpandRandomAsLVModel   Extend a user-specified model to include random
                        varibles
LinearOsc               Simulated time series data for a deterministic
                        linear damped oscillator model
LogisticSetPointSDE     Simulated time series data for a stochastic
                        linear damped oscillator model with logistic
                        time-varying setpoints
NonlinearDFAsim         Simulated multi-subject time series based on a
                        dynamic factor analysis model with nonlinear
                        relations at the latent level
Oscillator              Simulated time series data of a damped linear
                        oscillator
Outliers                Simulated time series data for detecting
                        outliers.
PFAsim                  Simulated time series data of a multisubject
                        process factor analysis
PPsim                   Simulated time series data for multiple
                        eco-systems based on a predator-and-prey model
RSPPsim                 Simulated time series data for multiple
                        eco-systems based on a regime-switching
                        predator-and-prey model
TrueInit_Y14            Simulated multilevel multi-subject time series
                        of a Van der Pol Oscillator
VARsim                  Simulated time series data for multiple
                        imputation in dynamic modeling.
autoplot.dynrTaste      The ggplot of the outliers estimates.
coef.dynrModel          Extract fitted parameters from a dynrCook
                        Object
confint.dynrCook        Confidence Intervals for Model Parameters
diag,character-method   Create a diagonal matrix from a character
                        vector
dynr-package            Dynamic Models with Regime-Switching
dynr.config             Check that dynr in configured properly
dynr.cook               Cook a dynr model to estimate its free
                        parameters
dynr.data               Create a list of data for parameter estimation
                        (cooking dynr) using 'dynr.cook'
dynr.flowField          A Function to plot the flow or velocity field
                        for a one or two dimensional autonomous ODE
                        system from the phaseR package written by
                        Michael J. Grayling.
dynr.ggplot             The ggplot of the smoothed state estimates and
                        the most likely regimes
dynr.ldl                LDL Decomposition for Matrices
dynr.mi                 Multiple Imputation of dynrModel objects
dynr.model              Create a dynrModel object for parameter
                        estimation (cooking dynr) using 'dynr.cook'
dynr.plotFreq           Plot of the estimated frequencies of the
                        regimes across all individuals and time points
                        based on their smoothed regime probabilities
dynr.taste              Detect outliers in state space models.
dynr.taste2             Re-fit state-space model using the estimated
                        outliers.
dynr.trajectory         A Function to perform numerical integration of
                        the chosen ODE system, for a user-specified set
                        of initial conditions. Plots the resulting
                        solution(s) in the phase plane. This function
                        from the phaseR package written by Michael J.
                        Grayling.
dynr.version            Current Version String
dynrCook-class          The dynrCook Class
dynrDynamics-class      The dynrDynamics Class
dynrInitial-class       The dynrInitial Class
dynrMeasurement-class   The dynrMeasurement Class
dynrModel-class         The dynrModel Class
dynrNoise-class         The dynrNoise Class
dynrRecipe-class        The dynrRecipe Class
dynrRegimes-class       The dynrRegimes Class
dynrTrans-class         The dynrTrans Class
getdx                   A wrapper function to call functions in the fda
                        package to obtain smoothed estimated
                        derivatives at a specified order
internalModelPrep       Do internal model preparation for dynr
logLik.dynrCook         Extract the log likelihood from a dynrCook
                        Object
names,dynrCook-method   Extract the free parameter names of a dynrCook
                        object
names,dynrModel-method
                        Extract the free parameter names of a dynrModel
                        object
nobs.dynrCook           Extract the number of observations for a
                        dynrCook object
nobs.dynrModel          Extract the number of observations for a
                        dynrModel object
oscData                 Another simulated multilevel multi-subject time
                        series of a damped oscillator model
plot.dynrCook           Plot method for dynrCook objects
plotFormula             Plot the formula from a model
plotGCV                 A function to evaluate the generalized
                        cross-validation (GCV) values associated with
                        derivative estimates via Bsplines at a range of
                        specified smoothing parameter (lambda) values
predict.dynrModel       'predict' method for 'dynrModel' objects
prep.formulaDynamics    Recipe function for specifying dynamic
                        functions using formulas
prep.initial            Recipe function for preparing the initial
                        conditions for the model.
prep.loadings           Recipe function to quickly create factor
                        loadings
prep.matrixDynamics     Recipe function for creating Linear Dynamics
                        using matrices
prep.measurement        Prepare the measurement recipe
prep.noise              Recipe function for specifying the measurement
                        error and process noise covariance structures
prep.regimes            Recipe function for creating regime switching
                        (Markov transition) functions
prep.tfun               Create a dynrTrans object to handle the
                        transformations and inverse transformations of
                        model paramters
printex                 The printex Method
summary.dynrCook        Get the summary of a dynrCook object
theta_plot              A function to plot simple slopes and region of
                        significance.
vcov.dynrCook           Extract the Variance-Covariance Matrix of a
                        dynrCook object
vdpData                 Another simulated multilevel multi-subject time
                        series of a Van der Pol Oscillator

Further information is available in the following vignettes:

Because the dynr package compiles C code in response to user input, more setup is required for the dynr package than for many others. We acknowledge that this additional setup can be bothersome, but we believe the ease of use for the rest of the package and the wide variety of models it is possible to fit with it will compensate for this initial burden. Hopefully you will agree!

See the installation vignette referenced in the Examples section below for installation instructions.

The naming convention for dynr exploits the pronunciation of the package name, dynr, pronounced the same as “dinner”. That is, the names of functions and methods are specifically designed to relate to things done surrounding dinner, such as gathering ingredients (e.g., the data), preparing recipes, cooking, and serving the finished product. The general procedure for using the dynr package can be summarized in five steps as below.

1. Data are prepared using with the dynr.data() function.

2. Recipes are prepared. To each part of a model there is a corresponding prep.*() recipe function. Examples of such prep.*() functions include: prep.measurement(), prep.matrixDynamics(), prep.formulaDynamics(), prep.initial(), prep.noise(), and prep.regimes().

3. The function dynr.model() mixes the data and recipes together into a model object of class dynrModel.

4. The model is cooked with dynr.cook().

5. Results from model fitting and related estimation are served using functions such as summary(), plot(), dynr.ggplot() (or its alias autoplot()), plotFormula(), and printex().

Note

State-space modeling, dynamic model, differential equation, regime switching, nonlinear

Author(s)

Lu Ou [aut], Michael D. Hunter [aut, cre] (<https://orcid.org/0000-0002-3651-6709>), Sy-Miin Chow [aut] (<https://orcid.org/0000-0003-1938-027X>), Linying Ji [aut], Meng Chen [aut], Hui-Ju Hung [aut], Jungmin Lee [aut], Yanling Li [aut], Jonathan Park [aut], Massachusetts Institute of Technology [cph], S. G. Johnson [cph], Benoit Scherrer [cph], Dieter Kraft [cph]

Maintainer: Michael D. Hunter <[email protected]>

References

Chow S, Grimm KJ, Guillaume F, Dolan CV, McArdle JJ (2013). “Regime-switching bivariate dual change score model.” Multivariate Behavioral Research, 48(4), 463-502. doi:10.1080/00273171.2013.787870.

Chow S, Zhang G (2013). “Nonlinear Regime-Switching State-Space (RSSS) Models.” Psychometrika: Application Reviews and Case Studies, 78(4), 740-768. doi:10.1007/s11336-013-9330-8.

Ou L, Hunter MD, Chow S (2019). “What's for dynr: A package for linear and nonlinear dynamic modeling in R.” The R Journal, 11(1), 1-20.

Yang M, Chow S (2010). “Using state-space model with regime switching to represent the dynamics of Facial electromyography (EMG) data.” Psychometrika: Application and Case Studies, 74(4), 744-771. doi:10.1007/s11336-010-9176-2.

Chow S, Ou L, Ciptadi A, Prince E, You D, Hunter MD, Rehg JM, Rozga A, Messinger DS (2018). “Representing sudden shifts in intensive dyadic interaction data using differential equation models with regime switching.” Psychometrika, 83, 476-510. doi:10.1007/s11336-018-9605-1.

See Also

For other annotated tutorials using the dynr package see https://quantdev.ssri.psu.edu/resources/what%E2%80%99s-dynr-package-linear-and-nonlinear-dynamic-modeling-r

Examples

# For installation instructions see the package vignette below
## Not run: 
vignette(package='dynr', 'InstallationForUsers')

## End(Not run)
# This should open a pdf/html file to guide you through proper
#  installation and configuration.

#For illustrations of the functions in dynr, check out some of the demo examples in:
## Not run: 
demo(package='dynr')

## End(Not run)

#For example, to run the demo 'LinearSDE' type
# the following without the comment character (#) in front of it.
## Not run: 
demo('LinearSDE', package='dynr')

## End(Not run)

Description

The ggplot of the outliers estimates.

Usage

## S3 method for class 'dynrTaste'
autoplot(object, numSubjDemo = 2, idtoPlot = NULL,
  names.state = NULL, names.observed = NULL, ...)

Arguments

Value

a list of ggplot objects for each ID. The plots of chi-square statistics (joint and independent), and the plots of t statistic for names.state and names.observed will be included. Users can modify the ggplot objects using ggplot grammar. If a filename is provided, a pdf of plots will be saved additionally.

Description

aliases coef.dynrModel coef<- coef<-.dynrModel

Usage

## S3 method for class 'dynrModel'
coef(object, ...)

coef(object) <- value

## S3 replacement method for class 'dynrModel'
coef(object) <- value

## S3 method for class 'dynrCook'
coef(object, ...)

Arguments

Value

A numeric vector of the fitted parameters.

See Also

Other S3 methods logLik.dynrCook

Examples

# Create a minimal cooked model called 'cook'
require(dynr)

meas <- prep.measurement(
	values.load=matrix(c(1, 0), 1, 2),
	params.load=matrix(c('fixed', 'fixed'), 1, 2),
	state.names=c("Position","Velocity"),
	obs.names=c("y1"))

ecov <- prep.noise(
	values.latent=diag(c(0, 1), 2),
	params.latent=diag(c('fixed', 'dnoise'), 2),
	values.observed=diag(1.5, 1),
	params.observed=diag('mnoise', 1))

initial <- prep.initial(
	values.inistate=c(0, 1),
	params.inistate=c('inipos', 'fixed'),
	values.inicov=diag(1, 2),
	params.inicov=diag('fixed', 2))

dynamics <- prep.matrixDynamics(
	values.dyn=matrix(c(0, -0.1, 1, -0.2), 2, 2),
	params.dyn=matrix(c('fixed', 'spring', 'fixed', 'friction'), 2, 2),
	isContinuousTime=TRUE)

data(Oscillator)
data <- dynr.data(Oscillator, id="id", time="times", observed="y1")

model <- dynr.model(dynamics=dynamics, measurement=meas,
	noise=ecov, initial=initial, data=data)

## Not run: 
cook <- dynr.cook(model,
	verbose=FALSE, optimization_flag=FALSE, hessian_flag=FALSE)

# Now grab the coef!
coef(cook)

## End(Not run)

Description

Confidence Intervals for Model Parameters

Usage

## S3 method for class 'dynrCook'
confint(object, parm, level = 0.95,
  type = c("delta.method", "endpoint.transformation"),
  transformation = NULL, ...)

Arguments

Details

The parm argument can be a numeric vector or a vector of names. If it is missing then it defaults to using all the parameters.

These are Wald-type confidence intervals based on the standard errors of the (transformed) parameters. Wald-type confidence intervals are known to be inaccurate for variance parameters, particularly when the variance is near zero (See references for issues with Wald-type confidence intervals).

Value

A matrix with columns giving lower and upper confidence limits for each parameter. These will be labelled as (1-level)/2 and 1 - (1-level)/2 as a percentage (e.g. by default 2.5

References

Pritikin, J.N., Rappaport, L.M. & Neale, M.C. (In Press). Likelihood-Based Confidence Intervals for a Parameter With an Upper or Lower Bound. Structural Equation Modeling. DOI: 10.1080/10705511.2016.1275969

Neale, M. C. & Miller M. B. (1997). The use of likelihood based confidence intervals in genetic models. Behavior Genetics, 27(2), 113-120.

Pek, J. & Wu, H. (2015). Profile likelihood-based confidence intervals and regions for structural equation models. Psychometrica, 80(4), 1123-1145.

Wu, H. & Neale, M. C. (2012). Adjusted confidence intervals for a bounded parameter. Behavior genetics, 42(6), 886-898.

Examples

# Minimal model
require(dynr)

meas <- prep.measurement(
	values.load=matrix(c(1, 0), 1, 2),
	params.load=matrix(c('fixed', 'fixed'), 1, 2),
	state.names=c("Position","Velocity"),
	obs.names=c("y1"))

ecov <- prep.noise(
	values.latent=diag(c(0, 1), 2),
	params.latent=diag(c('fixed', 'dnoise'), 2),
	values.observed=diag(1.5, 1),
	params.observed=diag('mnoise', 1))

initial <- prep.initial(
	values.inistate=c(0, 1),
	params.inistate=c('inipos', 'fixed'),
	values.inicov=diag(1, 2),
	params.inicov=diag('fixed', 2))

dynamics <- prep.matrixDynamics(
	values.dyn=matrix(c(0, -0.1, 1, -0.2), 2, 2),
	params.dyn=matrix(c('fixed', 'spring', 'fixed', 'friction'), 2, 2),
	isContinuousTime=TRUE)

data(Oscillator)
data <- dynr.data(Oscillator, id="id", time="times", observed="y1")

model <- dynr.model(dynamics=dynamics, measurement=meas,
	noise=ecov, initial=initial, data=data)

## Not run: 
cook <- dynr.cook(model,
	verbose=FALSE, optimization_flag=FALSE, hessian_flag=FALSE)

# Now get the confidence intervals
# But note that they are nonsense because we set hessian_flag=FALSE !!!!
confint(cook)

## End(Not run)

Description

Create a diagonal matrix from a character vector

Usage

## S4 method for signature 'character'
diag(x = 1, nrow, ncol)

Arguments

Details

We create a new method for diag with character input. The default behavior for missing nrow and/or ncol arguments is the same as for the diag function in the base package. Off-diagonal entries are filled with "0".

Value

A matrix

Examples

diag(letters[1:3])

Description

Check that dynr in configured properly

Usage

dynr.config(verbose = FALSE)

Arguments

Details

The 'dynr' package requires additional set-up and configuration beyond just installing the package. In particular, it requires compiling C code along with GSL to run (cook) models. This function runs some basic checks of the configuration. We check that (1) R is on the PATH variable, (2) Rtools exists and is on the PATH variable for Windows, (3) a C compiler is available, and (4) GSL is available and on the PATH.

In general, see the 'Installation for Users' vignette for set-up and configuration instructions.

Value

No return value.

Examples

## Not run: dynr.config()

Description

Cook a dynr model to estimate its free parameters

Usage

dynr.cook(dynrModel, conf.level = 0.95, infile, optimization_flag = TRUE,
  hessian_flag = TRUE, verbose = TRUE, weight_flag = FALSE,
  debug_flag = FALSE, perturb_flag = FALSE)

Arguments

Details

Free parameter estimation uses the SLSQP routine from NLOPT.

The typical items returned in the cooked model are the filtered and smoothed latent variable estimates. eta_smooth_final, error_cov_smooth_final and pr_t_given_T are respectively time-varying smoothed latent variable mean estimates, smoothed error covariance estimates, and smoothed regime probability. eta_filtered, error_cov_filtered and pr_t_given_t are respectively time-varying filtered latent variable mean estimates, filtered error covariance matrix estimates, and filtered regime probability. Note that if theta.formula is provided in dynrModel@dynamics, this assumes that random effects are present in the dynamic equation. This would call an internal function to insert the random effect components as additional state variables. In this case, the last set of elements (rows) in eta_smooth_final would contain the estimated random effect components.

When debug_flag is TRUE, then additional information is passed into the cooked model. eta_predicted, error_cov_predicted, innov_vec, and residual_cov are respectively time-varying predicted latent variable mean estimates, predicted error covariance matrix estimates, the error/residual estimates (innovation vector), and the error/residual covariance matrix estimates.

The exit flag given after optimization has finished is from the SLSQP optimizer. Generally, error codes have negative values and successful codes have positive values. However, codes 5 and 6 do not indicate the model converged, but rather simply ran out of iterations or time, respectively. A more full description of each code is available at https://nlopt.readthedocs.io/en/latest/NLopt_Reference/#return-values and is also listed in the table below.

The last row of this table corresponding to an exit code of -6, is not from NLOPT, but rather is specific to the dynr package.

Value

Object of class dynrCook.

See Also

autoplot, coef, confint, deviance, initialize, logLik, names, nobs, plot, print, show, summary, vcov.

Examples

# Minimal model
require(dynr)

meas <- prep.measurement(
	values.load=matrix(c(1, 0), 1, 2),
	params.load=matrix(c('fixed', 'fixed'), 1, 2),
	state.names=c("Position","Velocity"),
	obs.names=c("y1"))

ecov <- prep.noise(
	values.latent=diag(c(0, 1), 2),
	params.latent=diag(c('fixed', 'dnoise'), 2),
	values.observed=diag(1.5, 1),
	params.observed=diag('mnoise', 1))

initial <- prep.initial(
	values.inistate=c(0, 1),
	params.inistate=c('inipos', 'fixed'),
	values.inicov=diag(1, 2),
	params.inicov=diag('fixed', 2))

dynamics <- prep.matrixDynamics(
	values.dyn=matrix(c(0, -0.1, 1, -0.2), 2, 2),
	params.dyn=matrix(c('fixed', 'spring', 'fixed', 'friction'), 2, 2),
	isContinuousTime=TRUE)

data(Oscillator)
data <- dynr.data(Oscillator, id="id", time="times", observed="y1")

model <- dynr.model(dynamics=dynamics, measurement=meas,
	noise=ecov, initial=initial, data=data)

## Not run: 
# Now cook the model!
cook <- dynr.cook(model,
	verbose=FALSE, optimization_flag=FALSE, hessian_flag=FALSE)

## End(Not run)

Create a list of data for parameter estimation (cooking dynr) using dynr.cook

Description

Create a list of data for parameter estimation (cooking dynr) using dynr.cook

Usage

dynr.data(dataframe, id = "id", time = "time", observed, covariates)

Arguments

Value

A list with components as needed for dynr.model

Examples

data(EMGsim)
dd <- dynr.data(EMGsim, id = 'id', time = 'time', observed = 'EMG', covariates = 'self')

z <- ts(matrix(rnorm(300), 100, 3), start = c(1961, 1), frequency = 12)
dz <- dynr.data(z)

Description

A Function to plot the flow or velocity field for a one or two dimensional autonomous ODE system from the phaseR package written by Michael J. Grayling.

Usage

dynr.flowField(deriv, xlim, ylim, parameters = NULL, system = "two.dim",
  points = 21, col = "gray", arrow.type = "equal", arrow.head = 0.05,
  frac = 1, add = TRUE, xlab = "x", ylab = "y", state.names = c("x",
  "y"), ...)

Arguments

Value

Returns a list with the following components: add, arrow.head, arrow.type, col, deriv, dx, dy, frac, parameters, points, system, x, xlab, xlim, y, ylab, ylim. Most of these components correspond simply to their original input values.

The only new elements are:

dx = A matrix. In the case of a two dimensional system, the values of the derivative of the first dependent derivative at all evaluated points.

dy = A matrix. In the case of a two dimensional system, the values of the derivative of the second dependent variable at all evaluated points. In the case of a one dimensional system, the values of the derivative of the dependent variable at all evaluated points.

x = A vector. In the case of a two dimensional system, the values of the first dependent variable at which the derivatives were computed. In the case of a one dimensional system, the values of the independent variable at which the derivatives were computed.

y = A vector. In the case of a two dimensional system, the values of the second dependent variable at which the derivatives were computed. In the case of a one dimensional system, the values of the dependent variable at which the derivatives were computed.

Note

The phaseR package was taken off cran as off 10/1/2019 so we are exporting some selected functions from phaseR_2.0 published on 8/20/2018. For details of these functions please see original documentations on the phaseR package.

References

Grayling, Michael J. (2014). phaseR: An R Package for Phase Plane Analysis of Autonomous ODE Systems. The R Journal, 6(2), 43-51. DOI: 10.32614/RJ-2014-023. Available at https://doi.org/10.32614/RJ-2014-023

Description

The ggplot of the smoothed state estimates and the most likely regimes

Usage

dynr.ggplot(res, dynrModel, style = 1, numSubjDemo = 2, idtoPlot = c(),
  names.state, names.observed, names.regime, shape.values, title, ylab,
  is.bw = FALSE, colorPalette = "Set2", fillPalette = "Set2",
  mancolorPalette, manfillPalette, ...)

## S3 method for class 'dynrCook'
autoplot(object, dynrModel, style = 1, numSubjDemo = 2,
  idtoPlot = c(), names.state, names.observed, names.regime, shape.values,
  title, ylab, is.bw = FALSE, colorPalette = "Set2", fillPalette = "Set2",
  mancolorPalette, manfillPalette, ...)

Arguments

Details

This function outputs a ggplot layer that can be modified using functions in the package ggplot2. That is, one can add layers, scales, coords and facets with the "+" sign. In an example below, the ggplot2::ylim() function is used to modify the limits of the y axis of the graph. More details can be found on https://ggplot2.tidyverse.org/ and https://ggplot2.tidyverse.org/reference/.

The two functions dynr.ggplot() and autoplot() as identical aliases of one another. The autoplot() function is an S3 method from the package ggplot2 that allows many objects to be plotted and works like the base plot() function.

Value

ggplot object

ggplot object

Examples

# The following code is part of a demo example in dynr
## Not run: 
demo(RSLinearDiscreteYang, package='dynr')
p <- dynr.ggplot(yum, dynrModel = rsmod, style = 1,
 	names.regime = c("Deactivated", "Activated"),
 	title = "(B) Results from RS-AR model", numSubjDemo = 1,
 	shape.values = c(1),
 	text = element_text(size = 16),
 	is.bw = TRUE)
# One can modify the limits on the y axis by using '+'
p + ggplot2::ylim(-2, 4)

autoplot(yum, dynrModel = rsmod, style = 1,
	names.regime = c("Deactivated", "Activated"),
	title = "(B) Results from RS-AR model", numSubjDemo = 1,
	shape.values = c(1),
	text = element_text(size = 16),
	is.bw = TRUE)

## End(Not run)

Description

LDL Decomposition for Matrices

Usage

dynr.ldl(x)

Arguments

Value

A matrix

Description

Multiple Imputation of dynrModel objects

Usage

dynr.mi(dynrModel, which.aux = NULL, which.lag = NULL, lag = 0,
  which.lead = NULL, lead = 0, m = 5, iter = 5, imp.obs = FALSE,
  imp.exo = TRUE, diag = TRUE, Rhat = 1.1, conf.level = 0.95,
  verbose = TRUE, seed = NA)

Arguments

Details

See the demo, demo(package='dynr', 'MILinearDiscrete'), for an illustrative example of using dynr.mi to implement multiple imputation with a vector autoregressive model.

Value

an object of ‘dynrMi’ class that is a list containing: 1. the imputation information, including a data set containing structured lagged and leading variables and a ‘mids’ object from mice() function; 2. the diagnostic information, including trace plots, an Rhat plot and a matrix containing Rhat values; 3. the estimation results, including parameter estimates, standard error estimates and confidence intervals.

References

Ji, L., Chow, S-M., Schermerhorn, A.C., Jacobson, N.C., & Cummings, E.M. (2018). Handling Missing Data in the Modeling of Intensive Longitudinal Data. Structural Equation Modeling: A Multidisciplinary Journal, 1-22.

Yanling Li, Linying Ji, Zita Oravecz, Timothy R. Brick, Michael D. Hunter, and Sy-Miin Chow. (2019). dynr.mi: An R Program for Multiple Imputation in Dynamic Modeling. International Journal of Computer, Electrical, Automation, Control and Information Engineering, 13, 302-311.

Create a dynrModel object for parameter estimation (cooking dynr) using dynr.cook

Description

Create a dynrModel object for parameter estimation (cooking dynr) using dynr.cook

Usage

dynr.model(dynamics, measurement, noise, initial, data, ...,
  outfile = tempfile())

Arguments

Details

A dynrModel is a collection of recipes. The recipes are constructed with the functions prep.measurement, prep.noise, prep.formulaDynamics, prep.matrixDynamics, prep.initial, and in the case of regime-switching models prep.regimes. Additionally, data must be prepared with dynr.data and added to the model.

Several named arguments can be passed into the ... section of the function. These include

• Argument regimes is for a dynrRegimes object prepared with prep.regimes

• Argument transform is for a dynrTrans object prepared with prep.tfun.

• Argument options a list of options. Check the NLopt website https://nlopt.readthedocs.io/en/latest/NLopt_Reference/#stopping-criteria for details. Available options for use with a dynrModel object include xtol_rel, stopval, ftol_rel, ftol_abs, maxeval, and maxtime, all of which control the termination conditions for parameter optimization. The examples below show a case where options were set.

There are several available methods for dynrModel objects.

• The dollar sign ($) can be used to both get objects out of a model and to set pieces of the model.

• names returns the names of the free parameters in a model.

• printex prints LaTeX expressions for the equations that compose a model. The output can then be readily typeset for inclusion in presentations and papers.

• nobs gives the total number of observations (e.g. all times across all people)

• coef gives the free parameter starting values. Free parameters can also be assigned with coef(model) <- aNamedVectorOfCoefficients

Value

Object of class 'dynrModel'

Examples

# Create a minimal model called 'model'
# without 'cooking' (i.e., estimating parameters)
require(dynr)

meas <- prep.measurement(
	values.load=matrix(c(1, 0), 1, 2),
	params.load=matrix(c('fixed', 'fixed'), 1, 2),
	state.names=c("Position","Velocity"),
	obs.names=c("y1"))

ecov <- prep.noise(
	values.latent=diag(c(0, 1), 2),
	params.latent=diag(c('fixed', 'dnoise'), 2),
	values.observed=diag(1.5, 1),
	params.observed=diag('mnoise', 1))

initial <- prep.initial(
	values.inistate=c(0, 1),
	params.inistate=c('inipos', 'fixed'),
	values.inicov=diag(1, 2),
	params.inicov=diag('fixed', 2))

dynamics <- prep.matrixDynamics(
	values.dyn=matrix(c(0, -0.1, 1, -0.2), 2, 2),
	params.dyn=matrix(c('fixed', 'spring', 'fixed', 'friction'), 2, 2),
	isContinuousTime=TRUE)

data(Oscillator)
data <- dynr.data(Oscillator, id="id", time="times", observed="y1")

# Now here's the model!
model <- dynr.model(dynamics=dynamics, measurement=meas,
	noise=ecov, initial=initial, data=data)

Description

Plot of the estimated frequencies of the regimes across all individuals and time points based on their smoothed regime probabilities

Usage

dynr.plotFreq(res, dynrModel, names.regime, title, xlab, ylab, textsize = 12,
  print = TRUE)

Arguments

Value

ggplot object

Description

Compute shocks and chi-squared diagnostics following Chow, Hamaker, and Allaire (2009). Using Innovative Outliers to Detect Discrete Shifts in Dynamics in Group-Based State-Space Models

Usage

dynr.taste(dynrModel, dynrCook = NULL, which.state, which.obs,
  conf.level = 0.99, alternative = c("two.sided", "less", "greater"),
  debug_flag = FALSE)

Arguments

Value

an object of ‘dynrTaste’ class that is a list containing lists of results from the outlier detection process. Vectors of ID and measured time points are included for later use, such as in dynr.taste2. The values, p-values, and shock points related to ‘joint’ chi-square, ‘independent’ chi-square, and t statistic for innovative and additive outliers are following in that order. The estimated delta for innovative and additive components are in the last. If debug_flag is TRUE, The by-products of the Kalman filter and smoother (Q, S, s, F_inv, N, u, r) would be added at the end. See the reference for definition of the notations. The t statistic (estimate of an outlier divided by standard error of the outlier) of the last time point is NA, because the Kalman smoothing process starts with setting r and N to zero for the last time point (core elements of calculating estimates and the standard errors of outliers) that lead to 0/0 of the t statistic of the last time point. For the time-varing models, more NAs would appear at the end of times because the Kalman smoother needs more time points to obtain all elements of r nad N from limited number of observed variables in the model.

The ‘delta_chi’ list comprises magnitude of innovative (Latent) and additive (Observed) outliers, ‘delta.L’ and ‘delta.O’, when chi-square statitics is used to detect outliers. The ‘delta_t’ list comprises magnitude of innovative (Latent) and additive (Observed) outliers, ‘delta.L’ and ‘delta.O’, when t statitics is used to detect outliers.

References

Chow, S.-M., Hamaker, E. L., & Allaire, J. C. (2009). Using innovative outliers to detect discrete shifts in dynamics in group-based state-space models. _Multivariate Behavioral Research_, 44, 465-496.

Examples

## Not run: 
# See the demo for outlier detection, OutlierDetection.R
dynrCook <- dynr.cook(dynrModel)
dynrTaste <- dynr.taste(dynrModel, dynrCook)

# Detect outliers related to 'eta1' out of, say, three latent
# variables c("eta1", "eta2", "eta3"), and all measured variables.
dynrTaste <- dynr.taste(dynrModel, dynrCook, which.state=c("eta1"))

## End(Not run)

Description

The function dynr.taste2{} update the dynrModel object applying outliers from the dynrTaste object, or outliers from users. The function then re-cook the model.

Usage

dynr.taste2(dynrModel, dynrCook, dynrTaste, delta_inn = c("t", "ind", "jnt",
  "null"), delta_add = c("t", "ind", "jnt", "null"), delta_L = NULL,
  delta_O = NULL, cook = TRUE, verbose = FALSE,
  newOutfile = "new_taste.c")

Arguments

Details

The argument dynrTaste should be the dynrTaste object that is output of the dynr.taste function the argument dynrModel is applied.

The argument dynrTaste can be NULL, if user-specified outliers are offered by the arguments delta_L and delta_O.

Value

a list with the two arguments; a new dynrModel object the outliers are applied, and a dynrCook object the new dynrModel object is cooked.

Examples

## Not run: 
# See the demo for outlier detection, OutlierDetection.R
dynrCook <- dynr.cook(dynrModel)
dynrTaste <- dynr.taste(dynrModel, dynrCook)

# Detect outliers related to 'eta1' out of, say, three latent
# variables c("eta1", "eta2", "eta3"), and all measured variables.
taste2 <- dynr.taste2(dynrModel, dynrCook, dynrTaste)

## End(Not run)

Description

A Function to perform numerical integration of the chosen ODE system, for a user-specified set of initial conditions. Plots the resulting solution(s) in the phase plane. This function from the phaseR package written by Michael J. Grayling.

Usage

dynr.trajectory(deriv, y0 = NULL, n = NULL, tlim, tstep = 0.01,
  parameters = NULL, system = "two.dim", col = "black", add = TRUE,
  state.names = c("x", "y"), ...)

Arguments

Value

Returns a list with the following components: add, col, deriv, n, parameters, system, tlim, tstep, t, x, y, ylab, y0. Most of these components correspond simply to their original input values.

The only new elements are: t = A vector containing the values of the independent variable at each integration step.

x = In the two dimensional system case, a matrix whose columns are the numerically computed values of the first dependent variable for each set of ICs.

y = In the two dimensional system case, a matrix whose columns are the numerically computed values of the second dependent variable for each initial condition. In the one dimensional system case, a matrix whose columns are the numerically computed values of the dependent variable for each initial condition.

y0 = As per input, but converted to a matrix if supplied as a vector initially.

Note

The phaseR package was taken off cran as off 10/1/2019 so we are exporting some selected functions from phaseR_2.0 published on 8/20/2018. For details of these functions please see original documentations on the phaseR package.

References

Grayling, Michael J. (2014). phaseR: An R Package for Phase Plane Analysis of Autonomous ODE Systems. The R Journal, 6(2), 43-51. DOI: 10.32614/RJ-2014-023. Available at https://doi.org/10.32614/RJ-2014-023

Description

Current Version String

Usage

dynr.version(verbose = TRUE)

Arguments

Value

A (length-one) object of class 'package_version'

Examples

dynr.version()
dynr.version(verbose=FALSE)
packageVersion("dynr")

Description

The dynrCook Class

Details

This is an internal class structure. You should not use it directly. Use dynr.cook instead.

Description

The dynrDynamics Class

Details

This is an internal class structure. The classes dynrDynamicsFormula-class and dynrDynamicsMatrix-class are subclasses of this. However, you should not use it directly. Use prep.matrixDynamics or prep.formulaDynamics instead.

Description

The dynrInitial Class

Details

This is an internal class structure. You should not use it directly. Use prep.initial instead.

Description

The dynrMeasurement Class

Details

This is an internal class structure. You should not use it directly. Use prep.measurement or prep.loadings instead.

Description

The dynrModel Class

Details

This is an internal class structure. You should not use it directly. Use dynr.model instead.

Description

The dynrNoise Class

Details

This is an internal class structure. You should not use it directly. Use prep.noise instead.

Description

The dynrRecipe Class

Details

This is an internal class structure. You should not use it directly. The following are all subclasses of this class: dynrMeasurement-class, dynrDynamics-class, dynrRegimes-class, dynrInitial-class, dynrNoise-class, and dynrTrans-class. Recipes are the things that go into a dynrModel-class using dynr.model. Use the recipe prep functions (prep.measurement, prep.formulaDynamics, prep.matrixDynamics, prep.regimes, prep.initial, prep.noise, or prep.tfun) to create these classes instead.

Description

The dynrRegimes Class

Details

This is an internal class structure. You should not use it directly. Use prep.regimes instead.

Description

The dynrTrans Class

Details

This is an internal class structure. You should not use it directly. Use prep.tfun instead.

Description

A dataset obtained and analyzed in Yang and Chow (2010).

Usage

data(EMG)

Format

A data frame with 695 rows and 4 variables

Details

Reference: Yang, M-S. & Chow, S-M. (2010). Using state-space models with regime switching to represent the dynamics of facial electromyography (EMG) data. Psychometrika, 74(4), 744-771

The variables are as follows:

• id. ID of the participant (= 1 in this case, over 695 time points)

• time Time in seconds

• iEMG. Observed integrated facial electromyograhy data

• SelfReport. Covariate - the individual's concurrent self-reports

Description

A dataset simulated using an autoregressive model of order (AR(1)) with regime-specific AR weight, intercept, and slope for a covariate. This model is a special case of Model 1 in Yang and Chow (2010) in which the moving average coefficient is set to zero.

Reference: Yang, M-S. & Chow, S-M. (2010). Using state-space models with regime switching to represent the dynamics of facial electromyography (EMG) data. Psychometrika, 74(4), 744-771

Usage

data(EMGsim)

Format

A data frame with 500 rows and 6 variables

Details

The variables are as follows:

• id. ID of the participant (= 1 in this case, over 500 time points)

• EMG. Hypothetical observed facial electromyograhy data

• self. Covariate - the individual's concurrent self-reports

• truestate. The true score of the individual's EMG at each time point

• trueregime. The true underlying regime for the individual at each time point

Description

Extend a user-specified model to include random varibles

Usage

ExpandRandomAsLVModel(dynrModel)

Arguments

Details

A dynrModel is a collection of recipes. The recipes are constructed with the functions unctions prep.formulaDynamics, prep.measurement, prep.noise, prep.initial. Additionally, data must be prepared with dynr.data and added to the model.

Value

an object of dynrModel that is the expanede model.

Examples

# model <- dynr.model(dynamics=dynm, measurement=meas, noise=mdcov,
#  initial=initial, data=data, outfile="osc.cpp")
# extended_model <- ExpandRandomAsLVModel(model)


# For full demo examples, see:
# demo(OscWithRand, package="dynr")
# demo(VDPwithRand, package="dynr")

Description

A wrapper function to call functions in the fda package to obtain smoothed estimated derivatives at a specified order

Usage

getdx(theTimes, norder, roughPenaltyMax, lambda, dataMatrix, derivOrder)

Arguments

Value

A list containing: 1. out (a matrix containing the derivative estimates at the specified order that matches the dimension of dataMatrix); 2. basisCoef (estimated basis coefficients); 3. basis2 (basis functions)

References

Chow, S-M. (2019). Practical Tools and Guidelines for Exploring and Fitting Linear and Nonlinear Dynamical Systems Models. Multivariate Behavioral Research. https://www.nihms.nih.gov/pmc/articlerender.fcgi?artid=1520409

Chow, S-M., *Bendezu, J. J., Cole, P. M., & Ram, N. (2016). A Comparison of Two- Stage Approaches for Fitting Nonlinear Ordinary Differential Equation (ODE) Models with Mixed Effects. Multivariate Behavioral Research, 51, 154-184. Doi: 10.1080/00273171.2015.1123138.

Examples

data("LinearOsc")
# Number of subjects is 10
numP <- length(unique(LinearOsc$ID))
# Number of time points is 100
numT <- max(table(LinearOsc$ID))
out2 <- matrix(LinearOsc$x, ncol=numP, byrow=FALSE)
theTimes <- LinearOsc$theTimes[1:numT]
# Order of Bsplines - usually 2 higher than roughPenaltyMax
norder <- 6
# Penalization order
roughPenaltyMax <- 4 
# Pick lambda value that gives the low GCV
# Could/should use plotGCV instead
sp <- 1/2
# Smoothed level
x <- getdx(theTimes, norder, roughPenaltyMax, sp, out2, 0)[[1]]
# Smoothed 1st derivs
dx <- getdx(theTimes, norder, roughPenaltyMax, sp, out2, 1)[[1]]
# Smoothed 2nd derivs
d2x = getdx(theTimes, norder, roughPenaltyMax, sp, out2, 2)[[1]]

Description

Principally, this function takes a host of arguments and gives back a list that importantly includes the function addresses.

Usage

internalModelPrep(num_regime, dim_latent_var, xstart, ub, lb,
  options = default.model.options, isContinuousTime, infile, outfile,
  compileLib, verbose)

Arguments

Value

A list of model statements to be passed to dynr.cook().

Description

The variables are as follows:

Usage

data(LinearOsc)

Format

A data frame with 1000 rows and 3 variables

Details

• ID. ID of the systems (1 to 10)

• x. Latent level variable

• theTimes. Measured time Points

Examples

# The following was used to generate the data
#--------------------------------------
## Not run: 
Osc <- function(t, prevState, parms) {
  x1 <- prevState[1] # x1[t]
  x2 <- prevState[2] # x2[t]
  eta1 = parms[1]
  zeta1 = parms[2]
  with(as.list(parms), {
   dx1 <- x2
    dx2 <- eta1*x1 + zeta1*x2 
    res<-c(dx1,dx2)
    list(res)
  }
  )
}
n = 10 #Number of subjects
T = 100 #Number of time points
deltaT = .1 #dt
lastT = deltaT*T #Value of t_{i,T}
theTimes  = seq(0, lastT, length=T)  #A list of time values

eta = -.8
zeta = -.1
out1 = matrix(NA,T*n,1)
trueOut = matrix(NA,T*n,1)
parms = c(eta, zeta)
  for (i in 1:n){
  xstart = c(rnorm(1,0,2),rnorm(1,0,.5))
  out <- lsoda(as.numeric(xstart), theTimes, Osc, parms)
  trueOut[(1+(i-1)*T):(i*T)] = out[,2]
  out1[(1+(i-1)*T):(i*T)] = out[,2]+rnorm(T,0,1)
  }

LinearOsc= data.frame(ID=rep(1:n,each=T),x=out1[,1],
                  theTimes=rep(theTimes,n))
save(LinearOsc,file="LinearOsc.rda")

## End(Not run)

Description

A dataset simulated using a continuous-time stochastic linear damped oscillator model. The variables are as follows:

Usage

data(LogisticSetPointSDE)

Format

A data frame with 2410 rows and 6 variables

Details

• id. ID of the systems (1 to 10)

• times. Time index (241 time points for each system)

• x. Latent level variable

• y. Latent first derivative variable

• z. True values of time-varying setpoints

• obsy. Observed level

Examples

# The following was used to generate the data
#---------------------------------------
## Not run: 
require(Sim.DiffProc)
freq <- -1
damp <- -.1
mu <- -2
r <- .5
b <- .1
sigma1 <- 0.1
sigma2 <- 0.1
fx <- expression(y, freq*(x-z) + damp*y, r*z*(1-b*z))
gx <- expression(0, sigma1, 0)
r3dall <- c()
for (j in 1:10){
  r3dtemp <- c(-5,0,.1)
  r3d <- r3dtemp
  for (i in seq(0.125, 30, by=0.125)){
    mod3dtemp <- snssde3d(drift=fx, diffusion=gx, M=1, t0=i-0.125,
        x0=as.numeric(r3dtemp), T=i, N=500, type="str",
        method="smilstein")
    r3dtemp <- rsde3d(mod3dtemp,at=i)
    r3d <-rbind(r3d,r3dtemp)
  }
  r3dall <- rbind(r3dall, cbind(r3d, id=j))
}

r3dall$obsy <- r3dall$x+rnorm(length(r3dall$x),0,1)
write.table(r3dall, file="LogisticSetPointSDE.txt")

## End(Not run)

Description

Extract the log likelihood from a dynrCook Object

Usage

## S3 method for class 'dynrCook'
logLik(object, ...)

## S3 method for class 'dynrCook'
deviance(object, ...)

Arguments

Details

The 'df' attribute for this object is the number of freely estimated parameters. The 'nobs' attribute is the total number of rows of data, adding up the number of time points for each person.

The deviance method returns minus two times the log likelihood.

Value

In the case of logLik, an object of class logLik.

See Also

Other S3 methods coef.dynrCook

Examples

# Minimal model
require(dynr)

meas <- prep.measurement(
	values.load=matrix(c(1, 0), 1, 2),
	params.load=matrix(c('fixed', 'fixed'), 1, 2),
	state.names=c("Position","Velocity"),
	obs.names=c("y1"))

ecov <- prep.noise(
	values.latent=diag(c(0, 1), 2),
	params.latent=diag(c('fixed', 'dnoise'), 2),
	values.observed=diag(1.5, 1),
	params.observed=diag('mnoise', 1))

initial <- prep.initial(
	values.inistate=c(0, 1),
	params.inistate=c('inipos', 'fixed'),
	values.inicov=diag(1, 2),
	params.inicov=diag('fixed', 2))

dynamics <- prep.matrixDynamics(
	values.dyn=matrix(c(0, -0.1, 1, -0.2), 2, 2),
	params.dyn=matrix(c('fixed', 'spring', 'fixed', 'friction'), 2, 2),
	isContinuousTime=TRUE)

data(Oscillator)
data <- dynr.data(Oscillator, id="id", time="times", observed="y1")

model <- dynr.model(dynamics=dynamics, measurement=meas,
	noise=ecov, initial=initial, data=data)

## Not run: 
cook <- dynr.cook(model,
	verbose=FALSE, optimization_flag=FALSE, hessian_flag=FALSE)

# Now get the log likelihood!
logLik(cook)

## End(Not run)

Description

Extract the free parameter names of a dynrCook object

Usage

## S4 method for signature 'dynrCook'
names(x)

Arguments

Description

Extract the free parameter names of a dynrModel object

Usage

## S4 method for signature 'dynrModel'
names(x)

Arguments

Description

Extract the number of observations for a dynrCook object

Usage

## S3 method for class 'dynrCook'
nobs(object, ...)

Arguments

Details

We return the total number of rows of data, adding up the number of time points for each person. For some purposes, you may want the mean number of observations per person or the number of people instead. These are not currently supported via nobs.

Value

A single number. The total number of observations across all IDs.

Examples

# Minimal model
require(dynr)

meas <- prep.measurement(
	values.load=matrix(c(1, 0), 1, 2),
	params.load=matrix(c('fixed', 'fixed'), 1, 2),
	state.names=c("Position","Velocity"),
	obs.names=c("y1"))

ecov <- prep.noise(
	values.latent=diag(c(0, 1), 2),
	params.latent=diag(c('fixed', 'dnoise'), 2),
	values.observed=diag(1.5, 1),
	params.observed=diag('mnoise', 1))

initial <- prep.initial(
	values.inistate=c(0, 1),
	params.inistate=c('inipos', 'fixed'),
	values.inicov=diag(1, 2),
	params.inicov=diag('fixed', 2))

dynamics <- prep.matrixDynamics(
	values.dyn=matrix(c(0, -0.1, 1, -0.2), 2, 2),
	params.dyn=matrix(c('fixed', 'spring', 'fixed', 'friction'), 2, 2),
	isContinuousTime=TRUE)

data(Oscillator)
data <- dynr.data(Oscillator, id="id", time="times", observed="y1")

model <- dynr.model(dynamics=dynamics, measurement=meas,
	noise=ecov, initial=initial, data=data)

## Not run: 
cook <- dynr.cook(model,
	verbose=FALSE, optimization_flag=FALSE, hessian_flag=FALSE)

# Now get the total number of observations
nobs(cook)

## End(Not run)

Description

Extract the number of observations for a dynrModel object

Usage

## S3 method for class 'dynrModel'
nobs(object, ...)

Arguments

Details

We return the total number of rows of data, adding up the number of time points for each person. For some purposes, you may want the mean number of observations per person or the number of people instead. These are not currently supported via nobs.

Value

A single number. The total number of observations across all IDs.

Examples

# Create a minimal uncooked model called 'model'
# That is, without esimating parameters
require(dynr)

meas <- prep.measurement(
	values.load=matrix(c(1, 0), 1, 2),
	params.load=matrix(c('fixed', 'fixed'), 1, 2),
	state.names=c("Position","Velocity"),
	obs.names=c("y1"))

ecov <- prep.noise(
	values.latent=diag(c(0, 1), 2),
	params.latent=diag(c('fixed', 'dnoise'), 2),
	values.observed=diag(1.5, 1),
	params.observed=diag('mnoise', 1))

initial <- prep.initial(
	values.inistate=c(0, 1),
	params.inistate=c('inipos', 'fixed'),
	values.inicov=diag(1, 2),
	params.inicov=diag('fixed', 2))

dynamics <- prep.matrixDynamics(
	values.dyn=matrix(c(0, -0.1, 1, -0.2), 2, 2),
	params.dyn=matrix(c('fixed', 'spring', 'fixed', 'friction'), 2, 2),
	isContinuousTime=TRUE)

data(Oscillator)
data <- dynr.data(Oscillator, id="id", time="times", observed="y1")

model <- dynr.model(dynamics=dynamics, measurement=meas,
	noise=ecov, initial=initial, data=data)

# Now get the total number of observations!
nobs(model)

Description

A dataset simulated using a discrete-time nonlinear dynamic factor analysis model with 6 observed indicators for identifying two latent factors: individuals' positive and negative emotions. Proposed by Chow and Zhang (2013), the model was inspired by models of affect and it posits that the two latent factors follow a vector autoregressive process of order 1 (VAR(1)) with parameters that vary between two possible regimes: (1) an "independent" regime in which the lagged influences between positive and negative emotions are zero; (2) a "high-activation" regime to capture instances on which the lagged influences between PA and NA intensify when an individual's previous levels of positive and negative emotions were unusually high or low (see Model 2 in Chow & Zhang).

Reference: Chow, S-M, & Zhang, G. (2013). Regime-switching nonlinear dynamic factor analysis models. Psychometrika, 78(4), 740-768.

Usage

data(NonlinearDFAsim)

Format

A data frame with 3000 rows and 8 variables

Details

• id. ID of the participant (1 to 10)

• time. Time index (300 time points from each subject)

• y1-y3. Observed indicators for positive emotion

• y4-y6. Observed indicators for negative emotion

Description

The variables are as follows:

Usage

data(oscData)

Format

A data frame with 1,800 rows and 6 variables

Details

• id. Person ID

• times. Continuous time of measurement

• y1. Observed score 1

• u1. Covariate 1

• u2. Covariate 2

• trueb. True value of person-specific random effect

Description

A dataset simulated using a damped linear oscillator model in continuous time with 1 observed indicator for identifying two latent factors (position and velocity). The variables are as follows:

Usage

data(Oscillator)

Format

A data frame with 1000 rows and 5 variables

Details

• id. ID of the systems (1 to 1 because this is a single person)

• y1. Noisy observed position

• times. Time index (1000 time points) spaced at one unit intervals

• x1. True latent position

• x2. True latent velocity

Examples

# The following was used to generate the data
#--------------------------------------
# Data Generation
## Not run: 
require(mvtnorm)
require(Matrix)

xdim <- 2
udim <- 1
ydim <- 1
tdim <- 1000
set.seed(315)
tA <- matrix(c(0, -.3, 1, -.7), xdim, xdim)
tB <- matrix(c(0), xdim, udim)
tC <- matrix(c(1, 0), ydim, xdim)
tD <- matrix(c(0), ydim, udim)
tQ <- matrix(c(0), xdim, xdim); diag(tQ) <- c(0, 2.2)
tR <- matrix(c(0), ydim, ydim); diag(tR) <- c(1.5)

x0 <- matrix(c(0, 1), xdim, 1)
P0 <- diag(c(1), xdim)
tdx <- matrix(0, xdim, tdim+1)
tx <- matrix(0, xdim, tdim+1)
tu <- matrix(0, udim, tdim)
ty <- matrix(0, ydim, tdim)

tT <- matrix(0:tdim, nrow=1, ncol=tdim+1)

tI <- diag(1, nrow=xdim)

tx[,1] <- x0
for(i in 2:(tdim+1)){
	q <- t(rmvnorm(1, rep(0, xdim), tQ))
	tdx[,i] <- tA %*% tx[,i-1] + tB %*% tu[,i-1] + q
	expA <- as.matrix(expm(tA * (tT[,i]-tT[,i-1])))
	intA <- solve(tA) %*% (expA - tI)
	tx[,i] <- expA %*% tx[, i-1] + intA %*% tB %*% tu[,i-1] + intA %*% q
	ty[,i-1] <- tC %*% tx[,i] + tD %*% tu[,i-1] + t(rmvnorm(1, rep(0, ydim), tR))
}



rownames(ty) <- paste('y', 1:ydim, sep='')
rownames(tx) <- paste('x', 1:xdim, sep='')
simdata <- cbind(id=rep(1, tdim), t(ty), times=tT[,-1], t(tx)[-1,])
 write.table(simdata, file='Oscillator.txt', row.names=FALSE, col.names=TRUE)

plot(tx[1,], type='l')
plot(tT[,-1], ty[1,], type='l')

## End(Not run)

Description

This is a list object containing true outliers, the dataset, and the saved result from running dynr.taste.

Usage

data(Outliers)

Format

A data frame with 6000 rows and 6 variables

Details

The true outliers for observed variables are saved in ‘Outliers$generated$shockO’.

• id. Six outliers were added for each ID.

• time_O. Time points where the outliers were added.

• obs. Variable indices where the outliers were added.

• shock.O. The magnitude of outliers.

The true outliers for state variables are saved in ‘Outliers$generated$shockL’.

• id. Three outliers were added for each ID.

• time_L. Time points where the outliers were added.

• lat. Variable indices where the outliers were added.

• shock.L. The magnitude of outliers.

A dataset simulated based on state-space model including the outliers. The data is saved in ‘Outliers$generated$y’. The variables are as follows:

• id. ID of the systems (1 to 100)

• times. Time indices (100 time points for each participant)

• V1 - V6. observed variables

The detected innovative outliers from dynr.taste for this dataset, which is used for testing whether the dynr.taste replicate the same result. The data is saved in ‘Outliers$detect_O’. The variables are as follows:

• id. IDs

• time_L. Time points where the outliers were detected

• obs. Variable indices for observed variables where the outliers were detected

The detected additive outliers from dynr.taste for this dataset, which is used for testing whether the dynr.taste replicate the same result. The data is saved in ‘Outliers$detect_L’. The variables are as follows:

• id. IDs

• time_L. Time points where the outliers were detected

• obs. Variable indices for latent variables where the outliers were detected

Examples

## Not run: 
 #The following was used to generate the data
 #---------------------------------------
 lambda <- matrix(c(1.0, 0.0,
 0.9, 0.0,
 0.8, 0.0,
 0.0, 1.0,
 0.0, 0.9,
 0.0, 0.8), ncol=2, byrow=TRUE)
 psi <- matrix(c(0.3, -0.1,
                 -0.1, 0.3), ncol=2, byrow=TRUE)
 beta <- matrix(c(0.8, -0.2,
                  -0.2,  0.7), ncol=2, byrow=TRUE)
 theta <- diag(c(0.2, 0.2, 0.2, 0.2, 0.2, 0.2), ncol=6, nrow=6)
 nlat <- 2; nobs <- 6
 mean_0 <- rep(0, nlat)
 psi_inf <- diag(1, 2*2) - kronecker(beta, beta)
 psi_inf_inv <- try(solve(psi_inf), silent=TRUE)
 if("try-error" %in% class(psi_inf_inv)) {
   psi_inf_inv <- MASS::ginv(psi_inf)}
 psi_0 <- psi_inf_inv %*% as.vector(psi)
 dim(psi_0) <- c(2, 2)
 # measurement error covariance matrix
 mea_cov <- lambda %*% psi_0 %*% t(lambda) + theta
 resL <- lapply(1:100, function(subj) {
   # initial state
   eta_0 <- mvtnorm::rmvnorm(1, mean=mean_0, sigma=psi_0)#[1,nlat]
   zeta_0 <- mvtnorm::rmvnorm(1, mean=rep(0, nlat), sigma=psi)
   eta <- matrix(0, nrow=time, ncol=nlat)
   eta[1, ] <- beta %*% t(eta_0) + t(zeta_0) 
   zeta <- mvtnorm::rmvnorm(time, mean=rep(0, nlat), sigma=psi)
   # random shock generation
   # to avoid shock appearing too early or late (first and last 3)
   shkLat_time <- sample(4:(time-3), nshockLat)
   shk_lat <- sample(1:nlat, nshockLat, replace=TRUE)
   shockLatIdx <- matrix(c(shkLat_time, shk_lat), ncol=2)
   shockSignL <- sample(c(1,-1), nshockLat, replace=TRUE)
   colnames(shockLatIdx) <- c("time_L","lat")
   shockLatV <- shockSignL*( shockMag*sqrt(diag(shockPsi)))[shockLatIdx[,"lat"]]
   shockLatM <- matrix(0, time, nlat)
   shockLatM[shockLatIdx] <- shockLatV
   shkObs_time <- sample(4:(time-3), nshockObs)
   shk_obs <- sample(1:nobs, nshockObs, replace=TRUE)
   shockObsIdx <- matrix(c(shkObs_time, shk_obs), ncol=2)
   shockSignO <- sample(c(1,-1), nshockObs, replace=TRUE)
   colnames(shockObsIdx) <- c("time_O","obs")
   shockObsV <- shockSignO*( shockMag*sqrt(diag(mea_cov)) )[shockObsIdx[,"obs"]]
   shockObsM <- matrix(0, time, nobs)
   shockObsM[shockObsIdx] <- shockObsV
   # generate state process WITH shock
   for (t in 1:(time-1)) {
     eta[t+1, ] <- shockLatM[t, ] + beta %*% eta[t, ] + zeta[t, ]
   }
   # generate observed process
   y <- shockObsM + eta %*% t(lambda) +
     mvtnorm::rmvnorm(time, mean=rep(0, nobs), sigma=theta)# epsilon
 }
 
## End(Not run)

Description

A multiple subject dataset simulated using a two factor process factor analysis model in discrete time with 6 observed indicators for identifying two latent factors. The variables are as follows:

Usage

data(PFAsim)

Format

A data frame with 2,500 rows and 10 variables

Details

• ID. Person ID variable (1 to 50 because there are 50 simulated people)

• Time. Time ID variable (1 to 50 because there are 50 time points)

• V1. Noisy observed variable 1

• V2. Noisy observed variable 2

• V3. Noisy observed variable 3

• V4. Noisy observed variable 4

• V5. Noisy observed variable 5

• V6. Noisy observed variable 6

• F1. True latent variable 1 scores

• F2. True latent variable 2 scores

Variables V1, V2, and V3 load on F1, whereas variables V4, V5, V6 load on F2. The true values of the factor loadings are 1, 2, 1, 1, 2, and 1, respectively. The true measurement error variance is 0.5 for all variables. The true dynamic noise covariance has F1 with a variance of 2.77, F2 with a variance of 8.40, and their covariance is 2.47. The across-time dynamics have autoregressive effects of 0.5 for both F1 and F2 with a cross-lagged effect from F1 to F2 at 0.4. The cross-lagged effect from F2 to F1 is zero. The true initial latent state distribution as mean zero and a diagonal covariance matrix with var(F1) = 2 and var(F2) = 1. The generating model is the same for all individuals.

Examples

# The following was used to generate the data
## Not run: 
set.seed(12345678)
library(mvtnorm)
# setting up matrices
time      <- 50
# Occasions to throw out to wash away the effects of initial condition
npad      <- 0
np        <- 50
ne        <- 2 #Number of latent variables
ny        <- 6 #Number of manifest variables
# Residual variance-covariance matrix
psi       <- matrix(c(2.77, 2.47,
                      2.47, 8.40),
                    ncol = ne, byrow = T)
# Lambda matrix containing contemporaneous relations among
# observed variables and 2 latent variables. 
lambda    <- matrix(c(1, 0,
                      2, 0,
                      1, 0,
                      0, 1,
                      0, 2,
                      0, 1),
                    ncol = ne, byrow = TRUE)
# Measurement error variances
theta     <- diag(.5, ncol = ny, nrow = ny)
# Lagged directed relations among variables
beta      <- matrix(c(0.5, 0,
                      0.4, 0.5), 
                    ncol = ne, byrow = TRUE)
a0        <- mvtnorm::rmvnorm(1, mean = c(0, 0),
                                 sigma = matrix(c(2,0,0,1),ncol=ne))
yall <- matrix(0,nrow = time*np, ncol = ny)
eall <- matrix(0,nrow = time*np, ncol = ne)
for (p in 1:np){
  # Latent variable residuals
  zeta      <- mvtnorm::rmvnorm(time+npad, mean = c(0, 0), sigma = psi)
  # Measurement errors
  epsilon   <- rmvnorm(time, mean = c(0, 0, 0, 0, 0, 0), sigma = theta)
  # Set up matrix for contemporaneous variables
  etaC      <- matrix(0, nrow = ne, ncol = time + npad)
  # Set up matrix for lagged variables
  etaL      <- matrix(0, nrow = ne, ncol = time + npad + 1)
  
  etaL[,1]  <- a0
  etaC[,1] <- a0
  # generate factors
  for (i in 2:(time+npad)){
    etaL[ ,i] <- etaC[,i-1]
    etaC[ ,i]   <- beta %*% etaL[ ,i] + zeta[i, ]
  }
  etaC <- etaC[,(npad+1):(npad+time)]
  eta <- t(etaC)
  
  # generate observed series
  y   <- matrix(0, nrow = time, ncol = ny)
  for (i in 1:nrow(y)){
    y[i, ] <- lambda %*% eta[i, ] + epsilon[i, ]
  }
  yall[(1+(p-1)*time):(p*time),] <- y
  eall[(1+(p-1)*time):(p*time),] <- eta
}
yall <- cbind(rep(1:np,each=time),rep(1:time,np),yall)
yeall <- cbind(yall,eall)
write.table(yeall,'PFAsim.txt',row.names=FALSE,
  col.names=c("ID", "Time", paste0("V", 1:ny), paste0("F", 1:ne)))

## End(Not run)

Description

Plot method for dynrCook objects

Usage

## S3 method for class 'dynrCook'
plot(x, dynrModel, style = 1, names.state, names.observed,
  printDyn = TRUE, printMeas = TRUE, textsize = 4, ...)

Arguments

Details

This is a wrapper around dynr.ggplot. A great benefit of it is that it shows the model equations in a plot.

Value

ggplot object.

Description

Plot the formula from a model

Usage

plotFormula(dynrModel, ParameterAs, printDyn = TRUE, printMeas = TRUE,
  printRS = FALSE, textsize = 4)

Arguments

Details

This function typesets a set of formulas that represent the model. Typical inputs to the ParameterAs argument are (1) the starting values for a model, (2) the final estimated values for a model, and (3) the parameter names. These are accessible with (1) model$xstart, (2) coef(cook), and (3) model$param.names or names(coef(cook)), respectively.

Value

ggplot object

Description

A function to evaluate the generalized cross-validation (GCV) values associated with derivative estimates via Bsplines at a range of specified smoothing parameter (lambda) values

Usage

plotGCV(theTimes, norder, roughPenaltyMax, dataMatrix, lowLambda, upLambda,
  lambdaInt, isPlot)

Arguments

Value

A data frame containing: 1. lambda values; 2. edf (effective degrees of freedom); 3. GCV (Generalized cross-validation value as averaged across units (e.g., subjects))

References

Chow, S-M. (2019). Practical Tools and Guidelines for Exploring and Fitting Linear and Nonlinear Dynamical Systems Models. Multivariate Behavioral Research. https://www.nihms.nih.gov/pmc/articlerender.fcgi?artid=1520409

Chow, S-M., *Bendezu, J. J., Cole, P. M., & Ram, N. (2016). A Comparison of Two- Stage Approaches for Fitting Nonlinear Ordinary Differential Equation (ODE) Models with Mixed Effects. Multivariate Behavioral Research, 51, 154-184. Doi: 10.1080/00273171.2015.1123138.

Description

A dataset simulated using a continuous-time nonlinear predator-and-prey model with 2 observed indicators for identifying two latent factors. The variables are as follows:

Usage

data(PPsim)

Format

A data frame with 1000 rows and 6 variables

Details

• id. ID of the systems (1 to 20)

• time. Time index (50 time points for each system)

• prey. The true population of the prey species

• predator. The true population of the predator species

• x. Observed indicator for the population of the prey species

• y. Observed indicator for the population of the predator species

Description

predict method for dynrModel objects

Usage

## S3 method for class 'dynrModel'
predict(object, newdata = NULL, interval = c("none",
  "confidence", "prediction"), method = c("kalman", "ensemble"),
  level = 0.95, type = c("latent", "observed"), ...)

Arguments

Details

The newdata argument is either a data.frame or ts object. It passed as the dataframe argument of dynr.data and must accept the same further arguments as the data in the model passed in the object argument (e.g., same id, time, observed, and covariates arguments).

The available methods for prediction are 'kalman' and 'ensemble'. The 'kalman' method uses the Kalman filter to create predictions. The 'ensemble' method simulates a set of initial conditions and lets those run forward in time. The distribution of this ensemble provides the predictions. The mean is the value predicted. The quantiles of the distribution provide the intervals.

Value

A list of the prediction estimates, intervals, and ensemble members.

Description

Recipe function for specifying dynamic functions using formulas

Usage

prep.formulaDynamics(formula, startval = numeric(0),
  isContinuousTime = FALSE, jacobian, ...)

Arguments

Details

This function defines the dynamic functions of the model either in discrete time or in continuous time. The function can be either linear or nonlinear, with free or fixed parameters, numerical constants, covariates, and other mathematical functions that define the dynamics of the latent variables. Every latent variable in the model needs to be defined by a differential (for continuous time model), or difference (for discrete time model) equation. The names of the latent variables should match the specification in prep.measurement(). For nonlinear models, the estimation algorithm generally needs a Jacobian matrix that contains elements of first differentiations of the dynamic functions with respect to the latent variables in the model. For most nonlinear models, such differentiations can be handled automatically by dynr. However, in some cases, such as when the absolute function (abs) is used, the automatic differentiation would fail and the user may need to provide his/her own Jacobian functions. When theta.formula and other accompanying elements in "..." are provided, the program automatically inserts the random effect components specified in random.names as additional latent (state) variables in the model, and estimate (cook) this expanded model. Do check that the expanded model satisfies conditions such as observability for the estimation to work.

Value

Object of class 'dynrDynamicsFormula'

Examples

# In this example, we present how to define the dynamics of a bivariate dual change score model
# (McArdle, 2009). This is a linear model and the user does not need to worry about 
# providing any jacobian function (the default). 
 
# We start by creating a list of formula that describes the model. In this model, we have four 
# latent variables, which are "readLevel", "readSlope", "mathLevel", and "math Slope".  The right-
# hand side of each formula gives a function that defines the dynamics.   
 
 formula <- list(
          list(readLevel~ (1+beta.read)*readLevel + readSlope + gamma.read*mathLevel,
          readSlope~ readSlope,
          mathLevel~ (1+beta.math)*mathLevel + mathSlope + gamma.math*readLevel, 
          mathSlope~ mathSlope
          ))

# Then we use prep.formulaDynamics() to define the formula, starting value of the parameters in
# the model, and state the model is in discrete time by setting isContinuousTime=FALSE.
 
dynm  <- prep.formulaDynamics(formula=formula,
                             startval=c(beta.read = -.5, beta.math = -.5, 
                                        gamma.read = .3, gamma.math = .03
                             ), isContinuousTime=FALSE)


# For a full demo example of regime switching nonlinear discrete time model, you
# may refer to a tutorial on 
# \url{https://quantdev.ssri.psu.edu/tutorials/dynr-rsnonlineardiscreteexample}

#Not run: 
#For a full demo example that uses user-supplied analytic jacobian functions see:
#demo(RSNonlinearDiscrete, package="dynr")
formula <- list(
    list(
      x1 ~ a1*x1,
      x2 ~ a2*x2),
    list(
      x1 ~ a1*x1 + c12*(exp(abs(x2)))/(1+exp(abs(x2)))*x2,
      x2 ~ a2*x2 + c21*(exp(abs(x1)))/(1+exp(abs(x1)))*x1)
  )
jacob <- list(
  list(x1~x1~a1,
      x2~x2~a2),
  list(x1~x1~a1,
      x1~x2~c12*(exp(abs(x2))/(exp(abs(x2))+1)+x2*sign(x2)*exp(abs(x2))/(1+exp(abs(x2))^2)),
      x2~x2~a2,
      x2~x1~c21*(exp(abs(x1))/(exp(abs(x1))+1)+x1*sign(x1)*exp(abs(x1))/(1+exp(abs(x1))^2))))
dynm <- prep.formulaDynamics(formula=formula, startval=c( a1=.3, a2=.4, c12=-.5, c21=-.5),
                             isContinuousTime=FALSE, jacobian=jacob)

#For a full demo example that uses automatic jacobian functions (the default) see:
#demo(RSNonlinearODE , package="dynr")
formula=list(prey ~ a*prey - b*prey*predator, predator ~ -c*predator + d*prey*predator)
dynm <- prep.formulaDynamics(formula=formula,
                          startval=c(a = 2.1, c = 0.8, b = 1.9, d = 1.1),
                          isContinuousTime=TRUE)

#For a full demo example that includes unit-specific random effects in theta.formula see:
#demo(OscWithRand, package="dynr")
formula <- list(x ~ dx,
               dx ~ eta_i * x + zeta*dx)
theta.formula  = list (eta_i ~ 1 * eta0  + u1 * eta1 + u2 * eta2 + 1 * b_eta)
dynm <- prep.formulaDynamics(formula=formula,
                           startval=c(eta0=-1, eta1=.1, eta2=-.1,zeta=-.02),
                           isContinuousTime=TRUE,
                           theta.formula=theta.formula,
                           random.names=c('b_eta'),
                           random.params.inicov=matrix(c('sigma2_b_eta'), ncol=1,byrow=TRUE),
                           random.values.inicov=matrix(c(0.1), ncol=1,byrow=TRUE))

Description

Recipe function for preparing the initial conditions for the model.

Usage

prep.initial(values.inistate, params.inistate, values.inicov, params.inicov,
  values.regimep = 1, params.regimep = 0, covariates, deviation = FALSE,
  refRow)

Arguments

Details

The initial condition model includes specifications for the intial state vector, initial error covariance matrix, initial probabilities of being in each regime and all associated parameter specifications. The initial probabilities are specified in multinomial logistic regression form. When there are no covariates, this implies multinomial logistic regression with intercepts only. In particular, the initial probabilities not not specified on a 0 to 1 probability scale, but rather a negative infinity to positive infinity log odds scale. Fixing an initial regime probability to zero does not mean zero probability. It translates to a comparison log odds scale against which other regimes will be judged.

The structure of the initial state vector and the initial probability vector depends on the presence of covariates. When there are no covariates these should be vectors, or equivalently single-column matrices. When there are covariates they should have $c+1$ columns for $c$ covariates. For values.regimep and params.regimep the number of rows should be the number of regimes. For inistate and inicov the number of rows should be the number of latent states. Of course, inicov is a square and symmetric so its number of rows should be the same as its number of columns.

When deviation=FALSE, the non-deviation form of the multinomial logistic regression is used. This form has a separate intercept term for each entry of the initial probability vector. When deviation=TRUE, the deviation form of the multinomial logistic regression is used. This form has an intercept term that is common to all rows of the initial probability vector. The rows are then distinguished by their own individual deviations from the common intercept. The deviation form requires the same reference row constraint as the non-deviation form (described below). By default the reference row is taken to be the row with all zero covariate effects. Of course, if there are no covariates and the deviation form is desired, then the user must provide the reference row.

The refRow argument determines which row is used as the intercept row. It is only used in the deviation form (i.e. deviation=TRUE). In the deviation form, one row of values.regimep and params.regimep contains the intercepts, other rows contain deviations from these intercepts. The refRow argument says which row contains the intercept terms. The default behavior for refRow is to detect the reference row automatically based on which parameters are fixed. If we have problems detecting which is the reference row, then we provide error messages that are as helpful as we can make them.

Value

Object of class 'dynrInitial'

See Also

Methods that can be used include: print, printex, show

Examples

#### No-covariates
# Single regime, no covariates
# latent states are position and velocity
# initial position is free and called 'inipos'
# initial slope is fixed at 1
# initial covariance is fixed to a diagonal matrix of 1s
initialNoC <- prep.initial(
	values.inistate=c(0, 1),
	params.inistate=c('inipos', 'fixed'),
	values.inicov=diag(1, 2),
	params.inicov=diag('fixed', 2))

#### One covariate
# Single regime, one covariate on the inital mean
# latent states are position and velocity
# initial covariance is fixed to a diagonal matrix of 1s
# initial latent means have
#   nrow = numLatentState, ncol = numCovariates + 1
# initial position has free intercept and free u1 effect
# initial slope is fixed at 1
initialOneC <- prep.initial(
	values.inistate=matrix(
		c(0, .5,
		  1,  0), byrow=TRUE,
		nrow=2, ncol=2),
	params.inistate=matrix(
		c('iniPosInt', 'iniPosSlopeU1',
		'fixed', 'fixed'), byrow=TRUE,
		nrow=2, ncol=2),
	values.inicov=diag(1, 2),
	params.inicov=diag('fixed', 2),
	covariates='u1')

#### Regime-switching, one covariate
# latent states are position and velocity
# initial covariance is fixed to a diagonal matrix of 1s
# initial latent means have
#   nrow = numLatentState, ncol = numCovariates + 1
# initial position has free intercept and free u1 effect
# initial slope is fixed at 1
# There are 3 regimes but the mean and covariance
#   are not regime-switching.
initialRSOneC <- prep.initial(
	values.regimep=matrix(
		c(1, 1,
		  0, 1,
		  0, 0), byrow=TRUE,
		nrow=3, ncol=2),
	params.regimep=matrix(
		c('r1int', 'r1slopeU1',
		  'r2int', 'r2slopeU2',
		  'fixed', 'fixed'), byrow=TRUE,
		nrow=3, ncol=2),
	values.inistate=matrix(
		c(0, .5,
		  1,  0), byrow=TRUE,
		nrow=2, ncol=2),
	params.inistate=matrix(
		c('iniPosInt', 'iniPosSlopeU1',
		'fixed', 'fixed'), byrow=TRUE,
		nrow=2, ncol=2),
	values.inicov=diag(1, 2),
	params.inicov=diag('fixed', 2),
	covariates='u1')

Description

Recipe function to quickly create factor loadings

Usage

prep.loadings(map, params = NULL, idvar, exo.names = character(0),
  intercept = FALSE)

Arguments

Details

The default pattern for 'idvar' is to fix the first factor loading for each factor to one. The variable names listed in 'idvar' have their factor loadings fixed to one. However, if the names of the latent variables are used for 'idvar', then all the factor loadings will be freely estimated and you should fix the factor variances in the noise part of the model (e.g. prep.noise).

This function does not have the full set of features possible in the dynr package. In particular, it does not have any regime-swtiching. Covariates can be included with the exo.names argument, but all covariate effects are freely estimated and the starting values are all zero. Likewise, intercepts can be included with the intercept logical argument, but all intercept terms are freely estimated with zero as the starting value. For complete functionality use prep.measurement.

Value

Object of class 'dynrMeasurement'

Examples

#Single factor model with one latent variable fixing first loading
prep.loadings(list(eta1=paste0('y', 1:4)), paste0("lambda_", 2:4))

#Single factor model with one latent variable fixing the fourth loading
prep.loadings(list(eta1=paste0('y', 1:4)), paste0("lambda_", 1:3), idvar='y4')

#Single factor model with one latent variable freeing all loadings
prep.loadings(list(eta1=paste0('y', 1:4)), paste0("lambda_", 1:4), idvar='eta1')

#Single factor model with one latent variable fixing first loading
# and freely estimated intercept
prep.loadings(list(eta1=paste0('y', 1:4)), paste0("lambda_", 2:4),
 intercept=TRUE)

#Single factor model with one latent variable fixing first loading
# and freely estimated covariate effects for u1 and u2
prep.loadings(list(eta1=paste0('y', 1:4)), paste0("lambda_", 2:4),
 exo.names=paste0('u', 1:2))

# Two factor model with simple structure
prep.loadings(list(eta1=paste0('y', 1:4), eta2=paste0('y', 5:7)), 
paste0("lambda_", c(2:4, 6:7)))

#Two factor model with repeated use of a free parameter
prep.loadings(list(eta1=paste0('y', 1:4), eta2=paste0('y', 5:8)), 
paste0("lambda_", c(2:4, 6:7, 4)))

#Two factor model with a cross loading
prep.loadings(list(eta1=paste0('y', 1:4), eta2=c('y5', 'y2', 'y6')), 
paste0("lambda_", c("21", "31", "41", "22", "62")))

Description

Recipe function for creating Linear Dynamics using matrices

Usage

prep.matrixDynamics(params.dyn = NULL, values.dyn, params.exo = NULL,
  values.exo = NULL, params.int = NULL, values.int = NULL, covariates,
  isContinuousTime)

Arguments

Details

A recipe function for specifying the deterministic portion of a set of linear dynamic functions as:

Discrete-time model: eta(t+1) = int + dyn*eta(t) + exo*x(t), where eta(t) is a vector of latent variables, x(t) is a vector of covariates, int, dyn, and exo are vectors and matrices specified via the arguments *.int, *.dyn, and *.exo.

Continuous-time model: d/dt eta(t) = int + dyn*eta(t) + exo*x(t), where eta(t) is a vector of latent variables, x(t) is a vector of covariates, int, dyn, and exo are vectors and matrices specified via the arguments *.int, *.dyn, and *.exo.

The left-hand side of the dynamic model consists of a vector of latent variables for the next time point in the discrete-time case, and the vector of derivatives for the latent variables at the current time point in the continuous-time case.

For models with regime-switching dynamic functions, the user will need to provide a list of the *.int, *.dyn, and *.exo arguments. (when they are specified to take on values other than the default of zero vectors and matrices), or if a single set of vectors/matrices are provided, the same vectors/matrices are assumed to hold across regimes.

prep.matrixDynamics serves as an alternative to prep.formulaDynamics.

Value

Object of class 'dynrDynamicsMatrix'

See Also

Methods that can be used include: print, show

Examples

#Single-regime, continuous-time model. For further details run: 
#demo(RSNonlinearDiscrete, package="dynr"))
dynamics <- prep.matrixDynamics(
    values.dyn=matrix(c(0, -0.1, 1, -0.2), 2, 2),
    params.dyn=matrix(c('fixed', 'spring', 'fixed', 'friction'), 2, 2),
    isContinuousTime=TRUE)

#Two-regime, continuous-time model. For further details run: 
#demo(RSNonlinearDiscrete, package="dynr"))
dynamics <- prep.matrixDynamics(
    values.dyn=list(matrix(c(0, -0.1, 1, -0.2), 2, 2),
                    matrix(c(0, -0.1, 1, 0), 2, 2)),
    params.dyn=list(matrix(c('fixed', 'spring', 'fixed', 'friction'), 2, 2),
                    matrix(c('fixed', 'spring', 'fixed', 'fixed'), 2, 2)),
    isContinuousTime=TRUE)

Description

Prepare the measurement recipe

Usage

prep.measurement(values.load, params.load = NULL, values.exo = NULL,
  params.exo = NULL, values.int = NULL, params.int = NULL, obs.names,
  state.names, exo.names)

Arguments

Details

The values.* arguments give the starting and fixed values for their respective matrices. The params.* arguments give the free parameter labels for their respective matrices. Numbers can be used as labels. The number 0 and the character 'fixed' are reserved for fixed parameters.

When a single matrix is given to values.*, that matrix is not regime-switching. Correspondingly, when a list of length r is given, that matrix is regime-switching with values and params for the r regimes in the elements of the list.

Value

Object of class 'dynrMeasurement'

See Also

Methods that can be used include: print, printex, show

Examples

prep.measurement(diag(1, 5), diag("lambda", 5))
prep.measurement(matrix(1, 5, 5), diag(paste0("lambda_", 1:5)))
prep.measurement(diag(1, 5), diag(0, 5)) #identity measurement model

#Regime-switching measurement model where the first latent variable is
# active for regime 1, and the second latent variable is active for regime 2
# No free parameters are present.
prep.measurement(values.load=list(matrix(c(1,0), 1, 2), matrix(c(0, 1), 1, 2)))

Description

Recipe function for specifying the measurement error and process noise covariance structures

Usage

prep.noise(values.latent, params.latent, values.observed, params.observed, ...)

Arguments

Details

The arguments of this function should generally be either matrices or lists of matrices. Lists of matrices are used for regime-switching models with each list element corresponding to a regime. Thus, a list of three matrices implies a three-regime model. Single matrices are for non-regime-switching models. Some checking is done to ensure that the number of regimes implied by one part of the model matches that implied by the others. For example, the noise model (prep.noise) cannot suggest three regimes when the measurement model (prep.measurement) suggests two regimes. An exception to this rule is single-regime (i.e. non-regime-switching) components. For instance, the noise model can have three regimes even though the measurement model implies one regime. The single-regime components are simply assumed to be invariant across regimes.

Care should be taken that the parameters names for the latent covariances do not overlap with the parameters in the observed covariances. Likewise, the parameter names for the latent covariances in each regime should either be identical or completely distinct. Because the LDL' transformation is applied to the covariances, sharing a parameter across regimes may cause problems with the parameter estimation.

Use $ to show specific arguments from a dynrNoise object (see examples).

Value

Object of class 'dynrNoise'

See Also

printex to show the covariance matrices in latex.

Examples

# Two latent variables and one observed variable in a one-regime model
Noise <- prep.noise(values.latent=diag(c(0.8, 1)),
    params.latent=diag(c('fixed', "e_x")), 
    values.observed=diag(1.5,1), params.observed=diag("e_y", 1))
# For matrices that can be import to latex:
printex(Noise, show=TRUE)
# If you want to check specific arguments you've specified, for example,
# values for variance structure of the latent variables
Noise$values.latent

# Two latent variables and one observed variable in a two-regime model
Noise <- prep.noise(values.latent=list(diag(c(0.8, 1)), diag(c(0.8, 1))),
    params.latent=list(diag(c('fixed', "e_x1")), diag(c('fixed', "e_x2"))),
    values.observed=list(diag(1.5,1), diag(0.5,1)),
    params.observed=list(diag("e_y1", 1), diag("e_y2",1)))
# If the error and noise structures are assumed to be the same across regimes,
#  it is okay to use matrices instead of lists.

Description

Recipe function for creating regime switching (Markov transition) functions

Usage