1 Load data

data = read.csv('data.csv')
data$Resampling9ms.Axel = as.numeric(data$Resampling9ms.Axel)
data$Resampling10ms.Axel = as.numeric(data$Resampling10ms.Axel)
data_aggr = data %>% group_by(freq) %>%
                     summarize(Baseline.AxelM = mean(Baseline.Axel), ciBaseline.Axel.lower = ci(Baseline.Axel)["CI lower"], ciBaseline.Axel.upper = ci(Baseline.Axel)["CI upper"],
                               Resampling0ms.AxelM = mean(Resampling0ms.Axel), ciResampling0ms.Axel.lower = ci(Resampling0ms.Axel)["CI lower"], ciResampling0ms.Axel.upper = ci(Resampling0ms.Axel)["CI upper"],
                               MA.AxelM = mean(MA.Axel), ciMA.Axel.lower = ci(MA.Axel)["CI lower"], ciMA.Axel.upper = ci(MA.Axel)["CI upper"],
                               Resampling1ms.AxelM = mean(Resampling1ms.Axel), ciResampling1ms.Axel.lower = ci(Resampling1ms.Axel)["CI lower"], ciResampling1ms.Axel.upper = ci(Resampling1ms.Axel)["CI upper"],
                               Resampling2ms.AxelM = mean(Resampling2ms.Axel), ciResampling2ms.Axel.lower = ci(Resampling2ms.Axel)["CI lower"], ciResampling2ms.Axel.upper = ci(Resampling2ms.Axel)["CI upper"],
                               Resampling3ms.AxelM = mean(Resampling3ms.Axel), ciResampling3ms.Axel.lower = ci(Resampling3ms.Axel)["CI lower"], ciResampling3ms.Axel.upper = ci(Resampling3ms.Axel)["CI upper"],
                               Resampling4ms.AxelM = mean(Resampling4ms.Axel), ciResampling4ms.Axel.lower = ci(Resampling4ms.Axel)["CI lower"], ciResampling4ms.Axel.upper = ci(Resampling4ms.Axel)["CI upper"],
                               Resampling5ms.AxelM = mean(Resampling5ms.Axel), ciResampling5ms.Axel.lower = ci(Resampling5ms.Axel)["CI lower"], ciResampling5ms.Axel.upper = ci(Resampling5ms.Axel)["CI upper"],
                               Resampling6ms.AxelM = mean(Resampling6ms.Axel), ciResampling6ms.Axel.lower = ci(Resampling6ms.Axel)["CI lower"], ciResampling6ms.Axel.upper = ci(Resampling6ms.Axel)["CI upper"],
                               Resampling7ms.AxelM = mean(Resampling7ms.Axel), ciResampling7ms.Axel.lower = ci(Resampling7ms.Axel)["CI lower"], ciResampling7ms.Axel.upper = ci(Resampling7ms.Axel)["CI upper"],
                               Resampling8ms.AxelM = mean(Resampling8ms.Axel), ciResampling8ms.Axel.lower = ci(Resampling8ms.Axel)["CI lower"], ciResampling8ms.Axel.upper = ci(Resampling8ms.Axel)["CI upper"],
                               Resampling9ms.AxelM = mean(Resampling9ms.Axel), ciResampling9ms.Axel.lower = ci(Resampling9ms.Axel)["CI lower"], ciResampling9ms.Axel.upper = ci(Resampling9ms.Axel)["CI upper"],
                               Resampling10ms.AxelM = mean(Resampling10ms.Axel), ciResampling10ms.Axel.lower = ci(Resampling10ms.Axel)["CI lower"], ciResampling10ms.Axel.upper = ci(Resampling10ms.Axel)["CI upper"],
                               
                               TimeToFT.baseline.M = mean(timeToFT.baseline), TimeToTF.baseline.ci.lower = ci(timeToFT.baseline)["CI lower"], TimeToTF.baseline.ci.upper = ci(timeToFT.baseline)["CI upper"],
                               TimeToFT.MA.M = mean(timeToFT.MA), TimeToFT.MA.ci.lower = ci(timeToFT.MA)["CI lower"], TimeToFT.MA.ci.upper = ci(timeToFT.MA)["CI upper"],
                               TimeToFT.resamp0ms.M = mean(timeToFT.resamp.delay.0ms), TimeToFT.resamp0ms.ci.lower = ci(timeToFT.resamp.delay.0ms)["CI lower"], TimeToFT.resamp0ms.ci.upper = ci(timeToFT.resamp.delay.0ms)["CI upper"],
                               TimeToFT.resamp1ms.M = mean(timeToFT.resamp.delay.1ms), TimeToFT.resamp1ms.ci.lower = ci(timeToFT.resamp.delay.1ms)["CI lower"], TimeToFT.resamp1ms.ci.upper = ci(timeToFT.resamp.delay.1ms)["CI upper"],
                               TimeToFT.resamp2ms.M = mean(timeToFT.resamp.delay.2ms), TimeToFT.resamp2ms.ci.lower = ci(timeToFT.resamp.delay.2ms)["CI lower"], TimeToFT.resamp2ms.ci.upper = ci(timeToFT.resamp.delay.2ms)["CI upper"],
                               TimeToFT.resamp3ms.M = mean(timeToFT.resamp.delay.3ms), TimeToFT.resamp3ms.ci.lower = ci(timeToFT.resamp.delay.3ms)["CI lower"], TimeToFT.resamp3ms.ci.upper = ci(timeToFT.resamp.delay.3ms)["CI upper"],
                               TimeToFT.resamp4ms.M = mean(timeToFT.resamp.delay.4ms), TimeToFT.resamp4ms.ci.lower = ci(timeToFT.resamp.delay.4ms)["CI lower"], TimeToFT.resamp4ms.ci.upper = ci(timeToFT.resamp.delay.4ms)["CI upper"],
                               TimeToFT.resamp5ms.M = mean(timeToFT.resamp.delay.5ms), TimeToFT.resamp5ms.ci.lower = ci(timeToFT.resamp.delay.5ms)["CI lower"], TimeToFT.resamp5ms.ci.upper = ci(timeToFT.resamp.delay.5ms)["CI upper"],
                               TimeToFT.resamp6ms.M = mean(timeToFT.resamp.delay.6ms), TimeToFT.resamp6ms.ci.lower = ci(timeToFT.resamp.delay.6ms)["CI lower"], TimeToFT.resamp6ms.ci.upper = ci(timeToFT.resamp.delay.6ms)["CI upper"],
                               TimeToFT.resamp7ms.M = mean(timeToFT.resamp.delay.7ms), TimeToFT.resamp7ms.ci.lower = ci(timeToFT.resamp.delay.7ms)["CI lower"], TimeToFT.resamp7ms.ci.upper = ci(timeToFT.resamp.delay.7ms)["CI upper"],
                               TimeToFT.resamp8ms.M = mean(timeToFT.resamp.delay.8ms), TimeToFT.resamp8ms.ci.lower = ci(timeToFT.resamp.delay.8ms)["CI lower"], TimeToFT.resamp8ms.ci.upper = ci(timeToFT.resamp.delay.8ms)["CI upper"],
                               TimeToFT.resamp9ms.M = mean(timeToFT.resamp.delay.9ms), TimeToFT.resamp9ms.ci.lower = ci(timeToFT.resamp.delay.9ms)["CI lower"], TimeToFT.resamp9ms.ci.upper = ci(timeToFT.resamp.delay.9ms)["CI upper"],
                               TimeToFT.resamp10ms.M = mean(timeToFT.resamp.delay.10ms), TimeToFT.resamp10ms.ci.lower = ci(timeToFT.resamp.delay.10ms)["CI lower"], TimeToFT.resamp10ms.ci.upper = ci(timeToFT.resamp.delay.10ms)["CI upper"], .groups = 'drop'
                               )
data2 = read.csv('dataComparisonPixel3and4.csv')
data2$Resampling9ms.Axel = as.numeric(data2$Resampling9ms.Axel)
data2$Resampling10ms.Axel = as.numeric(data2$Resampling10ms.Axel)
data_aggr2 = data2 %>% group_by(freq) %>%
                     summarize(Baseline.AxelM = mean(Baseline.Axel), ciBaseline.Axel.lower = ci(Baseline.Axel)["CI lower"], ciBaseline.Axel.upper = ci(Baseline.Axel)["CI upper"],
                               Resampling0ms.AxelM = mean(Resampling0ms.Axel), ciResampling0ms.Axel.lower = ci(Resampling0ms.Axel)["CI lower"], ciResampling0ms.Axel.upper = ci(Resampling0ms.Axel)["CI upper"],
                               Resampling1ms.AxelM = mean(Resampling1ms.Axel), ciResampling1ms.Axel.lower = ci(Resampling1ms.Axel)["CI lower"], ciResampling1ms.Axel.upper = ci(Resampling1ms.Axel)["CI upper"],
                               Resampling2ms.AxelM = mean(Resampling2ms.Axel), ciResampling2ms.Axel.lower = ci(Resampling2ms.Axel)["CI lower"], ciResampling2ms.Axel.upper = ci(Resampling2ms.Axel)["CI upper"],
                               Resampling3ms.AxelM = mean(Resampling3ms.Axel), ciResampling3ms.Axel.lower = ci(Resampling3ms.Axel)["CI lower"], ciResampling3ms.Axel.upper = ci(Resampling3ms.Axel)["CI upper"],
                               Resampling4ms.AxelM = mean(Resampling4ms.Axel), ciResampling4ms.Axel.lower = ci(Resampling4ms.Axel)["CI lower"], ciResampling4ms.Axel.upper = ci(Resampling4ms.Axel)["CI upper"],
                               Resampling5ms.AxelM = mean(Resampling5ms.Axel), ciResampling5ms.Axel.lower = ci(Resampling5ms.Axel)["CI lower"], ciResampling5ms.Axel.upper = ci(Resampling5ms.Axel)["CI upper"],
                               Resampling6ms.AxelM = mean(Resampling6ms.Axel), ciResampling6ms.Axel.lower = ci(Resampling6ms.Axel)["CI lower"], ciResampling6ms.Axel.upper = ci(Resampling6ms.Axel)["CI upper"],
                               Resampling7ms.AxelM = mean(Resampling7ms.Axel), ciResampling7ms.Axel.lower = ci(Resampling7ms.Axel)["CI lower"], ciResampling7ms.Axel.upper = ci(Resampling7ms.Axel)["CI upper"],
                               Resampling8ms.AxelM = mean(Resampling8ms.Axel), ciResampling8ms.Axel.lower = ci(Resampling8ms.Axel)["CI lower"], ciResampling8ms.Axel.upper = ci(Resampling8ms.Axel)["CI upper"],
                               Resampling9ms.AxelM = mean(Resampling9ms.Axel), ciResampling9ms.Axel.lower = ci(Resampling9ms.Axel)["CI lower"], ciResampling9ms.Axel.upper = ci(Resampling9ms.Axel)["CI upper"],
                               Resampling10ms.AxelM = mean(Resampling10ms.Axel), ciResampling10ms.Axel.lower = ci(Resampling10ms.Axel)["CI lower"], ciResampling10ms.Axel.upper = ci(Resampling10ms.Axel)["CI upper"],
                               .groups = 'drop'
                               )

2 Impact of the display rate on the spatial jitter (Figure 7)

t <- list(
  family = "sans serif",
  size = 20,
  color = 'black')

t2 <- list(
  family = "sans serif",
  size = 30,
  color = 'black')

p <- plot_ly(data = data_aggr, type = "scatter", mode = 'lines+markers', x = ~freq, y = ~Baseline.AxelM, name = 'Simulated',
             error_y = list(symmetric=F, array=~ciBaseline.Axel.upper-Baseline.AxelM, arrayminus=~Baseline.AxelM-ciBaseline.Axel.lower, color = '#000000')) %>%

             add_trace(data = data_aggr2, y=~Baseline.AxelM, name = 'Real',
                       error_y = list(symmetric=F, array=~ciBaseline.Axel.upper-Baseline.AxelM, arrayminus=~Baseline.AxelM-ciBaseline.Axel.lower, color = '#000000')) %>%  
        
             layout(yaxis = list(title = "Spatial jitter (pixels)", range = c(0,35)),
                    xaxis = list(title = 'Output frequency (Hz)', type = "category"),
                    legend = list(x = 0.8, y = 0.95, bordercolor = "#000", borderwidth=2), 
                    font = t, title = list(font = t2)
                    )

p
## Warning: `arrange_()` is deprecated as of dplyr 0.7.0.
## Please use `arrange()` instead.
## See vignette('programming') for more help
## This warning is displayed once every 8 hours.
## Call `lifecycle::last_warnings()` to see where this warning was generated.
orca(p, "ImpactOfTheDisplayRateOnJitter.pdf")
htmlwidgets::saveWidget(p, file = 'ImpactOfTheDisplayRateOnJitter.html')

3 Plot Linear resampling - jitter (Figure 8)

p <- plot_ly(data = data_aggr, type = "scatter", mode = 'lines+markers', x = ~freq, y = ~Baseline.AxelM, name = 'Baseline',
             error_y = list(symmetric=F, array=~ciBaseline.Axel.upper-Baseline.AxelM, arrayminus=~Baseline.AxelM-ciBaseline.Axel.lower, color = '#000000')) %>%
  
             add_trace(data = data_aggr, y=~Resampling0ms.AxelM, name = '0 ms',
                       error_y = list(symmetric=F, array=~ciResampling0ms.Axel.upper-Resampling0ms.AxelM, arrayminus=~Resampling0ms.AxelM-ciResampling0ms.Axel.lower, color = '#000000')) %>%
             add_trace(data = data_aggr, y=~Resampling1ms.AxelM, name = '1 ms',
                       error_y = list(symmetric=F, array=~ciResampling1ms.Axel.upper-Resampling1ms.AxelM, arrayminus=~Resampling1ms.AxelM-ciResampling1ms.Axel.lower, color = '#000000')) %>%  
             add_trace(data = data_aggr, y=~Resampling2ms.AxelM, name = '2 ms',
                       error_y = list(symmetric=F, array=~ciResampling2ms.Axel.upper-Resampling2ms.AxelM, arrayminus=~Resampling2ms.AxelM-ciResampling2ms.Axel.lower, color = '#000000')) %>%  
             add_trace(data = data_aggr, y=~Resampling3ms.AxelM, name = '3 ms',
                       error_y = list(symmetric=F, array=~ciResampling3ms.Axel.upper-Resampling3ms.AxelM, arrayminus=~Resampling3ms.AxelM-ciResampling3ms.Axel.lower, color = '#000000')) %>%  
             add_trace(data = data_aggr, y=~Resampling4ms.AxelM, name = '4 ms',
                       error_y = list(symmetric=F, array=~ciResampling4ms.Axel.upper-Resampling4ms.AxelM, arrayminus=~Resampling4ms.AxelM-ciResampling4ms.Axel.lower, color = '#000000')) %>%  
             add_trace(data = data_aggr, y=~Resampling5ms.AxelM, name = '5 ms',
                       error_y = list(symmetric=F, array=~ciResampling5ms.Axel.upper-Resampling5ms.AxelM, arrayminus=~Resampling5ms.AxelM-ciResampling5ms.Axel.lower, color = '#000000')) %>%  
             add_trace(data = data_aggr, y=~Resampling6ms.AxelM, name = '6 ms',
                       error_y = list(symmetric=F, array=~ciResampling6ms.Axel.upper-Resampling6ms.AxelM, arrayminus=~Resampling6ms.AxelM-ciResampling6ms.Axel.lower, color = '#000000')) %>%  
             add_trace(data = data_aggr, y=~Resampling7ms.AxelM, name = '7 ms',
                       error_y = list(symmetric=F, array=~ciResampling7ms.Axel.upper-Resampling7ms.AxelM, arrayminus=~Resampling7ms.AxelM-ciResampling7ms.Axel.lower, color = '#000000')) %>%  
             add_trace(data = data_aggr, y=~Resampling8ms.AxelM, name = '8 ms',
                       error_y = list(symmetric=F, array=~ciResampling8ms.Axel.upper-Resampling8ms.AxelM, arrayminus=~Resampling8ms.AxelM-ciResampling8ms.Axel.lower, color = '#000000')) %>%  
             add_trace(data = data_aggr, y=~Resampling9ms.AxelM, name = '9 ms',
                       error_y = list(symmetric=F, array=~ciResampling9ms.Axel.upper-Resampling9ms.AxelM, arrayminus=~Resampling9ms.AxelM-ciResampling9ms.Axel.lower, color = '#000000')) %>%  
             add_trace(data = data_aggr, y=~Resampling10ms.AxelM, name = '10 ms', line=list(color='#2C5A39'), marker=list(color='#2C5A39'),
                       error_y = list(symmetric=F, array=~ciResampling10ms.Axel.upper-Resampling10ms.AxelM, arrayminus=~Resampling10ms.AxelM-ciResampling10ms.Axel.lower, color = '#000000')) %>%    

             layout(yaxis = list(title = "Spatial jitter (pixels)", range = c(0,35)),
                    xaxis = list(title = 'Output frequency (Hz)', type = "category"),
                    legend = list(x = 1, y = 0.95, bordercolor = "#000", borderwidth=2), 
                    font = t, title = list(font = t2)
                    )


p
orca(p, "LinearResamplingJitterVSfreq.pdf")
#htmlwidgets::saveWidget(p, file = 'JitterByfrequencyAxel.html')

4 Plot jitter latency tradeoff - linear resampling (Figure 9)

## `summarise()` regrouping output by 'freq' (override with `.groups` argument)
data_aggrBaseline = data_aggr5 %>% filter(technique == 'Baseline')
df0 = r.data_aggrBaseline
print(df0)
##    freq technique  JitterScore  timeToFrame
## 0    30  Baseline    13.260125     7.669513
## 1    40  Baseline    13.390629     6.588503
## 2    50  Baseline    27.186446     5.773825
## 3    60  Baseline    12.132154     4.898809
## 4    70  Baseline    22.462792     4.626359
## 5    80  Baseline    28.723137     4.640046
## 6    90  Baseline    26.042060     4.228488
## 7   100  Baseline    20.096991     4.650800
## 8   110  Baseline    13.109486     4.093190
## 9   120  Baseline     9.269057     4.036021
baseline = {}
for i, row in df0.iterrows():
  baseline[row['freq']] = {}
  baseline[row['freq']]['timeToFrame'] = row['timeToFrame']
#print(baseline)
import pandas
df = r.data_aggr5
for i, row in df.iterrows():
  df.at[i,"timeToFrame"] = row["timeToFrame"] - baseline[row['freq']]['timeToFrame']
data_aggr7 = py$df
data_aggrLatency = data_aggr7 %>% select(-JitterScore) %>% spread(technique, timeToFrame)
data_aggrJitter = data_aggr7 %>% select(-timeToFrame) %>% spread(technique, JitterScore)
t <- list(
  family = "sans serif",
  size = 15,
  color = 'black')
t2 <- list(
  family = "sans serif",
  size = 30,
  color = 'black')
symbols = c('circle', 'triangle-up', 'square', 'cross', 'diamond', 'star-triangle-down', 'hourglass', 'star', 'bowtie', 'x')
p <- plot_ly(type = "scatter", symbols = symbols, mode = 'markers', symbol = data_aggrLatency$freq, marker = list(size = 9, line = list(color = 'rgb(50,50,50)',width = 1))) %>%
  
    # Comment this line to disable Frequencies legend
    add_trace(x = data_aggrLatency$Baseline, y = data_aggrJitter$Baseline, name = paste0(data_aggrLatency$freq," Hz"), marker = list(color = 'rgba(0,0,0,0.5)'), legendgroup = "freq") %>%
    
    add_trace(x = data_aggrLatency$Baseline, y = data_aggrJitter$Baseline, mode = 'markers', name = 'Baseline', marker = list(color = 'rgb(0,0,0)'),legendgroup = "lat") %>%
    add_trace(x = data_aggrLatency$MA, y=data_aggrJitter$MA, name = 'MA', marker = list(color = 'rgb(255,255,255)'),legendgroup = "lat") %>%
    add_trace(x = data_aggrLatency$Resampling0ms, y=data_aggrJitter$Resampling0ms, name = '0 ms', marker = list(color = 'rgb(69,117,180)'),legendgroup = "lat") %>%
    add_trace(x = data_aggrLatency$Resampling2ms, y=data_aggrJitter$Resampling2ms, name = '2 ms', marker = list(color = 'rgb(215,48,39)'),legendgroup = "lat") %>%
    add_trace(x = data_aggrLatency$Resampling4ms, y=data_aggrJitter$Resampling4ms, name = '4 ms', marker = list(color = 'rgb(252,141,89)'),legendgroup = "lat") %>%
    add_trace(x = data_aggrLatency$Resampling6ms, y=data_aggrJitter$Resampling6ms, name = '6 ms', marker = list(color = 'rgb(255,255,191)'),legendgroup = "lat") %>%
    add_trace(x = data_aggrLatency$Resampling8ms, y=data_aggrJitter$Resampling8ms, name = '8 ms', marker = list(color = 'rgb(145,207,96)'),legendgroup = "lat") %>%
    add_trace(x = data_aggrLatency$Resampling10ms, y=data_aggrJitter$Resampling10ms, name = '10 ms', marker = list(color = 'rgb(26,152,80)'),legendgroup = "lat") %>%
  
    layout(yaxis = list(title = "Spatial jitter (pixels)"),
          xaxis = list(title = 'Latency compared to baseline (ms)'),
          legend = list(x = 1, y = 0.5, bordercolor = "#000", borderwidth=2, orientation = 'v', font = t, title = list(font = t2)
          ))
p
## Warning: The shape palette can deal with a maximum of 6 discrete values because
## more than 6 becomes difficult to discriminate; you have 10. Consider
## specifying shapes manually if you must have them.
orca(p, "JitterLatencyTradeoffComparedToBaseline.pdf")
## Warning: The shape palette can deal with a maximum of 6 discrete values because
## more than 6 becomes difficult to discriminate; you have 10. Consider
## specifying shapes manually if you must have them.
#htmlwidgets::saveWidget(p, file = 'ImpactOfTheDisplayRateOnJitter.html')

5 Test of other extrapolation techniques

5.1 Linear predictor

data = read.csv('dataLinear.csv')
data$Resampling9ms.Axel = as.numeric(data$Resampling9ms.Axel)
data$Resampling10ms.Axel = as.numeric(data$Resampling10ms.Axel)
data_aggrLinear = data %>% group_by(freq) %>%
                     summarize(Baseline.AxelM = mean(Baseline.Axel), ciBaseline.Axel.lower = ci(Baseline.Axel)["CI lower"], ciBaseline.Axel.upper = ci(Baseline.Axel)["CI upper"],
                               Resampling0ms.AxelM = mean(Resampling0ms.Axel), ciResampling0ms.Axel.lower = ci(Resampling0ms.Axel)["CI lower"], ciResampling0ms.Axel.upper = ci(Resampling0ms.Axel)["CI upper"],
                               Resampling1ms.AxelM = mean(Resampling1ms.Axel), ciResampling1ms.Axel.lower = ci(Resampling1ms.Axel)["CI lower"], ciResampling1ms.Axel.upper = ci(Resampling1ms.Axel)["CI upper"],
                               Resampling2ms.AxelM = mean(Resampling2ms.Axel), ciResampling2ms.Axel.lower = ci(Resampling2ms.Axel)["CI lower"], ciResampling2ms.Axel.upper = ci(Resampling2ms.Axel)["CI upper"],
                               Resampling3ms.AxelM = mean(Resampling3ms.Axel), ciResampling3ms.Axel.lower = ci(Resampling3ms.Axel)["CI lower"], ciResampling3ms.Axel.upper = ci(Resampling3ms.Axel)["CI upper"],
                               Resampling4ms.AxelM = mean(Resampling4ms.Axel), ciResampling4ms.Axel.lower = ci(Resampling4ms.Axel)["CI lower"], ciResampling4ms.Axel.upper = ci(Resampling4ms.Axel)["CI upper"],
                               Resampling5ms.AxelM = mean(Resampling5ms.Axel), ciResampling5ms.Axel.lower = ci(Resampling5ms.Axel)["CI lower"], ciResampling5ms.Axel.upper = ci(Resampling5ms.Axel)["CI upper"],
                               Resampling6ms.AxelM = mean(Resampling6ms.Axel), ciResampling6ms.Axel.lower = ci(Resampling6ms.Axel)["CI lower"], ciResampling6ms.Axel.upper = ci(Resampling6ms.Axel)["CI upper"],
                               Resampling7ms.AxelM = mean(Resampling7ms.Axel), ciResampling7ms.Axel.lower = ci(Resampling7ms.Axel)["CI lower"], ciResampling7ms.Axel.upper = ci(Resampling7ms.Axel)["CI upper"],
                               Resampling8ms.AxelM = mean(Resampling8ms.Axel), ciResampling8ms.Axel.lower = ci(Resampling8ms.Axel)["CI lower"], ciResampling8ms.Axel.upper = ci(Resampling8ms.Axel)["CI upper"],
                               Resampling9ms.AxelM = mean(Resampling9ms.Axel), ciResampling9ms.Axel.lower = ci(Resampling9ms.Axel)["CI lower"], ciResampling9ms.Axel.upper = ci(Resampling9ms.Axel)["CI upper"],
                               Resampling10ms.AxelM = mean(Resampling10ms.Axel), ciResampling10ms.Axel.lower = ci(Resampling10ms.Axel)["CI lower"], ciResampling10ms.Axel.upper = ci(Resampling10ms.Axel)["CI upper"],
                               .groups = 'drop'
                               )

5.2 Quadratic predictor

data = read.csv('dataQuadratic.csv')
data$Resampling9ms.Axel = as.numeric(data$Resampling9ms.Axel)
data$Resampling10ms.Axel = as.numeric(data$Resampling10ms.Axel)
data_aggrQuadratic = data %>% group_by(freq) %>%
                     summarize(Baseline.AxelM = mean(Baseline.Axel), ciBaseline.Axel.lower = ci(Baseline.Axel)["CI lower"], ciBaseline.Axel.upper = ci(Baseline.Axel)["CI upper"],
                               Resampling0ms.AxelM = mean(Resampling0ms.Axel), ciResampling0ms.Axel.lower = ci(Resampling0ms.Axel)["CI lower"], ciResampling0ms.Axel.upper = ci(Resampling0ms.Axel)["CI upper"],
                               Resampling1ms.AxelM = mean(Resampling1ms.Axel), ciResampling1ms.Axel.lower = ci(Resampling1ms.Axel)["CI lower"], ciResampling1ms.Axel.upper = ci(Resampling1ms.Axel)["CI upper"],
                               Resampling2ms.AxelM = mean(Resampling2ms.Axel), ciResampling2ms.Axel.lower = ci(Resampling2ms.Axel)["CI lower"], ciResampling2ms.Axel.upper = ci(Resampling2ms.Axel)["CI upper"],
                               Resampling3ms.AxelM = mean(Resampling3ms.Axel), ciResampling3ms.Axel.lower = ci(Resampling3ms.Axel)["CI lower"], ciResampling3ms.Axel.upper = ci(Resampling3ms.Axel)["CI upper"],
                               Resampling4ms.AxelM = mean(Resampling4ms.Axel), ciResampling4ms.Axel.lower = ci(Resampling4ms.Axel)["CI lower"], ciResampling4ms.Axel.upper = ci(Resampling4ms.Axel)["CI upper"],
                               Resampling5ms.AxelM = mean(Resampling5ms.Axel), ciResampling5ms.Axel.lower = ci(Resampling5ms.Axel)["CI lower"], ciResampling5ms.Axel.upper = ci(Resampling5ms.Axel)["CI upper"],
                               Resampling6ms.AxelM = mean(Resampling6ms.Axel), ciResampling6ms.Axel.lower = ci(Resampling6ms.Axel)["CI lower"], ciResampling6ms.Axel.upper = ci(Resampling6ms.Axel)["CI upper"],
                               Resampling7ms.AxelM = mean(Resampling7ms.Axel), ciResampling7ms.Axel.lower = ci(Resampling7ms.Axel)["CI lower"], ciResampling7ms.Axel.upper = ci(Resampling7ms.Axel)["CI upper"],
                               Resampling8ms.AxelM = mean(Resampling8ms.Axel), ciResampling8ms.Axel.lower = ci(Resampling8ms.Axel)["CI lower"], ciResampling8ms.Axel.upper = ci(Resampling8ms.Axel)["CI upper"],
                               Resampling9ms.AxelM = mean(Resampling9ms.Axel), ciResampling9ms.Axel.lower = ci(Resampling9ms.Axel)["CI lower"], ciResampling9ms.Axel.upper = ci(Resampling9ms.Axel)["CI upper"],
                               Resampling10ms.AxelM = mean(Resampling10ms.Axel), ciResampling10ms.Axel.lower = ci(Resampling10ms.Axel)["CI lower"], ciResampling10ms.Axel.upper = ci(Resampling10ms.Axel)["CI upper"],
                               .groups = 'drop'
                               )

5.3 LaViola predictor

data = read.csv('dataLaViola.csv')
data$Resampling9ms.Axel = as.numeric(data$Resampling9ms.Axel)
data$Resampling10ms.Axel = as.numeric(data$Resampling10ms.Axel)
data_aggrLaViola = data %>% group_by(freq) %>%
                     summarize(Baseline.AxelM = mean(Baseline.Axel), ciBaseline.Axel.lower = ci(Baseline.Axel)["CI lower"], ciBaseline.Axel.upper = ci(Baseline.Axel)["CI upper"],
                               Resampling0ms.AxelM = mean(Resampling0ms.Axel), ciResampling0ms.Axel.lower = ci(Resampling0ms.Axel)["CI lower"], ciResampling0ms.Axel.upper = ci(Resampling0ms.Axel)["CI upper"],
                               Resampling1ms.AxelM = mean(Resampling1ms.Axel), ciResampling1ms.Axel.lower = ci(Resampling1ms.Axel)["CI lower"], ciResampling1ms.Axel.upper = ci(Resampling1ms.Axel)["CI upper"],
                               Resampling2ms.AxelM = mean(Resampling2ms.Axel), ciResampling2ms.Axel.lower = ci(Resampling2ms.Axel)["CI lower"], ciResampling2ms.Axel.upper = ci(Resampling2ms.Axel)["CI upper"],
                               Resampling3ms.AxelM = mean(Resampling3ms.Axel), ciResampling3ms.Axel.lower = ci(Resampling3ms.Axel)["CI lower"], ciResampling3ms.Axel.upper = ci(Resampling3ms.Axel)["CI upper"],
                               Resampling4ms.AxelM = mean(Resampling4ms.Axel), ciResampling4ms.Axel.lower = ci(Resampling4ms.Axel)["CI lower"], ciResampling4ms.Axel.upper = ci(Resampling4ms.Axel)["CI upper"],
                               Resampling5ms.AxelM = mean(Resampling5ms.Axel), ciResampling5ms.Axel.lower = ci(Resampling5ms.Axel)["CI lower"], ciResampling5ms.Axel.upper = ci(Resampling5ms.Axel)["CI upper"],
                               Resampling6ms.AxelM = mean(Resampling6ms.Axel), ciResampling6ms.Axel.lower = ci(Resampling6ms.Axel)["CI lower"], ciResampling6ms.Axel.upper = ci(Resampling6ms.Axel)["CI upper"],
                               Resampling7ms.AxelM = mean(Resampling7ms.Axel), ciResampling7ms.Axel.lower = ci(Resampling7ms.Axel)["CI lower"], ciResampling7ms.Axel.upper = ci(Resampling7ms.Axel)["CI upper"],
                               Resampling8ms.AxelM = mean(Resampling8ms.Axel), ciResampling8ms.Axel.lower = ci(Resampling8ms.Axel)["CI lower"], ciResampling8ms.Axel.upper = ci(Resampling8ms.Axel)["CI upper"],
                               Resampling9ms.AxelM = mean(Resampling9ms.Axel), ciResampling9ms.Axel.lower = ci(Resampling9ms.Axel)["CI lower"], ciResampling9ms.Axel.upper = ci(Resampling9ms.Axel)["CI upper"],
                               Resampling10ms.AxelM = mean(Resampling10ms.Axel), ciResampling10ms.Axel.lower = ci(Resampling10ms.Axel)["CI lower"], ciResampling10ms.Axel.upper = ci(Resampling10ms.Axel)["CI upper"],
                               .groups = 'drop'
                               )

5.4 Plot Linear resampling - jitter (Figure 10)

p <- plot_ly(data = data_aggrLinear, type = "scatter", mode = 'lines+markers', x = ~freq, y = ~Baseline.AxelM, name = 'Baseline', line=list(color='#FF0000'), marker=list(color='#FF0000'), 
             error_y = list(symmetric=F, array=~ciBaseline.Axel.upper-Baseline.AxelM, arrayminus=~Baseline.AxelM-ciBaseline.Axel.lower, color = '#000000')) %>%
  
             add_trace(data = data_aggrLinear, y=~Resampling0ms.AxelM, name = 'Linear 0 ms', line=list(color='#A63603'), marker=list(color='#A63603'), 
                       error_y = list(symmetric=F, array=~ciResampling0ms.Axel.upper-Resampling0ms.AxelM, arrayminus=~Resampling0ms.AxelM-ciResampling0ms.Axel.lower, color = '#000000')) %>%
             add_trace(data = data_aggrLinear, y=~Resampling2ms.AxelM, name = 'Linear 2 ms', line=list(color='#E6550D'), marker=list(color='#E6550D'), 
                       error_y = list(symmetric=F, array=~ciResampling2ms.Axel.upper-Resampling2ms.AxelM, arrayminus=~Resampling2ms.AxelM-ciResampling2ms.Axel.lower, color = '#000000')) %>%
             add_trace(data = data_aggrLinear, y=~Resampling4ms.AxelM, name = 'Linear 4 ms', line=list(color='#FDBE85'), marker=list(color='#FC8D59'), 
                       error_y = list(symmetric=F, array=~ciResampling4ms.Axel.upper-Resampling4ms.AxelM, arrayminus=~Resampling4ms.AxelM-ciResampling4ms.Axel.lower, color = '#000000')) %>%

             add_trace(data = data_aggrQuadratic, y=~Resampling0ms.AxelM, name = 'Curve 0 ms', line=list(color='#006D2C'), marker=list(color='#006D2C'),
                       error_y = list(symmetric=F, array=~ciResampling0ms.Axel.upper-Resampling0ms.AxelM, arrayminus=~Resampling0ms.AxelM-ciResampling0ms.Axel.lower, color = '#000000')) %>%
             add_trace(data = data_aggrQuadratic, y=~Resampling2ms.AxelM, name = 'Curve 2 ms', line=list(color='#31A354'), marker=list(color='#31A354'),
                       error_y = list(symmetric=F, array=~ciResampling2ms.Axel.upper-Resampling2ms.AxelM, arrayminus=~Resampling2ms.AxelM-ciResampling2ms.Axel.lower, color = '#000000')) %>%
             add_trace(data = data_aggrQuadratic, y=~Resampling4ms.AxelM, name = 'Curve 4 ms', line=list(color='#BAE4B3'), marker=list(color='#BAE4B3'),
                       error_y = list(symmetric=F, array=~ciResampling4ms.Axel.upper-Resampling4ms.AxelM, arrayminus=~Resampling4ms.AxelM-ciResampling4ms.Axel.lower, color = '#000000')) %>%

             add_trace(data = data_aggrLaViola, y=~Resampling0ms.AxelM, name = 'DESP 0 ms', line=list(color='#08519C'), marker=list(color='#08519C'),
                       error_y = list(symmetric=F, array=~ciResampling0ms.Axel.upper-Resampling0ms.AxelM, arrayminus=~Resampling0ms.AxelM-ciResampling0ms.Axel.lower, color = '#000000')) %>%
             add_trace(data = data_aggrLaViola, y=~Resampling2ms.AxelM, name = 'DESP 2 ms', line=list(color='#3182BD'), marker=list(color='#3182BD'),
                       error_y = list(symmetric=F, array=~ciResampling2ms.Axel.upper-Resampling2ms.AxelM, arrayminus=~Resampling2ms.AxelM-ciResampling2ms.Axel.lower, color = '#000000')) %>%
             add_trace(data = data_aggrLaViola, y=~Resampling4ms.AxelM, name = 'DESP 4 ms', line=list(color='#BDD7E7'), marker=list(color='#BDD7E7'),
                       error_y = list(symmetric=F, array=~ciResampling4ms.Axel.upper-Resampling4ms.AxelM, arrayminus=~Resampling4ms.AxelM-ciResampling4ms.Axel.lower, color = '#000000')) %>%
  
             layout(yaxis = list(title = "Spatial jitter (pixels)", range = c(0,35)),
                    xaxis = list(title = 'Output frequency (Hz)', type = "category"),
                    legend = list(x = 1, y = 0.95, bordercolor = "#000", borderwidth=2), 
                    font = t, title = list(font = t2)
                    )


p
orca(p, "ComparisonPredictors.pdf")
#htmlwidgets::saveWidget(p, file = 'JitterByfrequencyAxel.html')