R notebook for analyses of median saccade latency in the sacc-tDCS dataset. Previous processing:
The .csv file with the eye tracking data was created in MATLAB.
anodal or cathodal sessionpre), during (tDCS), or after (post) tDCSlateral - fixation in center of display, saccade made towards the peripherycenter - fixation in periphery, saccade made back towards the center of the displayleft for saccades towards the left of current fixation position; right for saccades to the rightkable(head(subjectData))| subject | session.order | gender | age | dominant.eye |
|---|---|---|---|---|
| S01 | first.cathodal | female | 23 | right |
| S02 | first.anodal | male | 27 | left |
| S03 | first.cathodal | male | 34 | right |
| S04 | first.anodal | male | NA | right |
| S05 | first.cathodal | male | 24 | right |
| S06 | first.cathodal | female | 29 | left |
first.anodal) or cathodal stimulation in the first session (first.cathodal)The main use is to see if the nuisance factor session.order covaries with the factors of interest in the design. This could indicate the presence of carryover effects between the stimulation, or a difference in subgroups within the sample (see http://www.jerrydallal.com/lhsp/crossovr.htm for an introduction to these kinds of analyses.).
tooFast <- 50
tooSlow <- 400
badFix <- 1.8
badSacc <- 8
subs2exclude <- c("S28","S16","S22","S21","S25")Criteria for outlier saccades:
In Kanai et al. (2012), this was:
# Remove outliers and subjects
groupData <- filter(groupData,
# outliers
latency >= tooFast,
latency <= tooSlow,
deviation.start <= badFix,
deviation.end.x <= badSacc,
# subjects
!(subject %in% subs2exclude),
# missing values
complete.cases(groupData)
)# Compute median in each condition
latencyMedian <- groupData %>%
group_by(subject,stimulation,leg,block,type,direction) %>%
summarise(latency = median(latency))# Plot out all the data
fullPlot <- ggplot(latencyMedian, aes(interaction(block,leg), latency, color = stimulation, shape = stimulation)) +
facet_grid(type ~ direction) +
stat_summary(fun.y = mean, geom = "point", size = 3) +
stat_summary(fun.y = mean, geom = "line", aes(group = stimulation), size = 1) +
theme(axis.text.x = element_text(angle = 22, vjust = .5))
fullPlot
For some reason cathodal is always a little faster than anodal. There is probably a few ms difference between them on average, also in the baseline block. This is clearly just random variation, but it could be a problem, since the effects we expect are not much bigger…
The above plot also shows a lot of variability, so it might be best to average over each 3 consecutive blocks so the data come in 15-minute intervals.
# Compute median per leg
latencyMedianLeg <- groupData %>%
group_by(subject,stimulation,direction,type) %>%
summarise(baseline = median(latency[leg == "pre"]), # take average of 3 blocks, make new column
tDCS = median(latency[leg == "tDCS"]),
post.1 = median(latency[leg == "post" & block <= 3]),
post.2 = median(latency[leg == "post" & block >= 4])) %>%
gather(leg,latency,baseline,tDCS,post.1,post.2) %>% # gather new columns to use as factor
mutate(leg = factor(leg, levels = c("baseline", "tDCS", "post.1", "post.2"))) # reorder factor levelsNow make the same plot, but for the data collapsed over 3 blocks. Also draw the plot for each individual subject, so we can see which subjects drive the baseline difference, and in which direction the stimulation effect goes for each subject (if there is any).
kanaiSubsPlot <- ggplot(latencyMedianLeg, aes(leg, latency, color = stimulation, shape = type)) +
facet_wrap(~subject, ncol = 5) +
stat_summary(fun.y = mean, geom = "point", size = 2) +
stat_summary(fun.y = mean, geom = "line", aes(group = stimulation)) +
theme(axis.text.x = element_text(angle = 22, vjust = .5))
kanaiSubsPlot
There are a couple of subjects that show large differences in the baseline already, especially S01. The stimulation effects seem to not be very apparent, but they could also drown out at this scale.
Let’s look at the size of the baseline difference per subject.
baselineDiff <- latencyMedianLeg %>%
filter(leg == "baseline") %>% # keep only baseline data
spread(stimulation, latency) %>% # make separate columns for anodal and cathodal
mutate(latency.diff = anodal - cathodal) %>% # subtract the difference
group_by(subject) %>%
summarise(latency.diff = mean(latency.diff))# keep the average difference per subject
kable(baselineDiff, caption = 'Difference between baseline saccade latencies in anodal and cathodal session')| subject | latency.diff |
|---|---|
| S01 | 39.125 |
| S02 | 18.625 |
| S03 | 2.375 |
| S04 | -6.125 |
| S05 | -17.500 |
| S06 | -3.875 |
| S07 | 5.625 |
| S08 | 11.125 |
| S09 | 23.875 |
| S10 | 14.750 |
| S11 | 1.000 |
| S12 | 8.375 |
| S13 | 4.000 |
| S14 | -10.375 |
| S15 | 1.875 |
| S17 | -5.500 |
| S18 | -2.875 |
| S19 | 1.000 |
| S20 | 4.750 |
| S24 | 3.625 |
| S26 | 0.000 |
| S27 | 0.500 |
| S29 | -13.625 |
| S30 | 7.125 |
| S32 | 10.875 |
| S33 | -1.625 |
There are a few subjects with substantial latency differences between sessions, so it’s probably a good idea to subtract the baseline for the statistical analyses. The mean difference is 3.7 ms. So on average, the cathodal session is a little faster, so the counterbalancing is not perfect, but the difference is not very high.
Here we simply extract median RTs for each condition and use a repeated measures ANOVA for statistical analysis, following Kanai et al. (2012).
Let’s look at the group average plot for the first time.
kanaiPlot <- ggplot(latencyMedianLeg, aes(leg, latency, color = stimulation, shape = stimulation)) +
facet_grid(type ~ direction) +
stat_summary(fun.y = mean, geom = "point", size = 3) +
stat_summary(fun.y = mean, geom = "line", aes(group = stimulation), size = 1) +
stat_summary(fun.data = mean_cl_normal, geom = "errorbar", width = 0.3)
kanaiPlot
All differences between anodal & cathodal seem to be < 5 ms. Any differences that are there also 1) do not seem to differ much between the time blocks, or 2) be in the opposite direction for anodal and cathodal, with the possible exception of the first post-block in the center-right condition.
Let’s look at the same plot for individual subjects, to see whether there are any consistent patterns or huge outliers there
kanaiPlotSubsAnodal <- ggplot(latencyMedianLeg[latencyMedianLeg$stimulation == "anodal", ], aes(leg, latency)) +
facet_grid(type ~ direction) +
geom_line(aes(group=subject,color=subject)) +
stat_summary(fun.y = mean, aes(group = stimulation), geom = "line") +
stat_summary(fun.y = mean, geom = "point") +
ggtitle("Anodal session")
kanaiPlotSubsAnodal
There do appear to be some “outliers”, but they are fairly well balanced (below or above the average), and >80% or so of subjects seem to cluster together. I don’t see much evidence for systematic individual differences (e.g. slow subjects get faster; fast subjects get slower.)
kanaiPlotSubsCathodal <- ggplot(latencyMedianLeg[latencyMedianLeg$stimulation == "cathodal", ], aes(leg, latency)) +
facet_grid(type ~ direction) +
geom_line(aes(group=subject,color=subject)) +
stat_summary(fun.y = mean, aes(group = stimulation), geom = "line") +
stat_summary(fun.y = mean, geom = "point") +
ggtitle("Cathodal session")
kanaiPlotSubsCathodal
There seem to be less clear individual outliers, but the spread also seems a bit bigger (e.g. right-lateral)
Baseline differences (while small here) are not informative and could obscure real changes from baseline within subjects. Also it’s hard to see and compare the magnitude of the effects, also for instance because center saccades are much faster than lateral saccades. Therefore, subtract the baseline from each subsequent measurement.
latencyMedianBaseline <- latencyMedianLeg %>%
group_by(subject,stimulation,direction,type) %>% # for each condition, subtract baseline scores and make new columns
summarise(tDCS = latency[leg == "tDCS"] - latency[leg == "baseline"],
post.1 = latency[leg == "post.1"] - latency[leg == "baseline"],
post.2 = latency[leg == "post.2"] - latency[leg == "baseline"]) %>%
gather(leg, latency, tDCS, post.1, post.2) %>% # gather new columns to use as factor
mutate(leg = factor(leg, levels = c("tDCS", "post.1", "post.2"))) # reorder factor levelsBefore we look at the change-from-baseline data, it would be good to check the baseline data itself.
Analyzing change-from-baseline only makes sense if the baseline data are not very noisy. Subtracting a noisy measure from the data will increase the noise in the data, and thereby reduce power. One way to assess this is to look at the correlation between the baseline blocks in both sessions. If there is a high correlation, the pro-saccade task produces comparable latencies each time it is performed, and thus the baseline data should on average not be very noisy.
Let’s compute the correlation between all baselines (i.e. for the factors DIRECTION and TYPE).
baselineCorr <- latencyMedianLeg %>%
filter(leg == "baseline") %>% # keep only baseline data
group_by(direction,type) %>% # for each of these 4 pairs
spread(stimulation,latency) %>% # make 2 columns with baseline latencies of anodal and cathodal session
nest() %>% # create a separate data frame of all remaining columns. These data frames are stored in a list-column named "data"
mutate(stats = map(data, ~cor.test(formula = ~ anodal + cathodal, data =.))) %>% # run correlation test on baselines from each condition
mutate(tidy_model = map(stats, tidy)) %>% # force the four test outputs to also be data frames in a list column ("tidy_model")
unnest(tidy_model, .drop = TRUE) # unpack the list-column with results, drop the list-column we used to organize the dataNow a scatter plot of the data, with the results of the test:
latencyMedianLeg %>%
filter(leg == "baseline") %>%
spread(stimulation,latency) %>%
inner_join(., subjectData[ ,c("subject","session.order")], by = c("subject")) %>% # add column on session order from other data frame
ggplot(aes(anodal,cathodal)) +
facet_grid(type ~ direction) +
geom_abline(intercept = 0, slope = 1, linetype = "dashed") +
geom_smooth(method = "lm") +
geom_point(aes(color=session.order)) +
xlim(80,220) + ylim(80,220) +
geom_text(data = baselineCorr, x = 100, y = 200, aes(label = paste("italic(r) == ", round(estimate,2))), parse = TRUE) +
labs(title = "Baseline in anodal and cathodal sessions", subtitle = "scatterplot of baseline latencies")
The correlations are quite high, at least high enough for the baseline subtraction to help power rather than hurt (r’s > 0.5).
There are a few outliers though: some subjects deviate quite a bit from the idealized line (dashed).
When the data is split for session order, a clear pattern becomes apparent: latencies are generally faster in the 2nd session (for those who got cathodal first, are above the diagonal, meaning anodal < cathodal. The reverse pattern is visible for those who received anodal first). So in that sense there is a difference between two groups of subjects.
If the prosaccade task is indeed reliable, there shouldn’t be any significant differences at the group-level between the baseline blocks of both sessions. If there are, this may be problematic, because all the change-scores are expressed relative to the baseline. A significant difference between two conditions in later blocks could then be driven by the differences that already occur in the baseline.
Let’s plot the baseline scores in each session for each subject:
latencyMedianLeg %>%
filter(leg == "baseline") %>%
ggplot(aes(stimulation, latency)) +
facet_grid(type ~ direction) +
stat_summary(fun.y = mean, geom = "line", aes(group = 1), size = 1.5) +
geom_line(aes(group = subject, colour = subject), alpha = 0.5, position = position_dodge(width = 0.1)) +
labs(title = "Baseline in anodal and cathodal sessions", subtitle = "linked stripchart of baseline latencies")
Most lines appear fairly flat. However, there is quite some variability around the mean, and some subjects show quite a steep difference. There appears to be some sort of regression to the mean also: extreme scores in either session show a steeper slope than those scores that are less extreme.
Let’s look at the pairwise differences in more detail:
latencyMedianLeg %>%
filter(leg == "baseline") %>%
inner_join(., subjectData[ ,c("subject","session.order")], by = c("subject")) %>% # add column on session order from other data frame
group_by(subject,direction,type,session.order) %>%
summarise(latency.diff = latency[stimulation == "anodal"] - latency[stimulation == "cathodal"]) %>%
ggplot(aes(factor(0), latency.diff)) +
facet_grid(type ~ direction) +
geom_hline(yintercept = 0, linetype = "dashed") +
stat_summary(fun.data = mean_cl_normal) +
stat_summary(fun.y = mean, aes(label=round(..y.., digits=2), x = 1.3), geom = "label", alpha = 0.5) +
geom_point(shape = 21, aes(colour = session.order), position = position_jitter(width=.1)) +
labs(title = "Baseline in anodal and cathodal sessions", subtitle = "anodal - cathodal")
The sequence effect we observed earlier is very clear when we split the data by session order again, particularly in the lateral condition (i.e. latencies are always faster in the 2nd session).
Somewhat worryingly, the baseline differences are in the range of the tDCS effect size we expect, particularly in the center condition. That said, the differences do cluster around 0, there are really only a few subjects that fall far outside the boundary (one is really extreme).
The CIs appear to overlap with 0 though, so the differences shouldn’t be significant. Let’s test that:
latencyMedianLeg %>%
filter(leg == "baseline") %>%
group_by(direction,type) %>% # for each of these 4 pairs
nest() %>% # create a separate data frame of all remaining columns. These data frames are stored in a list-column named "data"
mutate(stats = map(data, ~t.test(formula = latency~stimulation, paired = TRUE, data =.))) %>% # run t-test on the data frames
mutate(tidy_model = map(stats, tidy)) %>% # force the four test outputs to also be data frames in a list column ("tidy_model")
unnest(tidy_model, .drop = TRUE) %>% # unpack the list-column with results, drop the list-column we used to organize the data
kable(.)| direction | type | estimate | statistic | p.value | parameter | conf.low | conf.high | method | alternative |
|---|---|---|---|---|---|---|---|---|---|
| left | lateral | 2.403846 | 1.024429 | 0.3154437 | 25 | -2.4289077 | 7.236600 | Paired t-test | two.sided |
| left | center | 5.557692 | 1.910225 | 0.0676402 | 25 | -0.4344207 | 11.549805 | Paired t-test | two.sided |
| right | lateral | 2.769231 | 1.060619 | 0.2990028 | 25 | -2.6081375 | 8.146599 | Paired t-test | two.sided |
| right | center | 4.211538 | 1.657484 | 0.1099194 | 25 | -1.0215893 | 9.444666 | Paired t-test | two.sided |
Indeed, they are not, although the center conditions come close.
kanaiPlotBase <- ggplot(latencyMedianBaseline, aes(leg, latency, color = stimulation, shape = stimulation)) +
facet_grid(type ~ direction) +
geom_hline(yintercept = 0, linetype = "dashed") +
stat_summary(fun.y = mean, geom = "point", size = 3) +
stat_summary(fun.y = mean, geom = "line", aes(group = stimulation), size = 1) +
stat_summary(fun.data = mean_cl_normal, geom = "errorbar", width = 0.3) +
coord_cartesian(ylim = c(-10,10))
kanaiPlotBase
Generally the differences between the stimulation conditions are small (< 4 ms). Except for the center-left condition, but this condition also had the largest baseline difference…
kanaiPlotBaseSubsAnodal <- ggplot(latencyMedianBaseline[latencyMedianBaseline$stimulation == "anodal", ], aes(leg, latency)) +
facet_grid(type ~ direction) +
geom_hline(yintercept = 0, linetype = "dashed") +
geom_line(aes(group=subject,color=subject)) +
stat_summary(fun.y = mean, aes(group = stimulation), geom = "line") +
stat_summary(fun.y = mean, geom = "point") +
ggtitle("Anodal difference from baseline")
kanaiPlotBaseSubsAnodal
kanaiPlotBaseSubsCathodal <- ggplot(latencyMedianBaseline[latencyMedianBaseline$stimulation == "cathodal", ], aes(leg, latency)) +
facet_grid(type ~ direction) +
geom_hline(yintercept = 0, linetype = "dashed") +
geom_line(aes(group=subject,color=subject)) +
stat_summary(fun.y = mean, aes(group = stimulation), geom = "line") +
stat_summary(fun.y = mean, geom = "point") +
ggtitle("Cathodal difference from baseline")
kanaiPlotBaseSubsCathodal
All of this is split pretty much 50-50, hence the average difference hovering around 0.
Kanai et al. (2012) specifically found an effect of anodal tDCS on the latency changes in contralateral (here: left) lateral saccades. Let’s examine the raw numbers here:
latencyMedianBaseline %>%
filter(type == "lateral", stimulation == "anodal", direction == "left") %>%
group_by(leg) %>%
summarise_at(vars(latency), funs(mean, sd)) %>%
kable(.)| leg | mean | sd |
|---|---|---|
| tDCS | -0.1730769 | 5.342176 |
| post.1 | -0.6153846 | 7.134855 |
| post.2 | 0.9615385 | 9.651863 |
These changes are all < 1 ms and highly variable, so it seems there’s truly nothing there…
# Make "subject" a factor, so we can model the repeated measures
latencyMedianBaseline <- latencyMedianBaseline %>%
ungroup() %>% # remove grouping info, because we need to refactor
inner_join(., subjectData[ ,c("subject","session.order")], by = c("subject")) %>% # add column from on session order from other data frame
mutate(subject = factor(subject)) # refactor
latencyMedianLeg <- latencyMedianLeg %>%
ungroup() %>%
inner_join(., subjectData[ ,c("subject","session.order")], by = c("subject")) %>% # add column from on session order from other data frame
mutate(subject = factor(subject))Data: Outlier trials removed, collapsed into 15-minute intervals.
Dependent measure: saccadic latency
Factors:
modelOmni <- ezANOVA(data = data.frame(latencyMedianLeg), # Repeated over subjects; type 3 sums of squares (cf. SPSS)
dv = .(latency), wid = .(subject), within = .(stimulation, leg, type, direction), type = 3)
kable(modelOmni$ANOVA)
| Effect | DFn | DFd | F | p | p<.05 | ges | |
|---|---|---|---|---|---|---|---|
| 2 | stimulation | 1 | 25 | 1.8399816 | 0.1870776 | 0.0046328 | |
| 3 | leg | 3 | 75 | 0.4395422 | 0.7253793 | 0.0004232 | |
| 4 | type | 1 | 25 | 28.5410095 | 0.0000154 | * | 0.1295560 |
| 5 | direction | 1 | 25 | 0.5725232 | 0.4563300 | 0.0015576 | |
| 6 | stimulation:leg | 3 | 75 | 1.4343849 | 0.2394786 | 0.0007140 | |
| 7 | stimulation:type | 1 | 25 | 0.4929953 | 0.4890804 | 0.0003222 | |
| 8 | leg:type | 3 | 75 | 4.7880531 | 0.0041682 | * | 0.0015673 |
| 9 | stimulation:direction | 1 | 25 | 0.6423151 | 0.4304257 | 0.0001620 | |
| 10 | leg:direction | 3 | 75 | 1.6943654 | 0.1754445 | 0.0003035 | |
| 11 | type:direction | 1 | 25 | 0.9424122 | 0.3409576 | 0.0008141 | |
| 12 | stimulation:leg:type | 3 | 75 | 0.6925407 | 0.5594502 | 0.0000998 | |
| 13 | stimulation:leg:direction | 3 | 75 | 2.2130816 | 0.0935459 | 0.0003838 | |
| 14 | stimulation:type:direction | 1 | 25 | 0.3842442 | 0.5409498 | 0.0000109 | |
| 15 | leg:type:direction | 3 | 75 | 6.8295412 | 0.0003940 | * | 0.0004416 |
| 16 | stimulation:leg:type:direction | 3 | 75 | 1.0933932 | 0.3572648 | 0.0000761 |
kable(modelOmni$`Mauchly's Test for Sphericity`)
| Effect | W | p | p<.05 | |
|---|---|---|---|---|
| 3 | leg | 0.2579115 | 0.0000058 | * |
| 6 | stimulation:leg | 0.3762730 | 0.0003186 | * |
| 8 | leg:type | 0.1267739 | 0.0000000 | * |
| 10 | leg:direction | 0.3612309 | 0.0002082 | * |
| 12 | stimulation:leg:type | 0.5303766 | 0.0102889 | * |
| 13 | stimulation:leg:direction | 0.5729125 | 0.0216415 | * |
| 15 | leg:type:direction | 0.7933796 | 0.3595870 | |
| 16 | stimulation:leg:type:direction | 0.8238579 | 0.4676623 |
kable(modelOmni$`Sphericity Corrections`)| Effect | GGe | p[GG] | p[GG]<.05 | HFe | p[HF] | p[HF]<.05 | |
|---|---|---|---|---|---|---|---|
| 3 | leg | 0.5904749 | 0.6231869 | 0.6322248 | 0.6364737 | ||
| 6 | stimulation:leg | 0.5986297 | 0.2486500 | 0.6420287 | 0.2482301 | ||
| 8 | leg:type | 0.4789969 | 0.0234582 | * | 0.5002439 | 0.0218426 | * |
| 10 | leg:direction | 0.6068634 | 0.1971371 | 0.6519485 | 0.1948609 | ||
| 12 | stimulation:leg:type | 0.6922238 | 0.5100514 | 0.7560487 | 0.5218275 | ||
| 13 | stimulation:leg:direction | 0.7236356 | 0.1150365 | 0.7949444 | 0.1090952 | ||
| 15 | leg:type:direction | 0.8657955 | 0.0008106 | * | 0.9750655 | 0.0004504 | * |
| 16 | stimulation:leg:type:direction | 0.9008895 | 0.3541166 | 1.0205912 | 0.3572648 |
latencyMedianLeg %>%
group_by(subject,type) %>%
summarise(latency = mean(latency)) %>%
ggplot(aes(type, latency)) +
stat_summary(fun.data = mean_cl_normal, size = 1) +
geom_jitter(width = 0.25)
This simply reflects that center saccades are faster than lateral saccades, because the location of the target is known.
latencyMedianLeg %>%
group_by(subject,leg,type) %>%
summarise(latency = mean(latency)) %>%
ggplot(aes(leg, latency, shape = type)) +
stat_summary(fun.y = mean, geom = "point") +
stat_summary(fun.y = mean, geom = "line", aes(group = type, linetype = type))
The effect of Type changes over time: lateral saccades steadily become slower; center saccades vary more randomly.
latencyMedianLeg %>%
group_by(subject,leg,type,direction) %>%
summarise(latency = mean(latency)) %>%
ggplot(aes(leg, latency, shape = type)) +
facet_wrap(~direction) +
stat_summary(fun.y = mean, geom = "point") +
stat_summary(fun.y = mean, geom = "line", aes(group = type, linetype = type))
The effect over time looks non-linear for left-center saccades. The “tDCS” block is the deviant here, but there is no interaction with stimulation.
Also, for right-lateral saccades, there is a big change that only emerges in the 2nd post-block.
Differing from the previous omnibus analysis, Kanai et al. (2012) analysed shifts from baseline and only had lateral saccades.
Data:
Dependent measure: saccadic latency
Factors:
modelKanai <- ezANOVA(data = data.frame(filter(latencyMedianBaseline, type == "lateral")),
dv = .(latency), wid = .(subject), within = .(stimulation,leg,direction), type = 3)
# OR, without the EZ package:
# modelKanai=aov(latency~stimulation*leg*direction + Error(subject/(stimulation*leg*direction)),data=latencyMedianBaselineLateral)
# summary(modelKanai)
kable(modelKanai$ANOVA)
| Effect | DFn | DFd | F | p | p<.05 | ges | |
|---|---|---|---|---|---|---|---|
| 2 | stimulation | 1 | 25 | 0.8044339 | 0.3783271 | 0.0089853 | |
| 3 | leg | 2 | 50 | 3.4589912 | 0.0391753 | * | 0.0164093 |
| 4 | direction | 1 | 25 | 0.0864110 | 0.7712172 | 0.0003663 | |
| 5 | stimulation:leg | 2 | 50 | 0.4014143 | 0.6715104 | 0.0007505 | |
| 6 | stimulation:direction | 1 | 25 | 0.5221926 | 0.4766141 | 0.0012979 | |
| 7 | leg:direction | 2 | 50 | 3.6626747 | 0.0327780 | * | 0.0023536 |
| 8 | stimulation:leg:direction | 2 | 50 | 2.5875456 | 0.0852459 | 0.0029246 |
kable(modelKanai$`Mauchly's Test for Sphericity`)
| Effect | W | p | p<.05 | |
|---|---|---|---|---|
| 3 | leg | 0.5213492 | 0.0004032 | * |
| 5 | stimulation:leg | 0.6360361 | 0.0043831 | * |
| 7 | leg:direction | 0.7573179 | 0.0355909 | * |
| 8 | stimulation:leg:direction | 0.9123864 | 0.3327710 |
kable(modelKanai$`Sphericity Corrections`)| Effect | GGe | p[GG] | p[GG]<.05 | HFe | p[HF] | p[HF]<.05 | |
|---|---|---|---|---|---|---|---|
| 3 | leg | 0.6762922 | 0.0597258 | 0.7012858 | 0.0578252 | ||
| 5 | stimulation:leg | 0.7331572 | 0.6087290 | 0.7674993 | 0.6179798 | ||
| 7 | leg:direction | 0.8047111 | 0.0434681 | * | 0.8517313 | 0.0406111 | * |
| 8 | stimulation:leg:direction | 0.9194442 | 0.0903598 | 0.9889658 | 0.0859306 |
latencyMedianBaseline %>%
filter(type == "lateral") %>%
group_by(subject,leg) %>%
summarise(latency = mean(latency)) %>%
ggplot(aes(leg, latency)) +
geom_hline(yintercept = 0, linetype = "dashed") +
stat_summary(fun.data = mean_cl_normal, size = 1) +
stat_summary(fun.y = mean, geom = "line", aes(group = 1)) +
geom_jitter(width = 0.25)
Saccades become slower over time with respect to the baseline. Note that this effect becomes just non-significant when correcting for sphericity.
latencyMedianBaseline %>%
filter(type == "lateral") %>%
group_by(subject,leg,direction) %>%
summarise(latency = mean(latency)) %>%
ggplot(aes(leg, latency, shape = direction)) +
geom_hline(yintercept = 0, linetype = "dashed") +
stat_summary(fun.y = mean, geom = "point") +
stat_summary(fun.y = mean, geom = "line", aes(group = direction, linetype = direction))
This main effect of leg is stronger for right saccades, but does not occur until the 2nd post-block.
This interaction is not significant, but it is of primary interest, so let’s still look at it in more detail.
latencyMedianBaseline %>%
filter(type == "lateral") %>%
group_by(subject,stimulation,leg,direction) %>%
summarise(latency = mean(latency)) %>%
ggplot(aes(leg, latency, shape = stimulation, color = stimulation)) +
geom_hline(yintercept = 0, linetype = "dashed") +
facet_wrap(~direction) +
stat_summary(fun.y = mean, geom = "point") +
stat_summary(fun.y = mean, geom = "line", aes(group = stimulation))
Interestingly, subjects get slower over time in all conditions, but this trend is reversed for left-saccades in the anodal session, as predicted. However, the effect barely reaches 1 ms…
Add an additional factor SESSION ORDER, which creates two groups: those subjects who received anodal tDCS in the first session vs. those who received cathodal tDCS in the first session. Note that these groups are not exactly balanced, which might affect (correcting for) violations of sphericity:
latencyMedianBaseline %>%
group_by(session.order) %>%
summarize(count = n_distinct(subject)) %>%
kable(.)| session.order | count |
|---|---|
| first.anodal | 14 |
| first.cathodal | 12 |
Data:
Dependent measure: saccadic latency
Factors:
modelKanaiOrder <- ezANOVA(data = data.frame(filter(latencyMedianBaseline, type == "lateral")), dv = .(latency),
wid = .(subject), within = .(stimulation,leg,direction), between = session.order, type = 3)
kable(modelKanaiOrder$ANOVA)
| Effect | DFn | DFd | F | p | p<.05 | ges | |
|---|---|---|---|---|---|---|---|
| 2 | session.order | 1 | 24 | 0.1433284 | 0.7083227 | 0.0020996 | |
| 3 | stimulation | 1 | 24 | 0.6340723 | 0.4336719 | 0.0068970 | |
| 5 | leg | 2 | 48 | 3.1888328 | 0.0500837 | 0.0158928 | |
| 7 | direction | 1 | 24 | 0.1014256 | 0.7528804 | 0.0004663 | |
| 4 | session.order:stimulation | 1 | 24 | 2.9615300 | 0.0981361 | 0.0314179 | |
| 6 | session.order:leg | 2 | 48 | 0.9521914 | 0.3930621 | 0.0047991 | |
| 8 | session.order:direction | 1 | 24 | 0.1559521 | 0.6963989 | 0.0007169 | |
| 9 | stimulation:leg | 2 | 48 | 0.2675333 | 0.7663984 | 0.0004736 | |
| 11 | stimulation:direction | 1 | 24 | 0.5243058 | 0.4760104 | 0.0014194 | |
| 13 | leg:direction | 2 | 48 | 3.4028814 | 0.0414916 | * | 0.0023321 |
| 10 | session.order:stimulation:leg | 2 | 48 | 3.6820297 | 0.0325327 | * | 0.0064786 |
| 12 | session.order:stimulation:direction | 1 | 24 | 0.0512017 | 0.8229013 | 0.0001388 | |
| 14 | session.order:leg:direction | 2 | 48 | 0.5664097 | 0.5713066 | 0.0003889 | |
| 15 | stimulation:leg:direction | 2 | 48 | 2.5819394 | 0.0860974 | 0.0030905 | |
| 16 | session.order:stimulation:leg:direction | 2 | 48 | 0.7368106 | 0.4839727 | 0.0008839 |
kable(modelKanaiOrder$`Mauchly's Test for Sphericity`)
| Effect | W | p | p<.05 | |
|---|---|---|---|---|
| 5 | leg | 0.5041662 | 0.0003798 | * |
| 6 | session.order:leg | 0.5041662 | 0.0003798 | * |
| 9 | stimulation:leg | 0.7028372 | 0.0173312 | * |
| 10 | session.order:stimulation:leg | 0.7028372 | 0.0173312 | * |
| 13 | leg:direction | 0.7680908 | 0.0481110 | * |
| 14 | session.order:leg:direction | 0.7680908 | 0.0481110 | * |
| 15 | stimulation:leg:direction | 0.8975070 | 0.2883601 | |
| 16 | session.order:stimulation:leg:direction | 0.8975070 | 0.2883601 |
kable(modelKanaiOrder$`Sphericity Corrections`)| Effect | GGe | p[GG] | p[GG]<.05 | HFe | p[HF] | p[HF]<.05 | |
|---|---|---|---|---|---|---|---|
| 5 | leg | 0.6685235 | 0.0724403 | 0.6933380 | 0.0704904 | ||
| 6 | session.order:leg | 0.6685235 | 0.3628233 | 0.6933380 | 0.3657030 | ||
| 9 | stimulation:leg | 0.7709133 | 0.7088451 | 0.8136384 | 0.7209771 | ||
| 10 | session.order:stimulation:leg | 0.7709133 | 0.0453943 | * | 0.8136384 | 0.0426597 | * |
| 13 | leg:direction | 0.8117481 | 0.0526303 | 0.8622305 | 0.0493806 | * | |
| 14 | session.order:leg:direction | 0.8117481 | 0.5370589 | 0.8622305 | 0.5469276 | ||
| 15 | stimulation:leg:direction | 0.9070352 | 0.0920383 | 0.9770103 | 0.0875336 | ||
| 16 | session.order:stimulation:leg:direction | 0.9070352 | 0.4720869 | 0.9770103 | 0.4811399 |
latencyMedianBaseline %>%
filter(type == "lateral") %>%
group_by(subject,session.order,stimulation,direction) %>%
summarise(latency = mean(latency)) %>%
ggplot(aes(direction, latency, shape = stimulation, color = stimulation)) +
geom_hline(yintercept = 0, linetype = "dashed") +
facet_wrap(~session.order) +
stat_summary(fun.y = mean, geom = "point") +
stat_summary(fun.y = mean, geom = "line", aes(group = stimulation))
The session order by stimulation direction interaction can be construed as a main effect of session. It appears that in the 2nd session, people slow down more than in the 1st. The difference between left and right also seems more pronounced, especially when the 2nd session is anodal.
This latter three-way interaction is very difficult to interpret though. In addition, the stimulation by direction interaction is not significant without session order, so we do not have to worry that this effect is present but confounded by session order.
Repeat the same ANOVA, but now for center saccades (which Kanai did not have).
modelKanaiCenter <- ezANOVA(data = data.frame(filter(latencyMedianBaseline, type == "center")),
dv = .(latency), wid = .(subject), within = .(stimulation,leg,direction), type = 3)
kable(modelKanaiCenter$ANOVA)
| Effect | DFn | DFd | F | p | p<.05 | ges | |
|---|---|---|---|---|---|---|---|
| 2 | stimulation | 1 | 25 | 3.0937191 | 0.0908311 | 0.0230246 | |
| 3 | leg | 2 | 50 | 2.3789819 | 0.1030547 | 0.0062189 | |
| 4 | direction | 1 | 25 | 4.7775590 | 0.0384115 | * | 0.0158632 |
| 5 | stimulation:leg | 2 | 50 | 0.1055936 | 0.8999904 | 0.0001735 | |
| 6 | stimulation:direction | 1 | 25 | 1.8802459 | 0.1824903 | 0.0055307 | |
| 7 | leg:direction | 2 | 50 | 3.1018543 | 0.0537187 | 0.0016853 | |
| 8 | stimulation:leg:direction | 2 | 50 | 1.9636026 | 0.1510262 | 0.0012036 |
kable(modelKanaiCenter$`Mauchly's Test for Sphericity`)
| Effect | W | p | p<.05 | |
|---|---|---|---|---|
| 3 | leg | 0.5173954 | 0.0003680 | * |
| 5 | stimulation:leg | 0.8690648 | 0.1856204 | |
| 7 | leg:direction | 0.9230394 | 0.3825099 | |
| 8 | stimulation:leg:direction | 0.7823535 | 0.0525819 |
kable(modelKanaiCenter$`Sphericity Corrections`)| Effect | GGe | p[GG] | p[GG]<.05 | HFe | p[HF] | p[HF]<.05 | |
|---|---|---|---|---|---|---|---|
| 3 | leg | 0.6744887 | 0.1241893 | 0.6991962 | 0.1225297 | ||
| 5 | stimulation:leg | 0.8842239 | 0.8774896 | 0.9465498 | 0.8902341 | ||
| 7 | leg:direction | 0.9285390 | 0.0579604 | 0.9999608 | 0.0537210 | ||
| 8 | stimulation:leg:direction | 0.8212564 | 0.1600417 | 0.8713552 | 0.1575428 |
latencyMedianBaseline %>%
filter(type == "center") %>%
group_by(subject,direction) %>%
summarise(latency = mean(latency)) %>%
ggplot(aes(direction, latency)) +
geom_hline(yintercept = 0, linetype = "dashed") +
stat_summary(fun.y = mean, geom = "line", aes(group = 1), size = 2) +
geom_line(aes(colour = subject, group = subject))
This seems to be a really tiny and variable effect, but apparently left saccades get somewhat slower and right saccades somewhat faster compared to baseline.
modelKanaiCenterOrder <- ezANOVA(data = data.frame(filter(latencyMedianBaseline, type == "center")), dv = .(latency),
wid = .(subject), within = .(stimulation,leg,direction), between = session.order, type = 3)
kable(modelKanaiCenterOrder$ANOVA)
| Effect | DFn | DFd | F | p | p<.05 | ges | |
|---|---|---|---|---|---|---|---|
| 2 | session.order | 1 | 24 | 0.1166399 | 0.7356800 | 0.0025860 | |
| 3 | stimulation | 1 | 24 | 3.3377548 | 0.0801702 | 0.0258116 | |
| 5 | leg | 2 | 48 | 2.3943806 | 0.1020456 | 0.0066197 | |
| 7 | direction | 1 | 24 | 6.5039074 | 0.0175653 | * | 0.0189442 |
| 4 | session.order:stimulation | 1 | 24 | 0.9506165 | 0.3392880 | 0.0074896 | |
| 6 | session.order:leg | 2 | 48 | 0.5737494 | 0.5672253 | 0.0015943 | |
| 8 | session.order:direction | 1 | 24 | 5.5454061 | 0.0270448 | * | 0.0161976 |
| 9 | stimulation:leg | 2 | 48 | 0.0965518 | 0.9081388 | 0.0001683 | |
| 11 | stimulation:direction | 1 | 24 | 2.4540824 | 0.1303112 | 0.0068754 | |
| 13 | leg:direction | 2 | 48 | 2.8467499 | 0.0678686 | 0.0016019 | |
| 10 | session.order:stimulation:leg | 2 | 48 | 0.5117829 | 0.6026597 | 0.0008915 | |
| 12 | session.order:stimulation:direction | 1 | 24 | 3.2609785 | 0.0835059 | 0.0091154 | |
| 14 | session.order:leg:direction | 2 | 48 | 1.1056981 | 0.3392615 | 0.0006228 | |
| 15 | stimulation:leg:direction | 2 | 48 | 2.2553582 | 0.1158348 | 0.0013961 | |
| 16 | session.order:stimulation:leg:direction | 2 | 48 | 1.7398007 | 0.1864427 | 0.0010773 |
kable(modelKanaiCenterOrder$`Mauchly's Test for Sphericity`)
| Effect | W | p | p<.05 | |
|---|---|---|---|---|
| 5 | leg | 0.5132456 | 0.0004664 | * |
| 6 | session.order:leg | 0.5132456 | 0.0004664 | * |
| 9 | stimulation:leg | 0.8520126 | 0.1585365 | |
| 10 | session.order:stimulation:leg | 0.8520126 | 0.1585365 | |
| 13 | leg:direction | 0.8959726 | 0.2827415 | |
| 14 | session.order:leg:direction | 0.8959726 | 0.2827415 | |
| 15 | stimulation:leg:direction | 0.8087403 | 0.0870572 | |
| 16 | session.order:stimulation:leg:direction | 0.8087403 | 0.0870572 |
kable(modelKanaiCenterOrder$`Sphericity Corrections`)| Effect | GGe | p[GG] | p[GG]<.05 | HFe | p[HF] | p[HF]<.05 | |
|---|---|---|---|---|---|---|---|
| 5 | leg | 0.6726060 | 0.1234676 | 0.6980931 | 0.1217371 | ||
| 6 | session.order:leg | 0.6726060 | 0.5031425 | 0.6980931 | 0.5090948 | ||
| 9 | stimulation:leg | 0.8710897 | 0.8835886 | 0.9334806 | 0.8962764 | ||
| 10 | session.order:stimulation:leg | 0.8710897 | 0.5781127 | 0.9334806 | 0.5903977 | ||
| 13 | leg:direction | 0.9057746 | 0.0738118 | 0.9754789 | 0.0693693 | ||
| 14 | session.order:leg:direction | 0.9057746 | 0.3351508 | 0.9754789 | 0.3382625 | ||
| 15 | stimulation:leg:direction | 0.8394475 | 0.1256368 | 0.8953942 | 0.1221703 | ||
| 16 | session.order:stimulation:leg:direction | 0.8394475 | 0.1923135 | 0.8953942 | 0.1903485 |
There is a main effect of direction in the ANOVA without session order, so this effect might be confounded.
latencyMedianBaseline %>%
filter(type == "center") %>%
group_by(subject,session.order,direction) %>%
summarise(latency = mean(latency)) %>%
ggplot(aes(direction, latency, shape = session.order)) +
geom_hline(yintercept = 0, linetype = "dashed") +
stat_summary(fun.y = mean, geom = "point") +
stat_summary(fun.y = mean, geom = "line", aes(group = session.order, linetype = session.order))
The main effect of direction is present for both groups, but it is much stronger in the group that first receives cathodal stimulation. This could be construed as a carryover effect, i.e. that the effect of cathodal stimulation was not yet washed-out by the second session and modulated the effect of anodal (and vice versa). Or, it could simply be that there is a difference between the two groups (anodal first vs. cathodal first) that has nothing to do with stimulation.
bfKanaiLateral = anovaBF(latency~stimulation*leg*direction+subject, data = data.frame(filter(latencyMedianBaseline, type == "lateral")), whichModels = "withmain", whichRandom = "subject", progress = FALSE, iterations = 100000) # compute Bayes Factors
View(latencyMedianBaseline)
bfKanaiLateral = sort(bfKanaiLateral, decreasing = TRUE) # sort such that winning model is at the topFirst we compare all models to the most simple (null) model, which is the intercept only + random effect model: latency ~ subject. This does not test for effects of SUBJECT but models it as a nuisance factor (whichRandom = "subject"). In addition, to decrease the model space, we do not consider models that have an interaction without the corresponding main effects (whichModels = "withmain").
kable(select(extractBF(bfKanaiLateral), bf)) # show only the Bayes factors in a table| bf | |
|---|---|
| leg + subject | 0.8399089 |
| stimulation + subject | 0.7416393 |
| stimulation + leg + subject | 0.6465591 |
| direction + subject | 0.1357138 |
| direction + leg + subject | 0.1114138 |
| stimulation + direction + subject | 0.0995612 |
| stimulation + direction + leg + subject | 0.0872255 |
| stimulation + leg + stimulation:leg + subject | 0.0478623 |
| stimulation + direction + stimulation:direction + leg + subject | 0.0246398 |
| stimulation + direction + stimulation:direction + subject | 0.0215074 |
| direction + leg + direction:leg + subject | 0.0110546 |
| stimulation + direction + leg + direction:leg + subject | 0.0085640 |
| stimulation + direction + leg + stimulation:leg + subject | 0.0063807 |
| stimulation + direction + stimulation:direction + leg + direction:leg + subject | 0.0018722 |
| stimulation + direction + stimulation:direction + leg + stimulation:leg + subject | 0.0014185 |
| stimulation + direction + leg + stimulation:leg + direction:leg + subject | 0.0006242 |
| stimulation + direction + stimulation:direction + leg + stimulation:leg + direction:leg + subject | 0.0001322 |
| stimulation + direction + stimulation:direction + leg + stimulation:leg + direction:leg + stimulation:direction:leg + subject | 0.0000286 |
The winning model is the one with just a main effect of LEG, with a Bayes factor of 0.8. However, so even for this model there is less evidence than for the null model, which by definition has a Bayes Factor of 1.
The conventional interpretation for Bayes Factors is the following:
0.333 < BF10 < 1: anecdotal evidence for H0
BF10 = 1: equivalent evidence for H0 and H1
BF10 > 10: strong evidence H1
So to switch between expressing the Bayes Factor as evidence for H1 (BF10) or H0 (BF01), you simply invert it (divide by 1).
We can compute the evidence for a particular effect by comparing this winning model with the best-fitting model that does not contain the effect. We can compute the evidence for absence of a particular effect by comparing the winning model with the best-fitting model that does contain the effect.
The evidence for the effect of LEG can be quantified by comparing the Bayes factors of the first and the 3rd model (because that is the first one that does not contain LEG as a factor): 1.3. This constitues only anocdotal evidence for the presence of a leg effect, even though it was significant in the classical ANOVA. But, given that the leg model itself is such a poor fit, in the Bayesian analysis this weak evidence is no surprise.
The evidence for the absence of the STIMULATION by DIRECTION interaction can be quantified by comparing the Bayes factors of the first and the 9th model (because that is the first one that does contain the interaction): 34.1. This constitues strong evidence for the absence of the interaction.
Another option for quantifying evidence for a particular effect is to compare the full model to a model where that effect is omitted (whichModels = top"). The full model is stimulation + direction + stimulation:direction + leg + stimulation:leg + direction:leg + stimulation:direction:leg + subject.
bfKanaiLateralFull = anovaBF(latency~stimulation*leg*direction+subject, data = data.frame(filter(latencyMedianBaseline, type == "lateral")), whichModels = "top", whichRandom = "subject", progress = FALSE, iterations = 100000) # compute Bayes Factors
bfKanaiLateralFullBayes factor top-down analysis
--------------
When effect is omitted from stimulation + direction + leg + stimulation:direction + stimulation:leg + direction:leg + stimulation:direction:leg + subject , BF is...
[1] Omit direction:leg:stimulation : 5.384841 <U+00B1>5.62%
[2] Omit direction:leg : 9.833237 <U+00B1>5.51%
[3] Omit leg:stimulation : 13.49779 <U+00B1>5.66%
[4] Omit direction:stimulation : 4.282142 <U+00B1>5.45%
[5] Omit leg : 1.191349 <U+00B1>7.91%
[6] Omit direction : 7.358471 <U+00B1>5.63%
[7] Omit stimulation : 1.295222 <U+00B1>8.68%
Against denominator:
latency ~ stimulation + direction + leg + stimulation:direction + stimulation:leg + direction:leg + stimulation:direction:leg + subject
---
Bayes factor type: BFlinearModel, JZSRemoving the LEG effect from the model yields a higher Bayes Factor, so removing this effect actually improved the model, although only by a little bit. The evidence thus goes in the direction of the null; when expressed in favor of the alternative, the Bayes Factir becomes less than one: 1 1.191 = 0.8
Similarly, removing the DIRECTION by STIMULATION effect improves the model, and this time a bit more: in the range of moderate evidence for the null.
The problem with the 1st option (comparing single models, for example to the winning model) is that for designs with many factors (and therefore models), it becomes risky and a bit subjective to base conclusions on just two models. In this case, the winning model is even a bad fit, so it doesn’t seem appropriate to use it as a benchmark.
The problem with the 2nd option (comparing against the full model) is actually very apparent in this dataset: the full model is a terrible fit, as it comes in last. There is even a lot of evidence against it when compared to the null model!
One solution is to do Bayesian Model Averaging: compare multiple models and aggregate the Bayes Factors.
In the JASP stats package, this analysis is also called the “inclusion Bayes factor”. Briefly, it compares all models that include the effect of interest vs. all models that do not. For examples and a conceptual explanation, see example 5 in: Bayesian inference for psychology. Part II: Example applications with JASP.
There’s also a different way of calculating inclusion Bayes factors (conceptualized by Sebastiaan Mathot) that is currently being implemented in JASP. This is called the “inclusion Bayes factor across matched models”, and is more selective in the set of models that’s being compared than the standard inclusion BF. Briefly, the inclusion BF across matched models compares:
Let’s compare the default inclusion Bayes Factors for all effects:
# Default inclusion Bayes Factors
kable(inclusionBF(bfKanaiLateral))| effect | Bayes.factor |
|---|---|
| stimulation | 0.2873386 |
| direction | 0.0556146 |
| stimulation:direction | 0.0287605 |
| leg | 0.3194814 |
| stimulation:leg | 0.0327913 |
| direction:leg | 0.0128231 |
| stimulation:direction:leg | 0.0001360 |
Doing the analysis this way, we also find strong evidence against most of these effects; particularly the interactions.
My preferred approach is the inclusion Bayes factor across matched models. The evidence is qualitatively similar, but less strong, probably because we aren’t including as many poorly fitting models in the model comparison (which have a very large BF10 and therefore a lot of effect on the composite Bayes Factor):
kable(inclusionBF(bfKanaiLateral, models = "matched"))| effect | Bayes.factor |
|---|---|
| stimulation | 0.7547571 |
| direction | 0.1344014 |
| stimulation:direction | 0.2449654 |
| leg | 0.8555487 |
| stimulation:leg | 0.0733784 |
| direction:leg | 0.0962754 |
| stimulation:direction:leg | 0.2163735 |
bfKanaiCenter = anovaBF(latency~stimulation*leg*direction+subject, data = data.frame(filter(latencyMedianBaseline, type == "center")), whichModels="withmain", whichRandom = "subject", progress = FALSE, iterations = 100000) # compute Bayes Factors
bfKanaiCenter <- sort(bfKanaiCenter, decreasing = TRUE) # sort such that winning model is at the topkable(select(extractBF(bfKanaiCenter), bf)) # show only the Bayes factors in a table| bf | |
|---|---|
| stimulation + direction + subject | 555.2527534 |
| stimulation + direction + stimulation:direction + subject | 401.8383357 |
| stimulation + direction + leg + subject | 94.0231519 |
| stimulation + direction + stimulation:direction + leg + subject | 71.1333686 |
| stimulation + subject | 56.4457669 |
| stimulation + leg + subject | 9.9269257 |
| stimulation + direction + leg + direction:leg + subject | 9.1764175 |
| direction + subject | 8.1385456 |
| stimulation + direction + stimulation:direction + leg + direction:leg + subject | 6.8923098 |
| stimulation + direction + leg + stimulation:leg + subject | 6.3767641 |
| stimulation + direction + stimulation:direction + leg + stimulation:leg + subject | 4.7562332 |
| direction + leg + subject | 1.2906652 |
| stimulation + direction + leg + stimulation:leg + direction:leg + subject | 0.6672675 |
| stimulation + leg + stimulation:leg + subject | 0.6516657 |
| stimulation + direction + stimulation:direction + leg + stimulation:leg + direction:leg + subject | 0.5806634 |
| leg + subject | 0.1579592 |
| direction + leg + direction:leg + subject | 0.1210615 |
| stimulation + direction + stimulation:direction + leg + stimulation:leg + direction:leg + stimulation:direction:leg + subject | 0.0672617 |
Interestingly, there is a lot more evidence across the board, especially for the models that feature Stimulation and Direction.
kable(inclusionBF(bfKanaiCenter, models = "matched"))| effect | Bayes.factor |
|---|---|
| stimulation | 67.6885828 |
| direction | 9.7544628 |
| stimulation:direction | 0.7290812 |
| leg | 0.1726179 |
| stimulation:leg | 0.0681792 |
| direction:leg | 0.0981963 |
| stimulation:direction:leg | 0.1158360 |
There seems to be a mismatch with the classical ANOVA: An effect of stimulation receives strong support, whereas it was non-significant. The effect of direction was significant, but it is less strongly supported by the inclusion Bayes Factor (although the evidence still goes in the right direction).
latencyMedianBaseline %>%
filter(type == "center") %>%
group_by(subject,stimulation) %>%
summarise(latency = mean(latency)) %>%
ggplot(aes(stimulation, latency)) +
geom_hline(yintercept = 0, linetype = "dashed") +
stat_summary(fun.y = mean, geom = "line", aes(group = 1), size = 2) +
geom_line(aes(colour = subject, group = subject))
It is easy to see why this effect is non-significant: the average difference is tiny and there is a lot of variability. The plot shows one major outlier though in terms of the effect size (S01): let’s see what happens to the Bayes Factor if we take their data out.
latencyNoS01 <- latencyMedianBaseline %>%
filter(subject != "S01") %>%
mutate(subject = factor(subject))
bfKanaiCenterNoS01 = anovaBF(latency~stimulation*leg*direction+subject, data = data.frame(filter(latencyNoS01, type == "center")), whichModels="withmain", whichRandom = "subject", progress = FALSE, iterations = 100000) # compute Bayes Factors
bfKanaiCenterNoS01 <- sort(bfKanaiCenterNoS01, decreasing = TRUE) # sort such that winning model is at the top
kable(select(extractBF(bfKanaiCenterNoS01), bf)) # show only the Bayes factors in a table
| bf | |
|---|---|
| stimulation + direction + subject | 1370.5021891 |
| direction + subject | 610.7330230 |
| stimulation + direction + leg + subject | 531.6196659 |
| stimulation + direction + stimulation:direction + subject | 352.9869220 |
| direction + leg + subject | 200.7493579 |
| stimulation + direction + stimulation:direction + leg + subject | 123.1974214 |
| stimulation + direction + leg + direction:leg + subject | 63.0498822 |
| stimulation + direction + leg + stimulation:leg + subject | 42.8446155 |
| direction + leg + direction:leg + subject | 24.4897141 |
| stimulation + direction + stimulation:direction + leg + direction:leg + subject | 14.9256981 |
| stimulation + direction + stimulation:direction + leg + stimulation:leg + subject | 8.3443542 |
| stimulation + direction + leg + stimulation:leg + direction:leg + subject | 4.0678602 |
| stimulation + subject | 1.9848227 |
| stimulation + direction + stimulation:direction + leg + stimulation:leg + direction:leg + subject | 1.0168697 |
| stimulation + leg + subject | 0.6129997 |
| leg + subject | 0.3006916 |
| stimulation + direction + stimulation:direction + leg + stimulation:leg + direction:leg + stimulation:direction:leg + subject | 0.1633244 |
| stimulation + leg + stimulation:leg + subject | 0.0422355 |
kable(inclusionBF(bfKanaiCenterNoS01, models = "matched"))| effect | Bayes.factor |
|---|---|
| stimulation | 2.3502132 |
| direction | 699.4732539 |
| stimulation:direction | 0.2487328 |
| leg | 0.3664546 |
| stimulation:leg | 0.0767869 |
| direction:leg | 0.1186097 |
| stimulation:direction:leg | 0.1606149 |
This completely abolishes the strong support for Stimulation, and greatly enhances the support for Direction, bringing the Bayesian and the classical ANOVAs more in line. At present it is unclear why the classical and Bayesian analyses differ in this regard. When simulating this case for normally distributed data, the Bayes Factors and p-values track each other nicely (see discussion on JASP/BayesFactor forum). This suggests there must be some assumption that is not met in this particular dataset, which is causing the divergence between the analyses.
Still, let’s do some follow-up tests with all of the data (i.e. including this subject) to see whether the anodal or cathodal change scores are significantly different from 0 on their own.
Bayesian one-sample t-tests:
latencyMedianBaseline %>%
filter(type == "center") %>% # keep only center saccades
group_by(stimulation,subject) %>% # for each session and subject
summarise(deviation.end = mean(latency)) %>% # average over all other variables
spread(stimulation,deviation.end) %>% # make separate columns with test data
summarise_if(is.numeric, funs(extractBF(ttestBF(.), onlybf = TRUE))) %>% # run Bayesian t-test on each column, keeping only the BF
gather(stimulation,BF,anodal,cathodal) %>% # make row for each stimulation condition
kable(.)| stimulation | BF |
|---|---|
| anodal | 0.3218926 |
| cathodal | 0.2840670 |
Frequentist one-sample t-tests:
latencyMedianBaseline %>%
filter(type == "center") %>% # keep only center saccades
group_by(stimulation,subject) %>% # for each session and subject
summarise(deviation.end = mean(latency)) %>% # average over all other variables (df is now still grouped per stimulation)
summarise_if(is.numeric, funs(list(tidy(t.test(.))))) %>% # run one-sample t-test for each stimulation condition, return tidy data frames
unnest() %>% # unpack the list-column with data frame for each test
kable(.)| stimulation | estimate | statistic | p.value | parameter | conf.low | conf.high | method | alternative |
|---|---|---|---|---|---|---|---|---|
| anodal | -1.608974 | -0.9875463 | 0.3328370 | 25 | -4.964508 | 1.746559 | One Sample t-test | two.sided |
| cathodal | 1.096154 | 0.8335773 | 0.4124131 | 25 | -1.612139 | 3.804446 | One Sample t-test | two.sided |
So neither are actually significant or have a BF with evidence for the alternative.