Skip to content

Commit

Permalink
odds ratio scale for dictator game results
Browse files Browse the repository at this point in the history
  • Loading branch information
@phelps-sg
phelps-sg committed Jan 19, 2024
1 parent f1475a4 commit 45a7bf0
Showing 1 changed file with 104 additions and 38 deletions.
142 changes: 104 additions & 38 deletions jupyter-book/R-MixedModel.R
View file
Original file line number Diff line number Diff line change
Expand Up @@ -42,7 +42,8 @@ libs <- c(
"R.cache",
"patchwork",
"ordinal",
"cmdstanr"
"cmdstanr",
"checkmate"
)

lapply(libs, function(lib) library(lib, character.only = TRUE))
Expand Down Expand Up @@ -256,7 +257,18 @@ summary(results_dictator$Num_cooperates)
levels(results_dictator$Model)

# %%
hist(results_dictator$Num_cooperates)
dictator_counts <- results_dictator %>% count(Num_cooperates)

# %%
pdf("figs/dictator-num-cooperates-hist.pdf")
ggplot(dictator_counts, aes(x = Num_cooperates, y = n)) +
geom_bar(stat = "identity") +
labs(
title = "Histogram of Dictator Game responses",
x = "Donation amount",
y = "Frequency"
)
dev.off()

# %%
results_dictator$Response <- as.factor(results_dictator$Num_cooperates)
Expand All @@ -266,32 +278,51 @@ results_dictator$Response <- as.factor(results_dictator$Num_cooperates)
model_dictator <- clmm(
Response ~
Participant_group + Participant_group:Model + t + Model + Temperature +
(1 | Participant_id),
(1 | Participant_id),
link = "logit", threshold = "equidistant",
data = results_dictator
)
summary(model_dictator)

#f%%
# f%%

model_dictator_hess <- clmm(
Response ~
Participant_group + Participant_group:Model + t + Model + Temperature +
(1 | Participant_id),
(1 | Participant_id),
link = "logit", threshold = "equidistant",
data = results_dictator,
Hess = TRUE,
nAGQ = 10
)
summary(model_dictator_hess)

# %% [markdown]

# Remove t parameter from the model since it had insigicant estimate and is not expected to have any causal effect.

# %%
model_dictator_1 <- clmm(
Response ~
Participant_group + Participant_group:Model + Model + Temperature +
(1 | Participant_id),
link = "logit", threshold = "equidistant",
data = results_dictator
)

# %%
summary(model_dictator_1)

# %%
print(odds_ratios_significant(model_dictator_1), type = "latex")

# %%
model_pd_poisson <- glmmTMB(
Num_cooperates ~
Participant_group + Participant_group:Partner_condition + t +
Model + Temperature +
Participant_group:Partner_condition:Model +
(1 | Participant_id),
Model + Temperature +
Participant_group:Partner_condition:Model +
(1 | Participant_id),
data = results_pd,
family = poisson
)
Expand All @@ -301,8 +332,8 @@ summary(model_pd_poisson)
model_pd <- glmmTMB(
cbind(Num_cooperates, 6 - Num_cooperates) ~
Participant_group + Participant_group:Partner_condition +
t + Model + Temperature +
(1 | Participant_id),
t + Model + Temperature +
(1 | Participant_id),
data = results_pd,
family = betabinomial
)
Expand All @@ -314,14 +345,27 @@ summary(model_pd)
model_pd_1 <- glmmTMB(
cbind(Num_cooperates, 6 - Num_cooperates) ~
Participant_group + Participant_group:Partner_condition +
t + Model + Temperature +
Participant_group:Partner_condition:Model +
(1 | Participant_id),
t + Model + Temperature +
Participant_group:Partner_condition:Model +
(1 | Participant_id),
data = results_pd,
family = betabinomial
)
summary(model_pd_1)

# %%
model_pd_2 <- glmmTMB(
cbind(Num_cooperates, 6 - Num_cooperates) ~
Participant_group + Participant_group:Partner_condition +
Model + Temperature +
Participant_group:Partner_condition:Model +
(1 | Participant_id),
data = results_pd,
family = betabinomial
)
summary(model_pd_2)


# %%
texreg(

Expand Down Expand Up @@ -349,41 +393,55 @@ hist(ranef(model_pd_1)$cond$Participant_id[, 1])
dev.off()

# %%
hist(ranef(model_dictator_1)$Participant_id[, 1])
pdf("figs/ranef_hist_dictator.pdf")
hist(ranef(model_dictator)$Participant_id[, 1])
dev.off()
# %%
xtable(lme4::formatVC(summary(model_pd)$varcor$cond))

# %%
m <- coef(summary(model_pd_1))$cond
m[, 1] <- exp(m[, 1])
m[, 2] <- exp(m[, 2])
colnames(m)[1] <- "Estimate (Odds ratio)"
as_odds <- function(model) {
m <- coef(summary(model))
if (testClass(m, "list")) {
m <- m$cond
}
m[, 1] <- exp(m[, 1])
m[, 2] <- exp(m[, 2])
colnames(m)[1] <- "Estimate (Odds ratio)"
return(m)
}

# %%
m[m[, 4] <= 0.05, ]
as_odds(model_pd_1)

# %%
odds_ratios_significant <-
odds_ratios_significant <- function(model) {
m <- as_odds(model)
xtable(
m[m[, 4] < 0.05, ],
caption =
"Model estimates for significant coefficients only ($p<0.05$) on odds ratio scale"
)
print(odds_ratios_significant, type = "latex")
}

# %%

odds_ratios_significant(model_pd_1)

# %%
print(odds_ratios_significant(model_pd_1), type = "latex")

# %%
xtable(coef(summary(model_pd))$cond, digits = 3)

# %%
predictions_by_group <- function(model) {
ggpredict(model, c("Participant_group", "Model")) %>%
rename_at("group", ~"Model")
ggpredict(model, c("Participant_group", "Model")) %>%
rename_at("group", ~"Model")
}

# %%
pd_predictions <- predictions_by_group(model_pd_1)
pd_predictions <- predictions_by_group(model_pd_2)

# %%
theoretical_data <- function(group, frequencies) {
Expand Down Expand Up @@ -423,8 +481,7 @@ theoretical_dictator <- function(group) {
}

theoretical_df_for_group <- function(
group, fn = function(x) x, experiment = theoretical_dilemma
) {
group, fn = function(x) x, experiment = theoretical_dilemma) {
d <- experiment(group)
theoretical_data(group, fn(d))
}
Expand Down Expand Up @@ -470,9 +527,17 @@ pd_predictions
predictions_plot <- function(predictions, title) {
ggplot(predictions) +
aes(x = x, y = predicted, group = Model) +
geom_errorbar(aes(ymin = conf.low, ymax = conf.high), width = .1, position = position_dodge(0.06)) +
geom_errorbar(
aes(ymin = conf.low, ymax = conf.high),
width = .1, position = position_dodge(0.06)
) +
geom_ribbon(
aes(x = x, y = predicted, ymax = conf.high, ymin = conf.low, fill = Model),
alpha = 0.2
) +
geom_line(aes(color = Model), size = 1) + # scale_y_continuous(limits = c(0, 1)) +
scale_color_brewer(palette = "Dark2", direction = -1) +
scale_fill_brewer(palette = "Dark2", direction = -1, guide = "none") +
labs(
title = title,
x = "Participant group", y = "Level of cooperation"
Expand All @@ -488,14 +553,11 @@ pd_predictions_plot <-
actual_and_theoretical(pd_predictions),
"Prisoners Dilemma"
)
print(pd_predictions_plot)
dev.off()
pd_predictions_plot

# %%
summary(model_dictator)

# %% jupyter={"outputs_hidden": false} pycharm={"name": "#%%\n"}
predictions_dictator <- predictions_by_group(model_dictator)
predictions_dictator <- predictions_by_group(model_dictator_1)

# %%
weighted_dictator <- predictions_dictator %>%
Expand All @@ -514,7 +576,7 @@ predictions_dictator <- actual_and_theoretical(weighted_dictator, h = hypothesiz
predictions_dictator

# %%
dictator_predictions_plot <-
dictator_predictions_plot <-
predictions_plot(predictions_dictator, "Dictator Game")
dictator_predictions_plot

Expand Down Expand Up @@ -554,10 +616,15 @@ interaction_plots <- function(participant_group, legend) {

p <- ggplot(combined) +
aes(x = x, y = predicted, group = Model) +
geom_errorbar(aes(ymin = conf.low, ymax = conf.high), width = .1, position = position_dodge(0.06)) +
geom_ribbon(
aes(x = x, y = predicted, ymax = conf.high, ymin = conf.low, fill = Model),
alpha = 0.2
) +
geom_errorbar(aes(ymin = conf.low, ymax = conf.high), width = .1, position = position_dodge(0.16)) +
geom_line(aes(color = Model), size = 1) +
scale_y_continuous(limits = c(0, 1)) +
scale_color_brewer(palette = "Dark2", direction = -1) +
scale_fill_brewer(palette = "Dark2", direction = -1, guide = "none") +
labs(
title = sprintf("Group: %s", participant_group),
x = "Partner condition", y = "Probability of cooperation"
Expand Down Expand Up @@ -602,7 +669,7 @@ combined.plots
# %%
options(repr.plot.height = 12)
pdf("figs/residuals-plots.pdf")
simulationOutput <-
simulationOutput <-
simulateResiduals(fittedModel = model_pd_1, plot = TRUE, integerResponse = TRUE)
dev.off()

Expand Down Expand Up @@ -718,7 +785,7 @@ pr <- prior(normal(0, 1), class = "b")
brm_model <- brm(
Num_cooperates | trials(6) ~
Participant_group + Partner_condition + t + Model +
Model:Partner_condition:Participant_group,
Model:Partner_condition:Participant_group,
~ (1 | Participant_id),
data = results_pd,
prior = NULL,
Expand Down Expand Up @@ -764,7 +831,7 @@ estimate_glmm <- function() {
MCMCglmm(
fixed =
Num_cooperates ~ Participant_group + Partner_condition + t +
Model + Model:Partner_condition:Participant_group,
Model + Model:Partner_condition:Participant_group,
random = ~ us(1):Participant_id,
prior = prior,
data = results_pd,
Expand Down Expand Up @@ -821,4 +888,3 @@ summary(m6)

# options(repr.plot.height = 24, repr.plot.width = 20)
# plot_estimates(m6)

0 comments on commit 45a7bf0

Please sign in to comment.