tDCS associate list differences
Molly Lewis
2020-03-30
This is structured to include only the final analyses in the paper.
rm(list = ls())
# load packages
library(knitr)
library(tidyverse)
library(tidytext)
library(RCurl)
library(forcats)
library(langcog)
library(data.table)
library(igraph)
library(intergraph)
library(ggnetwork)
library(stringr)
library(lme4)
opts_chunk$set(echo = T, message = F, warning = F,
error = F, cache = F, tidy = F)
theme_set(theme_bw())Read in data
# tdcs data
d.all = read_csv("../data/labDataPlusNorms.csv") %>%
rename(condition = electrode,
association.num = mention) %>%
mutate(condition = fct_recode(condition, "sham" = "na")) %>%
filter(as.character(target) != as.character(cue)) %>%
mutate(target= str_trim(target)) %>%
select(labSubj, condition, cue, target, association.num) %>%
arrange(condition, cue, association.num)
d_all2 <- read_csv("../data/WAS mentions.csv") %>%
rename(condition = electrode,
association.num = mention,
target = associate) %>%
filter(as.character(target) != as.character(cue)) %>%
mutate(condition = fct_recode(condition, "sham" = "na")) %>%
mutate(target= str_trim(target)) %>%
arrange(condition, cue, association.num)Associate-co-occurences
Response networks
Modularity setup
# get all d_all2
all.pairs = d_all2 %>%
filter(as.character(target) != as.character(cue)) %>%
group_by(labSubj, cue, condition) %>%
do(m = data.table(t(combn(.$target, 2)))) %>%
ungroup() %>%
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
closeness = mean(estimate_closeness(graph, cutoff= 100))
betweeness = mean(estimate_betweenness(graph, cutoff= 100),
weights = counts$n)
degree = mean(igraph::degree(graph))
weighted_degree = mean(igraph::graph.strength(graph))
data.frame(closeness = round(closeness, 4),
betweeness = round(betweeness, 4),
weighted_degree = round(weighted_degree, 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_all2 %>%
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) %>%
mutate(type = "x") %>%
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_all2 %>%
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) %>%
mutate(type = "x") %>%
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")
ga
#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::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()## 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(col = "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 |
|---|---|---|---|---|---|---|---|
| 0.6036974 | 0.7846923 | 0.4374987 | 38 | -0.9537557 | 2.161151 | Paired t-test | two.sided |
# mean degree
modularity %>%
group_by(condition) %>%
multi_boot_standard(col = "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.6749349 | -4.16159 | 0.0001743 | 38 | -1.003255 | -0.3466147 | Paired t-test | two.sided |
# weighted node degree
modularity %>%
group_by(condition) %>%
multi_boot_standard(col = "weighted_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("weighted_degree") +
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 |
|---|---|---|---|---|---|---|---|
| 0.6036974 | 0.7846923 | 0.4374987 | 38 | -0.9537557 | 2.161151 | Paired t-test | two.sided |
In the paper, report weighted mean degree (in previous version, unweighted mean degree). Here are the means:
mean.w.degree = modularity %>%
select(cue, condition, weighted_degree)
mean.w.degree %>%
group_by(condition) %>%
multi_boot_standard(col = "weighted_degree", na.rm = T) %>%
kable()| condition | ci_lower | ci_upper | mean |
|---|---|---|---|
| anodal | 8.027909 | 8.986647 | 8.504613 |
| cathodal | 8.720102 | 9.626863 | 9.158826 |
| sham | 8.550067 | 9.343082 | 8.930749 |
And, the model (M15a in the paper):
contrasts(mean.w.degree$condition) <-c(.5, -.5, 0)
colnames(contrasts(mean.w.degree$condition)) <- c('anodalVScathodal','experimentalVSsham')
M15a <- lmer(scale(weighted_degree) ~ condition + (1|cue), data = mean.w.degree)
degree.base <- lmer(scale(weighted_degree) ~ 1 + (1|cue), data = mean.w.degree)
kable(tidy(anova(M15a, degree.base, test = "Chi")))| term | df | AIC | BIC | logLik | deviance | statistic | Chi.Df | p.value |
|---|---|---|---|---|---|---|---|---|
| degree.base | 3 | 293.4318 | 301.7183 | -143.7159 | 287.4318 | NA | NA | NA |
| M15a | 5 | 286.0245 | 299.8354 | -138.0123 | 276.0245 | 11.40724 | 2 | 0.0033339 |
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.45 -0.19 -0.7"Model coefficients: -0.45 -0.19 -0.7
Pairwise comparision models:
summary(M15b <- lmer(scale(weighted_degree) ~ condition + (1|cue), data = mean.w.degree[mean.w.degree$condition != "cathodal",]))## Linear mixed model fit by REML ['lmerMod']
## Formula: scale(weighted_degree) ~ condition + (1 | cue)
## Data: mean.w.degree[mean.w.degree$condition != "cathodal", ]
##
## REML criterion at convergence: 199.1
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -1.95579 -0.47702 -0.09557 0.47180 2.04183
##
## Random effects:
## Groups Name Variance Std.Dev.
## cue (Intercept) 0.6698 0.8184
## Residual 0.3207 0.5663
## Number of obs: 78, groups: cue, 39
##
## Fixed effects:
## Estimate Std. Error t value
## (Intercept) -0.1485 0.1594 -0.932
## conditionsham 0.2970 0.1282 2.316
##
## Correlation of Fixed Effects:
## (Intr)
## conditinshm -0.402summary(M15c <- lmer(scale(weighted_degree) ~ condition + (1|cue), data = mean.w.degree[mean.w.degree$condition != "anodal",]))## Linear mixed model fit by REML ['lmerMod']
## Formula: scale(weighted_degree) ~ condition + (1 | cue)
## Data: mean.w.degree[mean.w.degree$condition != "anodal", ]
##
## REML criterion at convergence: 199.3
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -1.51829 -0.54194 -0.08103 0.52346 1.41774
##
## Random effects:
## Groups Name Variance Std.Dev.
## cue (Intercept) 0.6910 0.8313
## Residual 0.3153 0.5615
## Number of obs: 78, groups: cue, 39
##
## Fixed effects:
## Estimate Std. Error t value
## (Intercept) 0.08174 0.16063 0.509
## conditionsham -0.16348 0.12716 -1.286
##
## Correlation of Fixed Effects:
## (Intr)
## conditinshm -0.396summary(M15d <- lmer(scale(weighted_degree) ~ condition + (1|cue), data = mean.w.degree[mean.w.degree$condition != "sham",]))## Linear mixed model fit by REML ['lmerMod']
## Formula: scale(weighted_degree) ~ condition + (1 | cue)
## Data: mean.w.degree[mean.w.degree$condition != "sham", ]
##
## REML criterion at convergence: 202.3
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -1.85975 -0.59500 -0.05659 0.54649 1.73682
##
## Random effects:
## Groups Name Variance Std.Dev.
## cue (Intercept) 0.5962 0.7721
## Residual 0.3708 0.6089
## Number of obs: 78, groups: cue, 39
##
## Fixed effects:
## Estimate Std. Error t value
## (Intercept) -0.2121 0.1575 -1.347
## conditioncathodal 0.4242 0.1379 3.076
##
## Correlation of Fixed Effects:
## (Intr)
## condtncthdl -0.438