Brisk_Experiment_4
us
2023-07-14
Our design allows to test the within subject effect of background uncertainty as well as its between subject effect. For the power analysis to determine the sample size we calculate the power to find the between subject effect, as this will be the lower end for the within subject effect.
We want to expose participants to three levels of background uncertainty (control, risk, ambiguity). Participants complete four risk elicitation tasks. We count the number of risky choices and divide it by the number of safe choices, this results in a ratio measure. We then take the log of the ratio as our dependent variable.
Power analysis for the test of control vs. risk and control vs. ambiguity
We base our assumptions about the effect size on experiment 1 and experiment 3. In experiment 1 we observed b_control = 0.64 b1_risk = -.75 b2_ambiguity = -1.02
In experiment 3, considering mpl responses we observed b0_mpl = -0.13 b1_ambiguity = -0.32
bret was overall lower in the control group b2_bret = -.66 and not affected by the treatment b3_bret X ambiguity = .36
The interaction of treatment and measure was f = 0.11, p = 0.13.
As MPL has been used widely in previous and our own study, we use the mpl task as the reference group. Using the MPL task, we have 11 possible choices between safe and risky payoff. The point when the participant switches from risky to safe option indicates the degree of risk aversion. For this reason, the log ratio for each switching point presents a linear measure of risk aversion. The resulting log ratios for an MPL with 11 decisions is shown below. A ratio of 0 indicates risk neutrality, ratios below 0 indicate risk aversion, ratios above 0 indicate risk seeking.
risky_count <- c(0:10)
max <- length(risky_count)-1
ratio <- risky_count/(max-risky_count)
log(ratio)## [1] -Inf -2.1972246 -1.3862944 -0.8472979 -0.4054651 0.0000000
## [7] 0.4054651 0.8472979 1.3862944 2.1972246 Infplot(log(ratio))
We aim to test two main effects hypotheses: H1: background risk results in higher risk aversion compared to a neutral background H2: backgound risk results in higher risk aversion compared to background ambiguity. H3: background risk effect is weaker than the effect of background ambiguity.
To test the effects we use an indexed model for each group and compare the risk aversion in part1 and part2. This results in 4 models:
\[ \text{[1] } y[control,control] = \alpha_{[1]}+ \beta1_{[1]} * part2 \] \[ \text{[2] } y[control,risk] = \alpha_{[2]} + \beta1_{[2]}* part2\]
\[ \text{[3] } y[control,ambiguity] = \alpha_{[3]} + \beta1_{[3]} * part2 \] \[ \text{[4] } y[risk,ambiguity] = \alpha_{[4]} + \beta1_{[4]} * part2 \]
To test our hypotheses h1 and h2, we compare the difference of the posteriors of beta2. To test h3, we only need to focus on the posterior of in [4]. If they are more than p > 95% outside from 0 on the left (i.e. negative) we have an effect.
\[ H1: \beta1[control] - \beta1[risk]\] \[ H2: \beta1[control] - \beta1[ambiguity]\] \[ H3: \beta1 [ambiguity,risk]\]
Treating effects as fixed, we assume based on our experiment 1 and experiment 3, that after being conservative and halfing the effect to be more conservative.
\[\alpha_{[1-3]} \sim normal(\mu = - 0.1,\sigma = 0) \text { slight risk aversion in all groups during risk neutral condition in part 1 }\] \[\alpha_{[4]} \sim normal(\mu = - 0.4,\sigma = 0) \text { higher risk aversion in background risk condition in part 1}\]
Part 1 and Part 2 Differences
\[\beta_1{[1]} \sim normal(\mu =0,\sigma = 0) \text { because we assume difference when comparing part 1 and part 2 in the control group}\]
\[\beta_1{[2]} \sim normal(\mu =-0.6,\sigma = 0) \text { because we assume an increase of risk aversion by .5 }\]
\[\beta_1{[3]} \sim normal(\mu =-.9,\sigma = 0) \text { because we assume an increase of risk aversion by 0.3 when moving from risk to ambiguity}\]
\[ \beta_1{[4]} \sim normal(\mu =-.3,\sigma = 0) \text { because we assume an increase of risk aversion by 0.3 when moving from risk to ambiguity}\]
Sim starts here…NOT ADJUSTED YET
Using just a simple regression model for one measure and comparing three groups, with background risk being the baseline we have:
\[ \text{Regression } y = \beta_{0}+ \beta_{1}*Control + \beta_{2}*Ambiguity + \beta_{3}part2 + \beta_{4}part2*Control + \beta_{5}part2*Ambiguity\]
beta0 = -.12
beta1 = 0
beta2 = 0
beta3 = -.5
beta4 = 0.5
beta5 = -0.3
sigma0 = .15 # based on previous sdsampling_data <- function(n,beta0,beta1,beta2,beta3,beta4,beta5,sigma0) {
participant_N = n # how many participants?
iter = 1 # how many iterations
task <- c("mpl")
group <- factor(c("control","risk","ambiguity"),
ordered = TRUE,
levels = c("risk","control","ambiguity"))
part <- 1:2
participant_intercept <- rnorm(participant_N,beta0,1)
sim_responses <- participant_intercept %>%
tibble(beta0 = participant_intercept) %>%
mutate (id = 1:n()) %>%
expand_grid(part) %>%
group_by(id) %>%
mutate(group = sample(group,size = 1)) %>%
expand_grid(task) %>%
#mutate(length.out = participant_N) %>%
#group_by(id) %>%
mutate(
## add the rooms..
dummyControl = case_when(group== "control" ~ 1,
TRUE ~ 0),
dummyAmbiguity = case_when(group == "ambiguity" ~ 1,
TRUE ~ 0),
dummyPart = part-1
) %>%
ungroup() %>%
rowwise() %>%
mutate(response = rnorm(1,beta0,sigma0) +
rnorm(1, beta1, sigma0) * dummyControl +
rnorm(1, beta2, sigma0) * dummyAmbiguity +
rnorm(1, beta3, sigma0) * dummyPart+
rnorm(1, beta4, sigma0) * dummyControl * dummyPart+
rnorm(1, beta5, sigma0) * dummyAmbiguity * dummyPart
)
}#Testing if assumptions work
data<-sampling_data(n = 1000,
beta0 = beta0,
beta1 = beta1,
beta2 = beta2,
beta3 = beta3,
beta4 = beta4,
beta5 = beta5,
sigma0 = sigma0)
label_measure <- c(cem = "CEM",
scl = "SCL",
bret = "BRET",
mpl = "MPL")
ggplot(data,
aes(x = group, y = response, color = as.factor(part)))+
geom_violin()+
#geom_jitter(alpha = .2, width = .05) +
geom_hline(yintercept = 0, color = "darkgrey")+
#facet_wrap(~task, labeller = labeller(task = label_measure)) +
scale_y_continuous(name = "Risk preference") +
scale_x_discrete(name = "Experimental Group",
breaks = c("control","risk","ambiguity"),
labels = c("Control","Risk","Ambiguity"))+
theme_classic()
lmer(response ~ dummyControl*dummyPart + dummyAmbiguity*dummyPart + (1|id),data =data) %>%
summary(.)## Linear mixed model fit by REML ['lmerMod']
## Formula: response ~ dummyControl * dummyPart + dummyAmbiguity * dummyPart +
## (1 | id)
## Data: data
##
## REML criterion at convergence: 3526.2
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -2.7214 -0.4709 0.0048 0.4628 2.7416
##
## Random effects:
## Groups Name Variance Std.Dev.
## id (Intercept) 1.04771 1.0236
## Residual 0.05294 0.2301
## Number of obs: 2000, groups: id, 1000
##
## Fixed effects:
## Estimate Std. Error t value
## (Intercept) -0.12399 0.05846 -2.121
## dummyControl 0.07532 0.08101 0.930
## dummyPart -0.50165 0.01813 -27.665
## dummyAmbiguity 0.01713 0.08230 0.208
## dummyControl:dummyPart 0.52709 0.02513 20.978
## dummyPart:dummyAmbiguity -0.28959 0.02553 -11.345
##
## Correlation of Fixed Effects:
## (Intr) dmmyCn dmmyPr dmmyAm dmmC:P
## dummyContrl -0.722
## dummyPart -0.155 0.112
## dummyAmbgty -0.710 0.513 0.110
## dmmyCntrl:P 0.112 -0.155 -0.722 -0.080
## dmmyPrt:dmA 0.110 -0.080 -0.710 -0.155 0.513run_models <- function(n,beta0,beta1,beta2,beta3,beta4,beta5,sigma0){
data = sampling_data(n = n,
beta0,
beta1,
beta2,
beta3,
beta4,
beta5,
sigma0)
within_tidy <- lmerTest::lmer(response ~ dummyControl*dummyPart + dummyAmbiguity*dummyPart + (1|id),
data = data) %>%
broom.mixed::tidy() %>%
mutate (model = "within")
within_tidy
}
resampling_function <- function(sims,n,beta0,beta1,beta2,beta3,beta4,beta5,sigma0) {
purrr::map_dfr(1:sims, ~ run_models(n= n,
beta0 = beta0,
beta1 = beta1,
beta2 = beta2,
beta3 = beta3,
beta4 = beta4,
beta5 = beta5,
sigma0 = sigma0))
} %>%
mutate(SampleSize = n) %>%
filter(str_detect(term,'dummyControl|dummyPart|dummyAmbiguity|dummyControl:dummyPart|dummyPart:dummyAmbiguity'))
# In this function I can only run one model at the time, but we need two.
# this function works
# resampling_function <- function(sims,n,beta0,beta1,beta2,beta3,beta4,beta5,sigma0) {
# purrr::map_dfr(1:sims,
# ~ broom.mixed::tidy(
# lmerTest::lmer(response ~ dummyControl*dummyBRET + dummyAmbiguity*dummyBRET + (1|id),
# data = sampling_data(n = n,
# beta0,
# beta1,
# beta2,
# beta3,
# beta4,
# beta5,
# sigma0)))
# )
# } %>%
# mutate(SampleSize = n) %>%
# filter(str_detect(term,'dummyControl|dummyBRET|dummyAmbiguity|dummyControl:dummyBRET|dummyBRET:dummyAmbiguity'))
#
# resampling_function(sims = 3,
# n= 20,
# beta0 = beta0,
# beta1 = beta1,
# beta2 = beta2,
# beta3 = beta3,
# beta4 = beta4,
# beta5 = beta5,
# sigma0 = sigma0)set.seed(123)
nsims = 1000
resampled_est<-purrr::map_dfr(seq.int(from = 60, to = 240, length =6),~resampling_function(sims = nsims,
n= .x, # these are the values provided by the first function in map_dfr
beta0 = beta0,
beta1 = beta1,
beta2 = beta2,
beta3 = beta3,
beta4 = beta4,
beta5 = beta5,
sigma0 = sigma0), progress = TRUE)saveRDS(resampled_est,"resampled_coefs.RDS")
resampled_est<-readRDS("resampled_coefs.RDS")power_out <- resampled_est %>%
group_by(model,term,SampleSize) %>%
summarise(`Lower 95% CI` = quantile(estimate, (1-0.95)/2,na.rm = TRUE),
`Upper 95% CI` = quantile(estimate, 1-(1-0.95)/2,na.rm = TRUE),
'Mean Est' = mean(estimate,na.rm=TRUE),
'p_smaller_05' = sum(p.value<.025,na.rm = TRUE)/n(),
upper = binom::binom.exact(sum(p.value<.025,na.rm = TRUE),n())$upper,
lower =binom::binom.exact(sum(p.value<.025,na.rm = TRUE),n())$lower,
.groups="drop") ggplot(power_out,aes(x = factor(SampleSize), y = `Mean Est`)) +
geom_point(stat = 'identity') +
geom_errorbar(aes(ymin = `Lower 95% CI`, ymax = `Upper 95% CI`))+
geom_line() +
theme_classic() +
ylab("Estimate") +
xlab("Sample sizes") +
geom_hline(yintercept=0, linetype="dashed", color = "red") +
facet_grid(model~term)## `geom_line()`: Each group consists of only one observation.
## ℹ Do you need to adjust the group aesthetic?
## `geom_line()`: Each group consists of only one observation.
## ℹ Do you need to adjust the group aesthetic?
## `geom_line()`: Each group consists of only one observation.
## ℹ Do you need to adjust the group aesthetic?
## `geom_line()`: Each group consists of only one observation.
## ℹ Do you need to adjust the group aesthetic?
## `geom_line()`: Each group consists of only one observation.
## ℹ Do you need to adjust the group aesthetic?
ggsave("estimates_curve.png", units = "cm", width = 15, height = 10)## `geom_line()`: Each group consists of only one observation.
## ℹ Do you need to adjust the group aesthetic?
## `geom_line()`: Each group consists of only one observation.
## ℹ Do you need to adjust the group aesthetic?
## `geom_line()`: Each group consists of only one observation.
## ℹ Do you need to adjust the group aesthetic?
## `geom_line()`: Each group consists of only one observation.
## ℹ Do you need to adjust the group aesthetic?
## `geom_line()`: Each group consists of only one observation.
## ℹ Do you need to adjust the group aesthetic?ggplot(power_out, aes(x = factor(SampleSize), y = `p_smaller_05`, group= term)) +
geom_point(stat = 'identity') +
geom_errorbar(aes(ymin = `lower`, ymax = `upper`))+
geom_line() +
theme_classic() +
ylab("Estimate") +
xlab("Sample sizes") +
geom_hline(yintercept=0.9, linetype="dashed", color = "red") +
geom_hline(yintercept=0.95, linetype="dashed", color = "blue") +
geom_hline(yintercept=0.99, linetype="dashed", color = "blue") +
facet_grid(model~term)
ggsave("power_curve.png", units = "cm", width = 15, height = 30)