rm(list=ls())
# Load packages
require(dynr)

0.1 A demo for how to specify dynr recipes for a process factor analysis model

#Prepare dynr recipes
#Define the dynamic model
 dynamics <- prep.matrixDynamics(
    values.dyn=matrix(c(.5, 0.1, .3, .5), ncol=2,byrow=TRUE),
    params.dyn=matrix(c('phi11', 'phi12', 'phi21', 'phi22'), ncol=2,byrow=TRUE), 
    isContinuousTime=FALSE)
 
meas <- prep.measurement(
    values.load=matrix(c(1,0,
                         2,0,
                         1,0,
                         0,1,
                         0,2,
                         0,1),ncol=2,byrow=TRUE), #Starting values Lambda
    params.load=matrix(c('fixed','fixed',
                         'lambda21','fixed',
                         'lambda31','fixed',
                         'fixed','fixed',
                         'fixed','lambda52',
                         'fixed','lambda62'),
                       ncol=2,byrow=TRUE), #Labels for fixed and freed parameters 
    values.int = matrix(rep(0.1,6),ncol=1),
    params.int = matrix(paste0('int',1:6),ncol=1), 
    state.names=c("eta1","eta2"), #Labels for latent variables in eta(t)
    obs.names=paste0('V',1:6) #Labels for observed variables in y(t)
  )

#Note that in dynr, prep.initial sets the structure of E(eta(1|0)) and Cov(eta(1|0))
#Here, initial condition covariance matrix is fixed to a diagonal matrix of 2s. 
#Could also be freely estimated with #multiple-subject data.
#Iinitial means are fixed to a vector of zeros.
  initial <- prep.initial(
    values.inistate=c(0, 0),
    params.inistate=c('fixed', 'fixed'),
    values.inicov=matrix(c(2,0,0,2),ncol=2),
    params.inicov=matrix(c('fixed','fixed','fixed','fixed'),ncol=2)) 

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

0.2 Read in data and set up data structure in dynr

ch5 = read.table('./Data/ch5_data.csv',header=TRUE,sep=",")
ch5$ID = rep(1,dim(ch5)[1]) #Add subject ID to the data set
# Data
ch52 <- dynr.data(ch5, id="ID", time="Time", observed=paste0("V",1:6))

0.3 Cook it!

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

#Use the `$' sign to set upper and lower boundaries for the parameters
#For parameters that are subjected to user-specified (e.g., via
#"prep.tfun) transformations or system transformations (e.g., variance-
#covariance parameters in the process and measurement error cov matrices),
#it may be easier to use the `@' sign to set upper and lower boundaries
#on the unconstrained (untransformed) scales - e.g., for the log of a variance
#parameter as opposed to the variance.

  model@ub[!model$param.names %in% c('psi_11','psi_12','psi_22')] = 
    c(rep(2,3), rep(5,4), rep(log(10),6))
## Warning in model@ub[!model$param.names %in% c("psi_11", "psi_12",
## "psi_22")] = c(rep(2, : number of items to replace is not a multiple of
## replacement length
  model@lb[!model$param.names %in% c('psi_11','psi_12','psi_22')] = 
    c(rep(-2,3), rep(-5,4), rep(log(1e-10),6))
## Warning in model@lb[!model$param.names %in% c("psi_11", "psi_12",
## "psi_22")] = c(rep(-2, : number of items to replace is not a multiple of
## replacement length
  res <- dynr.cook(model)
  coef(res)
##       phi11       phi21       phi12       phi22    lambda21    lambda31 
##  0.42045967  0.36392171 -0.03780944  0.45171152  2.04807601  1.02716255 
##    lambda52    lambda62        int1        int2        int3        int4 
##  2.01686329  1.01327205  0.07979076  0.24471570  0.22645828  0.20659854 
##        int5        int6      psi_11      psi_12      psi_22      var_e1 
##  0.18984754  0.10645779  2.48464489  2.02878690  8.35279101  0.49062543 
##      var_e2      var_e3      var_e4      var_e5      var_e6 
##  0.53900800  0.36952354  0.47519320  0.54501597  0.39554840
  summary(res)
## Coefficients:
##          Estimate Std. Error t value ci.lower ci.upper Pr(>|t|)    
## phi11     0.42046    0.09766   4.305  0.22905  0.61187   <2e-16 ***
## phi21     0.36392    0.17782   2.047  0.01540  0.71244   0.0216 *  
## phi12    -0.03781    0.04797  -0.788 -0.13183  0.05621   0.2162    
## phi22     0.45171    0.08692   5.197  0.28136  0.62207   <2e-16 ***
## lambda21  2.04808    0.08872  23.086  1.87420  2.22196   <2e-16 ***
## lambda31  1.02716    0.05048  20.346  0.92821  1.12611   <2e-16 ***
## lambda52  2.01686    0.04125  48.890  1.93601  2.09772   <2e-16 ***
## lambda62  1.01327    0.02441  41.515  0.96543  1.06111   <2e-16 ***
## int1      0.07979    0.22543   0.354 -0.36205  0.52163   0.3621    
## int2      0.24472    0.44833   0.546 -0.63399  1.12342   0.2932    
## int3      0.22646    0.22896   0.989 -0.22230  0.67522   0.1625    
## int4      0.20660    0.51201   0.404 -0.79693  1.21012   0.3437    
## int5      0.18985    1.02721   0.185 -1.82344  2.20313   0.4269    
## int6      0.10646    0.51805   0.205 -0.90890  1.12182   0.4188    
## psi_11    2.48464    0.36893   6.735  1.76155  3.20774   <2e-16 ***
## psi_12    2.02879    0.46464   4.366  1.11812  2.93946   <2e-16 ***
## psi_22    8.35279    1.10979   7.526  6.17764 10.52795   <2e-16 ***
## var_e1    0.49063    0.07695   6.376  0.33980  0.64145   <2e-16 ***
## var_e2    0.53901    0.19503   2.764  0.15675  0.92127   0.0034 ** 
## var_e3    0.36952    0.06509   5.677  0.24195  0.49710   <2e-16 ***
## var_e4    0.47519    0.07641   6.219  0.32543  0.62496   <2e-16 ***
## var_e5    0.54502    0.20586   2.648  0.14154  0.94849   0.0047 ** 
## var_e6    0.39555    0.07085   5.583  0.25668  0.53442   <2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## -2 log-likelihood value at convergence = 2539.05
## AIC = 2585.05
## BIC = 2650.10