---
output:
 # html_document
  pdf_document
title: "Ch 7: Fit state-space models with TVPs in dynr"
author: "Sy-Miin Chow"
date: "June 5, 2021"
toc: true
toc_depth: 2
number_sections: true
---
  
  ```{r setup, include=FALSE}
knitr::opts_chunk$set(echo = FALSE, message = FALSE)
```

```{r, echo = TRUE}
#Start fresh, set display options and load required libraries
rm(list=ls())
set.seed(12345678)
options(scipen = 999, digits=5)
require(mvtnorm)
require(dynr)
require(xtable)
require(ggplot2)
```

# Background
## Introduction
This example mirrors the illustrative simulated example 7.5.1 in Chapter 7 of Gates, Chow, & Molenaar. The data generation model features two latent factors, each identified with 3 manifest indicators. The two latent factors conform to a common trend (changing local level or set-point) that follows the dual change score model (McArdle, & Hamagami, 2001), which includes as a special form sigmoid-shaped curves with individual differences in initial level and final asymptote. For this illustrative example, we generate data from $n$ = 10 participants, each with 30 time points. To make estimation feasible given the relatively small $T$ available from each participant, all parameters, except for the time-varying ones, are assumed and constrained to be invariant across individuals. Compared to the simulation example used in 7.5.1, substantially fewer participants and slightly larger $T$ were used to reduce computational time.

We will illustrate how to use four models to represent the trend as a time-varying parameter (TVP) using an R package, dynr. The models include the:

(1) random walk model;

(2) integrated random walk model (with fixed and correctly specified initial condition or IC); 

(3) integrated random walk model with freely estimated IC means and covariance matrix; 

(4) smoothed integrated random walk model; and

(5) model ignoring time trends (TVPs constrained as invariant).

## For first-time users of dynr
The dynr package compiles C code in response to user input, so more setup is required for the dynr package than for many others, including installation of the GSL library and C compiler. Steps for installation and automated installer for Windows can be downloaded here:
\href{https://dynrr.github.io/predownload.html}{https://dynrr.github.io/predownload.html}

## Reference
Gates, K., Chow, S-M., & Molenaar, P.C.M. Analysis of Intraindividual Variation. Systems Approaches to Human Process Analysis. New York, NY: Taylor & Francis.

McArdle, J. J., & Hamagami, F. (2001). Latent difference score structural models for linear dynamic analysis with incomplete longitudinal data. In L. Collins & A. Sayer (Eds.), New methods for the analysis of change (p. 139-175). Washington, DC: American Psychological Association.

Ou, L., Hunter, M. D., & Chow, S.-M. (2019). What’s for dynr: A Package for Linear and Nonlinear Dynamic Modeling in R. The R Journal. https://journal.r-project.org/archive/2019/RJ-2019-012/index.html

Dynr development team. Dynamic modeling in R (dynr) website. \href{https://dynrr.github.io/}{https://dynrr.github.io/}

Dynr development team. User Installation Instructions Manual for dynr.

\href{https://cran.r-project.org/web/packages/dynr/vignettes/InstallationForUsers.pdf}{https://cran.r-project.org/web/packages/dynr/vignettes/InstallationForUsers.pdf}

# Simulated Data Generation
## Initial set-up
Here we define the number of time points, participants, true parameter values, other data simulation settings, and various place holders.
```{r, echo = TRUE}
time      <- 30 #Total number of time points
npad      <- 0 #Occasions to throw out to wash away the effects of initial condition
np        <- 10 #Total number of participants
ne        <- 4 #Number of latent variables
ny        <- 6 #Number of manifest variables
psi       <- matrix(c(2, .5, 0,0,  # Process noise variance-covariance matrix
                      .5, 1.5,0,0,
                      0,0,0,0,
                      0,0,0,0),   
                    ncol = ne, byrow = T)

#Factor loading matrix
lambda    <- matrix(c(1, 0,1,0,        
                      1.5,0,1,0,       
                      1, 0,1,0, 
                      0, 1, 1,0,
                      0, .8,1,0, 
                      0, 1,1,0), 
                    ncol = ne, byrow = TRUE)
theta     <- diag(.5, ncol = ny, nrow = ny)  # Measurement error variances. 
beta      <- matrix(c(0.8, 0,0,0,    # Lagged directed relations among latent variables. 
                      -.2, .7,0,0,
                      0,0,.8,1,
                      0,0,0,1), 
                    ncol = ne, byrow = TRUE)
a0        <- c(0,0,1,5)
P0        <- matrix(c(1,0,0,0,
                      0,1,0,0,
                      0,0,.5,0,
                      0,0,0,.2),
                    ncol=ne,byrow=TRUE) #Initial covariance matrix for the trend elements
yall <- matrix(0,nrow = time*np, ncol = ny)
etaall <- matrix(0,nrow = time*np, ncol = ne)
```
## Data generation
Data generation starts here. The observed data are written out to a file called PFALDS.txt.
```{r, echo = TRUE}
for (p in 1:np){
# Set up matrix for contemporaneous variables. 
  etaC      <- matrix(0, nrow = ne, ncol = time + npad) 
  etaC[,1]  <- a0 + chol(P0)%*%rnorm(ne)
  
# Latent variable residuals.
  zeta      <- mvtnorm::rmvnorm(time+npad, mean = c(0, 0,0,0), sigma = psi)
  
# Measurement errors.
  epsilon   <- rmvnorm(time+npad, mean = c(0, 0, 0, 0, 0, 0), sigma = theta)  

# Generate latent factor scores
  for (i in 2:(time+npad)){
    etaC[ ,i]   <- beta %*% etaC[ ,i-1] + zeta[i, ]
  }#End of loop over time

  eta <- t(etaC[,(npad+1):(npad+time)])

  # 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
  etaall[(1+(p-1)*time):(p*time),] = eta
} #End of p loop

yall = cbind(rep(1:np,each=time),rep(1:time,np),yall)
etaall = cbind(rep(1:np,each=time),rep(1:time,np),etaall)
colnames(etaall) = c("ID","time",paste0("LV",1:ne))
etaall = data.frame(etaall)
colnames(yall) = c("ID","time",paste0("y",1:ny))
yall = data.frame(yall)
file0 = paste0('PFALDS.txt')
write.table(yall,file0,row.names=FALSE,col.names=FALSE,sep=",")
```
## Plots of simulated data
This is the code for creating Figure 7.2(A)-(D).
```{r, echo = TRUE}
par(mfrow=c(2,2))
par(cex.axis=1.3,cex.lab=2,mgp=c(2.0,.5,0),cex.main=2)
interaction.plot(etaall$time,etaall$ID,etaall$LV1,legend=F,lty=1,
                 ylab=expression(f[i1](t)),xlab="Time",
                 main=("(A)"),xaxt = "n")

par(cex.axis=1.3,cex.lab=2,mgp=c(2.0,.5,0),cex.main=2)
interaction.plot(etaall$time,etaall$ID,etaall$LV3,legend=F,lty=1,
                 ylab=expression(setpoint[i](t)),xlab="Time",
                 main=("(B)"),xaxt = "n")

par(cex.axis=1.3,cex.lab=2,mgp=c(2.0,.5,0),cex.main=2)
interaction.plot(etaall$time,etaall$ID,etaall$LV4,legend=F,lty=1,
                 ylab=expression(paste(Delta,setpoint[i](t))),
                 xlab="Time",main=("(C)"),xaxt = "n")

par(cex.axis=1.3,cex.lab=2,mgp=c(2.0,.5,0),cex.main=2)
interaction.plot(yall$time,yall$ID,yall$y1,legend=F,lty=1,
                 ylab=expression(paste(y[i1](t))),xlab="Time",
                 main=("(D)"),xaxt = "n")
```

# Fit State-Space Models with TVPs
Now we illustrate how to fit a series of state-space models with TVPs using dynr.

## Read in data and set up data structure in dynr
```{r, echo = TRUE}
# Read in data - normally you would start right here.
ch7 <- read.table(file0,header=FALSE,sep=",")
colnames(ch7) <- c("ID","Time",paste0("V",1:ny))

# Set up dynr Data
ch72 <- dynr.data(ch7, id="ID", time="Time", 
                  observed=paste0("V",1:ny))
```

## Model 1: Random walk model for the TVP
Now we prepare the dynr recipes for a two-factor model that follows a vector-autoregressive model of order 1, as subjected to the influence of a common trend (time-varying set-point) that follows a random walk model. This model is written as:
$\vartheta_i(t) = \vartheta_i(t-1) + \zeta_{\vartheta i}(t)$

```{r, echo = TRUE}
ne        <- 3 #Number of latent variables

#Define the dynamic model
dynamics <- prep.matrixDynamics(
  values.dyn=matrix(c(.6, -0.1, 0,
                      -.1, .5,0,
                      0,0,1), ncol=ne,byrow=TRUE),
  params.dyn=matrix(c('phi11', 'phi12', 'fixed',
                      'phi21', 'phi22', 'fixed',
                      rep('fixed',ne)), ncol=ne,byrow=TRUE), 
  isContinuousTime=FALSE)

#Meausurement models
meas <- prep.measurement(
  values.load=matrix(c(1,0,1,
                       2,0,1,
                       1,0,1,
                       0,1,1,
                       0,2,1,
                       0,1,1),ncol=ne,byrow=TRUE), #Starting values for entries in Lambda
  params.load=matrix(c(rep('fixed',ne),
                       'lambda21',rep('fixed',ne-1),
                       'lambda31',rep('fixed',ne-1),
                       rep('fixed',ne),
                       'fixed','lambda52',rep('fixed',ne-2),
                       'fixed','lambda62',rep('fixed',ne-2)),
                     ncol=ne,byrow=TRUE), #Labels for fixed and freed parameters 
  state.names=c("eta1","eta2","mu1"), #Labels for latent variables in eta(t)
  obs.names=paste0('V',1:ny) #Labels for observed variables in y(t)
)

#Initial condition means and covariance matrix. Here we free them up.
initial <- prep.initial(
  values.inistate=a0[1:3],
  params.inistate=c('m1','m2','m3'),
  values.inicov=P0[1:3,1:3],
  params.inicov=diag(c('v1','v2','v3')))


#Process and measurement noise covariance matrices
mdcov <- prep.noise(
  values.latent=matrix(c(1,-.2,0,
                         -.2,1,0,
                         0,0,.5),ncol=ne,byrow=TRUE),
  params.latent=matrix(c('psi_11','psi_12','fixed',
                         'psi_12','psi_22','fixed',
                         'fixed','fixed','psi_mu1'),ncol=ne,byrow=TRUE), 
  values.observed=diag(rep(.5,ny),ny), 
  params.observed=diag(paste0('var_e',1:ny),ny)
)

#Put recipes and data together to prepare the full model
model <- dynr.model(dynamics=dynamics, measurement=meas,
                    noise=mdcov, initial=initial, data=ch72,
                    outfile="LDS1.c")

# Cook it!
res <- dynr.cook(model,debug_flag=TRUE,verbose = FALSE)
coef(res)
summary(res)
```

## Model 2: Integrated random walk (IRW) Model 
This model is written as:
$$\vartheta_i(t) = \vartheta_i(t-1) + \Delta_i(t-1)$$
$\Delta_i(t) = \Delta_i(t-1) + \zeta_{\Delta i}(t)$

```{r, echo = TRUE}
ne        <- 4 #Number of latent variables
#Define the dynamic model
dynamics.IRW <- prep.matrixDynamics(
  values.dyn=matrix(c(.6, -0.1, 0,0,
                      -.1, .5,0,0,
                      0,0,1,1,
                      0,0,0,1), ncol=ne,byrow=TRUE),
  params.dyn=matrix(c('phi11', 'phi12', 'fixed','fixed',
                      'phi21', 'phi22', 'fixed','fixed',
                      rep('fixed',ne),
                      rep('fixed',ne)), ncol=ne,byrow=TRUE), 
  isContinuousTime=FALSE)

#Measurement model
meas.IRW <- prep.measurement(
  values.load=matrix(c(1,0,1,0,
                       2,0,1,0,
                       1,0,1,0,
                       0,1,1,0,
                       0,2,1,0,
                       0,1,1,0),ncol=ne,byrow=TRUE), #Starting values for entries in Lambda
  params.load=matrix(c(rep('fixed',ne),
                       'lambda21',rep('fixed',ne-1),
                       'lambda31',rep('fixed',ne-1),
                       rep('fixed',ne),
                       'fixed','lambda52',rep('fixed',ne-2),
                       'fixed','lambda62',rep('fixed',ne-2)),
                     ncol=ne,byrow=TRUE), #Labels for fixed and freed parameters 
  state.names=c("eta1","eta2","mu1","slope1"), #Labels for latent variables in eta(t)
  obs.names=paste0('V',1:ny) #Labels for observed variables in y(t)
)

#Initial condition. Here set to the known and correctly specified
#mean vector and covariance matrix
initial.IRW <- prep.initial(
  values.inistate=a0,
  params.inistate=c('fixed','fixed','fixed','fixed'),
  values.inicov=P0,
  params.inicov=matrix(rep('fixed',ne*ne),ncol=ne,byrow=TRUE))


#Process and measurement noise covariance matrices
mdcov.IRW <- prep.noise(
  values.latent=matrix(c(.5,-.2,0,0,
                         -.2,.5,0,0,
                         0,0,0,0,
                         0,0,0,.05),ncol=ne,byrow=TRUE),
  params.latent=matrix(c('psi_11','psi_12','fixed','fixed',
                         'psi_12','psi_22','fixed','fixed',
                         'fixed','fixed','fixed','fixed',
                         'fixed','fixed','fixed','psi_sl1'),ncol=ne,byrow=TRUE), 
  values.observed=diag(rep(.5,ny),ny), 
  params.observed=diag(paste0('var_e',1:ny),ny)
)

#Put recipes and data together to prepare the full model
model2 <- dynr.model(dynamics=dynamics.IRW, measurement=meas.IRW,
                     noise=mdcov.IRW, initial=initial.IRW, data=ch72,
                     outfile="LDS2.c")

#Cook it!
res2 <- dynr.cook(model2,debug_flag=TRUE,verbose = FALSE)
coef(res2)
summary(res2)
```


## Model 3: IRW model with empirically determined initial conditions (ICs)
Here we use the empirical means and variances from the first time point and the difference between the second and first time point to set the initial conditions for the local level and slope in the IRW model. This is in contrast to Model 2, in which these IC-related parameters are freely estimated. The code model for this model is essentially identical to that for Model 2:
$$\vartheta_i(t) = \vartheta_i(t-1) + \Delta_i(t-1)$$
$$\Delta_i(t) = \Delta_i(t-1) + \zeta_{\Delta i}(t)$$,
except for the IC specification using prep.initial.
```{r, echo = TRUE}
ytemp <- yall[yall$time <= 2,]
m3<-mean(unlist(ytemp[ytemp$time==1,3:8]),na.rm=T)
m4<-mean(unlist(ytemp[ytemp$time==2,3:8]-ytemp[ytemp$time==1,3:8]),na.rm=T)
v3<-var(rowMeans(ytemp[ytemp$time==1,3:8],na.rm=TRUE),na.rm=TRUE)
v4<-var(rowMeans(ytemp[ytemp$time==2,3:8]- ytemp[ytemp$time==1,3:8],na.rm=TRUE),na.rm=TRUE)

initial.free <- prep.initial(
  values.inistate=c(0, 0,m3,m4),
  params.inistate=c('m1','m2','m3','fixed'),
  values.inicov=matrix(c(1,0,0,0,
                         0,1,0,0,
                         0,0,v3,0,
                         0,0,0,v4),ncol=ne,byrow=TRUE),
  params.inicov=diag(c('v1','v2','v3','fixed')))

#Put recipes and data together to prepare the full model
model3 <- dynr.model(dynamics=dynamics.IRW, measurement=meas.IRW,
                     noise=mdcov.IRW, initial=initial.free, data=ch72,
                     outfile="LDS3.c")
#Cook it!
res3 <- dynr.cook(model3,debug_flag=TRUE,verbose = FALSE)
coef(res3)
summary(res3)
```


##  Model 4: Smoothed Integrated Random Walk (SIRW) Model
Note that this model is very similar to the IRW, with the exception that a parameter, $\alpha_1$, is estimated. This would induce some smoothing effect on the time-varying set-point when the absolute value of alpha1 is $<$ 1.
This model is written as:
$$\vartheta_i(t) = \alpha\vartheta_i(t-1) + \Delta_i(t-1)$$
$$\Delta_i(t) = \Delta_i(t-1) + \zeta_{\Delta i}(t).$$

```{r, echo = TRUE}
ne        <- 4 #Number of latent variables
#Prepare dynr recipes
#Define the dynamic model
dynamics.SIRW <- prep.matrixDynamics(
  values.dyn=matrix(c(.6, -0.1, 0,0,
                      -.1, .5,0,0,
                      0,0,.9,1,
                      0,0,0,1), ncol=ne,byrow=TRUE),
  params.dyn=matrix(c('phi11', 'phi12', 'fixed','fixed',
                      'phi21', 'phi22', 'fixed','fixed',
                      rep('fixed',2),'alpha1','fixed',
                      rep('fixed',ne)), ncol=ne,byrow=TRUE), 
  isContinuousTime=FALSE)

## Cook it!
#Put recipes and data together to prepare the full model
#We use fixed and correctly specified IC here. So all 
#subrecipes from this model are identical to those from
#Model 2 (IRW with correctly specified IC), except for
#the dynamic model (dynamics.SIRW in this case).
model4 <- dynr.model(dynamics=dynamics.SIRW, measurement=meas.IRW,
                     noise=mdcov.IRW, initial=initial.IRW, data=ch72,
                     outfile="LDS4.c")
res4 <- dynr.cook(model4,debug_flag=TRUE,verbose = FALSE)
coef(res4)
summary(res4)
```


## Model 5: Time-invariant trend 
Here to specify a time-invariant model, we use Model 1, the random walk model. All subrecipes are identical to Model 1, except that the process noise variances for the TVPs are constrained to be zero.
```{r, echo = TRUE}
ne        <- 3 #Number of latent variables

#Process and measurement noise covariance matrices
mdcov.inv <- prep.noise(
  values.latent=matrix(c(.5,-.2,0,
                         -.2,.5,0,
                         0,0,0),ncol=ne,byrow=TRUE),
  params.latent=matrix(c('psi_11','psi_12','fixed',
                         'psi_12','psi_22','fixed',
                         'fixed','fixed','fixed'),ncol=ne,byrow=TRUE), 
  values.observed=diag(rep(.5,ny),ny), 
  params.observed=diag(paste0('var_e',1:ny),ny)
)

#Initial condition means and covariance matrix. 
initial.inv <- prep.initial(
  values.inistate=a0[1:3],
  params.inistate=c('m1','m2','m3'),
  values.inicov=P0[1:3,1:3],
  params.inicov=diag(c('v1','v2','v3')))

model5 <- dynr.model(dynamics=dynamics, measurement=meas,
                     noise=mdcov.inv, initial=initial.inv, data=ch72,
                     outfile="LDS5.c")
# Cook it!
res5 <- dynr.cook(model5,debug_flag=TRUE,verbose = FALSE)
round(coef(res5),5)
summary(res5)
```


# Compare Estimation Results
Here we show how estimation results returned by dynr's cooking can be compared with the true values used in the data generation. We use the R package, xtable, to compile table of results.

## Extract estimates from different models and combine them
```{r, echo = TRUE}
# Extract true values of the latent variables, including the time-varying set-point
true <- etaall[,c(3:5)]

# True values of all parameters that are reasonable to be compared
truepar <- c(.8,-.2,0,.7,1.5,1,.8,1,2,.5,1.5,rep(.5,ny)) 

# Extract the corresponding estimated parameter values from different models
allPar <- cbind(truepar,
coef(res)[!model$'param.names' %in% c(paste0('m',1:3),paste0('v',1:3),'psi_mu1')],
coef(res2)[!model2$'param.names' %in% c('psi_sl1')],
coef(res3)[!model3$'param.names' %in% c('psi_sl1',paste0('m',1:4),paste0('v',1:4))],
coef(res4)[!model4$'param.names' %in% c('psi_sl1','alpha1')],
coef(res5)[!model5$'param.names' %in% c(paste0('m',1:3),paste0('v',1:3),'psi_mu1')])

colnames(allPar) <- c("True","RW","IRW (fixed IC)","IRW (free IC)",
                   "SIRW","Invariant")
rownames(allPar ) <- c("$\\phi_{11}$","$\\phi_{21}$","$\\phi_{12}$","$\\phi_{22}$",
                       "$\\lambda_{21}$","$\\lambda_{31}$",
                       "$\\lambda_{52}$","$\\lambda_{62}$",
                       "$\\psi_{11}$","$\\psi_{12}$","$\\psi_{22}$",
                       "$\\Var(\\epsilon_{1})$","$\\Var(\\epsilon_{2})$",
                       "$\\Var(\\epsilon_{3})$","$\\Var(\\epsilon_{4})$",
                       "$\\Var(\\epsilon_{5})$","$\\Var(\\epsilon_{6})$")

sum = xtable(allPar,include.rownames=FALSE,
             caption="Parameter Estimates for Illustrative Example 1.",
             digits=rep(2,dim(allPar)[2]+1))
write(paste(print(sum,sanitize.text.function = function(x) x)),
      file = "demo1Est.txt")

# Extract smoothed latent variables estimates

#First check the dimension of the smoothed latent variable estimates 
#stored in the cooked dynr object. This is of dimension ne x np*time.
dim(res@eta_smooth_final)

# Column-bind estimates of latent variable scores from different models.
etaSmoothAll <- cbind(rep(1:time,np),
                      rep(1:np,each=time),
                      etaall[,c(3:5)],
                      t(res@eta_smooth_final[1:3,]),
                      t(res2@eta_smooth_final[c(1:3),]), 
                      t(res3@eta_smooth_final[c(1:3),]), 
                      t(res4@eta_smooth_final[c(1:3),]), 
                      t(res5@eta_smooth_final[1:3,]))

#Convert into a data frame and apply column names
etaSmoothAll <- data.frame(etaSmoothAll)
colnames(etaSmoothAll) <- c("Time","ID",
                            paste0("Model0LV",1:3),
                            paste0("Model1LV",1:3),
                            paste0("Model2LV",1:3),
                            paste0("Model3LV",1:3),
                            paste0("Model4LV",1:3),
                            paste0("Model5LV",1:3))


#Reshape from wide to long format
eta2 <- reshape(etaSmoothAll, direction='long', 
                varying=3:dim(etaSmoothAll)[2], 
                timevar='Model',
                times=c("True","RW","IRW fixed IC","IRW free IC",
                        "SIRW", "Invariant"),
                v.names=paste0('LV',1:3),
                idvar=c("Time","ID"))

#Compile a data frame of biases to make it easier to compare estimates
etabiasAll <- cbind(etaSmoothAll$Time,etaSmoothAll$ID,
  etaSmoothAll[,6:20] - do.call(cbind, 
                                replicate(5, etaSmoothAll[,3:5], 
                                          simplify=FALSE)))
#Convert into a data frame and apply column names
etabiasAll <- data.frame(etabiasAll)
colnames(etabiasAll) <- c("Time","ID",
                            paste0("Model1LV",1:3),
                            paste0("Model2LV",1:3),
                            paste0("Model3LV",1:3),
                            paste0("Model4LV",1:3),
                            paste0("Model5LV",1:3))


#Reshape from wide to long format
eta3 <- reshape(etabiasAll, direction='long', 
                varying=3:dim(etabiasAll)[2], 
                timevar='Model',
                times=c("RW","IRW fixed IC","IRW free IC",
                        "SIRW", "Invariant"),
                v.names=paste0('BiasLV',1:3),
                idvar=c("Time","ID"))
```

##Plot and summarize comparisons of results
Code for this section can be used to replicate the results summarized in Tables 7.2-7.3, and Figures 7.3(A)-(D)
```{r, echo = TRUE}
#Use ggplot to compare the latent variables estimates - only one ID is used here.
eta.temp <- eta2[eta2$ID==10,]
eta.temp2 <- eta3[eta3$ID==10,]

ggplot(data = eta.temp, aes(x = Time, y = LV1, colour = Model, shape=Model)) +       
  geom_line() + geom_point() + ggtitle("(A)")+
  ylab(expression(paste("Smoothed estimates of ",f[i1](t))))+
  theme(plot.title = element_text(hjust = 0.5))+
  theme(text = element_text(size=20))


ggplot(data = eta.temp, aes(x = Time, y = LV3, colour = Model, shape=Model)) +       
  geom_line() + geom_point() + ggtitle("(B)")+
  ylab("Smoothed estimates of setpoint")+
  theme(plot.title = element_text(hjust = 0.5))+
  theme(text = element_text(size=20))

ggplot(data = eta.temp2, aes(x = Time, y = BiasLV3, colour = Model, shape=Model)) +       
  geom_line() + geom_point() + ggtitle("(C)")+
  ylab("Biases in estimates of setpoint")+
  theme(plot.title = element_text(hjust = 0.5))+
  theme(text = element_text(size=20))


#Calculate average biases for the latent variables
eta.bias <- data.frame(
  biasLV1 = colMeans(etaSmoothAll[,grep("LV1",colnames(etaSmoothAll))] - true[,1]),
  biasLV2 = colMeans(etaSmoothAll[,grep("LV2",colnames(etaSmoothAll))] - true[,2]),
  biasLV3 = colMeans(etaSmoothAll[,grep("LV3",colnames(etaSmoothAll))] - true[,3])
)              
eta.bias <- eta.bias[2:6,]
row.names(eta.bias) <- c("RW","IRW fixed IC","IRW free IC","SIRW", "Invariant")

#Calculate biases in parameter estimates for a single replication
parEst.bias <- data.frame(
  RW = mean(coef(res)[!model$'param.names' %in%
       c(paste0('m',1:3),paste0('v',1:3),'psi_mu1')]-truepar),
  IRW_fixedIC = mean(coef(res2)[!model2$'param.names' %in% 
                c('psi_sl1')]-truepar),
  IRW_freeIC = mean(coef(res3)[!model3$'param.names' %in%
               c('psi_sl1',paste0('m',1:4),paste0('v',1:4))]-truepar),
  SIRW = mean(coef(res4)[!model4$'param.names' %in% 
         c('psi_sl1','alpha1')]-truepar),
  TimeInvariant = mean(coef(res5)[!model5$'param.names' %in% 
                  c(paste0('m',1:3),paste0('v',1:3),'psi_mu1')]-truepar)
)

sum2 = xtable(eta.bias,include.rownames=FALSE,
             caption="Biases in Latent Variable Biases for Illustrative Example 1.",
             digits=rep(2,dim(eta.bias)[2]+1))
write(paste(print(sum2,sanitize.text.function = function(x) x)),
      file = "demo1Est.txt",append=TRUE)
write(paste(print(parEst.bias,sanitize.text.function = function(x) x)),
      file = "demo1Est.txt",append=TRUE)
parEst.bias
eta.bias
sum2
```