No dilution of causal strength with truely binary variables (main study; light switches)
Simon Stephan
#tdata.wide <- read.delim("dummy_data.txt", header=TRUE, sep="\t", na.strings="NA", dec=".", strip.white=TRUE)
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] 77mean(tdata_age$Age)## [1] 33.4213sd(tdata_age$Age)## [1] 13.39814# 1 = male, 2 = female, 3 = other
table(tdata$Sex)##
## 1 2 3
## 62 151 31 = 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
## 1.00000000 0.95425926 0.01105705 0.02179411 0.02640782 0.16250483
## coef.var
## 0.17029422
## ------------------------------------------------------------
## : CC
## median mean SE.mean CI.mean.0.95 var std.dev
## 0.99500000 0.94777778 0.01078760 0.02126300 0.02513643 0.15854474
## coef.var
## 0.16728050
## ------------------------------------------------------------
## : NC
## median mean SE.mean CI.mean.0.95 var std.dev
## 0.00000000 0.04171296 0.01135507 0.02238152 0.02785054 0.16688481
## coef.var
## 4.00079017library(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.71, 368.31 0.04 1864.08 *** .875 <.001
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '+' 0.1 ' ' 1
##
## Sphericity correction method: GGMain effect probably due to the low ratings for the non-cause.
# same ANOVA as before
lmeModel <- lmer(value ~ variable + (1|sID), data=tdata_sub)## boundary (singular) fit: see help('isSingular')# 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.9543 0.0111 215 0.9325 0.9761
## CC 0.9478 0.0108 215 0.9265 0.9690
## NC 0.0417 0.0114 215 0.0193 0.0641
##
## 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.9543 0.0111 215 0.9325 0.9761
## CC 0.9478 0.0108 215 0.9265 0.9690
## NC 0.0417 0.0114 215 0.0193 0.0641
##
## 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.00648 0.0132 215 0.490 0.8761
## SC - NC 0.91255 0.0190 215 47.963 <.0001
## CC - NC 0.90606 0.0187 215 48.411 <.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.00648 0.0132 215 -0.0247 0.0377
## SC - NC 0.91255 0.0190 215 0.8676 0.9574
## CC - NC 0.90606 0.0187 215 0.8619 0.9502
##
## 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.00648
CI_low <- -0.034
CI_up <- 0.0469
Mdiff## [1] 0.00648CI_low## [1] -0.034CI_up## [1] 0.0469# 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.2672249# 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.0333
## Standard Error: 0.0962
## Variance: 0.0093
## Lower CI: -0.1553
## Upper CI: 0.2220
## Weight: 107.9850d$ci.lo## [1] -0.1552624d$ci.hi## [1] 0.2219591d_ci <- (d$ci.hi - d$ci.lo)/2
d_ci## [1] 0.1886107Cluster 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
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 19
## 2 Invariance 186
## 3 Strengthening 11# 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"))(dilut <- 19/216)## [1] 0.08796296(invar <- 186/216)## [1] 0.8611111(stren <- 11/216)## [1] 0.05092593# get CIs for the proportions
prop.test(19,216,correct=FALSE)##
## 1-sample proportions test without continuity correction
##
## data: 19 out of 216, null probability 0.5
## X-squared = 146.69, df = 1, p-value < 2.2e-16
## alternative hypothesis: true p is not equal to 0.5
## 95 percent confidence interval:
## 0.05703566 0.13328995
## sample estimates:
## p
## 0.08796296prop.test(186,216,correct=FALSE)##
## 1-sample proportions test without continuity correction
##
## data: 186 out of 216, null probability 0.5
## X-squared = 112.67, df = 1, p-value < 2.2e-16
## alternative hypothesis: true p is not equal to 0.5
## 95 percent confidence interval:
## 0.8086529 0.9009493
## sample estimates:
## p
## 0.8611111prop.test(11,216,correct=FALSE)##
## 1-sample proportions test without continuity correction
##
## data: 11 out of 216, null probability 0.5
## X-squared = 174.24, df = 1, p-value < 2.2e-16
## alternative hypothesis: true p is not equal to 0.5
## 95 percent confidence interval:
## 0.02867101 0.08887487
## sample estimates:
## p
## 0.05092593## 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.