Simulated LMER
Dan Feeney
8/19/2021
setting up enviornment
# https://www.tjmahr.com/plotting-partial-pooling-in-mixed-effects-models/ provided a lot o
# basic code and inspiration for this. While journals won't allow citing a blog, this was a great
# resource
rm(list=ls())
library(tidyverse)## Warning: replacing previous import 'lifecycle::last_warnings' by
## 'rlang::last_warnings' when loading 'pillar'## Warning: replacing previous import 'lifecycle::last_warnings' by
## 'rlang::last_warnings' when loading 'tibble'## Warning: replacing previous import 'lifecycle::last_warnings' by
## 'rlang::last_warnings' when loading 'hms'## -- Attaching packages --------------------------------------- tidyverse 1.3.1 --## v ggplot2 3.3.5 v purrr 0.3.4
## v tibble 3.1.2 v dplyr 1.0.7
## v tidyr 1.1.3 v stringr 1.4.0
## v readr 1.4.0 v forcats 0.5.1## -- Conflicts ------------------------------------------ tidyverse_conflicts() --
## x dplyr::filter() masks stats::filter()
## x dplyr::lag() masks stats::lag()library(lme4)## Loading required package: Matrix##
## Attaching package: 'Matrix'## The following objects are masked from 'package:tidyr':
##
## expand, pack, unpacklibrary(emmeans)
library(MuMIn)## Warning: package 'MuMIn' was built under R version 4.1.3library(patchwork)
set.seed(5280)Simulating data
# set fixed effect parameters. https://journals.sagepub.com/doi/full/10.1177/2515245920965119
# set seed
set.seed(5280)
# set fixed effect parameters. https://journals.sagepub.com/doi/full/10.1177/2515245920965119
beta_0 <- 800 # intercept; i.e., the grand mean. this would be avg pwr
beta_1 <- 50 # slope; the effect of category this would be difference in pwr
# set random effect parameters
tau_0 <- 50 # by-subject random intercept sd. We assume this comes from a normal dist with mean 0 and unknown SD. SD for random intercept by sub
omega_0 <- 1 # by-item random intercept sd for seated and standing. by condition random intercept
# set more random effect and error parameters
tau_1 <- 25 # by-subject random slope sd for random slopes
rho <- .3 # correlation between intercept and slope
sigma <- 30 # residual (error) sd
#But note that we are sampling two random effects for each subject s, a random intercept T0s and a random slope T1s. It is possible for these values to be positively or negatively correlated, in which case we should not sample them independently. For instance, perhaps people who are faster than average overall (negative random intercept) also show a smaller than average effect of the in-group/out-group manipulation (negative random slope) because they allocate less attention to the task. We can capture this by allowing for a small positive correlation between the two factors, rho, which we assign to be .2.
# overall model: RTsi=β0+T0s+O0i+(β1+T1s)Xi+esi. The response time for subject s on item i, RTsi, is decomposed into a population grand mean, β0; a by-subject random intercept, T0s; a by-item random interce
# set number of subjects and items
n_subj <- 16 # number of subjects
# pt, O0i; a fixed slope, β1; a by-subject random slope, T1s; and a trial-level residual, esi. Our data-generating process is fully determined by seven population parameters, all denoted by Greek letters: β0, β1, τ0, τ1, ρ, ω0, and σ (see Table 2). In the next section, we apply this data-generating process to simulate the sampling of subjects, items, and trials (encounters).
n_seated <- 6 # number of seated observations
n_standing <- 6 # number of standing observations
#We need to create a table listing each item i, which category it is in, and its random effect, O0i:
# simulate a sample of items
# total number of items = n_ingroup +n_outgroup
items <- data.frame(
item_id = seq_len(n_seated + n_standing),
category = rep(c("seated", "standing"), c(n_seated, n_standing)),
O_0i = rnorm(n = n_seated + n_seated, mean = 0, sd = omega_0)
)
# effect-code category. this encodes a predictor as to which category each sim belongs to. seated should be less pwr, standing higher
items$X_i <- recode(items$category, "seated" = -0.5, "standing" = +0.5)
#We will later multiply this effect-coded factor by the fixed effect of category (beta_1 = 50) to simulate data in which the powers differ by postrue
# simulate a sample of subjects
# calculate random intercept / random slope covariance
covar <- rho * tau_0 * tau_1
# put values into variance-covariance matrix
cov_mx <- matrix(c(tau_0^2, covar, covar, tau_1^2),
nrow = 2, byrow = TRUE)
# generate the by-subject random effects
subject_rfx <- MASS::mvrnorm(n = n_subj,
mu = c(T_0s = 0, T_1s = 0),
Sigma = cov_mx)
# combine with subject IDs
subjects <- data.frame(subj_id = seq_len(n_subj),subject_rfx)
# cross subject and item IDs; add an error term
# nrow(.) is the number of rows in the table
trials <- crossing(subjects, items) %>%
mutate(e_si = rnorm(nrow(.), mean = 0, sd = sigma)) %>%
select(subj_id, item_id, category, X_i, everything())
# calculate the response variable
dat_sim <- trials %>%
mutate(pwr = beta_0 + T_0s + O_0i + (beta_1 + T_1s) * X_i + e_si) %>%
select(subj_id, item_id, category, X_i, pwr)Data visualization & Figure 1
# Visualize the data
p1 <- ggplot(data = dat_sim, aes(x = pwr, fill = category)) + geom_density() +
facet_wrap(~subj_id) + xlab('Power (W)') + ylab('Density') + ggtitle('LMEM Approach') +
scale_fill_grey()
ggplot(data = dat_sim, aes(x = category, y = pwr, fill = category)) + geom_boxplot() +
geom_point(position=position_jitterdodge())
dat_sim%>%
group_by(subj_id, category)%>%
summarize(
meanPow = mean(pwr)
) %>%
ggplot(aes(x = category, y = meanPow)) +
geom_point(aes(color=as.factor(subj_id))) +
geom_line(aes(group = subj_id, color = as.factor(subj_id) )) +
theme_bw()## `summarise()` has grouped output by 'subj_id'. You can override using the
## `.groups` argument.
p2 <- dat_sim%>%
group_by(subj_id, category)%>%
summarize(
meanPow = mean(pwr)
) %>%
ggplot(aes(x = category, y = meanPow)) +
geom_point(aes(color=category)) +
facet_wrap(~subj_id) + ylab('Mean Power (W)') +
xlab('Posture') + ggtitle('Averaged Data Approach') +
scale_color_grey() ## `summarise()` has grouped output by 'subj_id'. You can override using the
## `.groups` argument.# Figure 1
p1 | p2
## loading unbalanced dat from python creation
## to remove: subj 1 half of trials, subj
dat_mess <- read.csv('C:/Users/daniel.feeney/OneDrive - Boa Technology Inc/Desktop/Rethinking Neuro Data Manuscript/Code/Data/messDat.csv')
dat_unbal <- read.csv('C:/Users/daniel.feeney/OneDrive - Boa Technology Inc/Desktop/Rethinking Neuro Data Manuscript/Code/Data/unbalDat.csv')create a LMER model
dat_sim$Subject <- dat_sim$subj_id
# fit a linear mixed-effects model to data
mod_sim <- lmer(pwr ~ category + (1 + category | subj_id),
data = dat_sim)
# derive output from lmer
a <- summary(mod_sim, corr = FALSE)
r_sq <- r.squaredGLMM(mod_sim)## Warning: 'r.squaredGLMM' now calculates a revised statistic. See the help page.nullMod <- lmer(pwr ~ (category|subj_id),
data = dat_sim)
aovRes <- anova(mod_sim, nullMod)## refitting model(s) with ML (instead of REML)aovRes$`Pr(>Chisq)`[2]## [1] 2.617158e-06better RM model taken into account
avgDat <- dat_sim %>%
group_by(subj_id, category) %>%
summarize( meanPwr = mean(pwr))## `summarise()` has grouped output by 'subj_id'. You can override using the
## `.groups` argument.t.test(meanPwr ~ category, data = avgDat, paired = TRUE)##
## Paired t-test
##
## data: meanPwr by category
## t = -6.6797, df = 15, p-value = 7.357e-06
## alternative hypothesis: true difference in means is not equal to 0
## 95 percent confidence interval:
## -68.01426 -35.10840
## sample estimates:
## mean of the differences
## -51.56133Visualizations of no pooling, complete pooling, and partial pooling
ggplot(data = dat_sim, mapping=aes(x=category, y=pwr, color = category))+
geom_point() + facet_wrap(~subj_id)
# fit a no pooling model
no_pooling_mod <- lmList(pwr ~ category | subj_id, dat_sim)## Warning: Unknown or uninitialised column: `(weights)`.## Warning: Unknown or uninitialised column: `(offset)`.## Warning: Unknown or uninitialised column: `(weights)`.## Warning: Unknown or uninitialised column: `(offset)`.## Warning: Unknown or uninitialised column: `(weights)`.## Warning: Unknown or uninitialised column: `(offset)`.## Warning: Unknown or uninitialised column: `(weights)`.## Warning: Unknown or uninitialised column: `(offset)`.## Warning: Unknown or uninitialised column: `(weights)`.## Warning: Unknown or uninitialised column: `(offset)`.## Warning: Unknown or uninitialised column: `(weights)`.## Warning: Unknown or uninitialised column: `(offset)`.## Warning: Unknown or uninitialised column: `(weights)`.## Warning: Unknown or uninitialised column: `(offset)`.## Warning: Unknown or uninitialised column: `(weights)`.## Warning: Unknown or uninitialised column: `(offset)`.## Warning: Unknown or uninitialised column: `(weights)`.## Warning: Unknown or uninitialised column: `(offset)`.## Warning: Unknown or uninitialised column: `(weights)`.## Warning: Unknown or uninitialised column: `(offset)`.## Warning: Unknown or uninitialised column: `(weights)`.## Warning: Unknown or uninitialised column: `(offset)`.## Warning: Unknown or uninitialised column: `(weights)`.## Warning: Unknown or uninitialised column: `(offset)`.## Warning: Unknown or uninitialised column: `(weights)`.## Warning: Unknown or uninitialised column: `(offset)`.## Warning: Unknown or uninitialised column: `(weights)`.## Warning: Unknown or uninitialised column: `(offset)`.## Warning: Unknown or uninitialised column: `(weights)`.## Warning: Unknown or uninitialised column: `(offset)`.## Warning: Unknown or uninitialised column: `(weights)`.## Warning: Unknown or uninitialised column: `(offset)`.df_no_pooling <- lmList(pwr ~ category | subj_id, dat_sim) %>%
coef() %>%
# Subject IDs are stored as row-names. Make them an explicit column
rownames_to_column("Subject") %>%
rename(Intercept = `(Intercept)`, Slope_Cond = categorystanding) %>%
add_column(Model = "No pooling")## Warning: Unknown or uninitialised column: `(weights)`.
## Warning: Unknown or uninitialised column: `(offset)`.## Warning: Unknown or uninitialised column: `(weights)`.## Warning: Unknown or uninitialised column: `(offset)`.## Warning: Unknown or uninitialised column: `(weights)`.## Warning: Unknown or uninitialised column: `(offset)`.## Warning: Unknown or uninitialised column: `(weights)`.## Warning: Unknown or uninitialised column: `(offset)`.## Warning: Unknown or uninitialised column: `(weights)`.## Warning: Unknown or uninitialised column: `(offset)`.## Warning: Unknown or uninitialised column: `(weights)`.## Warning: Unknown or uninitialised column: `(offset)`.## Warning: Unknown or uninitialised column: `(weights)`.## Warning: Unknown or uninitialised column: `(offset)`.## Warning: Unknown or uninitialised column: `(weights)`.## Warning: Unknown or uninitialised column: `(offset)`.## Warning: Unknown or uninitialised column: `(weights)`.## Warning: Unknown or uninitialised column: `(offset)`.## Warning: Unknown or uninitialised column: `(weights)`.## Warning: Unknown or uninitialised column: `(offset)`.## Warning: Unknown or uninitialised column: `(weights)`.## Warning: Unknown or uninitialised column: `(offset)`.## Warning: Unknown or uninitialised column: `(weights)`.## Warning: Unknown or uninitialised column: `(offset)`.## Warning: Unknown or uninitialised column: `(weights)`.## Warning: Unknown or uninitialised column: `(offset)`.## Warning: Unknown or uninitialised column: `(weights)`.## Warning: Unknown or uninitialised column: `(offset)`.## Warning: Unknown or uninitialised column: `(weights)`.## Warning: Unknown or uninitialised column: `(offset)`.## Warning: Unknown or uninitialised column: `(weights)`.## Warning: Unknown or uninitialised column: `(offset)`.df_no_pooling$Subject <- as.integer(df_no_pooling$Subject)
#fit a complete pooling mode
m_pooled <- lm(pwr ~ category, dat_sim)
# Repeat the intercept and slope terms for each participant
df_pooled <- tibble(
Model = "Complete pooling",
Subject = unique(dat_sim$subj_id),
Intercept = coef(m_pooled)[1],
Slope_Cond = coef(m_pooled)[2]
)
head(df_pooled)## # A tibble: 6 x 4
## Model Subject Intercept Slope_Cond
## <chr> <int> <dbl> <dbl>
## 1 Complete pooling 1 778. 51.6
## 2 Complete pooling 2 778. 51.6
## 3 Complete pooling 3 778. 51.6
## 4 Complete pooling 4 778. 51.6
## 5 Complete pooling 5 778. 51.6
## 6 Complete pooling 6 778. 51.6compare models. No pooling is RM ANOVA, complete pooling is LM
# https://www.linguisticsociety.org/sites/default/files/e-learning/class_3_slides.pdf
# Join the raw data so we can use plot the points and the lines.
df_models <- bind_rows(df_pooled, df_no_pooling) %>%
left_join(dat_sim, by = "Subject")
p_model_comparison <- ggplot(df_models) +
aes(x = category, y = pwr) + ylab('Power (W)') +
# Set the color mapping in this layer so the points don't get a color
geom_abline(
aes(intercept = Intercept, slope = Slope_Cond, color = Model, linetype = Model),
size = .75
) +
geom_point() +
facet_wrap("Subject") +
theme(legend.position = "top", legend.justification = "left") +
scale_color_grey()
p_model_comparison 
#mod_sim
df_partial_pooling <- coef(mod_sim)[["subj_id"]] %>%
rownames_to_column("Subject") %>%
as_tibble() %>%
rename(Intercept = `(Intercept)`, Slope_Cond = categorystanding) %>%
add_column(Model = "Partial pooling")
df_partial_pooling$Subject <- as.integer(df_partial_pooling$Subject)
df_models <- bind_rows(df_pooled, df_no_pooling, df_partial_pooling) %>%
left_join(dat_sim, by = "Subject")
# Replace the data-set of the last plot
part1 <- p_model_comparison %+% df_modelsFigure 3 and visualizing all 3 models to show shrinkage
df_fixef <- tibble(
Model = "Partial pooling (average)",
Intercept = fixef(mod_sim)[1],
Slope_Cond = fixef(mod_sim)[2]
)
df_shrink <- df_pooled %>%
distinct(Model, Intercept, Slope_Cond) %>%
bind_rows(df_fixef)
df_shrink## # A tibble: 2 x 3
## Model Intercept Slope_Cond
## <chr> <dbl> <dbl>
## 1 Complete pooling 778. 51.6
## 2 Partial pooling (average) 778. 51.6df_pulled <- bind_rows(df_no_pooling, df_partial_pooling)
part2 <- ggplot(df_pulled) +
aes(x = Intercept, y = Slope_Cond, color = Model, shape = Model) +
geom_point(size = 2) +
geom_point(
data = df_shrink,
size = 5,
show.legend = FALSE
) +
# Draw an arrow connecting the observations between models
geom_path(
aes(group = Subject, color = NULL),
arrow = arrow(length = unit(.02, "npc")),
show.legend = FALSE
) +
#ggtitle("Pooling of regression parameters") +
xlab("Intercept estimate") +
ylab("Slope estimate") +
scale_color_grey() +
theme(legend.position="top")
part1 | part2
Comparing effect size estimats
print("Estimated Effect Sizes")## [1] "Estimated Effect Sizes"est_slope = df_pooled$Slope_Cond
df_pooled$Subject## [1] 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16ggplot(df_pooled) +
aes(x=Intercept, y=Slope_Cond) +
geom_point()
print(paste("Pooled model = ", mean(est_slope) / sd(est_slope))) #inf because sd of a single number is 0 ## [1] "Pooled model = Inf"est_slope_noPool = df_no_pooling$Slope_Cond
df_no_pooling$Subject## [1] 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16ggplot(df_no_pooling) +
aes(x=Intercept, y=Slope_Cond, label = Subject) +
geom_point()+ geom_text(hjust = 1, vjust = 1)
print(paste("No-pooling model = ", mean(est_slope_noPool) / sd(est_slope_noPool)))## [1] "No-pooling model = 1.66992101399288"est_slope_partialPool = df_partial_pooling$Slope_Cond
df_partial_pooling$Subject## [1] 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16ggplot(df_partial_pooling) +
aes(x=Intercept, y=Slope_Cond, label = Subject) +
geom_point() + geom_text(hjust = 1, vjust = 1)
print(paste("Partial-pooling model = ", mean(est_slope_partialPool) / sd(est_slope_partialPool)))## [1] "Partial-pooling model = 2.2122094573699"a <- avgDat %>% filter(category=="seated")
a <- a$meanPwr
b <- avgDat %>% filter(category=="standing")
b <- b$meanPwr
print(paste("Paired t-test = ", mean(b-a) / sd(b-a)))## [1] "Paired t-test = 1.66992101399289"plotting in R
df_all <- bind_rows(df_pooled, df_no_pooling, df_partial_pooling)
ggplot(df_all) +
aes(x = Subject, y = Slope_Cond, color=Model, label = Subject) +
geom_point(size = 2) +
geom_point(
data = df_all,
size = 2,
# Prevent size-5 point from showing in legend keys
show.legend = FALSE
) +
geom_text(vjust=1, hjust=1)+
scale_x_continuous(breaks=seq(0,15,1)) +
theme(panel.grid.minor = element_blank()) +
facet_grid(cols = vars(Model))
Appendix
add another fixed effect to show the model’s ability to generalize
dat_sim$footwear <- rep(c(rep('shoe1',3), rep('shoe2',3)),32)
ggplot(dat_sim, aes(x = category, y = pwr, color = footwear, fill = footwear)) +
geom_point() + geom_jitter(width = 0.1, height = 0.1) + facet_wrap(~Subject) +
scale_color_manual(values=c("red","blue")) + ylab("Power (W)") + xlab("Posture")
twoMod <- lmer(pwr ~ category + footwear + (category|Subject), data = dat_sim)
summary(twoMod)## Linear mixed model fit by REML ['lmerMod']
## Formula: pwr ~ category + footwear + (category | Subject)
## Data: dat_sim
##
## REML criterion at convergence: 1867.2
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -2.5310 -0.6474 0.1187 0.5417 2.6481
##
## Random effects:
## Groups Name Variance Std.Dev. Corr
## Subject (Intercept) 2776.9 52.70
## categorystanding 718.3 26.80 -0.24
## Residual 705.3 26.56
## Number of obs: 192, groups: Subject, 16
##
## Fixed effects:
## Estimate Std. Error t value
## (Intercept) 778.0801 13.5859 57.271
## categorystanding 51.5613 7.7192 6.680
## footwearshoe2 -0.8119 3.8333 -0.212
##
## Correlation of Fixed Effects:
## (Intr) ctgrys
## ctgrystndng -0.275
## footwearsh2 -0.141 0.000conduct F-test between anova models. In this case, there is no estimated main effect of footwear but there is a main effect of category (seated or standing)
fwMod <- lmer(pwr ~ category + footwear + (category|Subject), data = dat_sim)
baseMod <- lmer(pwr ~ category + (category|Subject), data = dat_sim)
anova(fwMod, baseMod)## refitting model(s) with ML (instead of REML)## Data: dat_sim
## Models:
## baseMod: pwr ~ category + (category | Subject)
## fwMod: pwr ~ category + footwear + (category | Subject)
## npar AIC BIC logLik deviance Chisq Df Pr(>Chisq)
## baseMod 6 1896.6 1916.1 -942.29 1884.6
## fwMod 7 1898.5 1921.3 -942.27 1884.5 0.0451 1 0.8318catMod <- lmer(pwr ~ category + footwear + (category|Subject), data = dat_sim)
baseMod <- lmer(pwr ~ footwear + (category|Subject), data = dat_sim)
anova(fwMod, baseMod)## refitting model(s) with ML (instead of REML)## Data: dat_sim
## Models:
## baseMod: pwr ~ footwear + (category | Subject)
## fwMod: pwr ~ category + footwear + (category | Subject)
## npar AIC BIC logLik deviance Chisq Df Pr(>Chisq)
## baseMod 6 1918.6 1938.2 -953.31 1906.6
## fwMod 7 1898.5 1921.3 -942.27 1884.5 22.079 1 2.617e-06 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1add interaction effects and plot them
library(sjPlot)## Warning: package 'sjPlot' was built under R version 4.1.2intMod <- lmer(pwr ~ category * footwear + (category|Subject), data = dat_sim)
plot_model(intMod) 
plot_model(intMod, type = 'pred', terms = c("category", "footwear"))
summary(intMod)## Linear mixed model fit by REML ['lmerMod']
## Formula: pwr ~ category * footwear + (category | Subject)
## Data: dat_sim
##
## REML criterion at convergence: 1858.8
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -2.4280 -0.6353 0.1279 0.5850 2.5472
##
## Random effects:
## Groups Name Variance Std.Dev. Corr
## Subject (Intercept) 2778.0 52.71
## categorystanding 720.5 26.84 -0.24
## Residual 698.7 26.43
## Number of obs: 192, groups: Subject, 16
##
## Fixed effects:
## Estimate Std. Error t value
## (Intercept) 781.106 13.718 56.941
## categorystanding 45.509 8.610 5.285
## footwearshoe2 -6.864 5.395 -1.272
## categorystanding:footwearshoe2 12.104 7.630 1.586
##
## Correlation of Fixed Effects:
## (Intr) ctgrys ftwrs2
## ctgrystndng -0.306
## footwearsh2 -0.197 0.313
## ctgrystnd:2 0.139 -0.443 -0.707Compare to only a linear model (worst choice)
# complete pooling
wrongMod <- lm(pwr ~ category, data = dat_sim)
wrongSum <- summary(wrongMod)
wrongSum$coefficients[8]## [1] 3.937432e-09wrongSum##
## Call:
## lm(formula = pwr ~ category, data = dat_sim)
##
## Residuals:
## Min 1Q Median 3Q Max
## -127.251 -46.277 1.735 40.313 128.029
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 777.674 5.905 131.700 < 2e-16 ***
## categorystanding 51.561 8.351 6.174 3.94e-09 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 57.86 on 190 degrees of freedom
## Multiple R-squared: 0.1671, Adjusted R-squared: 0.1627
## F-statistic: 38.12 on 1 and 190 DF, p-value: 3.937e-09# repeated measures error model is better #
aovMod <- aov( pwr ~ category + Error(factor(subj_id)), data = dat_sim )#
RMsum <- summary(aovMod)
RMsum$`Error: Within`## Df Sum Sq Mean Sq F value Pr(>F)
## category 1 127611 127611 144 <2e-16 ***
## Residuals 175 155079 886
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1summary(aovMod)##
## Error: factor(subj_id)
## Df Sum Sq Mean Sq F value Pr(>F)
## Residuals 15 480903 32060
##
## Error: Within
## Df Sum Sq Mean Sq F value Pr(>F)
## category 1 127611 127611 144 <2e-16 ***
## Residuals 175 155079 886
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1