In JoKra1/papss: What the package does (short line)

knitr::opts_chunk$set(
  collapse = TRUE,
  comment = "#>",
  fig.align = "center",
  fig.width = 6,
  fig.height = 4
)

library(papss)
set.seed(2022) # For replicability.

Simulate pupil data

Here we again start by simulating some pupil data. We only look at 1 subject here since we want to collect the history of coefficient updates this time, which quickly results in very large matrices.

n <- 1
sim_obj <- additive_pupil_sim(n_sub = n,
                              slope_dev = 1.5,
                              sub_dev = 0.15,
                              trial_dev = 0.25,
                              residual_dev=15.0,
                              should_plot=F)
dat <- sim_obj$data

Aggregate

We again aggregate the raw data, as is explained further in the artificial_data_analysis vignette.

aggr_dat <- aggregate(list("pupil"=dat$pupil),by=list("subject"=dat$subject,"time"=dat$time),FUN=mean)

aggr_dat <- aggr_dat[order(aggr_dat$subject),]

plot(aggr_dat$time[aggr_dat$subject == 1],
     aggr_dat$pupil[aggr_dat$subject == 1],
     type="l",ylim=c(min(aggr_dat$pupil),
                     max(aggr_dat$pupil)),
     xlab="Time",
     ylab="Pupil dilation",
     main= "Average pupil trajectories for all subjects",
     lwd=3)

Investigating the identifiability of the problems

How will the solver deal with different starting values?

iter <- 4 # More then 4 takes to long for vignette building
hist_collection <- list()

for(i in 1:iter){
  # Pick a new seed
  set.seed(i)

  # Solve the pupil spline
  solvedPupil <- papss::pupil_solve(pulse_spacing = 1,
                           data = aggr_dat,
                           maxiter_inner = 100000,
                           maxiter_outer = 30,
                           model="WIER_SHARED",
                           expand_by = 0,
                           drop_last = 300,
                           f=(1/(10^(24))),
                           should_collect_progress = T)

  # Coefficient history
  hist = solvedPupil$coefHistory
  hist_collection[[i]] <- hist
}

Coefficient changes over time

Evidently, the time it takes the algorithm to converge differs for different starting values. However, especially towards the end the trajectories become extremely similar, already suggesting that the algorithm ultimately approaches a similar solution independent of the exact starting value.

par(mfrow=c(1,2))
for(hi in 1:iter){
  chist <- hist_collection[[hi]]
  plot(1:ncol(chist),
       chist[1,],
       type="l",
       lwd=3,
       ylim=c(-1.75,1.25),
       xlab="Iterations",
       ylab="Coefficient value")

  for(ri in 2:nrow(chist)){
    lines(1:ncol(chist), chist[ri,],lwd=3,col=ri)
  }
}

This becomes especially evident if we only look at the final estimates.

par(mfrow=c(1,1))
# Visualize final coefficient set
chist <- hist_collection[[1]]
plot(1:nrow(chist),
       chist[,ncol(chist)],
       ylim=c(-3,3),
       xlab="Coefficient",
       ylab="Coefficient value",type="l",lwd=3)



for(hi in 2:iter){
  chist <- hist_collection[[hi]]
  lines(1:nrow(chist), chist[,ncol(chist)],
       col=hi,
       ylim=c(-3,3),
       xlab="Coefficient",
       ylab="Coefficient value",type="l",lwd=3,lty=3)
}


chist <- hist_collection[[1]]
# Visualize last 20 coefficient sets
heatmap(chist[,(ncol(chist) - 19):ncol(chist)],
        Colv = NA,
        Rowv = NA,
        scale = 'none',
        ylab="Coef",
        xlab="Iteration",
        main= "Run: 1")


for(hi in 2:iter){
  chist <- hist_collection[[hi]]
  heatmap(chist[,(ncol(chist) - 19):ncol(chist)],
          Colv = NA,
          Rowv = NA,
          scale = 'none',
          ylab="Coef",
          xlab="Iteration",main=paste0("Run: ",hi))


}