Introduction

This example illustrates how to fit a linear stochastic differential equation (SDE) model with a time-varying parameter in the dynr package.

Incorporating time-varying parameters is one way to model both intraindividual change (e.g. long-term trait development) and intraindividual variability (e.g. short-term state fluctuations; Nesselroade, 1991) at the same time.

This particular model was included as an illustrative example in the chapter Stochastic Differential Equation Models with Time-Varying Parameters, included in Continuous Time Modeling in the Behavioral and Related Sciences (van Montfort et al., 2018), as well as a dynr package demo in the name of "SDETVP".

We start with a stochastic damped linear oscillator model on the latent process \(\eta_1\) with frequency parameter \(\omega\), damping/amplification force \(\zeta\), and setpoint \(\mu\): \[d\eta_{1i}(t)=\eta_{2i}(t)dt\] \[d\eta_{2i}(t)=\bigg( \omega \Big( \eta_{1i}(t)-\mu_i(t) \Big) +\zeta \eta_{2i}(t) \bigg) dt+\sigma_p dw_i(t),\] where \(d\eta_1\) and \(d\eta_2\) are the first- and second-order differentials of \(\eta_1\), and \(dw\) is the first differential of a univariate standard Wiener process.

In this example, we let the setpoint \(\mu\) follow a logistic growth curve: \[d\mu_i(t) =r\mu_i(t)\Big(1-\frac{\mu_i(t)}{k}\Big)dt,\] and include it as an auxiliary latent state.

We assume the process is measured with errors, and the measurement model is:

\[y_i(t_{i,j}) = \begin{bmatrix} 1 & 0 & 0 \end{bmatrix}\begin{bmatrix} \eta_{1i}(t_{i,j})\\ \eta_{2i}(t_{i,j})\\ \mu_i(t_{i,j}) \end{bmatrix} + \epsilon_i(t_{i,j})\]

The process looks something like this for a single individual:

Data

We first load the dynr package.

The dynr.data function creates a dynr data object. The observed variable we are modeling is "obsy" . (The simulated data are store in "LogisticSetPointSDE" and loaded as part of the package.)

require(dynr)
# Data
data(LogisticSetPointSDE)
dataLog <- dynr.data(LogisticSetPointSDE, id="id", time="times", observed="obsy")

Dynamic Model

The prep.formulaDynamics function specifies the dynamic model of the latent states, entered as a list of formulae. In the list "formulaLogistic", the left-hand side corresponds to the first-order derivative of the element. x corresponds to \(\eta_1\); y corresponds to \(\eta_2\); and z corresponds to \(\mu\) in the equations under Introduction.

The startval argument takes a named vector as the starting point of the numerical optimization algorithm for maximum likelihood estimation of the parameters. The names in the vector need to match exactly with the parameter names in the list of formulae.

Since we are fitting a continuous-time differential equation model (as opposed to a difference equation model), we indicate isContinuousTime = TRUE that enables a filtering algorithm for continuous-time latent models (continuous-discrete extended Kalman filter).

The prep.noise function specifies the dynamic (*.latent) and observation (*.observed) noise components. The params.* arguments specify whether a component is fixed ('fixed') or freely estimated (a parameter name such as 'pnoise'). The values.* arguments specify the fixed value (if the corresponding position in params.* is "fixed") or the starting value for numerical optimization (otherwise).

formulaLogistic<-list(
  x~y,
  y~eta*(x-z)+damp*y,
  z~r*z*(1-b*z)
)


dynamicsLog<-prep.formulaDynamics(formula=formulaLogistic,
                                  startval=c(eta=-.1,damp=-.05,b=0.5,r=0.3),
                                  isContinuousTime=TRUE)


ecovLog <- prep.noise(
  values.latent=diag(c(0,0.1,0), 3), params.latent=diag(c('fixed', 'pnoise','fixed'), 3), 
  values.observed=diag(0.1, 1), params.observed=diag('mnoise', 1))

Measurement Model

The prep.measurement command is used to specify the factor loadings matrix. The values.load argument specifies that the factor loadings are fixed at 1 and 0s. In this model, all parameters are fixed, which is indicated by the params.load argument. The *.names argument gives names to the latent and observed variables. The obs.names should match the name specifed as observed in dynr.data.

measLog <- prep.measurement(
  values.load=matrix(c(1, 0, 0), 1, 3), 
  params.load=matrix("fixed", 1, 3),
  state.names=c("x","y","z"),
  obs.names=c("obsy"))

Initial Values

This step specifies the covariances and latent state values of the dynamical system at t=0. These initialize the recursive filtering algorithm (continuous-discrete extended Kalman filter) that dynr uses.

initialLog <- prep.initial(
  values.inistate=c(-5, 0 ,.1),
  params.inistate=c('fixed', 'fixed','fixed'), 
  values.inicov=diag(0, 3),
  params.inicov=diag('fixed', 3))

Dynr Model

Now we put together everything we've previously specified in dynr.model. This code connects the recipes we've written up with our data and writes a c file in our working directory. We can inspect c functions that go with each recipe in the c file.

modelLog <- dynr.model(dynamics=dynamicsLog, measurement=measLog, 
                       noise=ecovLog, initial=initialLog,
                       data=dataLog, outfile="setpointLog.c")

For bounded optimization, we can set upper and lower bounds easily with the $ub and $lb attributes.

modelLog$ub[c("eta",'r', 'pnoise')] <- c(0,10,10)

Tex Options

We can check our model specifications in a neatly printed pdf file using the following code.

The printex command is used to write the model into a Latex file, with a name given by the outFile argument. Then, the tools::texi2pdf command generates a pdf file from the latex file we just created. The system command prints out the pdf file:

printex(modelLog,ParameterAs=modelLog$param.names,show=FALSE,printInit=TRUE,
        outFile="setpointLog.tex")
tools::texi2pdf("setpointLog.tex")
system(paste(getOption("pdfviewer"), "setpointLog.pdf"))

We can also pull up the model specifications in the R plotting region, instead of generating a Latex file, using the command plotFormula.

plotFormula(modelLog,ParameterAs=modelLog$param.names)

Optimization Step and Results

Now the model is specified and we are ready to cook dynr! The summary command gives the model fitting results.

resLog <- dynr.cook(modelLog, verbose=FALSE)
summary(resLog)
## Coefficients:
##          Estimate Std. Error t value   ci.lower   ci.upper Pr(>|t|)    
## eta    -0.9966760  0.0032812 -303.75 -1.0031071 -0.9902449   <2e-16 ***
## damp   -0.1055465  0.0032592  -32.38 -0.1119344 -0.0991585   <2e-16 ***
## b       0.0999665  0.0002888  346.14  0.0994005  0.1005326   <2e-16 ***
## r       0.5011447  0.0021457  233.55  0.4969391  0.5053503   <2e-16 ***
## pnoise  0.0072619  0.0036314    2.00  0.0001446  0.0143793   0.0228 *  
## mnoise  0.9628490  0.0284114   33.89  0.9071637  1.0185343   <2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## -2 log-likelihood value at convergence = 6790.45
## AIC = 6802.45
## BIC = 6837.17

We can inspect the model again, with estimated parameter values instead of parameter names.

plotFormula(modelLog,ParameterAs=coef(resLog))

(Questions regarding this tutorial can be directed to Meng Chen: mengc2013@gmail.com)

References

Chen M., Chow S.-M., Hunter M.D. (2018) Stochastic Differential Equation Models with Time-Varying Parameters. In: van Montfort K., Oud J.H.L., Voelkle M.C. (eds) Continuous Time Modeling in the Behavioral and Related Sciences. Springer, Cham. https://doi.org/10.1007/978-3-319-77219-6_9

Nesselroade, J. R. (1991). The warp and the woof of the developmental fabric. In R. M. Downs, L. S. Liben, & D. S. Palermo (Eds.), Visions of aesthetics, the environment & development: The legacy of Joachim F. Wohlwill (p. 213–240). Lawrence Erlbaum Associates, Inc.