--- 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 ```