• Description
• Data
• Data analysis

Reproducing data analysis and figures reported in “Low intensity TMS enhances perception of visual stimuli” by Abrahamyan, Clifford, Arabzadeh and Harris (2015) published in Brain Stimulation. You can accesss single-subject data and reproduce the analysis by downloading files from github.com/armanabraham. If you use RStudio, load the downloaded project and reproduce the analysis by running Data_analysis_and_figures.Rmd file.

Subject data are provided in RData (Data_on_Github.RData) and Excel (Data_On_Github.xslx) formats.

library(plyr)
library(ggplot2)
library(ggthemes)
library(gridExtra)
library(effsize) # for cohen.d
library(pastecs) # For stat.desc
library(afex)
library(nlme) # For linear mixed-effects modelling, lme
library(BayesFactor) # Bayes factor analaysis of different experimental designs
## Load Winston Chang's functions to compute within-subjects error bars
source('WithinSbj_StdErrors.R')
# Load Exp 1 and 2 data
load("Data_On_Github.RData")

Show how data are formatted for the analysis

head(exp1_dataMean, 4)

##   PPT Condition   Th Sensitivity
## 1 S01       100 0.14608210    7.956563
## 2 S01        60 0.16881954    6.003824
## 3 S01        70 0.09309656   11.156462
## 4 S01        80 0.12383785    8.291708

Here we test if visual sensitivity changes as a function of TMS stimulation intensity when subjects discriminated gratings oriented ±45 deg. Due to repeated-measures design, we conducted mixed-modelling analysis wherein subject was a random effect and TMS stimulation intensity was fixed effect. We fitted linear, cubic and quadratic models to the data and found that quadratic fit explains data better than linear or cubic fits (χ2(1)=5.22, p=0.022).

################## --- Data analysis --- ###################
# Order Conditions for plotting convenience
exp1_dataMean$Condition <- factor(exp1_dataMean$Condition,
   levels=c("NoTMS", "ip80", "ip90", "60", "70", "80", "90", "100"))
# Find out if there is a quadratic trend in the data as a function of TMS intensity
# Select only experimental conditions (exclude NoTMS and ipsilatoral stimulation)
exp1_quadraticFitData <- droplevels(subset(exp1_dataMean, !(Condition %in% c("ip80", "ip90", "NoTMS"))))
exp1_quadraticFitData <- within(exp1_quadraticFitData, {CondAsNum <- as.numeric(as.character(Condition)) })
baseline <- lme(data=exp1_quadraticFitData, Sensitivity~1, random=~1|PPT/CondAsNum, method="ML")
linear <- update(baseline, .~. + CondAsNum)
quadratic <- update(baseline, .~. + CondAsNum + I(CondAsNum^2))
cubic <- update(baseline, .~. + CondAsNum + I(CondAsNum^2)+ I(CondAsNum^3))

anova(baseline, linear, quadratic, cubic)

##   Model df      AIC      BIC    logLik   Test  L.Ratio p-value
## baseline      1  4 248.5679 256.5973 -120.2840
## linear        2  5 250.1203 260.1570 -120.0602 1 vs 2 0.447629  0.5035
## quadratic     3  6 246.9039 258.9479 -117.4519 2 vs 3 5.216463  0.0224
## cubic         4  7 248.2253 262.2766 -117.1126 3 vs 4 0.678586  0.4101

This corresponds to Figure 1 in the manuscript. The dashed line shows quadratic trend in the data.

# Compute high-res quadratic trend for the plot
fitX <- data.frame(CondAsNum=seq(60,100, 0.1))
predVals <- data.frame(x=fitX$CondAsNum, y=predict(quadratic, fitX, level=0))
# Plot sensitivity from Exp 1 together with quadratic trend
# Summarise data for plotting
exp1_dataForPlot <- droplevels(subset(exp1_dataMean, Condition!="NoTMS"))
exp1_dataForPlot$CoilPosition <- "Expt"
exp1_dataForPlot$CoilPosition[exp1_dataForPlot$Condition %in% c("ip80", "ip90")] <- "Cntr"
exp1_dataForPlot$CoilPosition <- as.factor(exp1_dataForPlot$CoilPosition)
## Change the name of the column called "Condition" into "Intensity" cause
## it is more intuitive
colnames(exp1_dataForPlot)[which(names(exp1_dataForPlot) == "Condition")] <- "Intensity"
## Now change ip80 and ip90 into 80 and 90
exp1_dataForPlot$Intensity <- mapvalues(exp1_dataForPlot$Intensity,
   from = c("ip80", "ip90"), to = c("80", "90"))
## Change order of columns (not sure if this is needed, but looks clearer)
exp1_dataForPlot <- droplevels(exp1_dataForPlot[,c("PPT", "CoilPosition",
                                                   "Intensity", "Th", "Sensitivity")])
# Compute means and within-subjects std errors
exp1_qFitForPlot <- summarySEwithin(exp1_dataForPlot,
                                    measurevar="Sensitivity",
                                    withinvars=c("Intensity", "CoilPosition"),
                                    idvar="PPT")
gExp1_SensTrend <- ggplot(exp1_qFitForPlot, aes(x=as.numeric(as.character(Intensity)),
                                                y=Sensitivity, color=CoilPosition)) +
  theme_few() +
  theme(legend.position=c(0.19,0.87)) +
  xlab("Stimulation intensity\n (% phosphene threshold)") +
  ylab("Contrast sensitivity") +
  geom_line(data=predVals, aes(x=x, y=y), linetype="dashed", colour="grey70", size=.3) +
  geom_line(subset=.(CoilPosition=="Expt"), aes(group=1), color="grey50", size=1) +
  geom_errorbar(aes(ymin = Sensitivity - se, ymax = Sensitivity + se),
                width = 0.00,
                size = 0.2,
                color="grey50") +   # Change "se" to "ci" to plot 95% conf intervals
  geom_point(size=7, color="white") +
  geom_point(aes(shape=CoilPosition), size=4, color="grey20") +
  scale_shape_manual(values=c(1,19), name="Coil position",
                     labels=c("Control (Ipsi)","Experimental (Contra)"))
print(gExp1_SensTrend)

To confirm enhancement in visual sensitivity, here we compare experimental (contralateral) stimulation with control (ipsilateral) stimulation. Because ipsilateral control condition was only tested at two TMS intensities (80 and 90% of PT), we performed a 2×2 repeated measures ANOVA comparing the location of the TMS coil (experimental vs control) and stimulation intensities (80 and 90%).

# Compute 2x2 ANOVA comparing 80% and 90% ipsi and contra conditions
exp1_data8090 <- droplevels(subset(exp1_dataForPlot,  Intensity %in% c("80", "90")))
aov8090 <- aov.car(data=exp1_data8090,
   Sensitivity~CoilPosition*Intensity+Error(PPT/(CoilPosition*Intensity)),
                   return="Anova", args.return=list(es="pes"))

## Warning in stri_c(..., sep = sep, collapse = collapse, ignore_null = TRUE):
## longer object length is not a multiple of shorter object length

nice.anova(aov8090, es="pes")

##   Effect    df  MSE         F  pes     p
## 1           CoilPosition 1, 10 0.76 33.45 ***  .77 .0002
## 2              Intensity 1, 10 1.67      0.10 .009   .76
## 3 CoilPosition:Intensity 1, 10 1.67      0.50  .05   .49

This analysis confirmed that there was a significant main effect for the location of the TMS coil (F(1,10)=33.45, p<0.001, partial η2=0.77). However, the main effect of stimulation intensity was not significant, indicating no difference between 80 and 90% intensities (F(1,10)=0.10, p=0.76, partial η2=0.009). Further, there was no interaction between TMS location and stimulation intensity (F(1,10)=0.50, p=0.49, partial η2=0.05).

We followed up with paired t-tests, which showed significant differences between experimental and control conditions at both 80% and 90% stimulation (80%: Mdiff =1.24, t(10)=2.65, p=0.02, Cohen’s d=0.27; 90%: Mdiff =1.80, t(10)=3.80, p=0.003, Cohen’s d=0.47).

80% contra vs ipsi

# t-tests
# Compare 80% experimental and control conditions for sensitivity
exp1_data80 <- droplevels(subset(exp1_data8090, Intensity=="80"))
exp1_ttest80 <- t.test(Sensitivity~CoilPosition, data=exp1_data80, paired=TRUE)
# Compute effect size r using formula from DSUR (sqrt(t^2/(t^2+df)))
paste('r =',
   sqrt(exp1_ttest80$statistic[[1]]^2/(exp1_ttest80$statistic[[1]]^2+exp1_ttest80$parameter[[1]])))

## [1] "r = 0.641549923796786"

# Cohen's d is another estimate of effect size.
cohen.d(data=exp1_data80, Sensitivity~CoilPosition)

##
## Cohen's d
##
## d estimate: -0.268937 (small)
## 95 percent confidence interval:
##   inf        sup
## -1.2064354  0.6685615

90% contra vs ipsi

# 90% experimental vs control
exp1_data90 <- droplevels(subset(exp1_data8090, Intensity=="90"))
exp1_ttest90 <- t.test(Sensitivity~CoilPosition, data=exp1_data90, paired=TRUE)
print(exp1_ttest90)

##
##  Paired t-test
##
## data:  Sensitivity by CoilPosition
## t = -3.8084, df = 10, p-value = 0.003438
## alternative hypothesis: true difference in means is not equal to 0
## 95 percent confidence interval:
##  -2.8447441 -0.7446994
## sample estimates:
## mean of the differences
##   -1.794722

paste('r =',
   sqrt(exp1_ttest90$statistic[[1]]^2/(exp1_ttest90$statistic[[1]]^2+exp1_ttest90$parameter[[1]])))

## [1] "r = 0.769350361865236"

cohen.d(data=exp1_data90, Sensitivity~CoilPosition)
</code rsplus>
 
<code rsplus>
##
## Cohen's d
##
## d estimate: -0.472718 (small)
## 95 percent confidence interval:
##   inf        sup
## -1.4198125  0.4743765

Control (ipsi) vs no-TMS

## Compare ipsi and no-TMS for
exp1_ipsiNoTMS <- droplevels(subset(exp1_dataMean, Condition %in% c("ip80", "ip90", "NoTMS")))
exp1_ipsiNoTMSSens <- cast(exp1_ipsiNoTMS, PPT~Condition, value=.(Sensitivity))
exp1_ipsiNoTMSSens <- ddply(exp1_ipsiNoTMSSens, .(PPT, NoTMS, ip80, ip90), mutate, ipMean=(ip80+ip90)/2)
exp1_ttestNoTMS <- t.test(exp1_ipsiNoTMSSens$NoTMS, exp1_ipsiNoTMSSens$ipMean, paired=TRUE)
print(exp1_ttestNoTMS)

##
##  Paired t-test
##
## data:  exp1_ipsiNoTMSSens$NoTMS and exp1_ipsiNoTMSSens$ipMean
## t = 2.6245, df = 10, p-value = 0.0254
## alternative hypothesis: true difference in means is not equal to 0
## 95 percent confidence interval:
##  0.1650503 2.0208311
## sample estimates:
## mean of the differences
##   1.092941

# Compute effect size r
paste('r =',
   sqrt(exp1_ttestNoTMS$statistic[[1]]^2/(exp1_ttestNoTMS$statistic[[1]]^2+exp1_ttestNoTMS$parameter[[1]])))

## [1] "r = 0.638638128620297"

# Cohen's d
cohen.d(exp1_ipsiNoTMSSens$NoTMS, exp1_ipsiNoTMSSens$ipMean)

##
## Cohen's d
##
## d estimate: 0.2624841 (small)
## 95 percent confidence interval:
##   inf        sup
## -0.6747955  1.1997637

Best TMS condition (90% of PT) vs no-TMS

exp1_NoTMSAnd90 <- droplevels(subset(exp1_dataMean, Condition %in% c("90", "NoTMS")))
exp1_ttestNoTMSvs90 <- t.test(data=exp1_NoTMSAnd90, Sensitivity~Condition, paired=TRUE)
print(exp1_ttestNoTMSvs90)

##
##  Paired t-test
##
## data:  Sensitivity by Condition
## t = -0.82135, df = 10, p-value = 0.4306
## alternative hypothesis: true difference in means is not equal to 0
## 95 percent confidence interval:
##  -1.8679120  0.8617083
## sample estimates:
## mean of the differences
##   -0.5031019

# Effect size as r
paste('r =',
   sqrt(exp1_ttestNoTMSvs90$statistic[[1]]^2/(exp1_ttestNoTMSvs90$statistic[[1]]^2+exp1_ttestNoTMSvs90$parameter[[1]])))

## [1] "r = 0.251391125531333"

# Cohen's d
exp1_NoTMSAnd90Wide <- cast(exp1_NoTMSAnd90, PPT~Condition, value='Sensitivity')
cohen.d(exp1_NoTMSAnd90Wide$`90`, exp1_NoTMSAnd90Wide$NoTMS)

##
## Cohen's d
##
## d estimate: 0.1276009 (negligible)
## 95 percent confidence interval:
##   inf       sup
## -0.806314  1.061516