In this practical you will practice with fitting multiple (Bayesian) single subject models for intensive longitudinal data, specifically AR(1) and AR(2) models, an AR(1) model with a covariate, a VAR(1) model and two models that account for measurement error. You will fit the models for the data of two participants to get an impression of individual differences in dynamics.

Installing and loading packages

The analyses depend on four R-packages: Hmisc, for loading the spss data file; rjags, for fitting the Bayesian models; foreign, for loading an SPSS data file; coda, for running summary stats and plots on the Bayesian analysis results. In order for these analyses to run on your computer, you will also need to download and install JAGS, because rjags depends on this program.

#install(foreign)
#install(rjags)
#install(coda)
library(foreign)
library(rjags)
library(coda)

Loading the data

We will work with the two spss files for two female participants Jane and Mary (Names ficional, data courtesy of Schuurman, Houtveen and Hamaker, 2015). The data consists of daily diary measurements, for 91 days for Jane and 105 days for Mary. The data contains various measures on mood, positive and negative affect and symptoms surrounding the women’s periods. Load the packages with help from package foreign:

mydata_Jane= read.spss("dat2.sav", use.value.labels = FALSE, to.data.frame = FALSE, trim.factor.names = FALSE,
                  trim_values = TRUE, reencode = NA, use.missings = TRUE)  # This may return a warning about a "record type 7", but this warning can be ignored.
## re-encoding from CP1252
mydata_Mary= read.spss("dat4.sav", use.value.labels = FALSE, to.data.frame = FALSE, trim.factor.names = FALSE,
                  trim_values = TRUE, reencode = NA, use.missings = TRUE)  # This may return a warning about a "record type 7", but this warning can be ignored.
## re-encoding from CP1252

Note that the variable names in the files are in Dutch, but we will translate the relevant variables on the go when we use them.

We will start by looking at both women’s rating (scale 1-100 low-high) of their daily overall mood. Make a time series plot for the mood of both participants, as wel as inspect the autocorrelations and partial autocorrelations for the data of each women. You can use functions ts.plot(), acf(), and pacf() for this respectively. We have prepared this code for you for Jane, but not for Mary. Compare what you see for the two women.

###Extracting mood from Jane's dataset##
dat_Jane=data.frame(mydata_Jane$scale)
dat_Jane=as.matrix(dat_Jane)
colnames(dat_Jane) = c("mood")

###Extracting mood from Mary's dataset##
dat_Mary=data.frame(mydata_Mary$scale)
dat_Mary=as.matrix(dat_Jane)
colnames(dat_Mary) = c("mood")

###Exploring the data for Jane: Time series plot, autocorrelation of the data, and partial autocorrelations of the data.
ts.plot(dat_Jane,col="mediumseagreen",lwd=5, ylab= "mood")

acf(dat_Jane,na.action=na.pass)

pacf(dat_Jane,na.action=na.pass)

mean(dat_Jane, na.rm=TRUE)
## [1] 74.93421
var(dat_Jane, na.rm=TRUE)
##          mood
## mood 173.6623
###Exploring the data for Mary: Time series plot, autocorrelation of the data, and partial autocorrelations of the data.

Basic AR models

Based on the autocorrelations and partial autocorrelations for the different lags of the women’s data, what order of autoregressive model would you fit for each of them (how many lags would you like to include)?

We can easily fit n=1 time series models with package arima which is native in R. Fit at least an order 1, and an order 2 autoregressive model for both women using package arima. ARIMA uses a state space modeling framework which combines maximum likelihood estimation and the kalman filter. Use ?arima in R to learn more. Later we will fit the same models with Bayesian estimation techniques. Interpret and compare the results for the two women.

###AR(1) and AR(2) model fitted for Jane in arima##
arima(dat_Jane, c(1,0,0))
## 
## Call:
## arima(x = dat_Jane, order = c(1, 0, 0))
## 
## Coefficients:
##          ar1  intercept
##       0.0752    74.8807
## s.e.  0.1321     1.6063
## 
## sigma^2 estimated as 170.5:  log likelihood = -303.15,  aic = 612.29
arima(dat_Jane, c(2,0,0))
## 
## Call:
## arima(x = dat_Jane, order = c(2, 0, 0))
## 
## Coefficients:
##          ar1     ar2  intercept
##       0.0840  0.4070    74.7038
## s.e.  0.1053  0.1089     2.4736
## 
## sigma^2 estimated as 141.8:  log likelihood = -297.22,  aic = 602.43
###AR(1) and AR(2) model fitted for Mary in arima##

Bayesian n=1 AR1 Analysis

Now we will fit the same model in a Bayesian way in JAGS. This is the JAGS code for the model (everything in between parentheses), which we will supply to jags for the analysis. Take a look at the model specification, including the priors that were specified. The priors specified are intended to be uninformative.

AR1model <-"
model{
   

#####model######
    for (t in 2:nt) ###model for timepoint 2 onwards
    {
        y[t]~dnorm(muy[t],Ipre)  ###Ipre is innovation precision
        muy[t]<- mu + phi*z[t-1] ## mu is mean, phi is autoregressive parameters
        z[t] <- y[t]-mu     
        
    }

####timepoint 1#######

    z[1] <- y[1]-mu  ###model for time point one.
    
    
            
####priors###########
       
    Ipre <- inverse(Ivar) ## innovation precision
    Ivar ~dunif(0,300)   ##innovation variance

    phi~dunif(-1,1) ##autoregression coefficient
    mu~dnorm(50,.001) ##mean
        
              

}

"
AR1model.txt<-textConnection(AR1model)

We will fit the model in JAGS via package rjags with three chains. For each chain, initial values for each parameter can be specified in lists as follows

#specify initial values for each chain
ins = list(list( phi=.6, Ivar = 100,  mu=50),
           list( phi=.3, Ivar = 60,  mu=60),
           list( phi=-.2, Ivar = 80, mu = 70))  ###initial values for three chains

#specify the number of time points in the data (needed as data for the jags model)
nt=length(dat_Jane)

It is also possible to not specify any starting values, in which case the will be generated by jags based on the specified priors.

Next, we specify the model file, data to be used with the model file (for Jane), initial values, and the number of chains, for compiling the model in JAGS and rjags. We use a burnin of 20000 samples per chain. We take 20000 samples on which we will base our results.

###Specify modelfile, data, initial values, and number of chains for compiling the rjags model:
ins = list(list( phi=.6, Ivar = 100,  mu=50),
           list( phi=.3, Ivar = 60,  mu=60),
           list( phi=-.2, Ivar = 80, mu = 70))  ###initial values

jags_AR1 <- jags.model('AR1model.txt',
                   data = list('nt' = nt,
                               'y' = as.numeric(dat_Jane)
                               
                   ),
                   inits=ins,
                   n.chains = 3) ###compile and initialize jags AR model
## Compiling model graph
##    Resolving undeclared variables
##    Allocating nodes
## Graph information:
##    Observed stochastic nodes: 75
##    Unobserved stochastic nodes: 18
##    Total graph size: 184
## 
## Initializing model
update(jags_AR1, 2000) ####burn-in samples/iterations, to be discarded

samps_AR1=coda.samples(jags_AR1,
                   c('phi','mu','Ivar'),
                   2000) ### take the samplesiterations we will use for the results

The model has been fitted. Before showing the results, we’ll do some convergence checks, specifically, we check the density plots, trace plots, and the gelman-rubin statistics (plots and summary stats):

plot(samps_AR1,trace=FALSE) ##plot densities, should look smooth, unimodal, usually normal

gelman.diag(samps_AR1) ###plot gelman rubin stats Should be 1 or very very near 1 (smaller than 1.05)
## Potential scale reduction factors:
## 
##      Point est. Upper C.I.
## Ivar          1       1.01
## mu            1       1.00
## phi           1       1.00
## 
## Multivariate psrf
## 
## 1
plot(samps_AR1,density=FALSE) ##traceplots should look like a pile of different colored hairy catterpillars

The density plots are nice and smooth, trace plots look like fat catterpillars, gelman-rubin statistics are good. Convergence seems in order. The resulting parameter estimates:

summary(samps_AR1) ##results parameter estimates
## 
## Iterations = 3001:5000
## Thinning interval = 1 
## Number of chains = 3 
## Sample size per chain = 2000 
## 
## 1. Empirical mean and standard deviation for each variable,
##    plus standard error of the mean:
## 
##           Mean      SD Naive SE Time-series SE
## Ivar 170.29178 30.0324 0.387717       0.604890
## mu    75.35426  1.6707 0.021568       0.025582
## phi    0.08058  0.1237 0.001597       0.002551
## 
## 2. Quantiles for each variable:
## 
##          2.5%        25%       50%      75%    97.5%
## Ivar 121.2057 148.879378 167.04179 187.8023 239.3726
## mu    72.0328  74.275470  75.37984  76.4741  78.5375
## phi   -0.1725  -0.001797   0.08144   0.1648   0.3148

Rerun the model with a larger number of burnin and iterations. Do the results change (if the model has converged, the changes should be very small). Are the results similar to those from arima? Now fit the Bayesian AR(1) model for Mary.

Bayesian n=1 AR2 Analysis

Next, fit the AR2 model for Jane and Mary with the following model code:

AR2model <-"model{
   

#####model######
    for (t in 3:nt) ###model for timepoint 3 onwards
    {
        y[t]~dnorm(muy[t],Ipre)  ###Ipre is innovation precision
        muy[t]<- mu + phi*z[t-1] + phi2*z[t-2] ## mu is mean, phi is autoregressive parameter for lag 1, phi2 for lag 2
        z[t] <- y[t]-mu     
        
    }

####timepoint 1#######

    z[2] <- y[2]-mu  ###model for time point two.
    z[1] <- y[1]-mu ###model for time point one.
    
            
####priors###########
       
    Ipre <- inverse(Ivar) ## innovation precision
    Ivar ~dunif(0,300)   ##innovation variance

    phi~dunif(-1,1) ##autoregression coefficient
    phi2~dunif(-1,1) ##autoregression coefficient
    mu~dnorm(50,.001) ##mean
        
              

}

"
AR2model.txt<-textConnection(AR2model)
ins = list(list( phi=.6,phi2=.2, Ivar = 100,  mu=50),
           list( phi=.3,phi2=-.5, Ivar = 60,  mu=60),
           list( phi=-.2,phi2=0, Ivar = 80, mu = 70))  

jags_AR2 <- jags.model('AR2model.txt',
                   data = list('nt' = nt,
                               'y' = as.numeric(dat_Jane)
                               
                   ),
                   inits=ins,
                   n.chains = 3) 
## Compiling model graph
##    Resolving undeclared variables
##    Allocating nodes
## Graph information:
##    Observed stochastic nodes: 74
##    Unobserved stochastic nodes: 19
##    Total graph size: 246
## 
## Initializing model
update(jags_AR2, 3000) 

samps_AR2=coda.samples(jags_AR2,
                   c('phi','phi2', 'mu','Ivar'),
                   3000) 

plot(samps_AR2,trace=FALSE) 

gelman.diag(samps_AR2) 
## Potential scale reduction factors:
## 
##      Point est. Upper C.I.
## Ivar       1.00       1.00
## mu         1.02       1.02
## phi        1.00       1.01
## phi2       1.00       1.00
## 
## Multivariate psrf
## 
## 1
plot(samps_AR2,density=FALSE) 

summary(samps_AR2) 
## 
## Iterations = 4001:7000
## Thinning interval = 1 
## Number of chains = 3 
## Sample size per chain = 3000 
## 
## 1. Empirical mean and standard deviation for each variable,
##    plus standard error of the mean:
## 
##           Mean      SD Naive SE Time-series SE
## Ivar 143.38553 25.1256 0.264847       0.424922
## mu    76.10917  2.9208 0.030788       0.034063
## phi    0.04578  0.1114 0.001175       0.001805
## phi2   0.38787  0.1096 0.001155       0.001791
## 
## 2. Quantiles for each variable:
## 
##          2.5%       25%       50%      75%    97.5%
## Ivar 102.4454 125.68322 140.92909 157.6916 199.1517
## mu    71.0420  74.49574  76.05178  77.6841  81.7098
## phi   -0.1744  -0.02686   0.04585   0.1210   0.2650
## phi2   0.1667   0.31302   0.38810   0.4631   0.5996

Has the model converged? If so, interpret the results; in the summary results you find the mean estimate for each parameter and the posterior standard deviation, and the quantiles for the posterior distribution. The first and last quantile in the results can be used to form the lower and upper bound of a 95% credible interval for each parameter.

Does an AR(1) model seem to suit Jane best, or an AR(2) model, based on these results? What about for Mary?

AR1 model with a covariate

Next, lets add an exogenous covariate to the AR(1) model. We have information on when each women had her period - we will use this to predict each woman’s mood. For this we will load an external data file with the same mood data, and data on each women’s period (0 of not on period, 1 if on period). For the latter variable, for this example we have filled in scores where there were missing observations (because exogenous variables are not allowed to have missing observations in the analysis). An alternative way to deal with this is to make the covariate an endogenous variable (or impute for these observations in another way), but this is beyond the scope of this workshop.

Load the data for Mary and Jane, and explore their data with (for instance) a time series plot:

dat_predJane = read.table("Ar1pred_Jane.dat", header=TRUE, na.strings="*")
dat_predMary = read.table("Ar1pred_Mary.dat", header=TRUE, na.strings="*")
mood_J = dat_predJane[,1]
period_J =dat_predJane[,2]
mood_M = dat_predMary[,1]
period_M =dat_predMary[,2]

ts.plot(mood_J/10,col="mediumseagreen",lwd=5, ylab= "scores", ylim=c(0,10))
lines(period_J,col="pink",lwd=5)
lines(period_M,col="orange",lwd=5)
lines(mood_M/10,col="blue",lwd=5)

Take a look at the JAGS model code for the AR(1) model with an exogenous predictor. What are coefficients b1 and b2? Further, note that we are estimating an intercept rather than a mean for mood.

AR1predmodel <-"model{
   

#####model######
    for (t in 2:nt) ###model for timepoint 2 onwards
    {
        y[t]~dnorm(muy[t],Ipre)  ###Ipre is innovation precision
        muy[t]<- int + phi*y[t-1] + b1*x[t] + b2*x[t-1] ## mu is mean, phi is autoregressive parameters
        
    }

    #####priors###########
       
    Ipre <- inverse(Ivar) ## innovation precision
    Ivar ~dunif(0,500)   ##innovation variance

    phi~dunif(-1,1) ##autoregression coefficient
    b1~dnorm(0,.001) ##autoregression coefficient
    b2~dnorm(0,.001) ##autoregression coefficient   
    int~dnorm(50,.001) ##mean
        
              

}
"
AR1predmodel.txt<-textConnection(AR1predmodel)

Fit the model to the data for Jane and Mary, check the convergence, and interpret and compare the results.

####Jane#####
jags <- jags.model('AR1predmodel.txt',
                   data = list('nt' = length(mood_J),
                               'y' = as.numeric(mood_J),
                               'x' = as.numeric(period_J)
                               
                   ),
                   n.chains = 3) ###compile and initialize jags AR model
## Compiling model graph
##    Resolving undeclared variables
##    Allocating nodes
## Graph information:
##    Observed stochastic nodes: 74
##    Unobserved stochastic nodes: 20
##    Total graph size: 267
## 
## Initializing model
update(jags, 4000) ####burnin

samps=coda.samples(jags,
                   c('phi','b1', 'b2','int', 'Ivar'),
                   10000) ### take samples next to burnin


plot(samps,trace=FALSE) ##plot densities

gelman.diag(samps) ###plot gelman rubin stats
## Potential scale reduction factors:
## 
##      Point est. Upper C.I.
## Ivar       1.00       1.00
## b1         1.00       1.00
## b2         1.00       1.00
## int        1.01       1.02
## phi        1.01       1.02
## 
## Multivariate psrf
## 
## 1
plot(samps,density=FALSE) ##traceplots

summary(samps) ##results parameters estimates 
## 
## Iterations = 5001:15000
## Thinning interval = 1 
## Number of chains = 3 
## Sample size per chain = 10000 
## 
## 1. Empirical mean and standard deviation for each variable,
##    plus standard error of the mean:
## 
##           Mean      SD  Naive SE Time-series SE
## Ivar 169.59087 29.9784 0.1730802       0.287338
## b1     8.37124  7.4803 0.0431875       0.166928
## b2    -9.91473  7.9058 0.0456439       0.172831
## int   73.49214  9.5006 0.0548518       0.588733
## phi    0.02695  0.1266 0.0007312       0.008161
## 
## 2. Quantiles for each variable:
## 
##          2.5%       25%       50%      75%   97.5%
## Ivar 120.3801 148.60522 166.15851 186.7446 238.950
## b1    -6.6544   3.41544   8.47017  13.3909  22.899
## b2   -25.1778 -15.27575 -10.08725  -4.6628   5.836
## int   55.3833  67.07975  73.36937  79.9400  92.172
## phi   -0.2222  -0.05901   0.02951   0.1128   0.269

Basic Bivariate VAR(1) Model

Next we will fit a multivariate, Vector Autoregressive model (of order 1); a VAR(1) model. We will fit the model for two variables, overall mood like before, and nervous tension.

Start by loading the data for both woman, and take a look at their time series data. Do the variables seem to sync up over time? What kind of effects would you expect nervous tension and mood to have on each other over time? What about concurrently?

datvar_Jane=data.frame(mydata_Jane$scale,mydata_Jane$nerveuzegespannenheid)
colnames(datvar_Jane) = c("mood", "nervoustension")

datvar_Mary=data.frame(mydata_Mary$scale,mydata_Mary$nerveuzegespannenheid)
colnames(datvar_Mary) = c("mood", "nervoustension")

mood_J = datvar_Jane[,1]
tense_J = datvar_Jane[,2]
ts.plot(mood_J/10,col="mediumseagreen",lwd=5, ylab= "scores", ylim=c(1,10))
lines(tense_J,col="royalblue3",lwd=5)

Now, take a look at the JAGS model code for the bivariate VAR(1) model.

VAR1model <-"
model{
   

#####model######
    for (t in 2:nt) ###model for timepoint 2 onwards
    {
        y[t,1:2]~dmnorm(muy[t,1:2],Ipre[1:2,1:2])  ###Ipre is innovation precision
        muy[t,1]<- mu1 + phi11*z[t-1,1] + phi12*z[t-1,2] ## mu is mean, phi is autoregressive parameters
        muy[t,2]<- mu2 + phi22*z[t-1,2] + phi21*z[t-1,1] ## mu is mean, phi is autoregressive parameters
        z[t,1] <- y[t,1]-mu1    
        z[t,2] <- y[t,2]-mu2
    }

####timepoint 1#######

    z[1,1] <- y[1,1]-mu1  ###model for time point two.
    z[1,2] <- y[1,2]-mu2 ###model for time point one.
    
            
####priors###########
       
    Ipre <- inverse(Ivar) ## innovation precision
    Ivar[1,1] ~dunif(0,300)   ##innovation variance
    Ivar[2,2] ~dunif(0,300)   ##innovation variance
    Ivar[1,2] <- Icor*sqrt(Ivar[1,1])*sqrt(Ivar[2,2])
    Ivar[2,1] <- Ivar[1,2]
    Icor~dunif(-1,1)

    phi11~dnorm(0,.001) ##autoregression coefficient
    phi12~dnorm(0,.001) ##autoregression coefficient
    phi21~dnorm(0,.001) ##autoregression coefficient
    phi22~dnorm(0,.001) ##autoregression coefficient
    
    mu1~dnorm(50,.001) ##mean
    mu2~dnorm(3,.001) ##mean
        
              

}
"
VAR1model.txt<-textConnection(VAR1model)

Fit the models for Mary and Jane.

jags <- jags.model('VAR1model.txt',
                   data = list('nt' = nt,
                               'y' = datvar_Jane
                               
                   ),
                   n.chains = 3) ###compile and initialize jags AR model
## Compiling model graph
##    Resolving undeclared variables
##    Allocating nodes
## Graph information:
##    Observed stochastic nodes: 75
##    Unobserved stochastic nodes: 24
##    Total graph size: 627
## 
## Initializing model
update(jags, 2000) ####burnin

samps=coda.samples(jags,
                   c('phi11','phi22', 'mu1','phi12','phi21', 'mu2', 'Ivar', 'Icor'),
                   3000) ### take samples next to burnin

plot(samps,trace=FALSE) ##plot densities

#gelman.diag(samps) ###plot gelman rubin stats

plot(samps,density=FALSE) ##traceplots

summary(samps) ##results parameters estimates 
## 
## Iterations = 3001:6000
## Thinning interval = 1 
## Number of chains = 3 
## Sample size per chain = 3000 
## 
## 1. Empirical mean and standard deviation for each variable,
##    plus standard error of the mean:
## 
##                Mean       SD  Naive SE Time-series SE
## Icor       -0.26220  0.11183 0.0011788      0.0019965
## Ivar[1,1] 157.68290 27.98272 0.2949638      0.5003780
## Ivar[2,1]  -4.21739  2.04457 0.0215516      0.0377438
## Ivar[1,2]  -4.21739  2.04457 0.0215516      0.0377438
## Ivar[2,2]   1.61957  0.28400 0.0029936      0.0049049
## mu1        75.80379  2.05568 0.0216688      0.0482047
## mu2         2.37521  0.31047 0.0032727      0.0071924
## phi11      -0.02531  0.12888 0.0013585      0.0021968
## phi12      -2.99398  1.16236 0.0122524      0.0191290
## phi21      -0.01562  0.01266 0.0001335      0.0002056
## phi22       0.43249  0.11765 0.0012401      0.0019703
## 
## 2. Quantiles for each variable:
## 
##                2.5%       25%       50%        75%      97.5%
## Icor       -0.46743  -0.34179  -0.26622  -0.188162  -0.036774
## Ivar[1,1] 112.41183 137.67751 154.55336 174.260736 219.148990
## Ivar[2,1]  -8.55930  -5.50088  -4.09276  -2.831859  -0.563556
## Ivar[1,2]  -8.55930  -5.50088  -4.09276  -2.831859  -0.563556
## Ivar[2,2]   1.15765   1.41522   1.58672   1.791997   2.242123
## mu1        71.73091  74.51644  75.80350  77.131441  79.851295
## mu2         1.76201   2.18587   2.36992   2.566867   2.990892
## phi11      -0.28290  -0.11166  -0.02218   0.060877   0.226732
## phi12      -5.26378  -3.77231  -2.98953  -2.214692  -0.734975
## phi21      -0.04065  -0.02398  -0.01554  -0.007294   0.009178
## phi22       0.19801   0.35394   0.43320   0.511325   0.664079

Evaluate the convergence of the model. Draw a path model based on the results for Jane and for Mary, leaving out paths for which the sign of the effect is unclear based on the 95% credible intervals (see quantiles).

AR(1) Model that accounts for measurement error

The previous models did not account for measurement error, but it is possible to account for measurement error in dynamic models, even with only one indicator variable (see e.g., Schuurman, Houtveen & Hamaker, 2015; Schuurman & Hamaker, in press). The dynamic errors (innovations) and measurement errors can be distinguished because the dynamic errors are carried over through the dynamic process, while the measurement errors are unique to a single time point. The model AR(1) model with measurement error can be specified in JAGS as follows:

MEAR1model <-"
model{
   

#####model######
    for (t in 2:nt)  ###model for timepoint 2 onwards
    {
        y[t]~dnorm(muy[t],Epre)  ###Epre is measurement precision
        muy[t]<- mu + ytilde[t]  ##mu is mean
        ytilde[t] ~ dnorm(muytilde[t], Ipre) ### Ipre is innovation precision
        muytilde[t] <- phi*ytilde[t-1]       ##phi is autoregression coefficient
        
    }

####timepoint 1#######

    ytilde[1] <- y[1]-mu
            
####priors###########
       
    Epre <- inverse(Evar) ##Measurement error precision
    Evar ~dunif(0,300)   ## measurement error variance
    Ipre <- inverse(Ivar) ## innovation precision
    Ivar ~dunif(0,300) ##innovation variance
    
    phi~dunif(-1,1) ##autoregression coefficient
    mu~dnorm(50,.001) ##mean
        
              

}
"
MEAR1model.txt<-textConnection(MEAR1model)

Fit the model for Jane and Mary.

ins = list(list( phi=.6, Ivar = 30, Evar = 20, mu=50),
           list( phi=.3, Ivar = 60, Evar = 40, mu=70),
           list( phi=-.2, Ivar = 20, Evar = 30, mu = 80))  ###initial values ARWN Jags

nt=length(scale)
jags <- jags.model('MEAR1model.txt',
                   data = list('nt' = length(as.numeric(mydata_Jane$scale)),
                               'y' = as.numeric(mydata_Jane$scale)
                   ),
                   inits=ins,
                   n.chains = 3) ###compile and initialize jags ARWN model
## Compiling model graph
##    Resolving undeclared variables
##    Allocating nodes
## Graph information:
##    Observed stochastic nodes: 75
##    Unobserved stochastic nodes: 109
##    Total graph size: 378
## 
## Initializing model
update(jags, 4000) ####burnin

samps=coda.samples(jags,
                   c('phi','mu','Ivar', 'Evar'),
                   4000) ###samples taken next to burnin

##dic.samps=dic.samples(jags,
##                  20000) take samples fort the DIC
plot(samps,trace=FALSE) ###plot densities

#gelman.plot(samps) ##plot gelman rubin statistics
#gelman.diag(samps) ## table gelman rubin statistics
plot(samps,density=FALSE) ###traceplots

summary(samps) ##results parameter estimates
## 
## Iterations = 5001:9000
## Thinning interval = 1 
## Number of chains = 3 
## Sample size per chain = 4000 
## 
## 1. Empirical mean and standard deviation for each variable,
##    plus standard error of the mean:
## 
##          Mean      SD Naive SE Time-series SE
## Evar 112.0188 45.1830 0.412463        3.63664
## Ivar  50.8020 44.9931 0.410729        4.94248
## mu    75.7455  1.8468 0.016859        0.05320
## phi    0.3656  0.2638 0.002408        0.01572
## 
## 2. Quantiles for each variable:
## 
##         2.5%    25%      50%     75%    97.5%
## Evar 15.6824 84.896 115.5195 142.460 195.7378
## Ivar  2.1680 16.677  36.8827  71.311 164.0716
## mu   72.2344 74.505  75.7134  76.946  79.4605
## phi  -0.2075  0.184   0.4063   0.574   0.7607

Be sure to check the convergence of the model; the MEAR model can be tricky to estimate, and needs plenty of iterations. If the model has not converged yet, increase the number of burning and samples to keep.

Compare the results for Jane and Mary’s MEAR(1) model to the results for their AR(1) model. What are the main differences?

Latent AR(1) model or Dynamic Factor Model

Finally, we will model a latent variable, Postive Affect, with an AR(1) process. Both Jane and Mary 8 provided daily evaluations of 8 positive affect items, namely “attentive”, “strong”, “interested”, “enthausiastic”, “proud”, “determined”, “inspired”, and “exstatic”. We will model a common positive affect factor for these items, and see if there is an autoregressive effect for latent Positive Affect.

DFAR1model <-"model{
   

#####model######
    for (t in 2:nt) ###model for timepoint 2 onwards
    {
        y[t,1]~dnorm(muy[t,1],Ipre[1])  ###Ipre is residual/measurement precision
        muy[t,1]<- mu[1] + l[1]*f[t]
        y[t,2]~dnorm(muy[t,2],Ipre[2])  ###Ipre is residual precision
        muy[t,2]<- mu[2] + l[2]*f[t]
        y[t,3]~dnorm(muy[t,3],Ipre[3])  ###Ipre is residual precision
        muy[t,3]<- mu[3] + l[3]*f[t]
        y[t,4]~dnorm(muy[t,4],Ipre[4])  ###Ipre is residual/measurement  precision
        muy[t,4]<- mu[4] + l[4]*f[t]
        y[t,5]~dnorm(muy[t,5],Ipre[5])  ###Ipre is residual/measurement  precision
        muy[t,5]<- mu[5] + l[5]*f[t]
        y[t,6]~dnorm(muy[t,6],Ipre[6])  ###Ipre is residual/measurement  precision
        muy[t,6]<- mu[6] + l[6]*f[t]
        y[t,7]~dnorm(muy[t,7],Ipre[7])  ###Ipre is residual/measurement  precision
        muy[t,7]<- mu[7] + l[7]*f[t]
        y[t,8]~dnorm(muy[t,8],Ipre[8])  ###Ipre is residual/measurement  precision
        muy[t,8]<- mu[8] + l[8]*f[t]

        f[t]~dnorm(muf[t],Fpre)
        muf[t]<- phi*f[t-1]
        
        
    }

   
    #####priors###########
    f[1]~dnorm(0,1)
    Fpre[1] <- inverse(Ivar[1]) ## innovation precision latent variable 
    Fvar[1] ~dunif(0,50)  
    Ipre[1] <- inverse(Ivar[1]) 
    Ivar[1] ~dunif(0,50)   ##residual/measurement variance
    Ipre[2] <- inverse(Ivar[2])
    Ivar[2] ~dunif(0,50)   ##residual/measurement variance
    Ipre[3] <- inverse(Ivar[3]) 
    Ivar[3] ~dunif(0,50)  ###residual/measurement variance
    Ipre[4] <- inverse(Ivar[4]) 
    Ivar[4] ~dunif(0,50)   ##residual/measurement variance
    Ipre[5] <- inverse(Ivar[5]) 
    Ivar[5] ~dunif(0,50)   ##residual/measurement variance
    Ipre[6] <- inverse(Ivar[6]) 
    Ivar[6] ~dunif(0,50)  #residual/measurement variance
    Ipre[7] <- inverse(Ivar[7]) #
    Ivar[7] ~dunif(0,50)   ##residual/measurement variance
    Ipre[8] <- inverse(Ivar[8]) 
    Ivar[8] ~dunif(0,50)   ##residual/measurement variance

    phi~dunif(-1,1) ##autoregression coefficient
    mu[1]~dnorm(0,.001) ##mean
    mu[2]~dnorm(0,.001) ##mean
    mu[3]~dnorm(0,.001) ##mean
    mu[4]~dnorm(0,.001) ##mean
    mu[5]~dnorm(0,.001) ##mean
    mu[6]~dnorm(0,.001) ##mean
    mu[7]~dnorm(0,.001) ##mean
    mu[8]~dnorm(0,.001) ##mean
    l[1]<-1
    l[2]~dnorm(0,.001) ##loading
    l[3]~dnorm(0,.001) ##loading
    l[4]~dnorm(0,.001)  ##loading
    l[5]~dnorm(0,.001) ##loading
    l[6]~dnorm(0,.001) ##loading
    l[7]~dnorm(0,.001) ##loading
    l[8]~dnorm(0,.001) ##loading
              

}
"
DFAR1model.txt<-textConnection(DFAR1model)

Load the data for Jane and Mary, and look at the time series for the different items. Fit the model for Jane and Mary, and draw a path model for both of them based on the results.

datser=data.frame(mydata_Jane$aandachtig, mydata_Jane$sterk,mydata_Jane$geinteresseerd,mydata_Jane$enthousiast, mydata_Jane$trots,mydata_Jane$vastberaden,mydata_Jane$geinspireerd, mydata_Jane$uitgelaten)
datser=as.matrix(datser)
colnames(datser) = c("attentive", "strong", "interested", "enthausiastic", "proud", "determined", "inspired", "exstatic")

ts.plot(datser[,1],col="mediumseagreen",lwd=5, ylab= "scores", ylim=c(1,6))
lines(datser[,2],col="pink",lwd=5)
lines(datser[,3],col="blue",lwd=5)
lines(datser[,4],col="green",lwd=5)
lines(datser[,5],col="orange",lwd=5)
lines(datser[,6],col="red",lwd=5)
lines(datser[,7],col="purple",lwd=5)
lines(datser[,8],col="yellow",lwd=5)

jags <- jags.model('DFAR1model.txt',
                   data = list('nt' = length(datser[,1]),
                               'y' = datser
                               
                               
                   ),
                   n.chains = 3) ###compile and initialize jags AR model
## Compiling model graph
##    Resolving undeclared variables
##    Allocating nodes
## Graph information:
##    Observed stochastic nodes: 597
##    Unobserved stochastic nodes: 239
##    Total graph size: 2437
## 
## Initializing model
update(jags, 5000) ####burnin

samps=coda.samples(jags,
                   c('phi','l', 'mu','Ivar'),
                   10000) ### take samples next to burnin

#dic.samps=dic.samples(jags,
#                     20000) ##take samples for DIC


plot(samps,trace=FALSE) ##plot densities

#gelman.plot(samps) ###plot gelman rubin stats

plot(samps,density=FALSE) ##traceplots

summary(samps) ##results parameters estimates 
## 
## Iterations = 6001:16000
## Thinning interval = 1 
## Number of chains = 3 
## Sample size per chain = 10000 
## 
## 1. Empirical mean and standard deviation for each variable,
##    plus standard error of the mean:
## 
##           Mean      SD  Naive SE Time-series SE
## Ivar[1] 0.3539 0.05839 0.0003371      0.0008399
## Ivar[2] 0.2586 0.06538 0.0003775      0.0010154
## Ivar[3] 0.1496 0.03437 0.0001985      0.0004750
## Ivar[4] 0.2574 0.05581 0.0003222      0.0006872
## Ivar[5] 0.5778 0.11905 0.0006873      0.0013835
## Ivar[6] 0.3960 0.06967 0.0004022      0.0005860
## Ivar[7] 0.6742 0.11955 0.0006902      0.0010282
## Ivar[8] 0.7695 0.13963 0.0008061      0.0012054
## l[1]    1.0000 0.00000 0.0000000      0.0000000
## l[2]    0.9878 0.17925 0.0010349      0.0026836
## l[3]    0.6655 0.12738 0.0007355      0.0016753
## l[4]    0.8218 0.16399 0.0009468      0.0020959
## l[5]    1.0077 0.23452 0.0013540      0.0028704
## l[6]    0.3894 0.16989 0.0009808      0.0014569
## l[7]    0.4635 0.22178 0.0012805      0.0019125
## l[8]    0.5408 0.23548 0.0013596      0.0019373
## mu[1]   3.4330 0.11870 0.0006853      0.0029147
## mu[2]   3.5852 0.11187 0.0006459      0.0028272
## mu[3]   3.8658 0.07925 0.0004575      0.0019282
## mu[4]   4.4388 0.09859 0.0005692      0.0023717
## mu[5]   3.5414 0.13183 0.0007611      0.0028907
## mu[6]   3.9061 0.08297 0.0004790      0.0011754
## mu[7]   3.5053 0.10656 0.0006152      0.0014256
## mu[8]   4.0510 0.11587 0.0006690      0.0015938
## phi     0.2709 0.15376 0.0008877      0.0018677
## 
## 2. Quantiles for each variable:
## 
##             2.5%    25%    50%    75%  97.5%
## Ivar[1]  0.25602 0.3122 0.3480 0.3898 0.4822
## Ivar[2]  0.14900 0.2131 0.2519 0.2977 0.4048
## Ivar[3]  0.09085 0.1259 0.1467 0.1703 0.2262
## Ivar[4]  0.16460 0.2178 0.2522 0.2908 0.3831
## Ivar[5]  0.37956 0.4934 0.5658 0.6487 0.8443
## Ivar[6]  0.28200 0.3468 0.3881 0.4370 0.5528
## Ivar[7]  0.47817 0.5898 0.6611 0.7439 0.9460
## Ivar[8]  0.54376 0.6700 0.7539 0.8507 1.0905
## l[1]     1.00000 1.0000 1.0000 1.0000 1.0000
## l[2]     0.64651 0.8671 0.9843 1.1031 1.3542
## l[3]     0.42549 0.5781 0.6608 0.7488 0.9270
## l[4]     0.51108 0.7099 0.8184 0.9288 1.1536
## l[5]     0.55889 0.8509 1.0040 1.1587 1.4837
## l[6]     0.06348 0.2741 0.3862 0.5008 0.7320
## l[7]     0.03306 0.3146 0.4619 0.6087 0.9089
## l[8]     0.08849 0.3836 0.5379 0.6949 1.0137
## mu[1]    3.19897 3.3555 3.4325 3.5112 3.6690
## mu[2]    3.36504 3.5119 3.5849 3.6577 3.8104
## mu[3]    3.71084 3.8144 3.8655 3.9166 4.0233
## mu[4]    4.24413 4.3746 4.4387 4.5025 4.6343
## mu[5]    3.27972 3.4546 3.5403 3.6283 3.8029
## mu[6]    3.74215 3.8508 3.9057 3.9617 4.0690
## mu[7]    3.29524 3.4344 3.5053 3.5765 3.7143
## mu[8]    3.82420 3.9741 4.0507 4.1280 4.2808
## phi     -0.02990 0.1675 0.2710 0.3750 0.5723

Compare the models for Mary and Jane. Would you say based on these results that Positive Affect is measurement invariant for Jane and Mary (that the latent variable represent the same thing for Jane and Mary)? If not, how do they differ?

End of n=1 practical.

This is the end of the n=1 exercises. I hope you’ve been able to get an impression of n=1 modeling of intensive longitudinal data, and the individual differences in dynamic processes. In the multilevel extensions of these models, these individual differences are explicitly modeled. You can practice with such a model in the multilevel practical.