Created
October 8, 2021 18:06
-
-
Save arcaldwell49/fe54ba23c0e08f030d2b9949edfe2400 to your computer and use it in GitHub Desktop.
Simulate multilevel brain data
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
# 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