mxExpectationStateSpaceContinuousTime: Create an MxExpectationStateSpace Object

Description

This function creates a new MxExpectationStateSpace object.

Usage

mxExpectationStateSpaceContinuousTime(A, B, C, D, Q, R, x0, P0, u, t = NA,
                        dimnames = NA, thresholds = NA, threshnames = dimnames,
                        ..., scores=FALSE)
mxExpectationSSCT(A, B, C, D, Q, R, x0, P0, u, t = NA,
                        dimnames = NA, thresholds = NA, threshnames = dimnames,
                        ..., scores=FALSE)

Arguments

A character string indicating the name of the 'A' matrix.

A character string indicating the name of the 'B' matrix.

A character string indicating the name of the 'C' matrix.

A character string indicating the name of the 'D' matrix.

A character string indicating the name of the 'Q' matrix.

A character string indicating the name of the 'R' matrix.

A character string indicating the name of the 'x0' matrix.

A character string indicating the name of the 'P0' matrix.

A character string indicating the name of the 'u' matrix.

A character string indicating the name of the 't' matrix.

dimnames

An optional character vector to be assigned to the row names of the 'C' matrix.

thresholds

Not Yet Implemented. An optional character string indicating the name of the thresholds matrix.

threshnames

Not Yet Implemented. An optional character vector to be assigned to the column names of the thresholds matrix.

...

Unused. Requires further arguments to be named.

scores

Not to be used

Value

Returns a new MxExpectationStateSpace object. mxExpectationStateSpace objects should be included with models with referenced MxAlgebra, MxData and MxMatrix objects.

Details

The mxExpectationStateSpaceContinuousTime and mxExpectationSSCT functions are identical. The latter is simply an abbreviated name. When using the former, tab completion is strongly encouraged to save tedious typing. Both of these functions are wrappers for the mxExpectationStateSpace function, which could be used for both discrete and continuous time modeling. However, there is a strong possibility of misunderstanding the model parameters when switching between discrete time and continuous time. The expectation matrices have the same names, but mean importantly different things so caution is warranted. The best practice is to use mxExpectationStateSpace for discrete time models, and mxExpectationStateSpaceContinuousTime for continuous time models. Expectation functions define the way that model expectations are calculated. That is to say, expectation functions define how a set of model matrices get turned into expectations for the data. When used in conjunction with the mxFitFunctionML, the mxExpectationStateSpace uses maximum likelihood prediction error decomposition (PED) to obtain estimates of free parameters in a model of the raw MxData object. Continuous time state space expectations treat the raw data as a multivariate time series of possibly unevenly spaced times with each row corresponding to a single occasion. Continuous time state space expectations implement a hybrid Kalman filter to produce expectations. The hybrid Kalman filter uses a Kalman-Bucy filter for the prediction step and the classical Kalman filter for the update step. It is a hybrid between the classical Kalman filter used for the discrete (but possibly unequally spaced) measurement occastions and the continous time Kalman-Bucy filter for latent variable predictions. Missing data handling is implemented in the same fashion as full information maximum likelihood for partially missing rows of data. Additionally, completely missing rows of data are handled by only using the prediction step from the Kalman-Bucy filter and omitting the update step. This model uses notation for the model matrices commonly found in engineering and control theory. The 'A', 'B', 'C', 'D', 'Q', 'R', 'x0', and 'P0' arguments must be the names of MxMatrix or MxAlgebraobjects with the associated properties of the A, B, C, D, Q, R, x0, and P0 matrices in the state space modeling approach. The 't' matrix must be a 1x1 matrix using definition variables that gives the times at which measurements occurred. The state space expectation is defined by the following model equations. $$\frac{d}{dt} x(t) = A x(t) + B u_t + q(t)$$ $$y_t = C x_t + D u_t + r_t$$ with $q(t)$ and $r_t$ both independently and identically distributed random Gaussian (normal) variables with mean zero and covariance matrices $Q$ and $R$, respectively. Subscripts or square brackets indicate discrete indices; parentheses indicate continuous indices. The derivative of $x(t)$ with respect to $t$ is $\frac{d}{dt} x(t)$. The first equation is called the state equation. It describes how the latent states change over time with a first-order linear differential equation. Unlike some other programs, we do not require that the continuous time $A$ matrix has an inverse. This allows zero dynamics (i.e. no growth models) and many other important kinds of processes. The second equation is called the output equation. It describes how the latent states relate to the observed states at a single point in time. The output equation shows how the observed output is produced by the latent states. Also, the output equation in state space modeling is directly analogous to the measurement model in LISREL structural equation modeling. Note that the covariates, $u$, have "instantaneous" effects on both the state and output equations. If lagged effects are desired, then the user must create a lagged covariate by shifting their observed variable to the desired lag. The state and output equations, together with some minimal assumptions and the Kalman filter, imply a new expected covariance matrix and means vector for every row of data. The expected covariance matrix of row $t$ is $$S_{t} = C ( A P_{t-1} A^{\sf T} + Q ) C^{\sf T} + R$$ The expected means vector of row $t$ is $$\hat{y}_{t} = C x_{t} + D u_{t}$$ The 'dimnames' arguments takes an optional character vector. The 'A' argument refers to the $A$ matrix in the State Space approach. This matrix gives the dynamics. Entries in the diagonal give the strength of the influence of a variable's position on its slope. Entries in the off-diagonal give the coupling strength from one variable to another. The $A$ matrix is sometimes called the state-transition model. The 'B' argument refers to the $B$ matrix in the State Space approach. This matrix consists of exogenous forces that influence the dynamics. Note that the covariate effect is contemporaneous: the covariate at time $t$ has influence on the slope of the latent state also at time $t$. A lagged effect can be created by lagged the observed variable. The $B$ matrix is sometimes called the control-input model. The 'C' argument refers to the $C$ matrix in the State Space approach. This matrix consists of contemporaneous regression coefficients from the latent variable in column $j$ to the observed variable in row $i$. This matrix is directly analogous to the factor loadings matrix in LISREL and Mplus models. The $C$ matrix is sometimes called the observation model. The 'D' argument refers to the $D$ matrix in the State Space approach. This matrix consists of contemporaneous regressive coefficients from the input (manifest covariate) variable $j$ to the observed variable in row $i$. The $D$ matrix is sometimes called the feedthrough or feedforward matrix. The 'Q' argument refers to the $Q$ matrix in the State Space approach. This matrix gives the covariance of the dynamic noise. The dynamic noise can be thought of as unmeasured covariate inputs active at all times. This matrix must be symmetric, diagonal, or zero. As a special case, it is often diagonal. The $Q$ matrix is the covariance of the process noise. Just as in factor analysis and general structural equation modeling, the scale of the latent variables is usually set by fixing some factor loadings in the $C$ matrix, or fixing some factor variances in the $Q$ matrix. The 'R' argument refers to the $R$ matrix in the State Space approach. This matrix consists of residual covariances among the observed (manifest) variables. This matrix must be symmetric As a special case, it is often diagonal. The $R$ matrix is the covariance of the observation noise. The 'x0' argument refers to the $x_0$ matrix in the State Space approach. This matrix consists of the column vector of the initial values for the latent variables. The state space expectation uses the $x_0$ matrix as the starting point to recursively estimate the latent variables' values at each time. These starting values can be difficult to pick, however, for sufficiently long time series they often do not greatly impact the estimation. The 'P0' argument refers to the $P_0$ matrix in the State Space approach. This matrix consists of the initial values of the covariances of the error in the initial latent variable estimates given in $x_0$. That is, the $P_0$ matrix gives the covariance of $x_0 - xtrue_0$ where $xtrue_0$ is the vector of true initial values. $P_0$ is a measure of the accuracy of the intial latent state estimates. The Kalman filter uses this initial covariance to recursively generated a new covariance for each time point based on the previous time point. The Kalman filter updates this covariance so that it is as small as possible (minimum trace). Similar to the $x_0$ matrix, these starting values are often difficult to choose. The 'u' argument refers to the $u$ matrix in the State Space approach. This matrix consists of the inputs or manifest covariates of the state space expectation. The $u$ matrix must be a column vector with the same number of rows as the $B$ and $D$ matrices have columns. If no inputs are desired, $u$ can be a zero matrix. If time-varying inputs are desired, then they should be included as columns in the MxData object and referred to in the labels of the $u$ matrix as definition variables. There is an example of this below. The 't' argument refers to the $t$ matrix in the State Space approach. This matrix should be 1x1 (1 row and 1 column) and not free. The label for the element of this matrix should be 'data.YourTimeVariable'. The 'data' part does not change, but 'YourTimeVariable' should be a name in your data set that gives the times at which measurement happened. The units of time are up to you. Your choice of time units will influence of the values of the parameters you estimate. Also, recall that the model is given $x_0$ and $P_0$. These always happen at $t=0$. So the first row of data happens some amount of time after zero. The MxMatrix objects included as arguments may be of any type, but should have the properties described above. The mxExpectationStateSpace will not return an error for incorrect specification, but incorrect specification will likely lead to estimation problems or errors in the mxRun function. mxExpectationStateSpaceContinuousTime evaluates with respect to an MxData object. The MxData object need not be referenced in the mxExpectationStateSpace function, but must be included in the MxModel object. mxExpectationStateSpace requires that the 'type' argument in the associated MxData object be equal to 'raw'. Neighboring rows of the MxData object are treated as adjacent, equidistant time points increasing from the first to the last row. To evaluate, place an mxExpectationStateSpaceContinuousTime object, the mxData object for which the expected covariance approximates, referenced MxAlgebra and MxMatrix objects, and optional MxBounds and MxConstraint objects in an MxModel object. This model may then be evaluated using the mxRun function. The results of the optimization can be found in the 'output' slot of the resulting model, and may be obtained using the mxEval function.

References

K.J. and R.M. Murray (2010). Feedback Systems: An Introduction for Scientists and Engineers. Princeton University Press. J. Durbin and S.J. Koopman. (2001). Time Series Analysis by State Space Methods. Oxford University Press. R.E. Kalman (1960). A New Approach to Linear Filtering and Prediction Problems. Basic Engineering, 82, 35-45. R.E. Kalman and R.S. Bucy (1961). New Results in Linear Filtering and Prediction Theory. Transactions of the ASME, Series D, Journal of Basic Engineering, 83, 95-108. G. Petris (2010). An R Package for Dynamic Linear Models. Journal of Statistical Software, 36, 1-16. The OpenMx User's guide can be found at http://openmx.ssri.psu.edu/documentation.

Examples

Run this code

#------------------------------------------------------------------------------
# Example 1
# Undamped linear oscillator, i.e. a noisy sine wave.
#	Measurement error, but no dynamic error, single indicator.
#	This example works great.

#--------------------------------------
# Data Generation

require(OpenMx)

set.seed(405)
tlen <- 200
t <- seq(1.2, 50, length.out=tlen)

freqParam <- .5
initialCond <- matrix(c(2.5, 0))
x <- initialCond[1,1]*cos(freqParam*t)
plot(t, x, type='l')

measVar <- 1.5
y <- cbind(obs=x+rnorm(tlen, sd=sqrt(measVar)), tim=t)

plot(t, y[,1], type='l')

#--------------------------------------
# Model Specification

#Note: the bounds are here only to keep SLSQP from
# stepping too far off a cliff.  With the bounds in
# place, SLSQP finds the right solution.  Without
# the bounds, SLSQP goes crazy.


cdim <- list('obs', c('ksi', 'ksiDot'))

amat <- mxMatrix('Full', 2, 2, c(FALSE, TRUE, FALSE, TRUE), c(0, -.1, 1, -.2),
	name='A', lbound=-10)
bmat <- mxMatrix('Zero', 2, 1, name='B')
cmat <- mxMatrix('Full', 1, 2, FALSE, c(1, 0), name='C', dimnames=cdim)
dmat <- mxMatrix('Zero', 1, 1, name='D')
qmat <- mxMatrix('Zero', 2, 2, name='Q')
rmat <- mxMatrix('Diag', 1, 1, TRUE, .4, name='R', lbound=1e-6)
xmat <- mxMatrix('Full', 2, 1, TRUE, c(0, 0), name='x0', lbound=-10, ubound=10)
pmat <- mxMatrix('Diag', 2, 2, FALSE, 1, name='P0')
umat <- mxMatrix('Zero', 1, 1, name='u')
tmat <- mxMatrix('Full', 1, 1, name='time', labels='data.tim')

osc <- mxModel("LinearOscillator", 
	amat, bmat, cmat, dmat, qmat, rmat, xmat, pmat, umat, tmat,
	mxExpectationSSCT('A', 'B', 'C', 'D', 'Q', 'R', 'x0', 'P0', 'u', 'time'),
	mxFitFunctionML(),
	mxData(y, 'raw'))


oscr <- mxRun(osc)


#--------------------------------------
# Results Examination

summary(oscr)


(ssFreqParam <- mxEval(sqrt(-A[2,1]), oscr))
freqParam


(ssMeasVar <- mxEval(R, oscr))
measVar


dampingParam <- 0
(ssDampingParam <- mxEval(-A[2,2], oscr))
dampingParam

Run the code above in your browser using DataLab