%\VignetteIndexEntry{Linear discrete-time regime-switching models}
\documentclass{article}
%\VignetteKeyword{linear dynamic models}
\usepackage[utf8]{inputenc}
\usepackage[T1]{fontenc}

\usepackage{amsmath}%align
\usepackage{arydshln}%dashline in array
\usepackage{bm}%bm

\newcommand\pkg[1]{{\bf #1}}
\newcommand{\CRANpkg}[1]{\href{https://CRAN.R-project.org/package=#1}{\pkg{#1}}}%
\newcommand\proglang[1]{{\bf #1}}
\newcommand\code[1]{{\textit{#1}}}

\newcommand\refequ[1]{Equation~\ref{#1}}
\newcommand\reffig[1]{Figure~\ref{#1}}
\newcommand\reftb[1]{Table~\ref{#1}}

\newcommand{\nbf} [1] {\ensuremath{\boldsymbol{#1}}}
\newcommand{\nbff} [1] {\ensuremath{\boldsymbol{\mathrm{#1}}}}
\newcommand{\x}{\nbf{x}_i}
\newcommand{\bb}{{\bf b}}
\newcommand{\bd}{\bf d}
\newcommand{\btau}{\bm{\tau}}
\newcommand{\bLambda}{{\bm{\Lambda}}}
\newcommand{\bepsilon}{\mbox{\boldmath $\epsilon$}}
\newcommand{\bR}{{\bf R}}
\newcommand{\bQ}{{\bf Q}}
\newcommand{\St}{S_i(t_{i,j})}
\newcommand{\Stb}{S_i(t_{i,j-1})}
\newcommand{\definedas}{\stackrel{\Delta}{=}}

\usepackage[american]{babel}
\usepackage{csquotes}
\usepackage[backend=bibtex,sorting=none]{biblatex}
\DeclareLanguageMapping{american}{american-apa}
\addbibresource{ou-hunter-chow.bib}


\title{Example: Regime-switching Linear Discrete-time Model}
\author{Lu Ou, Michael D. Hunter, Sy-Miin Chow}
%\date{\today}

\begin{document}
%\SweaveOpts{concordance=TRUE}

\maketitle

This file demonstrates the utilization of \pkg{dynr} in fitting a regime-switching linear dynamic models. A complete modeling script for this example is available as one of the demo examples in \pkg{dynr} and can be found using 
<<include = TRUE, cache = FALSE, echo = TRUE,eval=FALSE,size='footnotesize'>>=
file.edit(system.file("demo", "RSLinearDiscreteYang.R", package = "dynr"))
@

\section{Regime-switching State Space Model}
Facial electromyography (EMG) has been used in the behavioral sciences as one possible indicator of
human emotions \cite[e.g.,][]{Schwartz75a,Cacioppo81a,Cacioppo86a,Dimberg00a}. When human subjects are exposed to emotion induction procedures, researchers have detected changes in individuals' facial EMG recordings even when the corresponding changes in facial expression are too subtle to be detected by human raters \cite{Schwartz75a,Dimberg90a,Cacioppo81a}.

A time series of EMG data contains bursts of electrical activity that are typically magnified when an individual is under emotion induction. \cite{Yang10a} proposed using a regime-switching linear state-space model in which the individual may transition between regimes with and without facial EMG activation. As such, heterogeneities in the dynamic patterns and variance of EMG data are also accounted for through the incorporation of these latent regimes. Model fitting was previously performed at the individual level. Data from the participant shown in Figure \ref{EMGdata}(A) are made available as part of the demonstrative examples in \pkg{dynr}.  

The model of interest is the final model selected for this participant by \cite{Yang10a}:
\begin{eqnarray}
&&{y}_i(t_{i,j}) = \mu_{y\St} + \beta_{\St} \text{Self-report}(t_{i,j}) + \eta_i (t_{i,j}), \label{eq:emgMeas}\\
&&\eta_i (t_{i,j+1}) = \phi_{\St} \eta_i (t_{i,j}) + \zeta_i (t_{i,j+1}) , \label{eq:emgDyn}
\end{eqnarray}\label{SSTranARMA}in which we allowed  the intercept, $\mu_{y\St}$, the regression slope, $\beta_{\St}$, and the autoregression coefficient, $\phi_{\St}$, to be regime-dependent. By allowing $\phi_{\St}$ to be regime-specific, we indirectly allowed the total variance of the latent component, $\eta_i (t_{i,j+1})$, to be heterogeneous across the deactivation and activation stages, in spite of requiring the dynamic noise variance, $E(\zeta_i(t)^2)$, to be constant across regimes.

\begin{figure}[htbp]
\begin{center}
\includegraphics[width=.45\textwidth]{./Figures/NewplotEMG1Subj.pdf}
\includegraphics[width=.45\textwidth]{./Figures/plotRSGG.pdf}
\end{center}
\caption{(A) A plot of integrated electromyography (iEMG) and self--report affect ratings for one participant with a time interval of 0.2 seconds between two adjacent observations. Self--report = self--report affect ratings; iEMG = integrated EMG signals. (B) An automatic plot of the smoothed state estimates for the regime-switching linear state-space model.}
\label{EMGdata}
\end{figure}
\subsection{Prepare the data}
The first step in \pkg{dynr} modeling is to structure the data. This is done with the \code{dynr.data()} function.
<<include = TRUE, cache = FALSE, echo = TRUE,eval=FALSE,size='footnotesize'>>=
#---- Load packages ----
require("dynr")
#---- Read in data and create dynr data object----
data("EMG")
EMGdata <- dynr.data(EMG, id = 'id', time = 'time', 
  observed = 'iEMG', covariates = 'SelfReport')
@

The first argument of this function is either a \code{ts} class object of single-subject time series or a \code{data.frame} structured in a long (relational) format (i.e., with different measurement occasions from the same subject appearing as different rows in the data frame). Missing values in the observed variables should be indicated by \code{NA}. When a \code{ts} class object is passed to \code{dynr.data()}, no other inputs are needed. Otherwise, the \code{id} argument needs the name of the ID variable as input, and allows multiple people to be estimated in a single model by distinguishing different individuals with the ID variable. That is, it indicates which rows should be modeled together as a time series. Thus, multi-subject modeling is as easy as single-subject modeling; only the data differ. The \code{time} argument needs the name of the TIME variable that indicates subject-specific measurement occasions. If a discrete-time model is desired, the TIME variable should contain subject-specific sequences of (subsets of) consecutively equally spaced numbers (e.g, 1, 2, 3, $\cdots$). In other words, the program assumes that the input \code{data.frame} is equally spaced with potential missingness. If the measurement occasions for a subject are a subset of an arithmetic sequence but are not consecutive, \code{NA}s will be inserted automatically to create an equally spaced data set before estimation. If a continuous-time model is being specified, the TIME variable can contain subject-specific increasing sequences of irregularly spaced real numbers. That is, the data may be input at their original, irregularly spaced intervals without the need to insert missingness. In this particular example, a discrete time model is used. 

The \code{observed} and \code{covariates} arguments are used to indicate the names of the observed variables and covariates in the data. Covariates are defined as fixed predictors that are hypothesized to affect the modeling functions in one or more ways, but are otherwise not of interest (i.e., not modeled as dependent variables) to the user. Missing values in covariates are not allowed. That is, missing values in the covariates, if there are any, should be imputed first. The \code{dynr.data()} function lets users include data sets with many variables, but only use a few. The output of the function combines with the model recipe information later to map the model onto the data. 

\subsection{Prepare the recipes}
The next step in \pkg{dynr} modeling is to build the recipes for the various parts of a model. The recipes are created with \code{prep.*()} functions. 

\subsubsection{Model specification: the dynamic functions}
The dynamic model can take on the form of continuous-time models as 
\begin{equation}
d \nbf{\eta}_i (t) = \nbf{f}_{S_i(t)}\left(\nbf{\eta}_i(t), t, \x(t) \right)dt + d\nbf{w}_i(t),\label{RS-SDE}
\end{equation} or the form of discrete-time state-space models \cite{Durbin01a} as 
\begin{equation}
\nbf{\eta_i} (t_{i,j+1}) = \nbf{f}_{S_i(t)}\left(\nbf{\eta_i} (t_{i,j}), t_{i,j}, \x(t_{i,j}) \right) + \nbf{w}_i(t_{i,j+1}), \label{RS-dyn}
\end{equation} where $i$ indexes person, $t$ indexes time, $\nbf{\eta}_i(t)$ is the $r\times 1$ vector of latent variables at time $t$, $\nbf{x}_i(t)$ is the vector of covariates at time $t$, and $\nbf{f}_{S_i(t)}(.)$ is the vector of (possibly nonlinear) dynamic functions. $\nbf{f}_{S_i(t)}(.)$ depends on the latent regime indicator, $S_i(t)$, the discrete-valued latent variable that indexes the operating regime at time $t$.

The dynamic functions, $\nbf{f}_{S_i(t)}()$ in Equations \ref{RS-SDE} and \ref{RS-dyn}, can be specified using one of two possible functions in \pkg{dynr}: \code{prep.formulaDynamics()} and \code{prep.matrixDynamics()}. The dynamic model in this particular example consists only of linear functions, although the parameters that appear in these linear functions are regime-dependent. In this special case, the dynamic model in Equation \ref{RS-dyn} reduces to: 
\begin{equation}
\nbf{\eta_i} (t_{i,j+1}) = \nbf{\alpha}_{S_i(t_{i,j})} + \nbf{F}_{S_i(t_{i,j})} \nbf{\eta_i} (t_{i,j}) + \nbf{B}_{S_i(t_{i,j})} \x(t_{i,j}) + \nbf{w}_i(t_{i,j+1}),\label{eq:linDyn}
\end{equation}
where the general, possibly nonlinear function $\nbf{f}_{S_i(t)}()$ is replaced with a linear function consisting of (1) an intercept term $\nbf{\alpha}_{S_i(t_{i,j})}$, (2) linear dynamics instantiated as an $r \times r$ matrix $\nbf{F}_{S_i(t_{i,j})}$, (3) linear covariate regression effects $\nbf{B}_{S_i(t_{i,j})}$, and the same additive noise term $\nbf{w}_i(t_{i,j+1})$. As indicated by the subscript $S_i(t_{i,j})$, all of these can also be regime-dependent.  Of course, the same structure is possible in continuous time as the linear analog of Equation \ref{RS-SDE}.
\begin{equation}
d \nbf{\eta_i} (t) = \left( \nbf{\alpha}_{S_i(t)} + \nbf{F}_{S_i(t)} \nbf{\eta_i} (t) + \nbf{B}_{S_i(t)} \x(t) \right) dt + d \nbf{w}_i(t),\label{eq:linDynC}
\end{equation}

In this example, the dynamics as in Equation \ref{eq:emgDyn} are linear and discrete-time, so we can describe the dynamics in terms of Equation \ref{eq:linDyn} as
\begin{equation}
\label{eq:ex1dyn}
\eta_i (t_{i,j+1}) = \underbrace{0}_{\nbf{\alpha}_{S_i(t_{i,j})}} + \underbrace{\phi_{\St}}_{\nbf{F}_{S_i(t_{i,j})}} \eta_i (t_{i,j}) + \underbrace{0}_{\nbf{B}_{S_i(t_{i,j})}} \x(t_{i,j}) + \underbrace{\zeta_i (t_{i,j+1})}_{\nbf{w}_i(t_{i,j+1})}. 
\end{equation}
The \code{prep.matrixDynamics()} function allows the user to specify the structures of the intercept vector $\nbf{\alpha}_{S_i(t_{i,j})}$, through \code{values.int} and \code{params.int}, the covariate regression matrix $\nbf{B}_{S_i(t_{i,j})}$, through \code{values.exo} and \code{params.exo}, and the one-step-ahead transition matrix $\nbf{F}_{S_i(t_{i,j})}$, through \code{values.dyn} and \code{params.dyn}, in the linear special case for those who prefer to work in such a matrix algebraic framework. We illustrate this function in the current example below. The \code{values.dyn} argument gives a list of matrices for the starting values of $\nbf{F}_{S_i(t_{i,j})}$. The \code{params.dyn} argument names the free parameters. These are the $\phi_{S_t}$ in Equation \ref{eq:emgDyn}. The \code{isContinuousTime} argument switches between continuous-time modeling (when true) and discrete-time modeling (when false). Because this argument is false, the dynamics are in a discrete-time form that matches Equation \ref{eq:emgDyn}. The arguments corresponding to the intercepts (\code{values.int} and \code{params.int}) and the covariate effects (\code{values.exo} and \code{params.exo}) are omitted to leave these matrices as zeros. 
<<include = TRUE, cache = FALSE, echo = TRUE,eval=FALSE,size='footnotesize'>>=
#---- Dynamic functions ----
recDyn <- prep.matrixDynamics(
  values.dyn = list(matrix(.1, 1, 1), matrix(.5, 1, 1)),
  params.dyn = list(matrix('phi_1', 1, 1), matrix('phi_2', 1, 1)),
  isContinuousTime = FALSE)
@

\subsubsection{Model specification: the linear measurement function}
For both discrete- and continuous-time models, we assume that we have a discrete-time measurement model in which $\nbf{\eta}_i(t_{i,j})$ at discrete time point $t_{i,j}$ is indicated by a $p$ $\times$ 1 vector of manifest observations, $\nbf{y}_i(t_{i,j})$ as
\begin{eqnarray}
\nbf{y}_i(t_{i,j})=\btau_{\St}+ \bLambda_{\St} \nbf{\eta}_i(t_{i,j})+ \nbff{A}_{\St} \nbf{x}_i(t_{i,j}) + \nbf{\epsilon}_i(t_{i,j}),
\text{\hskip .1in} \bepsilon_i(t_{i,j})
\sim N\left(\mathbf{0},
             \bR_{\St}\right), \label{SSMeaSt}
\end{eqnarray}
where $\btau_{\St}$ is a $p \times 1$ vector of intercepts, $\nbff{A}_{\St}$ is a matrix of regression weights for the covariates observed at time $t_{i,j}$, $\bLambda_{\St}$ is a $p \times r$ factor loadings matrix that links the observed variables to the latent variables, and $\bepsilon_i(t_{i,j})$ is a $p \times 1$ vector of measurement errors assumed to be serially uncorrelated over time and normally distributed with zero means and (possibly) regime-specific covariance matrix, $\bR_{\St}$.

<<include = TRUE, cache = FALSE, echo = TRUE,eval=FALSE,size='footnotesize'>>=
#---- Measurement ----
recMeas <- prep.measurement(
  values.load = rep(list(matrix(1, 1, 1)), 2),
  values.int = list(matrix(4, 1, 1), matrix(3, 1, 1)),
  params.int = list(matrix('mu_1', 1, 1), matrix('mu_2', 1, 1)),
  values.exo = list(matrix(0, 1, 1), matrix(1, 1, 1)),
  params.exo = list(matrix('fixed', 1, 1), matrix('beta_2', 1, 1)),
  obs.names = c('iEMG'),
  state.names = c('eta'),
  exo.names = c("SelfReport"))
@


\subsubsection{Model specification: the latent and observed noise covariance structures}
The noise recipe is created with \code{prep.noise()}. $\nbf{w}_i(t)$ in \refequ{RS-SDE} is an $r$-dimensional Wiener process (i.e., continuous-time analog of a random walk process). The differentials of the Wiener processes have zero means and regime-specific covariance matrix, $\bQ_{S_i(t)}$, often called the {\em diffusion} matrix. In \refequ{RS-dyn}, however, $\nbf{w}_i(t)$ denotes a vector of Gaussian distributed process noise with regime-specific covariance matrix, $\bQ_{S_i(t)}$. In both continuous- and discrete-time models, $\bQ_{S_i(t)}$ can be specified by the \code{*.latent} arguments in \code{prep.noise()}. The \code{*.observed} arguments are for $\bR_{\St}$ in \refequ{SSMeaSt}.

The code below creates the noise recipe by calling the \code{prep.noise()} function. The noise recipe is stored in the \code{recNoise} object, an abbreviation for ``recipe noise''. The latent noise covariance matrix is a $1\times 1$ matrix with a free parameter called \code{dynNoise}, short for ``dynamic noise.'' The observed noise covariance matrix is also a $1 \times 1$ matrix, but has the measurement noise variance fixed to zero. These covariance matrices need to be positive definite. The zero's in the diagonal of a covariance matrix are internally replaced by a small positive number automatically before estimation. To ensure the matrix stays positive definite in estimation, we will apply a set of transformations to the matrix in each iteration of the optimization, so the starting or fixed values of these matrices are automatically adjusted for this purpose.
<<include = TRUE, cache = FALSE, echo = TRUE,eval=FALSE,size='footnotesize'>>=
#---- Dynamic and measurement noise cov structures----
recNoise <- prep.noise(
  values.latent = matrix(1, 1, 1),
  params.latent = matrix('dynNoise', 1, 1),
  values.observed = matrix(0, 1, 1),
  params.observed = matrix('fixed', 1, 1))
@

\subsubsection{Model specification: the initial condition}
In both the discrete- and continuous-time cases, the initial conditions for the dynamic functions are defined explicitly to be the latent variables at an individual-specific initial time point, $t_{i,1}$ (i.e., the first observed time point), denoted as $\nbf{\eta}_i(t_{i,1})$, and are specified to be normally distributed with means $\nbf{\mu_{\eta_1}}$ and covariance matrix, $\nbf{\Sigma_{\eta_1}}$:
  \begin{equation}
\nbf{\eta}_i(t_{i,1})\sim N\left(\nbf{\mu_{\eta_1}}, \nbf{\Sigma_{\eta_1}}\right). \label{initial}
\end{equation} 

The subscript $S_i(t)$ that appears in Equations \ref{RS-SDE}--\ref{SSMeaSt} indicates that the values of the parameters in these functions and matrices may depend on $S_i(t)$, the operating regime for individual $i$ at time point, $t$. Often in practice, only some of these elements are freed to vary by regime. To make inferences on $\St$, it is essential to specify a model or mechanism through which $\St$ changes over individuals and time. Just as with the continuous latent variables in $\nbf{\eta}_i(t)$, we initialize the categorical latent variable $\St$ on the first occasion and then provide a model for how $\St$ changes over time.  The initial class (or regime) probabilities for $S_i(t_{i,1})$ are represented using a multinomial regression model as
\begin{eqnarray}
\Pr\big(S_i(t_{i,1}) = m |\nbf{x}_i(t_{i,1})\big) &\definedas \pi_{m,i1} = \frac{\exp(a_{m} + \bb^T_{m}\nbf{x}_i(t_{i,1}))}{\sum_{k=1}^{M} \exp(a_{k} +\bb^T_{k} \nbf{x}_i(t_{i,1}))}, \label{MultinomialReg0}
\end{eqnarray}
where $M$ denotes the total number of regimes, $a_{m}$ denotes the logit intercept for the $m$th regime and $\bb_{m}$ is a $n_b \times 1$ vector of regression slopes linked to a vector of covariates used to explain possible interindividual differences in initial log-odds (LO) of being in a regime relative to the reference regime selected by the user, operationalized as the regime where $a_{m}$ and all entries in $\bb_{m}$ are set to zero. Setting these entries to be zero in at least the reference regime is necessary for identification purposes: this ensures that the initial regime probabilities across all the hypothesized regimes sum to 1.0. In the simplest case without covariates, Equation \ref{MultinomialReg0} reduces to a specification of $M$ initial regime prevalence parameters - either on a LO (ranging from $-\infty$ to $+\infty$) or a probability (ranging between 0 and 1) scale, as preferred by the user.

The \code{prep.initial()} function is used to specify the $\nbf{\mu_{\eta_1}}$ and $\nbf{\Sigma_{\eta_1}}$ in \refequ{initial}, as well as the model for the initial regime probabilities (i.e., \refequ{MultinomialReg0}). 

<<include = TRUE, cache = FALSE, echo = TRUE,eval=FALSE,size='footnotesize'>>=
#---- Initial condition specification ----
recIni <- prep.initial(
  values.inistate = matrix(0, 1, 1),
  params.inistate = matrix('fixed', 1, 1),
  values.inicov = matrix(1, 1, 1),
  params.inicov = matrix('fixed', 1, 1),
  values.regimep = c(1, 0),
  params.regimep = c('fixed', 'fixed'))
@

\subsubsection{Model specification: the regime-switching model}
With the initial class probabilities specified, it remains to create a model for how the classes change over time. A simple strategy for this is to create a first-order Markov model for the categorical latent variable, $S_i(t_{i,j})$, $j=2, \ldots, T$, for the remaining time span. Such a model assumes that the probability of entering the current regime depends only on the previous regime. All possible transitions from one regime to another can be arranged into a matrix of transition probabilities, in which the rows index the previous regime at time $t_{i,j-1}$ and the columns index the regime to which the system transitions at time $t_{i,j}$. The rows of this matrix sum to 1.0 because the probability of transitioning from a particular state to any other state must be 1.0.
Hence, we use a first-order Markov process to define how the classes change over time in a transition probability matrix. This transition matrix may also depend on covariates. Thus, a multinomial logistic regression equation is assumed to govern the probabilities of transitions between regimes as:
  \begin{eqnarray}
\Pr\big(\St = m | \Stb=l, \nbf{x}_i(t_{i,j})\big) &\definedas \pi_{lm,it} = \frac{\exp(c_{lm} +{\bd}^T_{lm}\nbf{x}_i(t_{i,j})}{\sum_{k
  =1}^{M} \exp(c_{lk} +{\bd}^T_{lk}\nbf{x}_i(t_{i,j}))}, \label{MultinomialRegNoDev}
\end{eqnarray}
where $\pi_{lm,it}$ denotes individual $i$'s probability of transitioning from class $l$ at time $t_{i,j-1}$ to class $m$ at time $t_{i,j}$ (i.e., the entry in the $l$th row and $m$th column of the transition probability matrix), $c_{lm}$ denotes the logit intercept for the transition probability, and ${\bd}_{lm}$ is a $n_d\times 1$ vector of logit slopes summarizing the effects of the covariates in $\nbf{x}_i(t_{i,j})$ on that transition probability. The coefficients in ${\bd}_{lm}$ are LO parameters representing the effects of the covariates on the LO of transitioning from the $l$th regime into the $m$th regime relative to transitioning into the reference regime - namely, the regime in which all LO parameters (including $c_{lM}$ and all elements in ${\bd}^T_{lM}$) are set to 0. One regime, again, has to be specified as the reference regime for identification purposes to ensure that conditional on being in a particular regime at time $t_{i,j-1}$, the probabilities of transitioning to each of the $M$ regimes sum to 1.0 (i.e.,  $\sum^M_{m=1}\pi_{lm}=1$).

The \code{prep.regimes()} function specifies the structure of the regime switching functions shown in Equation \ref{MultinomialRegNoDev}. Note that based on Equation \ref{MultinomialRegNoDev}, a total of $n_d+1$ parameters, including an intercept, $c_{lm}$, and $n_d$ regression slopes in ${\bd}_{lm}$, have to be defined for each of the functions governing the transition from the $l$th regime ($l = 1, \ldots, M$) to the $m$th regime ($m = 1,\ldots, M$). In total, there are $M \times M$ of such transition functions, corresponding to entries in an $M \times M$ transition probability matrix. The function \code{prep.regimes()} requires the user to provide the starting values (through the \code{values} argument) and names (through the \code{params} argument) for  these $M\times(n_d+1)$ parameters as a matrix whose number of rows equals to the number of regimes (i.e., $M$) and number of columns equals to the product of the number of regimes and the total number of parameters (i.e., $(n_d+1)M$) as:
\begin{equation}
\label{MultinomialRegMat}
\left[
\begin{array}{cc;{1pt/1pt}cc;{1pt/1pt}c;{1pt/1pt}cc}
c_{11} & {\bd}^\top_{11} & c_{12} & {\bd}^\top_{12} & \cdots &c_{1M} & {\bd}^\top_{1M}\\\hdashline[1pt/1pt]
c_{21} & {\bd}^\top_{21} & c_{22} & {\bd}^\top_{22} & \cdots &c_{2M} & {\bd}^\top_{2M}\\\hdashline[1pt/1pt]
\vdots & \vdots & \vdots & \vdots & \vdots & \vdots & \vdots \\\hdashline[1pt/1pt]
c_{M1} & {\bd}^\top_{M1} & c_{M2} & {\bd}^\top_{M2} & \cdots &c_{MM} & {\bd}^\top_{MM}
\end{array}
\right].
\end{equation}
In this example, we do not have any covariates in the regime-switching (RS) functions. Thus, all the cells corresponding to the entries in ${\bd}_{lm}$ drop out. The problem then reduces to the specification of a 2 $\times$ 2 transition log-odds (LO) matrix. Here, we are interested in specifying a RS model in which conditional on any of the two previous regimes, there is a higher probability of staying in the current regime than transitioning to a different regime. To accomplish this, we first note that we set the LO entries in second column of the 2 $\times$ 2 transition LO matrix (corresponding to the Activated Regime) to zero for identification purposes. Thus, the Activated Regime serves in this case as the reference regime. The first column of the transition LO matrix, which consists of freely estimated LO parameters named \code{c11} and \code{c21}, is populated with the starting values of: (1) \code{c11} $= 0.7$, corresponding to $exp(0.7)$ $=$ $2.01$ times greater LO of staying within the Deactivated Regime as transitioning from the Deactivated to the Activated Regime, the reference regime; and (2) \code{c21} $=-1$, corresponding to $exp(-1)$ $=$ $0.37$ times lower LO of transitioning from the Activated Regime to the Deactivated Regime relative to the LO of staying Activated.
<<include = TRUE, cache = FALSE, echo = TRUE,eval=FALSE,size='footnotesize'>>=
# ---- Regimes-switching model ----
recReg <- prep.regimes(
  values = matrix(c(.7, -1, 0, 0), 2, 2),
  params = matrix(c('c11', 'c21', 'fixed', 'fixed'), 2, 2))
@
In essence, the above code creates the following transition probability matrix:
\begin{equation}
\label{eq:regParam}
\bordermatrix{ & Deactivated_{t_{i,j+1}} & Activated_{t_{i,j+1}} \cr
Deactivated_{t_{i,j}} & \frac{exp(c_{11})}{exp(c_{11})+exp(0)} & \frac{exp(0)}{exp(c_{11})+exp(0)} \cr
Activated_{t_{i,j}} & \frac{exp(c_{21})}{exp(c_{21})+exp(0)} & \frac{exp(0)}{exp(c_{21})+exp(0)} \cr}
\end{equation}
with starting values of
\begin{equation}
\label{eq:refParamStart}
\bordermatrix{ & Deactivated_{t_{i,j+1}} & Activated_{t_{i,j+1}} \cr
Deactivated_{t_{i,j}} & .668 & .332 \cr
Activated_{t_{i,j}} & .269 & .731 \cr}.
\end{equation}

In many situations it is useful to specify the structure of the transition LO matrix in deviation form - that is, to express the LO intercepts in all but the reference regime as deviations from the LO intercept in the reference regime. This creates a comparison class with all other transition intercepts evaluated as compared to that class. Note that the deviation reference regime differs from that described previously. The former description had only a reference {\em column}, whereas the deviation reference regime adds to this a reference {\em row}. In the deviation case it is expedient to reformulate the intercept as the sum of a baseline and a deviation:
\begin{equation}
c_{lm} = c_{m} + c_{\Delta,lm}
\end{equation}
where $c_{m}$ denotes the logit intercept for the probability of switching into latent class $m$ from the reference {\em row} class, $c_{\Delta,lm}$ denotes the deviation in LO of switching into latent class $m$ at time $t_{i,j}$ from latent class $l$ (i.e., from $\Stb = l$ to $\St = m$), as compared to switching from the reference {\em row} class. In this case, the multinomial logistic regression equation in Equation \ref{MultinomialRegNoDev} now appears as:
\begin{eqnarray}
\Pr\big(\St = m | \Stb=l, \nbf{x}_i(t_{i,j})\big) &\definedas \pi_{lm,it} = \frac{\exp(c_{m} + c_{\Delta,lm} +{\bd}^T_{lm}\nbf{x}_i(t_{i,j})}{\sum_{k
=1}^{M} \exp(c_{k} + c_{\Delta, lk} +{\bd}^T_{lk}\nbf{x}_i(t_{i,j}))}, \label{MultinomialReg}
\end{eqnarray} and all parameters in \refequ{MultinomialReg} can be summarized into
% \begin{equation}
% \label{MultinomialRegNoDevMat}
% \left[
% \begin{array}{cc;{2pt/2pt}cc;{2pt/2pt}c;{2pt/2pt}cc}
% c_{1}+c_{\Delta,11} & {\bd}^T_{11} & c_{2}+c_{\Delta,12} & {\bd}^T_{12} & \cdots &c_{M}+c_{\Delta,1M} & {\bd}^T_{1M}\\\hdashline[2pt/2pt]
% c_{1}+c_{\Delta,21} & {\bd}^T_{21} & c_{2}+c_{\Delta,22} & {\bd}^T_{22} & \cdots &c_{M}+c_{\Delta,2M} & {\bd}^T_{2M}\\\hdashline[2pt/2pt]
% \vdots & \vdots & \vdots & \vdots & \vdots & \vdots & \vdots \\\hdashline[2pt/2pt]
% c_{1}+ c_{\Delta,M1} & {\bd}^T_{M1} & c_{2}+c_{\Delta,M2} & {\bd}^T_{M2} & \cdots &c_{M}+c_{\Delta,MM} & {\bd}^T_{MM}
% \end{array}
% \right].
% \end{equation}
\begin{equation}
\label{MultinomialRegNoDevMat}
\left[
\begin{array}{cc;{2pt/2pt}cc;{2pt/2pt}c;{2pt/2pt}cc}
c_{\Delta,11} & {\bd}^T_{11} & c_{\Delta,12} & {\bd}^T_{12} & \cdots & c_{\Delta,1M} & {\bd}^T_{1M}\\\hdashline[2pt/2pt]
c_{\Delta,21} & {\bd}^T_{21} & c_{\Delta,22} & {\bd}^T_{22} & \cdots  &c_{\Delta,2M} & {\bd}^T_{2M}\\\hdashline[2pt/2pt]
\vdots & \vdots & \vdots & \vdots & \vdots & \vdots & \vdots \\\hdashline[2pt/2pt]
c_{1} & {\bd}^T_{M1} & c_{2} & {\bd}^T_{M2} & \cdots &c_{M} & {\bd}^T_{MM}
\end{array}
\right],
\end{equation}if we use regime M for the reference row and hence have $c_{\Delta,Ml}$ $=$ 0 for $l$ $=$ 1, $\ldots$, $M$. This allows the same parameter matrix structure for both the deviation form (Equation \ref{MultinomialRegNoDevMat}) and non-deviation form (Equation \ref{MultinomialRegMat}) of the regime switching probabilities by dropping the $c_{\Delta,Ml}$ terms in the reference row and replacing them with the $c_{m}$ logit intercepts. For convenience, both the deviation and the non-deviation form (Equations \ref{MultinomialReg} and \ref{MultinomialRegNoDev}) are available in \pkg{dynr}. For identification purposes, we can again choose regime M as the reference column and impose the constraints that $c_M$ $=$ $c_{\Delta,lM}$ $=$ 0 and ${\bd}_{lM}$ $=$ $\nbf{0}$ for $l$ $=$ 1, $\ldots$, $M$ to ensure that $\sum^M_{m=1}\pi_{lm}=1$.  Likewise, above we have shown an example that uses regime $M$ for the reference row and column, but these are independent choices that can be made by the user.  Users can specify regime-switching log-odds in deviation form by invoking the optional argument, \code{deviation=TRUE} in the call to \code{prep.regimes}.

<<include = TRUE, cache = FALSE, echo = TRUE,eval=FALSE,size='footnotesize'>>=
recReg2 <- prep.regimes(
  values = matrix(c(.8, -1, 0, 0), 2, 2),
  params = matrix(c('c_Delta11', 'c1', 'fixed', 'fixed'), 2, 2),
  deviation = TRUE, refRow = 2)
@
By default the reference row is set to the automatically detected reference column, but the code makes this choice explicit. Importantly, this code creates the same starting values as seen in Equation \ref{eq:refParamStart} but parameterized in the form of Equation \ref{MultinomialRegNoDevMat}.  The deviation form can be extremely useful for testing hypotheses about the relationships between LO intercepts and for making constraints across regimes.


\subsection{Create and cook the model}

After the recipes for all parts of the model are defined, the \code{dynr.model()} function creates the model and stores it in the \code{dynrModel} object. Each recipe (i.e., objects of class \code{dynrRecipe} created by \code{prep.*()}) and the data prepared by \code{dynr.data()} are given to this function. The function requires \code{dynamics}, \code{measurement}, \code{noise}, \code{initial}, and \code{data} as mandatory inputs for all models. When there are multiple regimes in the model, the \code{regimes} argument should be provided as shown below. When parameters are subject to transformation functions, a \code{transform} argument can be added, which will be discussed in the second example. The \code{dynr.model()} function takes the recipes and the data and combines information from both. In doing so, this function uses the information from each recipe to write the text for a \proglang{C} function. Optionally, the \proglang{C} functions can be written to a file named by the \code{outfile} argument so that the user can inspect the automatically generated \proglang{C} code. Ideally of course, there is no need to ever examine this file; however, it is sometimes useful for debugging purposes and may be helpful for specifying models that extend those supported by the \proglang{R} interface functions. More frequently, inspecting the \code{dynrModel} object and ``serving it'' will provide the needed information.
<<include = TRUE, cache = FALSE, echo = TRUE, eval=FALSE, size='footnotesize'>>=
#---- Create model  ----

rsmod <- dynr.model(
  dynamics = recDyn,
  measurement = recMeas,
  noise = recNoise,
  initial = recIni,
  regimes = recReg,
  data = EMGdata)

#---- Create model and cook it all up  ----

yum <- dynr.cook(rsmod)
@
In the last line above, the model is ``cooked'' with the \code{dynr.cook()} function to estimate the free parameters and their standard errors. When cooking, the \proglang{C} code that was written by \code{dynr.model()} is compiled and dynamically linked to the rest of the compiled \pkg{dynr} code. Then the \proglang{C} is executed to optimize the free parameters while calling the dynamically linked \proglang{C} functions that were created from the user-specified recipes. There are two points worth emphasizing in this regard. First, the user never has to write \proglang{C} functions. Second, the user benefits from the \proglang{C} functions because of their speed. In this way, \pkg{dynr} provides an \proglang{R} interface for dynamical systems modeling while maintaining much of the speed associated with \proglang{C}.

\subsection{Serve the results}
The final step associated with \pkg{dynr} modeling is serving results (a \code{dynrCook} object) after the model has been cooked. To this end, several standard, popular S3 methods are defined for the \code{dynrCook} class, including \code{coef()}, \code{confint()}, \code{deviance()}, \code{logLik()} (and thus implicitly \code{AIC()} and \code{BIC()}), \code{names()}, \code{nobs()}, \code{summary()}, and \code{vcov()}. These methods perform the same tasks as their counterparts for regression models (i.e., \code{lm} class objects). Besides, \pkg{dynr} also provides a few other model-serving functions. Here we illustrate in turn: \code{summary()}, \code{plot()}, \code{dynr.ggplot()} (or \code{autoplot()}), \code{plotFormula()}, and \code{printex()}. The \code{summary()} method provides a table of free parameter names, estimates, standard errors, t-values, and Wald-type confidence intervals. 
<<include = TRUE, cache = FALSE, echo = TRUE,eval=FALSE,size='footnotesize'>>=
#---- Serve it! ----
summary(yum)
@
These parameter estimates, standard errors, and likelihood values closely mirror those reported in \cite[p. 755-756]{Yang10a}.  In the Deactivated Regime, the autoregressive parameter (\code{phi\_1}) and the intercept (\code{mu\_1}) are lower than in the Activated Regime. So, neighboring EMG measurements are more closely related in the Activated Regime and the overall level is slightly higher. This matches very well with the idea that the Activated Regime consists of bursts of facial muscular activities and an elevated emotional state. Similarly, the effect of the self-reported emotional level is positive in the Activated Regime and fixed to zero in the Deactivated Regime. In the nested model that freely estimated this covariate effect in the Deactivated Regime, it was estimated at -0.00258 with a $t$-value of -0.097 and thus was subsequently fixed at zero. So, in the Deactivated Regime there is no relationship between the self-reported emotional level and the facial muscular activity. Essentially, in the Activated Regime the facial EMG and the self-reported emotions become coupled, but in the Deactivated Regime they are unrelated. The dynamic noise parameter gives a sense of the size of the intrinsic unmeasured disturbances that act on the system. These forces perturb the system with a typical magnitude (i.e., standard deviation) of a little less than half a point on the EMG scale seen in Figure \ref{EMGdata}(A).  Lastly, the log-odds parameters (\code{c11} and \code{c21}) can be turned into the transition probability matrix yielding
\begin{equation}
\bordermatrix{ & Deactivated_{t_{i,j+1}} & Activated_{t_{i,j+1}} \cr
Deactivated_{t_{i,j}} & .9959 & .0041 \cr
Activated_{t_{i,j}} & .0057 & .9943 \cr}
\end{equation}
which implies that both the Deactivated and the Activated Regimes are strongly persistent with high self-transistion probabilities. Next we consider some of the visualization options for serving a model.

The default \code{plot()} method is used to visualize the time series in a collection of plots: (1) a plot of time series created by \code{dynr.ggplot()} (or \code{autoplot()}), (2) a histogram of predicted regimes, and (3) a plot of equations created by \code{plotFormula()}. 
<<include = TRUE, cache = FALSE, echo = TRUE,eval=FALSE,size='footnotesize'>>=
plot(yum, dynrModel = rsmod, style = 1, textsize = 5)
@
The \code{dynr.ggplot()} (or \code{autoplot()}) method creates a plot of the smoothed state estimates overlaying the predicted regimes. It needs the result object and model object as inputs, and allows for plotting (1) user-selected smoothed state variables by default and (2) user-selected observed-versus-predicted values by setting a \code{style} to 2. An illustrative plot is created from the code below and shown in \reffig{EMGdata}(B).
<<include = TRUE, cache = FALSE, echo = TRUE,eval=FALSE,size='footnotesize'>>=
#pdf('./Figures/plotRSGG.pdf', height=7, width=12)
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 = 24),
  is.bw = TRUE)
#dev.off()

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)
@
This shows that for the first 99 seconds the participant is in the Deactivated Regime, with their latent state $\eta_i (t_{i,j+1})$ varying according to the lower autocorrelation model and having no relation to the variation in the self-reported emotional data in Figure \ref{EMGdata}(A). Then the participant switches to the Activated Regime and their latent state becomes more strongly autocorrelated and coupled to the self-report data. There follows a brief period in the Deactivated Regime around time=130 seconds with a subsequent return to the Activated Regime for the remainder of the observation.  Of course, note that Figure \ref{EMGdata}(A) shows the observed EMG data whereas Figure \ref{EMGdata}(B) shows the latent state which is related to the observed data by Equation \ref{eq:emgMeas}.

For all users,  the \code{plotFormula()} method can be used to display equations on \proglang{R} plots. Equations can be viewed in several ways after the model is specified: (1) with free parameter names and fixed values, as illustrated here in \reffig{plotformula}(A), (2) with parameter starting values, or (3) after estimation with fitted parameter values as in \reffig{plotformula}(B). Each of these desired characteristics can be embedded in the neatly typeset equations. The \code{ParameterAs} argument changes which of these characteristics is used in the equations.  Here the user-supplied parameter names and estimated parameters are typeset in \reffig{plotformula} because \code{ParameterAs} was respectively given \code{names(rsmod)}, namely, the parameter names stored in the \code{dynrModel} object, \code{rsmod}, and \code{coef(yum)}, namely, the estimated free parameter values stored in the \code{dynrCook} object, \code{yum}. Starting values for parameters are also possible values for this argument. The \code{plotFormula()} method does not require the user to install \LaTeX ~facilities and compile \LaTeX ~code in a separate step, and hence are convenient to use. To maximize the readability of the equations, it is only shown here using equations for the dynamic model and measurement model, which can be obtained by respectively setting the \code{printDyn} and \code{printMeas} arguments to true.
<<include = TRUE, cache = FALSE, echo = TRUE,eval=FALSE,size='footnotesize'>>=
plotFormula(dynrModel = rsmod, ParameterAs = names(rsmod),
  printDyn = TRUE, printMeas = TRUE) +
  ggtitle("(A)") +
  theme(plot.title = element_text(hjust = 0.5, vjust = 0.01, size = 16))

plotFormula(dynrModel = rsmod, ParameterAs = coef(yum),
  printDyn = TRUE, printMeas = TRUE) +
  ggtitle("(B)") +
  theme(plot.title = element_text(hjust = 0.5, vjust = 0.01, size = 16))
@
We can see that the equations in \reffig{plotformula}(A) are precisely those from Equations \ref{eq:emgMeas} and \ref{eq:emgDyn} which we used to define the model except that we have fixed $\beta_1$ to zero. If these equations did not match, it may indicate that we made a mistake in our model specification.

\begin{figure}[htbp]
\begin{center}
\includegraphics [width=.45\textwidth,page=1]{./Figures/Example1FormulaNames.pdf}
\includegraphics [width=.45\textwidth]{./Figures/Example1FormulaEstimates.pdf}
\end{center}
\caption{Automatic plots of model equations with (A) parameter names and (B) estimated parameters for the regime-switching linear state-space model.}
\label{plotformula}
\end{figure}

Finally, for \LaTeX ~users, the \code{printex()} method helps generate equations for the model in \LaTeX ~form.
<<include = TRUE, cache = FALSE, echo = TRUE,eval=FALSE,size='footnotesize'>>=
printex(rsmod,
  ParameterAs = names(rsmod),
  printInit = TRUE, printRS = TRUE,
  outFile = "RSLinearDiscreteYang.tex")
@
The \code{ParameterAs} argument functions the same as that in the \code{plotFormula()} method. Here we have selected to use the names of the free parameters as evidenced by giving \code{names(rsmod)} to the \code{ParameterAs} argument. In this case the initial conditions and regime-switching functions are included in the equations, as indicated by the \code{printInit} and \code{printRS} arguments being set to true. The \LaTeX ~code for the equations is written to the file specified, ``RSLinearDiscreteYang.tex'', which the user can then work with and modify as he/she wishes. Of course, this function is designed more as a convenience feature for users who are already using \LaTeX ~as a writing tool and requires all the \LaTeX-related facilities already in place on the user's computer. If so desired, the tex file can also be compiled within \proglang{R} and viewed as a pdf via the \code{texi2pdf()} function in the \pkg{tools} library:
<<include = TRUE, cache = FALSE, echo = TRUE,eval=FALSE,size='footnotesize'>>=
tools::texi2pdf("RSLinearDiscreteYang.tex")
system(paste(getOption("pdfviewer"), "RSLinearDiscreteYang.pdf"))
@

This example has used real EMG data from a previous study \cite{Yang10a} to illustrate many parts of the user-interface for \pkg{dynr}. Of particular note are the various ``serving'' functions which allow users to both verify their model and examine their results in presentation-ready formats. In the next example, we will use simulated data to further illustrate features of \pkg{dynr}, especially the nonlinear formula interface for dynamics. 

\printbibliography
\end{document}