Skip to content

Instantly share code, notes, and snippets.

@arcaldwell49
Created October 8, 2021 18:06
Show Gist options
  • Star 0 You must be signed in to star a gist
  • Fork 0 You must be signed in to fork a gist
  • Save arcaldwell49/fe54ba23c0e08f030d2b9949edfe2400 to your computer and use it in GitHub Desktop.
Save arcaldwell49/fe54ba23c0e08f030d2b9949edfe2400 to your computer and use it in GitHub Desktop.
Simulate multilevel brain data
# Generate data with simstudy ------
library(simstudy)
library(tidyverse)
set.seed(10082021)
# Brain level effects ----
## generate random intercept at brain -----
gen.brain <- defData(varname = "bi0",
dist = "normal",
formula = 0, variance = 3,
id = "idBrain")
## brain random slope at brian -------
gen.brain <- defData(gen.brain, varname = "bs0",
dist = "normal",
formula = 0, variance = 1)
## generate number of wells (n = 8) -----
gen.brain <- defData(gen.brain, varname = "nWell",
dist = "nonrandom", formula = 8)
## generate data with 12 brains ----
dtbrain <- genData(12, gen.brain)
# check -----
head(dtbrain, 3)
# Generate data at well level ------
## random intercept at well ------
gen.well <- defDataAdd(varname = "wi0",
dist = "normal",
formula = 0, variance = 2)
# Number of observations per well (n=6) ------
gen.well <- defDataAdd(gen.well, varname = "nObs",
dist = "nonrandom",
formula = 6)
# Generate well/clusters ----
dtWell <- genCluster(dtbrain, "idBrain",
numIndsVar = "nWell",
level1ID = "idWell")
dtWell <- addColumns(gen.well, dtWell)
## check again -----
head(dtWell, 5)
# Assign 2 treatments -----
dtWell <- trtAssign(dtWell,
nTrt=2,
balanced = TRUE,
strata = "idBrain")
# Outcome at cell level ------
gen.cell <- defDataAdd(varname = "y",
dist = "normal",
# formula equal 50 + random intercepts - (1.5 + random slopes) times treat group
formula = "(50+bi0+wi0) - (1.5+bs0) * trtGrp",
variance = 2)
## Generate cell/clusters ----
dtCell <- genCluster(dtWell,
cLevelVar = "idWell",
numIndsVar = "nObs",
level1ID = "idCell")
# One data set ----
dtCell <- addColumns(gen.cell, dtCell)
# Modify dataset ------
brain2 = dtCell %>%
select(idBrain,idWell, trtGrp, idCell, y) %>%
mutate(idBrain = as.factor(idBrain),
treat = factor(trtGrp,
level = c(0,1),
labels = c("Control", "Treatment")),
idWell = as.factor(idWell))
# Check data again ----
head(brain2)
# ANALYSIS OPTIONS ----
## Linear Mixed Model with lmerTest/lme4 -----
library(lmerTest)
fit1 = lme4::lmer(y ~ trtGrp + (1+trtGrp|idBrain) + (1|idBrain:idWell),
brain2)
anova(fit1)
# emmeans::emmeans(fit1, pairwise ~ trtGrp) # pairwise
### Variance components ----
library(mixedup)
extract_vc(fit1) # trtGrp v. small prop of variance
## Analysis of Variance -----
# Fails to compute....
fit2 = aov(y ~ treat + Error(idBrain/idWell), # nesting structure
data = brain2)
# emmeans::emmeans(fit2, pairwise ~ trtGrp) # pairwise
# Fit with nlme as an alt. to fit1 ----
## To make it different let's drop the well component...
library(nlme)
fit3 = lme(y ~ trtGrp,
random = list(idBrain = ~1+trtGrp),
data= brain2, method = "REML")
anova(fit3)
extract_vc(fit3)
# Visualize -----
library(ggbeeswarm)
# Get brain level averages
ReplicateAverages <- brain2 %>%
select(idBrain, treat, y) %>%
group_by(treat, idBrain) %>%
summarise_each(list(mean))
ggplot(brain2,
aes(x=1,
y=y,
color=idBrain)) +
facet_wrap(~treat) +
geom_beeswarm(cex=3, alpha = .5) +
scale_color_viridis_d() +
geom_beeswarm(data=ReplicateAverages, size=5, dodge.width = .1) +
theme_bw() +
theme(legend.position="none",
axis.text.x=element_blank(),
axis.ticks.x=element_blank(),
axis.title.x=element_blank())
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment