Module 7: Modeling Trajectories and Time Series

library("markovchain")
mcStats <- new("markovchain",
state = c("bored", "horrified", "sleeping", "somewhat awake"),
transitionMatrix = matrix(c(0.60, 0.05, 0.25, 0.10,
0.15, 0.30, 0.10, 0.45,
0.10, 0.10, 0.50, 0.30,
0.30, 0.10, 0.50, 0.10),
byrow = TRUE, nrow = 4),
name = "StatsClass")
print(mcStats)

               bored horrified sleeping somewhat awake
bored           0.60      0.05     0.25           0.10
horrified       0.15      0.30     0.10           0.45
sleeping        0.10      0.10     0.50           0.30
somewhat awake  0.30      0.10     0.50           0.10

#The rows above all sum to 1.

library("qgraph")
A <- mcStats@transitionMatrix
rowSums(A)

         bored      horrified       sleeping somewhat awake 
             1              1              1              1 

qgraph(A, edge.labels = TRUE, edge.color = "black")

#Here we are setting the state as “horrified”

pi <- c(0, 1, 0, 0)

#After 10 minutes

pi * mcStats^2 

     bored horrified sleeping somewhat awake
[1,]  0.28      0.15     0.34           0.22

After 15 minutes

pi * mcStats^3 

     bored horrified sleeping somewhat awake
[1,]  0.29      0.12     0.37           0.22

#Long-run

p <- steadyStates(mcStats)
p

     bored horrified sleeping somewhat awake
[1,]   0.3      0.11     0.38           0.21

library("MPsychoR")

package 㤼㸱MPsychoR㤼㸲 was built under R version 4.0.3

library("depmixS4")

package 㤼㸱depmixS4㤼㸲 was built under R version 4.0.5Loading required package: nnet
Loading required package: MASS
package 㤼㸱MASS㤼㸲 was built under R version 4.0.5Loading required package: Rsolnp
package 㤼㸱Rsolnp㤼㸲 was built under R version 4.0.4Loading required package: nlme

data("iatfaces")
p1dat <- subset(iatfaces, id == 1)

set.seed(123)
p1obj1 <- depmix(log(latency) ~ 1, data = p1dat, nstates = 1)
p1fit1 <- fit(p1obj1)

converged at iteration 1 with logLik: -26 

p1obj2 <- depmix(log(latency) ~ 1, data = p1dat, nstates = 2)
p1fit2 <- fit(p1obj2)

converged at iteration 23 with logLik: 1.5 

p1obj3 <- depmix(log(latency) ~ 1, data = p1dat, nstates = 3)
p1fit3 <- fit(p1obj3)

converged at iteration 192 with logLik: 8.7 

p1obj4 <- depmix(log(latency) ~ 1, data = p1dat, nstates = 4)
p1fit4 <- fit(p1obj4)

converged at iteration 107 with logLik: 12 

c(BIC(p1fit1), BIC(p1fit2), BIC(p1fit3), BIC(p1fit4))

[1] 60 28 44 77

#If the participant is in state 1, he/she tends to stay there (a11 = 0.975). There is a low probability for switching to state 2 (a12 = 0.025). Once switched to state 2, he/she remains there (a22 = 1), and there is virtually no chance to go back to state 1 (a21 = 0).

summary(p1fit2)

Initial state probabilities model 
pr1 pr2 
  1   0 

Transition matrix 
       toS1  toS2
fromS1 0.97 0.025
fromS2 0.00 1.000

Response parameters 
Resp 1 : gaussian 
    Re1.(Intercept) Re1.sd
St1             6.5   0.17
St2             6.9   0.29

#Time Series was really difficult for me. I was able to run one model. #In this section we are interested whether a specific event, occurred at a known point in time, has an influence on the time series.

aca <- c(rep(0, 39), rep(1, 69))

library("MPsychoR")
library(tseries)
library(ggplot2)

package 㤼㸱ggplot2㤼㸲 was built under R version 4.0.5Need help? Try Stackoverflow:
https://stackoverflow.com/tags/ggplot2

data("ageiat")
yts <- ts(ageiat, start = c(2007, 1), frequency = 12)
plot(yts, ylab = "d-measure", main = "Age IAT Time Series")

preaca <- window(yts, end = c(2010, 3)) ## pre-event series
tspre <- arima(preaca, order = c(0,1,1), include.drift = TRUE)

Error in arima(preaca, order = c(0, 1, 1), include.drift = TRUE) : 
  unused argument (include.drift = TRUE)

library("MPsychoR")
library("fda.usc")

package 㤼㸱fda.usc㤼㸲 was built under R version 4.0.5Loading required package: fda
package 㤼㸱fda㤼㸲 was built under R version 4.0.3Loading required package: splines
Loading required package: Matrix
package 㤼㸱Matrix㤼㸲 was built under R version 4.0.5Loading required package: fds
package 㤼㸱fds㤼㸲 was built under R version 4.0.3Loading required package: rainbow
package 㤼㸱rainbow㤼㸲 was built under R version 4.0.3Loading required package: pcaPP
package 㤼㸱pcaPP㤼㸲 was built under R version 4.0.3Loading required package: RCurl
package 㤼㸱RCurl㤼㸲 was built under R version 4.0.3
Attaching package: 㤼㸱fda㤼㸲

The following object is masked from 㤼㸱package:graphics㤼㸲:

    matplot

Loading required package: mgcv
This is mgcv 1.8-31. For overview type 'help("mgcv-package")'.

Attaching package: 㤼㸱mgcv㤼㸲

The following object is masked from 㤼㸱package:nnet㤼㸲:

    multinom

----------------------------------------------------------------------------------
 Functional Data Analysis and Utilities for Statistical Computing
 fda.usc version 2.0.2 (built on 2020-02-17) is now loaded
 fda.usc is running sequentially usign foreach package
 Please, execute ops.fda.usc() once to run in local parallel mode
 Deprecated functions: min.basis, min.np, anova.hetero, anova.onefactor, anova.RPm
 New functions: optim.basis, optim.np, fanova.hetero, fanova.onefactor, fanova.RPm
----------------------------------------------------------------------------------

data("tension")
tension1 <- as.matrix(tension[,1:800]) ## tension time series
cond <- tension$cond ## condition
ftension <- fdata(tension1,
argvals = seq(1, 80, length.out = 800),
names = list(main = "Music tension", xlab = "Time (sec)",
ylab = "Tension"))