No dilution of causal strength with truely binary variables (Pilot; 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] 71mean(tdata_age$Age)## [1] 34.31034sd(tdata_age$Age)## [1] 14.14997# 1 = male, 2 = female, 3 = other
table(tdata$Sex)##
## 1 2 4
## 14 14 11 = 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)Descriptive Stats
## : SC
## median mean SE.mean CI.mean.0.95 var std.dev
## 1.00000000 0.95586207 0.02562287 0.05248606 0.01903941 0.13798336
## coef.var
## 0.14435489
## ------------------------------------------------------------
## : CC
## median mean SE.mean CI.mean.0.95 var std.dev
## 0.99000000 0.94586207 0.02619798 0.05366413 0.01990369 0.14108045
## coef.var
## 0.14915542
## ------------------------------------------------------------
## : NC
## median mean SE.mean CI.mean.0.95 var std.dev
## 0.0000000000 0.0048275862 0.0025120738 0.0051457499 0.0001830049 0.0135279313
## coef.var
## 2.8022143472library(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.68, 47.09 0.02 664.72 *** .940 <.001
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '+' 0.1 ' ' 1
##
## Sphericity correction method: GGSign. Effect of “Variable” most likely due to the “non-cause” condition (see graph).
# same ANOVA as before
lmeModel <- lmer(value ~ variable + (1|sID), data=tdata_sub)
# follow-up analysis to get means and CIs
ls1 <- lsmeans(a1, c("variable")) # joint evaluation (basically gives the same table)
ls1## variable lsmean SE df lower.CL upper.CL
## SC 0.95586 0.02562 28 0.903376 1.00835
## CC 0.94586 0.02620 28 0.892198 0.99953
## NC 0.00483 0.00251 28 -0.000318 0.00997
##
## 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.95586 0.02562 28 0.903376 1.00835
## CC 0.94586 0.02620 28 0.892198 0.99953
## NC 0.00483 0.00251 28 -0.000318 0.00997
##
## 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.010 0.0359 28 0.279 0.9582
## SC - NC 0.951 0.0265 28 35.906 <.0001
## CC - NC 0.941 0.0265 28 35.485 <.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.010 0.0359 28 -0.0788 0.0988
## SC - NC 0.951 0.0265 28 0.8855 1.0166
## CC - NC 0.941 0.0265 28 0.8754 1.0067
##
## 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.010
CI_low <- -0.0621
CI_up <- 0.0821
Mdiff## [1] 0.01CI_low## [1] -0.0621CI_up## [1] 0.0821# 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.0407351# 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.0517
## Standard Error: 0.2627
## Variance: 0.0690
## Lower CI: -0.4631
## Upper CI: 0.5665
## Weight: 14.4951d$ci.lo## [1] -0.4630594d$ci.hi## [1] 0.5665363d_ci <- (d$ci.hi - d$ci.lo)/2
d_ci## [1] 0.5147979# How many subjects would we need to make the ci-width small enough?
d <- esc_mean_sd(
grp1m = m1[,1], grp1sd = sd1[,1], grp1n = 216,
grp2m = m2[,1], grp2sd = sd2[,1], grp2n = 216,
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.0517
## Standard Error: 0.0962
## Variance: 0.0093
## Lower CI: -0.1369
## Upper CI: 0.2404
## Weight: 107.9639d_ci <- (d$ci.hi - d$ci.lo)/2
d_ci## [1] 0.1886292