This is structured to include only the final analyses in the paper.

Read in data

# tdcs data
d.all = read.csv("../data/labDataPlusNorms.csv") %>%
  rename(condition = electrode,
         association.num = mention) %>%
  select(1:4,7,9:19,21,22) %>%
  mutate(condition = fct_recode(condition, "sham" = "na")) %>%
  #filter(as.character(target) != as.character(cue)) %>%
  mutate(target= str_trim(target)) 

# turker accuracy data
d.bw.raw = read.csv("../data/backALLfinal_withALLSimilaritiesPlusNorms.csv")  
d.bw = d.bw.raw %>%
          select(cue, labSubj, subjCode, isRight)

Verify analytical approach

Try to reproduce number of associates analysis in manuscript

num.associates = count(d.all, condition, labSubj, cue) 

num.associates %>%
  group_by(condition) %>%
  summarize(mean = mean(n)) %>%
  kable()

condition	mean
anodal	5.655046
cathodal	5.648598
sham	4.849910

This is almost identical to manuscript (sham differs slightly).

Fit models

contrasts(num.associates$condition) <-c(.5, -.5, 0)
colnames(contrasts(num.associates$condition)) <- c('anodalVScathodal','experimentalVSsham')
  
m1a = lmer(n~condition+(1|labSubj)+ (1|cue), data=num.associates)
ma = lmer(n~1+(1|labSubj)+ (1|cue), data=num.associates)
anova(m1a, ma, test = "Chi")

## Data: num.associates
## Models:
## ma: n ~ 1 + (1 | labSubj) + (1 | cue)
## m1a: n ~ condition + (1 | labSubj) + (1 | cue)
##     Df    AIC    BIC  logLik deviance  Chisq Chi Df Pr(>Chisq)  
## ma   4 4159.9 4181.5 -2075.9   4151.9                           
## m1a  6 4158.0 4190.4 -2073.0   4146.0 5.8496      2    0.05367 .
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

“X2(2)=5.87, p=.053”

Exploring random effects

m1a = lmer(n~condition+(1|labSubj)+ (condition|cue), data=num.associates)
ma = lmer(n~1+(1|labSubj)+ (condition|cue), data=num.associates)
anova(m1a, ma, test = "Chi")

## Data: num.associates
## Models:
## ma: n ~ 1 + (1 | labSubj) + (condition | cue)
## m1a: n ~ condition + (1 | labSubj) + (condition | cue)
##     Df    AIC    BIC  logLik deviance  Chisq Chi Df Pr(>Chisq)  
## ma   9 4168.0 4216.6 -2075.0   4150.0                           
## m1a 11 4166.1 4225.5 -2072.1   4144.1 5.8425      2    0.05387 .
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

m1b = lmer(n~condition+(1|labSubj)+ (1|cue), 
           data=filter(num.associates, condition !="cathodal"))
mb = lmer(n~1+(1|labSubj)+ (1|cue), 
          data=filter(num.associates, condition !="cathodal"))
lmerTest::difflsmeans(m1b)

## Differences of LSMEANS:
##                         Estimate Standard Error   DF t-value Lower CI
## condition anodal - sham      0.8          0.365 27.0    2.26   0.0753
##                         Upper CI p-value  
## condition anodal - sham     1.57    0.03 *
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

anova(m1b, mb, test = "Chi")

## Data: filter(num.associates, condition != "cathodal")
## Models:
## mb: n ~ 1 + (1 | labSubj) + (1 | cue)
## m1b: n ~ condition + (1 | labSubj) + (1 | cue)
##     Df  AIC  BIC  logLik deviance  Chisq Chi Df Pr(>Chisq)  
## mb   4 2760 2780 -1376.0     2752                           
## m1b  5 2757 2782 -1373.5     2747 5.0127      1    0.02516 *
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

….for this m1a model, the model converges when condition is included as a random slope. Check this.

Pairwise comparisions:

m1b = lmer(n~condition+(1|labSubj)+ (1|cue), 
           data=filter(num.associates, condition !="cathodal"))
mb = lmer(n~1+(1|labSubj)+ (1|cue), 
          data=filter(num.associates, condition !="cathodal"))
lmerTest::difflsmeans(m1b)

## Differences of LSMEANS:
##                         Estimate Standard Error   DF t-value Lower CI
## condition anodal - sham      0.8          0.365 27.0    2.26   0.0753
##                         Upper CI p-value  
## condition anodal - sham     1.57    0.03 *
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

anova(m1b, mb, test = "Chi")

## Data: filter(num.associates, condition != "cathodal")
## Models:
## mb: n ~ 1 + (1 | labSubj) + (1 | cue)
## m1b: n ~ condition + (1 | labSubj) + (1 | cue)
##     Df  AIC  BIC  logLik deviance  Chisq Chi Df Pr(>Chisq)  
## mb   4 2760 2780 -1376.0     2752                           
## m1b  5 2757 2782 -1373.5     2747 5.0127      1    0.02516 *
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

“b=-.83, 95% CI [-1.54, -.11]; X2(1)=5.03, p=.025”

Near identical, but slightly different confidence intervals – what package is being used here to calculate confidence intervals? I’m using lmerTest.

m1c = lmer(n~condition+(1|labSubj)+ (1|cue), 
           data=filter(num.associates, condition !="anodal"))
mc = lmer(n~1+(1|labSubj)+ (1|cue), 
          data=filter(num.associates, condition !="anodal"))

lmerTest::difflsmeans(m1c)

## Differences of LSMEANS:
##                           Estimate Standard Error   DF t-value Lower CI
## condition cathodal - sham      0.8          0.346 27.0    2.37    0.109
##                           Upper CI p-value  
## condition cathodal - sham     1.53    0.03 *
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

anova(m1c, mc, test = "Chi")

## Data: filter(num.associates, condition != "anodal")
## Models:
## mc: n ~ 1 + (1 | labSubj) + (1 | cue)
## m1c: n ~ condition + (1 | labSubj) + (1 | cue)
##     Df    AIC    BIC  logLik deviance  Chisq Chi Df Pr(>Chisq)  
## mc   4 2774.1 2794.0 -1383.0   2766.1                           
## m1c  5 2770.6 2795.6 -1380.3   2760.6 5.4625      1    0.01943 *
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

“b=-.82, 95% CI [-1.50, -.14]; X2(1)=5.49, p=.019”

m1d = lmer(n~condition+(1|labSubj)+ (1|cue), data=filter(num.associates, condition !="sham"))
md = lmer(n~1+(1|labSubj)+ (1|cue), data=filter(num.associates, condition !="sham"))
anova(m1d, md, test = "Chi")

## Data: filter(num.associates, condition != "sham")
## Models:
## md: n ~ 1 + (1 | labSubj) + (1 | cue)
## m1d: n ~ condition + (1 | labSubj) + (1 | cue)
##     Df    AIC    BIC  logLik deviance Chisq Chi Df Pr(>Chisq)
## md   4 2804.9 2824.8 -1398.4   2796.9                        
## m1d  5 2806.9 2831.8 -1398.4   2796.9 2e-04      1     0.9876

“p=.99”

Overall, this looks reasonable.

Associate types

Spoken frequency and concreteness

Descriptive statistics

# read in frequency and concreteness
subtlexus.url <- getURL("https://raw.githubusercontent.com/mllewis/RC/master/data/corpus/SUBTLEXus_corpus.txt")
freqs <- read.table(text = subtlexus.url, header=TRUE) %>%
          select(Word,Lg10WF)
concreteness <-read.csv("brysbaert_corpus.csv", header=TRUE) %>%
  select(Word, Conc.M)

associate.types = d.all %>%
  left_join(freqs, by=c("target"= "Word")) %>%
  left_join(concreteness, by=c("target"= "Word")) %>%
  select(condition,labSubj, cue, target, Lg10WF, Conc.M) 

subjs.means = associate.types %>%
  group_by(condition, labSubj, cue) %>%
  summarize(log.freq = mean(Lg10WF, na.rm = T),
            conc = mean(Conc.M, na.rm = T)) %>%
  group_by(condition, labSubj) %>%
  summarize(log.frequency = mean(log.freq, na.rm = T),
            concreteness = mean(conc, na.rm = T)) 

cond.means = subjs.means %>%
  gather(variable, value, -1:-2) %>%
  group_by(variable, condition) %>%
  multi_boot_standard(column = "value", na.rm = T) 

kable(cond.means)

variable	condition	mean	ci_lower	ci_upper
concreteness	anodal	4.060420	3.924228	4.178522
concreteness	cathodal	4.108435	4.012533	4.198627
concreteness	sham	4.080662	3.956819	4.189281
log.frequency	anodal	3.214721	3.172736	3.254843
log.frequency	cathodal	3.237422	3.182005	3.293175
log.frequency	sham	3.208827	3.115074	3.278188

ggplot(cond.means, aes(x = condition, y = mean, fill = condition)) +
    facet_wrap(~variable, scales = "free") +
    geom_bar(position = "dodge", stat = "identity") +
    geom_linerange(aes(ymax = ci_upper, ymin = ci_lower)) +
    ggtitle("Associate-type measures") +
    theme_bw() +
    theme(legend.position = "none")

Models

Fit mixed-effects models.

First, setup contrasts.

contrasts(associate.types$condition) <- c(.5, -.5, 0)
colnames(contrasts(associate.types$condition)) <- c('anodalVScathodal',
'experimentalVSsham')
kable(contrasts(associate.types$condition))

	anodalVScathodal	experimentalVSsham
anodal	0.5	-0.4082483
cathodal	-0.5	-0.4082483
sham	0.0	0.8164966

Frequency

Model M6a in paper

M6a <-lmer(scale(Lg10WF) ~ condition +  (1|labSubj) + (condition|cue), data=associate.types)
kable(tidy(M6a))

term	estimate	std.error	statistic	group
(Intercept)	0.0003009	0.0495119	0.0060771	fixed
conditionanodalVScathodal	-0.0271363	0.0560200	-0.4844040	fixed
conditionexperimentalVSsham	-0.0063399	0.0395229	-0.1604101	fixed
sd_(Intercept).labSubj	0.1310104	NA	NA	labSubj
sd_(Intercept).cue	0.2743833	NA	NA	cue
sd_conditionanodalVScathodal.cue	0.0117969	NA	NA	cue
sd_conditionexperimentalVSsham.cue	0.0146148	NA	NA	cue
cor_(Intercept).conditionanodalVScathodal.cue	-1.0000000	NA	NA	cue
cor_(Intercept).conditionexperimentalVSsham.cue	-1.0000000	NA	NA	cue
cor_conditionanodalVScathodal.conditionexperimentalVSsham.cue	1.0000000	NA	NA	cue
sd_Observation.Residual	0.9552074	NA	NA	Residual

freq.base<-lmer(scale(Lg10WF) ~ 1 +  (1|labSubj) + (condition|cue) , data=associate.types)

kable(tidy(anova(M6a, freq.base, test = "Chi")))

term	df	AIC	BIC	logLik	deviance	statistic	Chi.Df	p.value
freq.base	9	22077.38	22140.23	-11029.69	22059.38	NA	NA	NA
M6a	11	22081.11	22157.93	-11029.55	22059.11	0.2700137	2	0.8737099

Concreteness

Model M6b in paper.

M6b <- lmer(scale(Conc.M) ~ condition+ (1|labSubj) + (1|cue), data=associate.types)
kable(tidy(M6b))

term	estimate	std.error	statistic	group
(Intercept)	0.0005930	0.0587381	0.0100959	fixed
conditionanodalVScathodal	-0.0497765	0.1050779	-0.4737107	fixed
conditionexperimentalVSsham	-0.0030181	0.0730804	-0.0412988	fixed
sd_(Intercept).labSubj	0.2698494	NA	NA	labSubj
sd_(Intercept).cue	0.2528691	NA	NA	cue
sd_Observation.Residual	0.9288013	NA	NA	Residual

conc.base <- lmer(scale(Conc.M) ~ 1 + (1|labSubj) + (1|cue), data=associate.types) # does not converge with random slopes by condition|cue

kable(tidy(anova(M6b, conc.base, test = "Chi")))

term	df	AIC	BIC	logLik	deviance	statistic	Chi.Df	p.value
conc.base	4	21927.94	21955.92	-10959.97	21919.94	NA	NA	NA
M6b	6	21931.70	21973.67	-10959.85	21919.70	0.2406801	2	0.8866189

Cue-associate relationships

Conditional probabilitiy

Merge in Nelson norms

# nelson norms
cp.url <- getURL("https://raw.githubusercontent.com/wjhopper/USF-WAN-SS/master/norms.csv")
cps <- read.csv(text = cp.url, header=FALSE) 

names(cps) = c("cue", "target", "FS", "BS")

cps = cps %>%
      mutate(cue = tolower(as.character(cue)),
             target = tolower(as.character(target)))

dcp =  d.all %>%
    filter(as.character(target) != as.character(cue)) %>%
    left_join(cps) %>%
    select(labSubj, cue, target, condition, association.num, FS, BS)

All cues are present in Nelson norms

cuetdcs = unique(d.all$cue)
cueNelson = unique(cps$cue)

length(intersect(cuetdcs, cueNelson)) # all cues are present in nelson

## [1] 39

Lots of NAs but no difference between critical conditions

dcp %>%
  group_by(condition) %>%
  summarize(na = length(which(is.na(FS))),
            notna = length(which(!is.na(FS)))) %>%
  kable()

condition	na	notna
anodal	1883	1110
cathodal	1887	1040
sham	1668	924

Replacing missing associate values with zero. Note that BS and FS are skewed. Take the log.

dcp <- dcp %>%
       mutate(FS = ifelse(is.na(FS), 0, FS),
              BS = ifelse(is.na(BS), 0, BS),
              log.FS = log(FS),
              log.BS = log(BS))  

# non-finite values are treated as NA

dcp %>%
  gather(variable, value, 6:9) %>%
  ggplot(aes(x = value)) +
  facet_wrap(~variable, scales = "free")+
  geom_histogram() +
  theme_bw()

Do conditional probabilities differ by condition?

dcp.item = dcp %>%
  gather(variable, value, 6:9) %>%
  mutate(value = ifelse(is.infinite(value), NA, value)) %>%
  group_by(condition, labSubj, cue, variable) %>%
  summarize(value = mean(value, na.rm = T)) %>%
  group_by(labSubj, condition, variable) %>%
  summarize(value = mean(value, na.rm = T))

dcp.means = dcp.item %>%
  filter(is.finite(value)) %>%
  group_by(condition, variable) %>%
  multi_boot_standard(column = "value", na.rm = T)

ggplot(dcp.means, aes(x = condition, y = mean, fill = condition)) +
    facet_wrap(~variable, scales = "free") +
    geom_bar(position = "dodge", stat = "identity") +
    geom_linerange(aes(ymax = ci_upper, ymin = ci_lower)) +
    ggtitle("Measures of forward/backward probability") +
    theme_bw() +
    theme(legend.position = "none")

No difference here for either backwards or forward probability.

Models

First, setup contrasts.

contrasts(dcp$condition) <- c(.5, -.5, 0)
colnames(contrasts(dcp$condition)) <- c('anodalVScathodal',
'experimentalVSsham')
kable(contrasts(dcp$condition))

	anodalVScathodal	experimentalVSsham
anodal	0.5	-0.4082483
cathodal	-0.5	-0.4082483
sham	0.0	0.8164966

FW probability

Descriptive statistics

dcp.means %>%
  filter(variable == "FS") %>%
  mutate(mean = round(mean,2))

## Source: local data frame [3 x 5]
## Groups: condition [3]
## 
##   condition variable  mean   ci_lower   ci_upper
##      <fctr>    <chr> <dbl>      <dbl>      <dbl>
## 1    anodal       FS  0.05 0.04277467 0.05660685
## 2  cathodal       FS  0.05 0.04149103 0.05685115
## 3      sham       FS  0.05 0.04152063 0.05637608

Model

M9a <- lmer(scale(FS) ~ condition+ (1|labSubj) + (condition|cue), data = dcp)

FS.base <- lmer(scale(FS) ~ 1 + (1|labSubj) + (condition|cue), data = dcp)

kable(tidy(anova(M9a, FS.base, test = "Chi")))

term	df	AIC	BIC	logLik	deviance	statistic	Chi.Df	p.value
FS.base	9	23657.39	23720.83	-11819.69	23639.39	NA	NA	NA
M9a	11	23661.37	23738.91	-11819.68	23639.37	0.0188213	2	0.9906335

BW probability

Descriptive statistics

dcp.means %>%
  filter(variable == "BS") %>%
  mutate(mean = round(mean,2))

## Source: local data frame [3 x 5]
## Groups: condition [3]
## 
##   condition variable  mean   ci_lower   ci_upper
##      <fctr>    <chr> <dbl>      <dbl>      <dbl>
## 1    anodal       BS  0.05 0.03802430 0.05355545
## 2  cathodal       BS  0.04 0.03594373 0.05066615
## 3      sham       BS  0.04 0.03342589 0.04779287

Model

M9b <- lmer(scale(BS) ~ condition+ (1|labSubj) + (association.num*condition|cue), data = dcp)

BS.base <- lmer(scale(BS) ~ 1 + (1|labSubj) + (association.num*condition|cue), data = dcp)

kable(tidy(anova(M9b, BS.base, test = "Chi")))

term	df	AIC	BIC	logLik	deviance	statistic	Chi.Df	p.value
BS.base	24	23367.64	23536.82	-11659.82	23319.64	NA	NA	NA
M9b	26	23369.41	23552.69	-11658.71	23317.41	2.228804	2	0.3281114

cue-associate similarity

This is the existing plot in the paper

d.all %>%
  select(condition, taxInclusive, partsInclusive, thematic, 
                       action, descriptor, holonym, metaphor) %>%
  gather(word_type, value, -1) %>%
  group_by(condition, word_type) %>%
  summarise(n= sum(value)) %>%
  mutate(freq = n / sum(n)) %>%
  ungroup() %>%
  mutate(word_type = as.factor(word_type),
         word_type = fct_relevel(word_type, "thematic", "taxInclusive", "partsInclusive",
                                 "holonym", "action", "descriptor", "metaphor"),
         condition = fct_relevel(condition, "anodal", "sham", "cathodal")) %>%
  ggplot(aes(y = freq, fill = condition, x = word_type)) +
  ylab("Proportion of listed associates")+
  ylim(0,.4)+
  geom_bar(stat = "identity", position = "dodge") +
    theme(axis.text.x = element_text(angle = 60, hjust = 1))

Associate-co-occurences

Response networks

A third approach: I constructed a network for each cue for each condition. Nodes are the associates given by participants. A link is added between two nodes if a participant listed both associates. Weights are determined by the number of people who listed two associates. If participants are giving more idiosyncratic responses in the cathodal group, we should expect less centrality in the cathodal networks. In the limit, if participants give ONLY idiosnycratic responses, then the network should look like a bunch of separate clusters. Conversely, if everyone gives the same responses, it should look like a hairball.

First, here are the graphs for “bicycle” in the cathodal and anodal conditions.

It looks like participants in the cathodal condition are doing more idiosyncratic responses in general. In principle, it seems like forward and backward probabilites should capture this.

Modularity setup

# get all pairs
all.pairs = d.all %>%
  #filter(association.num < 6) %>%
  filter(as.character(target) != as.character(cue)) %>%
  group_by(labSubj, cue, condition) %>%
  do(m = data.table(t(combn(.$target, 2)))) %>%
  unnest() %>%
  rename(a1 = V1,
         a2 = V2)

# modularity function
get_modularity <- function(df){
    counts <- count(df, a1, a2) 
    graph <- counts %>%
      select(a1, a2) %>%
      graph_from_data_frame(directed = FALSE) %>%
      simplify() 
    
    E(graph)$weight=counts$n # add in weights
    
    clustering = cluster_leading_eigen(graph, 
                              options = list(maxiter=1000000), 
                              weights = counts$n)
    
    closeness = mean(estimate_closeness(graph, cutoff= 100))
    betweeness = mean(estimate_betweenness(graph, cutoff= 100),  
                      weights = counts$n)
    degree = mean(igraph::degree(graph))

    data.frame(Q = round(clustering$modularity,4), 
               n.groups = round(length(clustering),4),
               closeness  = round(closeness, 4),
               betweeness  = round(betweeness, 4),
               degree = mean(degree))
}

getCathodal <- function() {    
  # single item example
  test = all.pairs %>%
    filter(cue == "bicycle" & condition == "cathodal") %>%
    filter(a1 != a2)  %>%
      mutate(a1 = str_trim(a1),
             a2 = str_trim(a2)) 
  
  counts.ex <- count(test, a1, a2) 
  
  v.attr = d.all %>%
    filter(as.character(target) != as.character(cue)) %>%
    mutate(target = str_trim(target)) %>%
    filter(cue == "bicycle" & condition == "cathodal") %>%
    distinct(target, .keep_all = T) %>%
    gather(type, value, c(10, 12,14,11,16, 17,18)) %>%
    filter(value == 1) %>%
    select(target, type) 
  
  graph <- counts.ex %>%
        select(a1, a2) %>%
        graph_from_data_frame(directed = FALSE, vertices = v.attr) %>%
        simplify() 
      
  E(graph)$weight=counts.ex$n # add in weights
  
  gc = ggplot(ggnetwork(asNetwork(graph)), 
            aes(x = x, y = y, xend = xend, yend = yend)) +
         geom_edges(aes(size = weight), color = "grey60", 
                    curvature = .15, 
                    show.legend = FALSE) +
         geom_nodes(aes(color = type), show.legend = FALSE) +
         geom_nodelabel(aes(label =  vertex.names,
                                  fill = type), 
                        show.legend = FALSE, 
                        label.padding = unit(0.2,"lines"), 
                        label.size = 0.4, size = 3,
                        force = .1) +
        ggtitle("Cathodal Stimulation") +
        theme_blank()
  
  gt <- ggplot_gtable(ggplot_build(gc))
  gt$layout$clip[gt$layout$name == "panel"] <- "off"
  gt
}
 
getAnodal <- function() {    
  test = all.pairs %>%
    filter(cue == "bicycle" & condition == "anodal") %>%
    filter(a1 != a2)  %>%
      mutate(a1 = str_trim(a1),
              a2 = str_trim(a2)) 
  
  counts.ex <- count(test, a1, a2) 
  
  v.attr = d.all %>%
    filter(as.character(target) != as.character(cue)) %>%
    mutate(target = str_trim(target)) %>%
    filter(cue == "bicycle" & condition == "anodal") %>%
    distinct(target, .keep_all = T) %>%
    gather(type, value, c(10, 12, 14, 11, 16, 17,18)) %>%
    filter(value == 1) %>%
    select(target, type)  %>%
    mutate(type =  fct_recode(type, 
                              "Part" = "partsInclusive",
                              "Taxonomic" = "taxInclusive",
                              "Action" = "action",
                              "Thematic" = "thematic",
                              "Descriptor" = "descriptor")) 
  
  
  graph <- counts.ex %>%
        select(a1, a2) %>%
        graph_from_data_frame(directed = FALSE, vertices = v.attr) %>%
        simplify() 
      
  E(graph)$weight=counts.ex$n # add in weights
  
  ga = ggplot(ggnetwork(asNetwork(graph)), 
            aes(x = x, y = y, xend = xend, yend = yend)) +
         geom_edges(aes(size = weight), color = "grey60", 
                    curvature = .15, 
                    show.legend = FALSE) +
         geom_nodes(aes(color = type), size = 3) +
         geom_nodelabel(aes(label =  vertex.names , fill = type), 
                        show.legend = FALSE, 
                        label.padding = unit(0.2,"lines"), 
                        label.size = 0.4, size = 3) +
        ggtitle("Anodal Stimulation") +
        theme_blank()  +
        scale_fill_discrete(name = "New Legend Title")
  
      
  
  gt <- ggplot_gtable(ggplot_build(ga))
  gt$layout$clip[gt$layout$name == "panel"] <- "off"
  gt
}

#multiplot(grid.draw(getCathodal()), grid.draw(getAnodal()), cols=2)
cath = getCathodal()
anodal = getAnodal()

#png("legend.png",width = 6, height = 6, units = 'in', res = 500) 
#grid.draw(anodal) 
#dev.off()

#png("legend.png",width = 6, height = 6, units = 'in', res = 500) 

#v.attr %>%
 # count(type) %>%
#  ggplot( aes(x = type, y = n )) +
#  geom_bar(stat= "identity", aes(fill = type)) +
#  scale_fill_discrete(name = "Similarity Relation")
#dev.off()

…it looks like the cathodal network is more modular.

Now, let’s do this across all items and quantify modularity, using three measures: Betweenness, mean degree, and modularity.

## get all modularity
modularity = all.pairs %>%
      filter(a1 != a2)  %>%
      mutate(a1 = str_trim(a1),
             a2 = str_trim(a2)) %>%
     group_by(cue, condition) %>%
     do(get_modularity(.))

#  mean betweeness
modularity %>%
    group_by(condition) %>%
    multi_boot_standard(column = "betweeness", na.rm = T) %>%
    ggplot(aes(x = condition, y = mean, fill = condition)) +
    geom_bar(position = "dodge", stat = "identity") +
    geom_linerange(aes(ymax = ci_upper, ymin = ci_lower)) +
    theme_bw()  +
    ggtitle("betweeness") +
    theme(legend.position = "none")

modularity %>%
  filter(condition != "sham") %>%
  select(cue, condition, betweeness) %>%
  ungroup() %>%
  spread(condition, betweeness) %>%
  do(tidy(t.test(.$anodal, .$cathodal, paired=TRUE)))   %>%
  kable()

estimate	statistic	p.value	parameter	conf.low	conf.high	method	alternative
-2.714844	-3.421774	0.0015021	38	-4.321003	-1.108685	Paired t-test	two.sided

#  mean degree
modularity %>%
    group_by(condition) %>%
    multi_boot_standard(column = "degree", na.rm = T) %>%
    ggplot(aes(x = condition, y = mean, fill = condition)) +
    geom_bar(position = "dodge", stat = "identity") +
    geom_linerange(aes(ymax = ci_upper, ymin = ci_lower)) +
    theme_bw()  +
    ggtitle("mean degree") +
    theme(legend.position = "none")

modularity %>%
  filter(condition != "sham") %>%
  select(cue, condition, degree) %>%
  ungroup() %>%
  spread(condition, degree) %>%
  do(tidy(t.test(.$anodal, .$cathodal, paired=TRUE)))   %>%
  kable()

estimate	statistic	p.value	parameter	conf.low	conf.high	method	alternative
0.3877107	2.541442	0.0152419	38	0.0788784	0.6965431	Paired t-test	two.sided

#  clustering
modularity %>%
    group_by(condition) %>%
    multi_boot_standard(column = "Q", na.rm = T) %>%
    ggplot(aes(x = condition, y = mean, fill = condition)) +
    geom_bar(position = "dodge", stat = "identity") +
    geom_linerange(aes(ymax = ci_upper, ymin = ci_lower)) +
    ggtitle("modularity") +
    theme_bw() +
    theme(legend.position = "none")

modularity %>%
  filter(condition != "sham") %>%
  select(cue, condition, Q) %>%
  ungroup() %>%
  spread(condition, Q) %>%
  do(tidy(t.test(.$anodal, .$cathodal, paired=TRUE)))   %>%
  kable()

estimate	statistic	p.value	parameter	conf.low	conf.high	method	alternative
-0.0350205	-2.870941	0.0066548	38	-0.0597146	-0.0103264	Paired t-test	two.sided

Across three different measures, anodal group is more clustered than cathodal group, suggesting they are all sampling from a shared “concept” of the cue.

In the paper, report mean degree. Here are the means:

mean.degree = modularity %>%
  select(cue, condition, degree) 

mean.degree %>%
    group_by(condition) %>%
    multi_boot_standard(column = "degree", na.rm = T) %>%
   kable()

condition	mean	ci_lower	ci_upper
anodal	8.255716	8.007472	8.495233
cathodal	7.868005	7.651740	8.111151
sham	6.255249	6.021826	6.506498

And, the model (M15a in the paper):

contrasts(mean.degree$condition) <-c(.5, -.5, 0)
colnames(contrasts(mean.degree$condition)) <- c('anodalVScathodal','experimentalVSsham')
  
M15a <- lmer(scale(degree) ~ condition + (1|cue), data = mean.degree)
degree.base <- lmer(scale(degree) ~ 1 + (1|cue), data = mean.degree)
kable(tidy(anova(M15a, degree.base, test = "Chi")))

term	df	AIC	BIC	logLik	deviance	statistic	Chi.Df	p.value
degree.base	3	337.0273	345.3138	-165.5137	331.0273	NA	NA	NA
M15a	5	229.7178	243.5286	-109.8589	219.7178	111.3096	2	0

names =  rownames(summary(M15a)$coefficients)
coefficients = summary(M15a)$coefficients %>%
     as.data.frame() %>%
     mutate(comparision = names) %>%
     rename(SE = `Std. Error`) %>%
     filter(comparision == "conditionanodalVScathodal")

print(paste(round(coefficients$Estimate, 2),
            round(coefficients$Estimate + (coefficients$SE * 1.96),2),
            round(coefficients$Estimate - (coefficients$SE * 1.96),2)))

## [1] "0.34 0.58 0.1"

Model coefficients: 0.34 0.58 0.1

Pairwise comparision models; all are significant.

summary(M15b <- lmer(scale(degree) ~ condition + (1|cue), data = mean.degree[mean.degree$condition != "cathodal",]))

## Linear mixed model fit by REML ['lmerMod']
## Formula: scale(degree) ~ condition + (1 | cue)
##    Data: mean.degree[mean.degree$condition != "cathodal", ]
## 
## REML criterion at convergence: 142.7
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -1.6901 -0.6604 -0.0120  0.5196  2.2336 
## 
## Random effects:
##  Groups   Name        Variance Std.Dev.
##  cue      (Intercept) 0.1324   0.3639  
##  Residual             0.2395   0.4894  
## Number of obs: 78, groups:  cue, 39
## 
## Fixed effects:
##               Estimate Std. Error t value
## (Intercept)    0.79046    0.09765   8.095
## conditionsham -1.58091    0.11082 -14.265
## 
## Correlation of Fixed Effects:
##             (Intr)
## conditinshm -0.567

summary(M15c <- lmer(scale(degree) ~ condition + (1|cue), data = mean.degree[mean.degree$condition != "anodal",]))

## Linear mixed model fit by REML ['lmerMod']
## Formula: scale(degree) ~ condition + (1 | cue)
##    Data: mean.degree[mean.degree$condition != "anodal", ]
## 
## REML criterion at convergence: 155.7
## 
## Scaled residuals: 
##      Min       1Q   Median       3Q      Max 
## -1.82389 -0.64042 -0.03304  0.52346  2.22542 
## 
## Random effects:
##  Groups   Name        Variance Std.Dev.
##  cue      (Intercept) 0.1732   0.4162  
##  Residual             0.2744   0.5238  
## Number of obs: 78, groups:  cue, 39
## 
## Fixed effects:
##               Estimate Std. Error t value
## (Intercept)     0.7423     0.1071   6.929
## conditionsham  -1.4847     0.1186 -12.516
## 
## Correlation of Fixed Effects:
##             (Intr)
## conditinshm -0.554

summary(M15d <- lmer(scale(degree) ~ condition + (1|cue), data = mean.degree[mean.degree$condition != "sham",]))

## Linear mixed model fit by REML ['lmerMod']
## Formula: scale(degree) ~ condition + (1 | cue)
##    Data: mean.degree[mean.degree$condition != "sham", ]
## 
## REML criterion at convergence: 217.6
## 
## Scaled residuals: 
##      Min       1Q   Median       3Q      Max 
## -1.80243 -0.62591 -0.00911  0.51941  2.12670 
## 
## Random effects:
##  Groups   Name        Variance Std.Dev.
##  cue      (Intercept) 0.1704   0.4128  
##  Residual             0.7768   0.8813  
## Number of obs: 78, groups:  cue, 39
## 
## Fixed effects:
##                   Estimate Std. Error t value
## (Intercept)         0.2536     0.1558   1.627
## conditioncathodal  -0.5072     0.1996  -2.541
## 
## Correlation of Fixed Effects:
##             (Intr)
## condtncthdl -0.640

Number of incorrect guesses

Here are the means:

d.bw2 = d.bw.raw %>%
       select(electrode,  labSubj,  subjCode, cue, 
             turkerGuess, isRight) %>%
        rename(condition = electrode)  %>%     
        mutate(condition = fct_recode(condition, "sham" = "na")) %>%
        mutate(turkerGuess = str_trim(turkerGuess))

cue.ns = d.bw2 %>%
  filter(isRight == 0) %>%
  group_by(condition, cue) %>%
  summarize(n = length(unique(turkerGuess))) 
  
cue.ms = cue.ns %>%
  group_by(condition) %>%
  multi_boot_standard(column = "n", na.rm = T)  

ggplot(cue.ms, aes(x = condition, y = mean, fill = condition)) +
    geom_bar(position = "dodge", stat = "identity") +
    geom_linerange(aes(ymax = ci_upper, ymin = ci_lower)) +
    ggtitle("mean number of unique wrong guesses") +
    theme_bw() +
    theme(legend.position = "none")

 kable(cue.ms)

condition	mean	ci_lower	ci_upper
anodal	46.33333	41.53718	50.59103
cathodal	49.30769	45.20449	53.00128
sham	52.02564	47.99936	56.43590

The model (M16a in paper)

contrasts(cue.ns$condition) <- c(.5, -.5, 0)
colnames(contrasts(cue.ns$condition)) <- c('anodalVScathodal','experimentalVSsham')
  
M16a <- lmer(scale(n) ~ condition + (1|cue), data = cue.ns)
wrong.base <- lmer(scale(n) ~ 1 + (1|cue), data = cue.ns)
kable(tidy(anova(M16a, wrong.base, test = "Chi")))

term	df	AIC	BIC	logLik	deviance	statistic	Chi.Df	p.value
wrong.base	3	278.2508	286.5373	-136.1254	272.2508	NA	NA	NA
M16a	5	270.5438	284.3546	-130.2719	260.5438	11.70703	2	0.0028698

names =  rownames(summary(M16a)$coefficients)
coefficients = summary(M16a)$coefficients %>%
     as.data.frame() %>%
     mutate(comparision = names) %>%
     rename(SE = `Std. Error`) %>%
     filter(comparision == "conditionexperimentalVSsham")

print(paste(round(coefficients$Estimate, 2),
            round(coefficients$Estimate + (coefficients$SE * 1.96),2),
            round(coefficients$Estimate - (coefficients$SE * 1.96),2)))

## [1] "0.25 0.41 0.09"

Model coefficients: 0.25 0.41 0.09

M16b in paper

M16b<- lmer(scale(n) ~ condition + (1|cue), data = cue.ns[cue.ns$condition != "cathodal",])

wrong.base <- lmer(scale(n) ~ 1 + (1|cue), data = cue.ns[cue.ns$condition != "cathodal",])
kable(tidy(anova(M16b, wrong.base, test = "Chi")))

term	df	AIC	BIC	logLik	deviance	statistic	Chi.Df	p.value
wrong.base	3	209.9556	217.0257	-101.97780	203.9556	NA	NA	NA
M16b	4	203.6639	213.0907	-97.83195	195.6639	8.291691	1	0.0039827

names =  rownames(summary(M16b)$coefficients)
coefficients = summary(M16b)$coefficients %>%
     as.data.frame() %>%
     mutate(comparision = names) %>%
     rename(SE = `Std. Error`) %>%
    slice(2)

print(paste(round(coefficients$Estimate, 2),
            round(coefficients$Estimate + (coefficients$SE * 1.96),2),
            round(coefficients$Estimate - (coefficients$SE * 1.96),2)))

## [1] "0.4 0.65 0.14"

Model coefficients: 0.4 0.65 0.14

M16c in paper

M16c<- lmer(scale(n) ~ condition + (1|cue), data = cue.ns[cue.ns$condition != "anodal",])
wrong.base <- lmer(scale(n) ~ 1 + (1|cue), data = cue.ns[cue.ns$condition != "anodal",])
kable(tidy(anova(M16c, wrong.base, test = "Chi")))

term	df	AIC	BIC	logLik	deviance	statistic	Chi.Df	p.value
wrong.base	3	202.7162	209.7863	-98.35808	196.7162	NA	NA	NA
M16c	4	202.1227	211.5495	-97.06133	194.1227	2.59351	1	0.1073023

names =  rownames(summary(M16c)$coefficients)
coefficients = summary(M16c)$coefficients %>%
     as.data.frame() %>%
     mutate(comparision = names) %>%
     rename(SE = `Std. Error`) %>%
     slice(2)

print(paste(round(coefficients$Estimate, 2),
            round(coefficients$Estimate + (coefficients$SE * 1.96),2),
            round(coefficients$Estimate - (coefficients$SE * 1.96),2)))

## [1] "0.2 0.45 -0.04"

Model coefficients: 0.2 0.45 -0.04

M16d in paper

M16d<- lmer(scale(n) ~ condition + (1|cue), data = cue.ns[cue.ns$condition != "sham",])
wrong.base <- lmer(scale(n) ~ 1 + (1|cue), data = cue.ns[cue.ns$condition != "sham",])
kable(tidy(anova(M16d, wrong.base, test = "Chi")))

term	df	AIC	BIC	logLik	deviance	statistic	Chi.Df	p.value
wrong.base	3	181.2404	188.3105	-87.62021	175.2404	NA	NA	NA
M16d	4	177.5321	186.9590	-84.76607	169.5321	5.708273	1	0.0168851

names =  rownames(summary(M16d)$coefficients)
coefficients = summary(M16d)$coefficients %>%
     as.data.frame() %>%
     mutate(comparision = names) %>%
     rename(SE = `Std. Error`) %>%
    slice(2)

print(paste(round(coefficients$Estimate, 2),
            round(coefficients$Estimate + (coefficients$SE * 1.96),2),
            round(coefficients$Estimate - (coefficients$SE * 1.96),2)))

## [1] "0.22 0.39 0.04"

Model coefficients: 0.22 0.39 0.04

tDCS associate list differences

Molly Lewis

2017-04-26

Verify analytical approach

Associate types

Descriptive statistics

Models

Frequency

Concreteness

Cue-associate relationships

Conditional probabilitiy

Models

FW probability

BW probability

cue-associate similarity

Associate-co-occurences

Response networks

Number of incorrect guesses