From e483ad1c9e93243593b8922bb14b42bb2e20e4b8 Mon Sep 17 00:00:00 2001 From: Steve Phelps Date: Thu, 21 Dec 2023 10:06:54 +0000 Subject: [PATCH] interaction terms and effect sizes for GLMM --- environment.yml | 2 + jupyter-book/R-MixedModel.py | 154 ++++++++++++++++++++++++++++++++++- 2 files changed, 154 insertions(+), 2 deletions(-) diff --git a/environment.yml b/environment.yml index 647c324..f03592e 100644 --- a/environment.yml +++ b/environment.yml @@ -42,6 +42,8 @@ dependencies: - r-modelsummary - r-texreg - r-MCMCglmm + - r-effects + - r-sjPlot - pip - pip: - -e ./ diff --git a/jupyter-book/R-MixedModel.py b/jupyter-book/R-MixedModel.py index 284dd32..3ccedf2 100644 --- a/jupyter-book/R-MixedModel.py +++ b/jupyter-book/R-MixedModel.py @@ -25,6 +25,12 @@ library(dplyr) library(xtable) library(texreg) +library(geepack) +library(MCMCglmm) +library(parallel) +library(coda) +library(effects) +library(ggeffects) # %% options(repr.plot.width = 20, repr.plot.height = 10) @@ -161,28 +167,106 @@ # %% model.pd <- glmmTMB(cbind(Num_cooperates, 6 - Num_cooperates) ~ - Participant_group + Partner_condition + t + Model + Temperature + - Partner_condition:Model + Participant_group:Model + + Participant_group + Participant_group:Partner_condition + t + Model + Temperature + (1|Participant_id), data = results.pd, family = betabinomial) summary(model.pd) +# %% +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), + data = results.pd, + family = betabinomial) +summary(model.pd.1) + # %% texreg(model.pd, caption="Fitted model for Prisoners Dilemma", label="table:pd-estimates", file="pd-estimates.tex", single.row=TRUE, fontsize="small") +# %% + # %% 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)' +m + +# %% +m[m[, 4] <= 0.05, ] + +# %% + # %% xtable(coef(summary(model.pd))$cond, digits=3) +# %% +plot(ggpredict(model.pd, c("Participant_group", "Model"))) + +# %% +plot(ggpredict(model.pd, c("Partner_condition", "Model"))) + +# %% +predicted <- ggpredict(model.pd.1, c("Participant_group", "Model", "Partner_condition")) +predicted + +# %% +plot(predicted) + +# %% +predicted <- ggpredict(model.pd.1, "Temperature") +predicted + +# %% +plot(allEffects(model.pd)) + +# %% +options(repr.plot.height=12) +simulationOutput <- simulateResiduals(fittedModel = model.pd, plot = TRUE, integerResponse=TRUE) + # %% options(repr.plot.height=10) +pdf("residuals-plots.pdf") simulationOutput <- simulateResiduals(fittedModel = model.pd, plot = TRUE, integerResponse=TRUE) +dev.off() # %% +pdf("residuals-hist.pdf") hist(residuals(simulationOutput)) +dev.off() + +# %% +# Fit a GEE model with an interaction term +gee_model <- geeglm(cbind(Num_cooperates, 6 - Num_cooperates) ~ Participant_group + Partner_condition + t + Model + Temperature + + Partner_condition:Model + Participant_group:Model, + id = Participant_id, + data = results.pd, + family = binomial, + corstr = "exchangeable") + +# %% +summary(gee_model) + +# %% +residuals_gee <- residuals(gee_model, type = "pearson") +# Plotting Pearson residuals against fitted values +fitted_values <- fitted(gee_model) +plot(fitted_values, residuals_gee, xlab = "Fitted Values", ylab = "Pearson Residuals") +abline(h = 0, col = "red") + + +# %% +# Overdispersion check +overdispersion_check <- sum(residuals_gee^2) / gee_model$df.residual +overdispersion_check # %% model.pd.factorial <- glmmTMB(cbind(Num_cooperates, 6 - Num_cooperates) ~ @@ -201,3 +285,69 @@ # %% plotResiduals(simulationOutput.factorial, results.clean$Participant_group) + +# %% [markdown] +# ## Bayesian MCMC GLMM + +# %% +# Set up a list for priors +prior <- list( + R = list(V = 1, nu = 0.002), # Prior for the residual covariance matrix + G = list(G1 = list(V = 1, nu = 0.002)) # Prior for the random effects +) + +estimate_glmm <- function() { +# Fit the model with the specified priors + MCMCglmm(fixed=Num_cooperates ~ Participant_group + Partner_condition + t + Model + Model:Partner_condition:Participant_group, + random=~ us(1):Participant_id, + prior = prior, + data = results.pd, + verbose = FALSE) +} + +model.pd.mcmc <- estimate_glmm() + +# %% +plot(model.pd.mcmc$Sol) + +# %% +plot(model.pd.mcmc$VCV) + +# %% +m6 <- mclapply(1:12, function(i) { estimate_glmm() }, mc.cores=12) +m6 <- lapply(m6, function(m) m$Sol) +m6 <- do.call(mcmc.list, m6) + +# %% +par(mfrow=c(4,2), mar=c(2,2,1,2)) +gelman.plot(m6, auto.layout=F) + +# %% +gelman.diag(m6) + +# %% +par(mfrow=c(8,2), mar=c(2, 1, 1, 1)) +plot(m6, ask=F, auto.layout=F) + +# %% +summary(m6) + +# %% +plot.estimates <- function(x) { + if (class(x) != "summary.mcmc") + x <- summary(x) + n <- dim(x$statistics)[1] + par(mar=c(2, 7, 4, 1)) + plot(x$statistics[,1], n:1, + yaxt="n", ylab="", + xlim=range(x$quantiles)*1.2, + pch=19, + main="Posterior means and 95% credible intervals") + grid() + axis(2, at=n:1, rownames(x$statistics), las=2) + arrows(x$quantiles[,1], n:1, x$quantiles[,5], n:1, code=0) + abline(v=0, lty=2) +} + +options(repr.plot.height=24, repr.plot.width=20) +plot.estimates(m6)