No dilution of causal strength with truely binary variables (Probabilistic causes; switches)
Simon Stephan
tdata <- read.delim("data.txt", header=TRUE, sep="\t", na.strings="NA", dec=".", strip.white=TRUE)Results
Demographics
# demographics
tdata_age <- tdata
min(tdata_age$Age)## [1] 18max(tdata_age$Age)## [1] 75mean(tdata_age$Age)## [1] 34.94444sd(tdata_age$Age)## [1] 12.31927# 1 = male, 2 = female, 3 = other
table(tdata$Sex)##
## 1 2 3
## 89 125 21 = male, 2 = female, 3 = non-binary
Graphs
myTheme <- theme(plot.title = element_text(face="bold", size = 22),
axis.title.x = element_blank(),
axis.title.y = element_text(face = "bold", size = 20),
axis.text.x = element_text(size = 18, angle = 0),
axis.text.y = element_text(size = 16, angle = 0),
legend.text = element_text(size = 18),
legend.title = element_text(face = "bold", size = 18),
strip.text.x = element_text(size = 18),
#panel.grid.major = element_blank(),
panel.grid.minor = element_blank(),
panel.background = element_blank(),
axis.line.x = element_line(colour = "black"),
axis.line.y = element_line(colour = "black"),
axis.text = element_text(colour ="black"),
axis.ticks = element_line(colour ="black"))
tdata_sub <- tdata_long
library(see)
## first, turn sID into a factor
tdata_sub$sID <- factor(tdata_sub$sID)
pd <- position_dodge(width = 0.3)
tdata_sub$valueJitter <- jitter(tdata_sub$value, factor = 1, amount = 0.04)
theme_set(theme_light(base_size = 20, base_family = "Poppins"))
# new labes for the facets
g <- ggplot(tdata_sub, aes(x=variable, y=valueJitter, group = sID)) +
guides(fill=FALSE)+
#facet_grid(Query_order ~ Cause_order)+
#ggtitle("Subjects' causal srength ratings") +
scale_y_continuous(limits = c(-0.05, 1.05), breaks=seq(0, 1, 0.1), expand = c(0,0)) +
scale_x_discrete(labels=c("Single-effect \ncause", "Common \ncause", "No \ncause")) +
#stat_summary(fun.y = mean, geom = "bar", position = "dodge", colour = "black", alpha =0.5) +
geom_violinhalf(aes(y = value, group = variable, fill = variable), color = NA, position=position_dodge(1), alpha = 0.2)+
geom_line(position = pd, color = "black", size = 1, alpha=0.04) +
geom_point(aes(color = variable), position = pd, alpha = 0.2) +
stat_summary(aes(y = value,group=1), fun.data = mean_cl_boot, geom = "errorbar", width = 0, size = 1) +
stat_summary(aes(y = value,group=1), fun.y=mean, colour="black", geom="line",group=1, size = 1.5, linetype = "solid", alpha = 1)+
stat_summary(aes(y = value,group=1, fill = variable), fun.y=mean, geom="point", color = "black", shape = 22, size = 5, group=1, alpha = 1)+
stat_summary(aes(y = value,group=1), fun.y=median, geom="point", color = "black", shape = 3, size = 4, group=1, alpha = 1, position = position_dodge(width = 0.5))+
labs(x = "Number Cause's Effects", y = "Causal Strength Rating") +
#scale_color_manual(name = "Entity",values=c("#fc9272", "#3182bd"))+
#scale_fill_manual(name = "Entity",values=c("#fc9272", "#3182bd"))+
theme(legend.position = "none")+
myTheme## Warning: `guides(<scale> = FALSE)` is deprecated. Please use `guides(<scale> =
## "none")` instead.## Warning: `fun.y` is deprecated. Use `fun` instead.
## `fun.y` is deprecated. Use `fun` instead.
## `fun.y` is deprecated. Use `fun` instead.g
#ggsave("results_lines.svg",width=6,height=4.3)
#ggsave("results_lines.pdf",width=6,height=4.3)myTheme <- theme(plot.title = element_text(face="bold", size = 22),
axis.title.x = element_blank(),
axis.title.y = element_text(face = "bold", size = 20),
axis.text.x = element_text(size = 18, angle = 0),
axis.text.y = element_text(size = 16, angle = 0),
legend.text = element_text(size = 18),
legend.title = element_text(face = "bold", size = 18),
strip.text.x = element_text(size = 18),
#panel.grid.major = element_blank(),
panel.grid.minor = element_blank(),
panel.background = element_blank(),
axis.line.x = element_line(colour = "black"),
axis.line.y = element_line(colour = "black"),
axis.text = element_text(colour ="black"),
axis.ticks = element_line(colour ="black"))
tdata_sub <- tdata_long
tdata_sub_graph <- subset(tdata_sub, variable == "CC" | variable == "SC")
library(see)
## first, turn sID into a factor
tdata_sub$sID <- factor(tdata_sub$sID)
pd <- position_dodge(width = 0.3)
tdata_sub_graph$valueJitter <- jitter(tdata_sub_graph$value, factor = 1, amount = 0.04)
theme_set(theme_light(base_size = 20, base_family = "Poppins"))
# new labes for the facets
g <- ggplot(tdata_sub_graph, aes(x=variable, y=valueJitter, group = sID)) +
guides(fill=FALSE)+
#facet_grid(Query_order ~ Cause_order)+
#ggtitle("Subjects' causal srength ratings") +
scale_y_continuous(limits = c(-0.05, 1.05), breaks=seq(0, 1, 0.1), expand = c(0,0)) +
scale_x_discrete(labels=c("single-effect \ncause", "common \ncause", "No \ncause")) +
#stat_summary(fun.y = mean, geom = "bar", position = "dodge", colour = "black", alpha =0.5) +
geom_violinhalf(aes(y = value, group = variable, fill = variable), color = NA, position=position_dodge(1), alpha = 0.2)+
geom_line(position = pd, color = "black", size = 1, alpha=0.04) +
geom_point(aes(color = variable), position = pd, alpha = 0.2) +
stat_summary(aes(y = value,group=1), fun.data = mean_cl_boot, geom = "errorbar", width = 0, size = 1) +
stat_summary(aes(y = value,group=1), fun.y=mean, colour="black", geom="line",group=1, size = 1.5, linetype = "solid", alpha = 1)+
stat_summary(aes(y = value,group=1, fill = variable), fun.y=mean, geom="point", color = "black", shape = 22, size = 5, group=1, alpha = 1)+
stat_summary(aes(y = value,group=1), fun.y=median, geom="point", color = "black", shape = 3, size = 4, group=1, alpha = 1, position = position_dodge(width = 0.5))+
labs(x = "Number Cause's Effects", y = "Causal Strength Rating") +
scale_color_manual(name = "Entity",values=c("#fc9272", "#3182bd"))+
scale_fill_manual(name = "Entity",values=c("#fc9272", "#3182bd"))+
theme(legend.position = "none")+
myTheme## Warning: `guides(<scale> = FALSE)` is deprecated. Please use `guides(<scale> =
## "none")` instead.## Warning: `fun.y` is deprecated. Use `fun` instead.
## `fun.y` is deprecated. Use `fun` instead.
## `fun.y` is deprecated. Use `fun` instead.g
#ggsave("results_lines_dilut.svg",width=6,height=4.3)
#ggsave("results_lines_dilut.pdf",width=4.5,height=4.3)Descriptive Stats
## : SC
## median mean SE.mean CI.mean.0.95 var std.dev
## 0.79500000 0.78384259 0.01497235 0.02951138 0.04842098 0.22004768
## coef.var
## 0.28072942
## ------------------------------------------------------------
## : CC
## median mean SE.mean CI.mean.0.95 var std.dev
## 0.77000000 0.74712963 0.01621575 0.03196220 0.05679730 0.23832185
## coef.var
## 0.31898327
## ------------------------------------------------------------
## : NC
## median mean SE.mean CI.mean.0.95 var std.dev
## 0.00000000 0.05787037 0.01303991 0.02570243 0.03672847 0.19164672
## coef.var
## 3.31165539library(afex)## ************
## Welcome to afex. For support visit: http://afex.singmann.science/## - Functions for ANOVAs: aov_car(), aov_ez(), and aov_4()
## - Methods for calculating p-values with mixed(): 'S', 'KR', 'LRT', and 'PB'
## - 'afex_aov' and 'mixed' objects can be passed to emmeans() for follow-up tests
## - NEWS: emmeans() for ANOVA models now uses model = 'multivariate' as default.
## - Get and set global package options with: afex_options()
## - Set orthogonal sum-to-zero contrasts globally: set_sum_contrasts()
## - For example analyses see: browseVignettes("afex")
## ************##
## Attache Paket: 'afex'## Das folgende Objekt ist maskiert 'package:lme4':
##
## lmerlibrary(emmeans)
a1 <- aov_car(value ~ variable + Error(sID/(variable)), tdata_sub)
a1## Anova Table (Type 3 tests)
##
## Response: value
## Effect df MSE F ges p.value
## 1 variable 1.69, 363.47 0.05 883.12 *** .703 <.001
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '+' 0.1 ' ' 1
##
## Sphericity correction method: GG# same ANOVA as before
lmeModel <- lmer(value ~ variable + (1|sID), data=tdata_sub)
# follow-up analysis
ls1 <- lsmeans(a1, c("variable")) # joint evaluation (basically gives the same table)
ls1## variable lsmean SE df lower.CL upper.CL
## SC 0.7838 0.0150 215 0.7543 0.8134
## CC 0.7471 0.0162 215 0.7152 0.7791
## NC 0.0579 0.0130 215 0.0322 0.0836
##
## Confidence level used: 0.95###############
# a conditional analysis
ls2 <- lsmeans(a1, c("variable")) # group means by between-condition
ls2## variable lsmean SE df lower.CL upper.CL
## SC 0.7838 0.0150 215 0.7543 0.8134
## CC 0.7471 0.0162 215 0.7152 0.7791
## NC 0.0579 0.0130 215 0.0322 0.0836
##
## Confidence level used: 0.95# simple main effects
t <- pairs(ls2) # compares rep-measure differences separately for each between-factor level
t## contrast estimate SE df t.ratio p.value
## SC - CC 0.0367 0.0147 215 2.490 0.0359
## SC - NC 0.7260 0.0211 215 34.342 <.0001
## CC - NC 0.6893 0.0217 215 31.725 <.0001
##
## P value adjustment: tukey method for comparing a family of 3 estimatesconfint(t, level = 0.95)## contrast estimate SE df lower.CL upper.CL
## SC - CC 0.0367 0.0147 215 0.00192 0.0715
## SC - NC 0.7260 0.0211 215 0.67608 0.7759
## CC - NC 0.6893 0.0217 215 0.63798 0.7405
##
## Confidence level used: 0.95
## Conf-level adjustment: tukey method for comparing a family of 3 estimatesNo dilution
Make a difference plot:
#t <- qt(0.975, 29, lower.tail = TRUE, log.p = FALSE)
#t
effect <- "Mdiff"
Mdiff <- 0.0367
CI_low <- -0.00906
CI_up <- 0.0825
Mdiff## [1] 0.0367CI_low## [1] -0.00906CI_up## [1] 0.0825# Plot
myTheme <- theme(plot.title = element_text(face="bold", size = 22),
axis.title.x = element_text(face = "bold", size = 20),
axis.title.y = element_blank(),
axis.text.x = element_text(size = 18, angle = 0),
axis.text.y = element_text(size = 40, angle = 0),
legend.text = element_text(size = 18),
legend.title = element_text(size = 22),
strip.text.x = element_text(size = 18),
#panel.grid.major = element_blank(),
#panel.grid.minor = element_blank(),
#panel.background = element_blank(),
axis.line.x = element_line(colour = "black"),
axis.line.y = element_line(colour = "black"),
axis.text = element_text(colour ="black"),
axis.ticks = element_line(colour ="black"))
theme_set(theme_light(base_size = 30, base_family = "Poppins"))
barchart <- ggplot()+
myTheme+
#guides(fill=FALSE)+
#facet_wrap(~Latency + SampleSize, ncol=2)+
#ggtitle("Mean difference (95% CI)") +
#coord_cartesian(ylim=c(-0.1,2)) +
scale_y_continuous(limits = c(-0.1, 0.5), breaks=seq(-0.1, 0.5, 0.1), expand = c(0,0)) +
scale_x_discrete(labels=c("r")) +
#annotate("rect", xmin=1.7, xmax=2.3, ymin=0.95, ymax=1.05, color="#31a354", fill = "white", size = 1) +
#stat_summary(fun.y=mean, colour="grey20", geom="point", shape = 21, size = 3)+
#stat_summary(fun.y = mean, geom = "bar", position = "dodge", colour = "black")+
#stat_summary(fun.data = mean_cl_boot, geom = "errorbar", position = position_dodge(width = 0.90), width = 0.2) +
#geom_jitter(width = 0.3, height = 0.02, alpha = 0.6, colour = "red") +
#ggtitle("Means (95% bootstr. CIs)") +
#theme(axis.text.x = element_text(size = 10, angle = 0, hjust = 0.5))+
labs(x= "", y = "Mean change") +
#scale_color_manual(values=c("#005083", "#f0b64d"))# +
#scale_fill_manual(values=c("#969696", "#969696"))
#annotate("point", x = 1, y = 100, colour = "firebrick", size = 2)+
#annotate(xmin = -Inf, xmax = Inf, ymin = 4.77-1.96*0.297, ymax = 4.77+1.96*0.297, geom = "rect", alpha = 0.2, fill = "firebrick")+
#annotate(xmin = -Inf, xmax = Inf, ymin = 5.02-1.96*0.372, ymax = 5.02+1.96*0.372, geom = "rect", alpha = 0.2, fill = "blue")+
#annotate(geom = "hline",yintercept = 100, y = 100, color = "red")+
annotate("pointrange", x = 1, y = Mdiff, ymin = CI_low, ymax = CI_up, colour = "black", size = 2, shape = 24, fill = "darkgrey")+
#annotate("pointrange", x = 2, y = 5.02, ymin = 5.02-1.96*0.372, ymax = 5.02+1.96*0.372, colour = "blue", size = 0.8, shape = 15)+
#annotate("text", x = 0.5, y = 2.6, family = "Poppins", size = 6, color = "gray20", label = "Impfeffekt")+
#geom_curve(aes(x = 0.5, y = 3, xend = 0.9, yend = 4),arrow = arrow(length = unit(0.03, "npc")),color = "gray20", curvature = +0.2)+
#annotate("text", x = 1.8, y = 2.6, family = "Poppins", size = 6, color = "gray20", label = "Dosierungseffekt")+
#geom_curve(aes(x = 1.8, y = 3, xend = 2, yend = 4),arrow = arrow(length = unit(0.03, "npc")),color = "gray20", curvature = +0.2)+
annotate(geom = "hline",yintercept = 0, y = 0, color = "red", size = 1.2)+
theme(plot.background = element_rect(
fill = "white",
colour = "white",
size = 1
))## Warning: Ignoring unknown aesthetics: ybarchart
#ggsave("delta.svg",width=2.5,height=4)
#ggsave("delta.pdf",width=2.5,height=4)What Cohen’s d is this?
dat <- tdata_sub
# since we have a repeated-meausres design, we now need the correlations of the ratings
library(dplyr) # for pipe operator
tdata -> t
r <- cor(t$single_strength_rating, t$multiple_strength_rating)
r## [1] 0.5555428# now compute ES and SE and CI of it
# using the esc package because it gives SE of the ES directly
library(esc)
# get means and sds
m1 <- dat %>%
filter(variable == "SC")%>%
summarize(Mean1 = mean(value))
sd1 <- dat %>%
filter(variable == "SC")%>%
summarize(SD1 = sd(value))
m2 <- dat %>%
filter(variable == "CC")%>%
summarize(Mean2 = mean(value))
sd2 <- dat %>%
filter(variable == "CC")%>%
summarize(SD2 = sd(value))
d <- esc_mean_sd(
grp1m = m1[,1], grp1sd = sd1[,1], grp1n = length(dat$sID)/3,
grp2m = m2[,1], grp2sd = sd2[,1], grp2n = length(dat$sID)/3,
r = r,
es.type = "d"
)
d##
## Effect Size Calculation for Meta Analysis
##
## Conversion: mean and sd (within-subject) to effect size d
## Effect Size: 0.1694
## Standard Error: 0.0964
## Variance: 0.0093
## Lower CI: -0.0195
## Upper CI: 0.3584
## Weight: 107.6138d$ci.lo## [1] -0.01950226d$ci.hi## [1] 0.3583692d_ci <- (d$ci.hi - d$ci.lo)/2
d_ci## [1] 0.1889357Cluster Analysis
data_cluster <- tdata[c(1:5,7)]
# append absolut deviations as new columns
data_cluster$Delta_Rating <- data_cluster$single_strength_rating - data_cluster$multiple_strength_rating
crit = 3 # define a threshold that needs to be exceeded to leave the "invariance" cluster
data_cluster$Change[data_cluster$Delta_Rating > crit] <- "Dilution"
data_cluster$Change[data_cluster$Delta_Rating >= -crit & data_cluster$Delta_Rating <= crit] <- "Invariance"
data_cluster$Change[data_cluster$Delta_Rating < -crit] <- "Strengthening"
data_cluster %>% count(Change)## Change n
## 1 Dilution 43
## 2 Invariance 151
## 3 Strengthening 22# now append the clustering results to the main data frame
data_cluster <- subset(data_cluster, select = c(1,8))
tdata_sub <- subset(tdata_sub, variable != "NC")
tdata_sub <- merge(tdata_sub, data_cluster, by = c("sID"))# get CIs for the proportions
prop.test(43,216,correct=FALSE)##
## 1-sample proportions test without continuity correction
##
## data: 43 out of 216, null probability 0.5
## X-squared = 78.241, df = 1, p-value < 2.2e-16
## alternative hypothesis: true p is not equal to 0.5
## 95 percent confidence interval:
## 0.1512878 0.2573770
## sample estimates:
## p
## 0.1990741prop.test(151,216,correct=FALSE)##
## 1-sample proportions test without continuity correction
##
## data: 151 out of 216, null probability 0.5
## X-squared = 34.241, df = 1, p-value = 4.87e-09
## alternative hypothesis: true p is not equal to 0.5
## 95 percent confidence interval:
## 0.6348662 0.7563248
## sample estimates:
## p
## 0.6990741prop.test(22,216,correct=FALSE)##
## 1-sample proportions test without continuity correction
##
## data: 22 out of 216, null probability 0.5
## X-squared = 136.96, df = 1, p-value < 2.2e-16
## alternative hypothesis: true p is not equal to 0.5
## 95 percent confidence interval:
## 0.06822738 0.14939063
## sample estimates:
## p
## 0.1018519## Warning: `guides(<scale> = FALSE)` is deprecated. Please use `guides(<scale> =
## "none")` instead.## Warning: `fun.y` is deprecated. Use `fun` instead.
## `fun.y` is deprecated. Use `fun` instead.