rm(list = ls())
suppressWarnings(suppressMessages( {
library(brms)
library(beeswarm)
library(vioplot)
library(glmmTMB)
library(car)
library(cowplot)
library(ggplot2)
library(tidyverse)
library(bayesplot)
library(gtsummary)
library(vegan)
library(knitr)
library(ggplot2)
} ) )Experimental treatment with edaravone in a mouse model of spinocerebellar ataxia 1
Analysis of motor deficit on rotarod
This page shows R code for the study Sucha et al. (2023, International Journal of Molecular Science).
Citation:
Sucha M, Benediktova S, Tichanek F, Jedlicka J, Kapl S, Jelinkova D, Purkartova Z, Tuma J, Kuncova J, Cendelin J. Experimental Treatment with Edaravone in a Mouse Model of Spinocerebellar Ataxia 1. International Journal of Molecular Sciences. 2023; 24(13):10689. https://doi.org/10.3390/ijms241310689
GitHub page: https://github.com/filip-tichanek/edaravonSCA1
1 Upload of packages
2 Data upload and wrangling
load("source_data/myEnvironment.RData")
## defining new variables
rotarod$group <- interaction(rotarod$treatment, rotarod$genotype)
rotarod$group <- factor(rotarod$group,levels=c('ctrl.wt','eda.wt','ctrl.sca','eda.sca'))
rotarod$sex_num <- ifelse(rotarod$sex == "f", -0.5, 0.5)
rotarod$group_sex <- interaction(rotarod$treatment, rotarod$genotype, rotarod$sex)
rotarod$group_sex <- factor(rotarod$group_sex, levels=c('ctrl.wt.f','ctrl.wt.m',
'eda.wt.f','eda.wt.m',
'ctrl.sca.f','ctrl.sca.m',
'eda.sca.f','eda.sca.m'))
## subject as factor
rotarod$id <- as.factor(rotarod$id)
## order of all roatrod tests
rotarod$order <- (rotarod$day-1) * 4 + rotarod$session
rotarod$order2 <- ifelse(rotarod$pre_post == 'pre', rotarod$order, rotarod$order+12)
## data summary
summary(rotarod) id treatment genotype sex day session
7 : 24 ctrl:960 sca:960 f:936 Min. :1 Min. :1.00
9 : 24 eda :960 wt :960 m:984 1st Qu.:1 1st Qu.:1.75
11 : 24 Median :2 Median :2.50
12 : 24 Mean :2 Mean :2.50
13 : 24 3rd Qu.:3 3rd Qu.:3.25
14 : 24 Max. :3 Max. :4.00
(Other):1776
rot_latency pre_post group sex_num group_sex
Min. : 2.30 post:960 ctrl.wt :480 Min. :-0.5000 ctrl.wt.m :264
1st Qu.: 94.97 pre :960 eda.wt :480 1st Qu.:-0.5000 eda.wt.f :240
Median :126.90 ctrl.sca:480 Median : 0.5000 eda.wt.m :240
Mean :127.79 eda.sca :480 Mean : 0.0125 ctrl.sca.f:240
3rd Qu.:157.60 3rd Qu.: 0.5000 ctrl.sca.m:240
Max. :299.20 Max. : 0.5000 eda.sca.f :240
(Other) :456
order order2
Min. : 1.00 Min. : 1.00
1st Qu.: 3.75 1st Qu.: 6.75
Median : 6.50 Median :12.50
Mean : 6.50 Mean :12.50
3rd Qu.: 9.25 3rd Qu.:18.25
Max. :12.00 Max. :24.00
Sub-setting pre-treatment and post-treatment rotarod data and data aggregation.
We will create 3 aggregated datasets, each with different pre-treatment rotarod latency (preTreat RotLat) variable. The reason is that our previous observations suggested that the genotype-related difference is less pronounced for the first 1-2 rotarod session(s) and do not predict later rotarod well. Models will be fitted using all three datasets and the model with the best predictive accuracy (as estimated via leave-one-out cross-validation - LOO-CV) will be used
## subseting post-treatment dataset
rot_post <- subset(rotarod, rotarod$pre_post == 'post')
## pre-treatment rotarod
rot_pre3 <- subset(rotarod, rotarod$pre_post == 'pre'& rotarod$day >2)
rot_pre23 <- subset(rotarod,rotarod$pre_post == 'pre'& rotarod$day >1)
rot_pre <- subset(rotarod, rotarod$pre_post == 'pre')
## aggregated data sets. Each has different definition of pre-treatment rotarod:
### mean from all pre-treatment rotarod runs
rot_pre_agg <- rot_pre %>%
group_by(id) %>%
summarise(mean_lat = mean(rot_latency))
### mean from 2nd and 3rd day
rot_pre_agg_23 <- rot_pre23 %>%
group_by(id) %>%
summarise(mean_lat = mean(rot_latency))
### mean from 3rd rotarod day only
rot_pre_agg_3 <- rot_pre3 %>%
group_by(id) %>%
summarise(mean_lat = mean(rot_latency))
## joining new variable to post-treatment rotarod data
rot_post_co <- rot_post %>%
left_join(rot_pre_agg, by ='id') %>%
left_join(rot_pre_agg_23, by ='id') %>%
left_join(rot_pre_agg_3, by ='id')
## centering the variables 'day' and 'session'
rot_post_co$day_cen <- rot_post_co$day-2
rot_post_co$session_cen <- rot_post_co$session -2.5
summary(rot_post_co) id treatment genotype sex day session
7 : 12 ctrl:480 sca:480 f:468 Min. :1 Min. :1.00
9 : 12 eda :480 wt :480 m:492 1st Qu.:1 1st Qu.:1.75
11 : 12 Median :2 Median :2.50
12 : 12 Mean :2 Mean :2.50
13 : 12 3rd Qu.:3 3rd Qu.:3.25
14 : 12 Max. :3 Max. :4.00
(Other):888
rot_latency pre_post group sex_num group_sex
Min. : 16.10 post:960 ctrl.wt :240 Min. :-0.5000 ctrl.wt.m :132
1st Qu.: 93.75 pre : 0 eda.wt :240 1st Qu.:-0.5000 eda.wt.f :120
Median :125.75 ctrl.sca:240 Median : 0.5000 eda.wt.m :120
Mean :126.67 eda.sca :240 Mean : 0.0125 ctrl.sca.f:120
3rd Qu.:155.03 3rd Qu.: 0.5000 ctrl.sca.m:120
Max. :291.20 Max. : 0.5000 eda.sca.f :120
(Other) :228
order order2 mean_lat.x mean_lat.y
Min. : 1.00 Min. :13.00 Min. : 55.98 Min. : 55.67
1st Qu.: 3.75 1st Qu.:15.75 1st Qu.:108.61 1st Qu.:117.49
Median : 6.50 Median :18.50 Median :127.36 Median :136.45
Mean : 6.50 Mean :18.50 Mean :128.90 Mean :140.28
3rd Qu.: 9.25 3rd Qu.:21.25 3rd Qu.:147.82 3rd Qu.:166.12
Max. :12.00 Max. :24.00 Max. :195.66 Max. :218.11
mean_lat day_cen session_cen
Min. : 60.2 Min. :-1 Min. :-1.50
1st Qu.:121.4 1st Qu.:-1 1st Qu.:-0.75
Median :141.0 Median : 0 Median : 0.00
Mean :146.2 Mean : 0 Mean : 0.00
3rd Qu.:173.2 3rd Qu.: 1 3rd Qu.: 0.75
Max. :255.6 Max. : 1 Max. : 1.50
New variables mean_lat.x, mean_lat.y, mean_lat indicate pre-treatment rotarod latency (preTreat RotLat) used as covariate in models to adjuste for baseline performance
- mean_lat.x: average latency from all pre-treatment measurements
- mean_lat.y: latency averaged from the days 2 and 3
- mean_lat latency averaged from the 3rd day
3 Data exploration
Figure showing individual trajectories or rotarod (all 24 sessions, blue vertical line shows transition from pre-treatment to post-treatment)
rotarod %>% ggplot(aes(x=order2, y=rot_latency, by=id, col=sex)) +
geom_line(alpha=0.8) +
geom_vline(xintercept = 12.5,
color = "blue", linewidth=0.8, linetype= 'dashed') +
facet_wrap(~group) 
Sex does not exhibit consistent effect, or only minor specifically in WT group. We will ignore the sex throughout the analysis
Figure showing data distribution per group
par(mfrow=c(2,2))
hist(rotarod[which(rotarod$group == 'ctrl.wt'),c('rot_latency')], 20, main="WT_0", xlim=c(0,300))
hist(rotarod[which(rotarod$group == 'eda.wt'),c('rot_latency')], 20, main="WT_E", xlim=c(0,300))
hist(rotarod[which(rotarod$group == 'ctrl.sca'),c('rot_latency')], 20,
main="SCA_0", xlim=c(0,300))
hist(rotarod[which(rotarod$group == 'eda.sca'),c('rot_latency')], 20, main="SCA_E", xlim=c(0,300))
Distributions seem more or less symmetric. I do not see marks of heteroscedasticity (width of distribution seems similar across groups, similar magnitude of changes in low- and high-performing individuals on the previous plot of individual trajectories).
4 Modelling
As rotarod latency in non-negative, and right tail-distribution is expected, gamma likelihood may be the first choice. On the other hand, the values of latency is restricted as the rotarod increases the speed; the latency is thus strongly restricted, resulting in normal-like distribution.
Also data exploration suggest that Gaussian likelihood distribution may be reasonable choice here. We will fit both and look on posterior predictive check (PPC) to choose between them.
Firstly, we will select likelihood used, with models including pre-treatment roator latency averaged over all pre-treatment rotarod measurements (mean_lat.x).
4.1 Prior specification
First, we look at data summary statistics and estimate from previous publication Tichanek et al.
## changing reference group
rot_post_co$group<-relevel(rot_post_co$group,ref='ctrl.sca')
## wrangling data from Tichanek et al.
young_scirep <- young_scirep %>% mutate(
rotarod_latency = mean(grepl('rotarod', young_scirep)))
young_scirep$rot_lat <- rowMeans(young_scirep[
,which(grepl("rotarod", colnames(young_scirep)))])
## prior for genotype effect, using Gaussian likelihood
mprio_norm <- glmmTMB(rot_lat~genotype, data = young_scirep)
summary(mprio_norm) Family: gaussian ( identity )
Formula: rot_lat ~ genotype
Data: young_scirep
AIC BIC logLik deviance df.resid
569.5 575.6 -281.8 563.5 52
Dispersion estimate for gaussian family (sigma^2): 1.65e+03
Conditional model:
Estimate Std. Error z value Pr(>|z|)
(Intercept) 134.771 7.416 18.174 <2e-16 ***
genotypew 24.259 10.999 2.206 0.0274 *
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1prios(mprio_norm) Estimate
2.701332 32.337246 ## prior for genotype effect with Gamma likelihood
mprio_norm <- glmmTMB(rot_lat~genotype, data = young_scirep,family=Gamma(link="log"))
summary(mprio_norm) Family: Gamma ( log )
Formula: rot_lat ~ genotype
Data: young_scirep
AIC BIC logLik deviance df.resid
567.2 573.2 -280.6 561.2 52
Dispersion estimate for Gamma family (sigma^2): 0.079
Conditional model:
Estimate Std. Error z value Pr(>|z|)
(Intercept) 4.90358 0.05133 95.53 <2e-16 ***
genotypew 0.16552 0.07613 2.17 0.0297 *
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1prios(mprio_norm) Estimate
0.01629459 0.22383596 ## summary statistics
mean(rot_post_co$rot_latency)[1] 126.674mean(rot_post_co[rot_post_co$genotype == 'sca', c('rot_latency')])[1] 109.2742sd(rot_post_co$rot_latency) *c(5, 2, 1.2)[1] 229.65135 91.86054 55.11632sd(rot_post_co[rot_post_co$genotype == 'sca', c('rot_latency')])*c(5, 2, 1.2)[1] 191.03904 76.41562 45.84937mean(log(rot_post_co$rot_latency))[1] 4.765482mean(log(rot_post_co[rot_post_co$genotype == 'sca', c('rot_latency')]))[1] 4.621306sd(log(rot_post_co$rot_latency))*c(5)[1] 2.071834.2 Setting priors, preTreat RotLat: all days
## priors for Gaussian model
prior01 <- c(set_prior("student_t(3, 127, 229)", class = "b", coef = "Intercept"),
set_prior("normal(0,92)", class = "b", coef = "groupctrl.wt"),
set_prior("normal(0,92)", class = "b", coef = "groupeda.wt"),
set_prior("normal(0,46)" , class = "b", coef = "groupeda.sca"),
set_prior("normal(10,20)", class = "b", coef = "day_cen"),
set_prior("normal(0,20)", class = "b", coef = "session_cen"),
set_prior("normal(0,2)", class = "b", coef = "session_cen:day_cen"),
set_prior("normal(0.5,2)", class = "b", coef = "mean_lat.x"))
## prior for Gamma (log-linked) model
prior02 <- c(set_prior("student_t(3, 4.8, 2.1)", class = "b", coef = "Intercept"),
set_prior("normal(0,2)", class = "b", coef = "groupctrl.wt"),
set_prior("normal(0,2)", class = "b", coef = "groupeda.wt"),
set_prior("normal(0,1.2)", class = "b", coef = "groupeda.sca"),
set_prior("normal(0.2,0.5)", class = "b", coef = "day_cen"),
set_prior("normal(0,0.5)", class = "b", coef = "session_cen"),
set_prior("normal(0,0.2)", class = "b", coef = "session_cen:day_cen"),
set_prior("normal(0.5,1)", class = "b", coef = "mean_lat.x"))4.3 Fitting and diagnsotics, preTreat RotLat: all days
## fitting gaussian model
model1 <- run_model(
brm(rot_latency ~ 0 + Intercept + mean_lat.x + session_cen +
day_cen + session_cen:day_cen + group + (1|id),
family=gaussian(),
data = rot_post_co,
prior= prior01,
save_pars = save_pars(all = TRUE),
iter=8000, warmup=2000,chains=2,cores=1,seed=17,
control = list(adapt_delta = 0.99, max_treedepth = 15)),
'output/brm/model1', reuse = TRUE)
### fitting gamma model
model2 <- run_model(
brm(rot_latency ~ 0 + Intercept + mean_lat.x + session_cen +
day_cen + session_cen:day_cen + group + (1|id),
family=Gamma(link = 'log'),
data = rot_post_co,
prior= prior02,
save_pars = save_pars(all = TRUE),
iter=8000, warmup=2000,chains=2,cores=1,seed=17,
control = list(adapt_delta = 0.99, max_treedepth = 15)),
'output/brm/model2', reuse = TRUE)
## models summary
summary(model1) Family: gaussian
Links: mu = identity; sigma = identity
Formula: rot_latency ~ 0 + Intercept + mean_lat.x + session_cen + day_cen + session_cen:day_cen + group + (1 | id)
Data: rot_post_co (Number of observations: 960)
Draws: 2 chains, each with iter = 8000; warmup = 2000; thin = 1;
total post-warmup draws = 12000
Group-Level Effects:
~id (Number of levels: 80)
Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
sd(Intercept) 24.22 2.29 20.09 29.10 1.00 2140 4087
Population-Level Effects:
Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
Intercept 39.03 14.23 10.97 66.87 1.00 1272 2805
mean_lat.x 0.56 0.11 0.35 0.77 1.00 1358 2896
session_cen 5.90 0.89 4.15 7.63 1.00 12064 9673
day_cen 8.99 1.23 6.56 11.37 1.00 10391 8646
groupctrl.wt 29.56 8.21 13.30 46.05 1.00 1211 2165
groupeda.wt 27.82 8.16 11.59 43.94 1.00 1536 3203
groupeda.sca 3.29 7.95 -12.34 18.69 1.00 1313 2785
session_cen:day_cen -0.26 0.96 -2.13 1.63 1.00 11634 9794
Family Specific Parameters:
Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
sigma 30.88 0.73 29.48 32.33 1.00 10860 9755
Draws were sampled using sampling(NUTS). For each parameter, Bulk_ESS
and Tail_ESS are effective sample size measures, and Rhat is the potential
scale reduction factor on split chains (at convergence, Rhat = 1).summary(model2) Family: gamma
Links: mu = log; shape = identity
Formula: rot_latency ~ 0 + Intercept + mean_lat.x + session_cen + day_cen + session_cen:day_cen + group + (1 | id)
Data: rot_post_co (Number of observations: 960)
Draws: 2 chains, each with iter = 8000; warmup = 2000; thin = 1;
total post-warmup draws = 12000
Group-Level Effects:
~id (Number of levels: 80)
Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
sd(Intercept) 0.20 0.02 0.17 0.24 1.00 3088 5140
Population-Level Effects:
Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
Intercept 4.11 0.12 3.88 4.35 1.00 3155 4839
mean_lat.x 0.00 0.00 0.00 0.01 1.00 3220 4962
session_cen 0.05 0.01 0.03 0.06 1.00 17402 8648
day_cen 0.07 0.01 0.05 0.09 1.00 16963 8730
groupctrl.wt 0.26 0.07 0.12 0.40 1.00 2524 4237
groupeda.wt 0.25 0.07 0.11 0.38 1.00 2755 4614
groupeda.sca 0.05 0.07 -0.08 0.19 1.00 2617 4652
session_cen:day_cen -0.01 0.01 -0.02 0.01 1.00 16026 8456
Family Specific Parameters:
Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
shape 12.82 0.61 11.66 14.04 1.00 14064 8030
Draws were sampled using sampling(NUTS). For each parameter, Bulk_ESS
and Tail_ESS are effective sample size measures, and Rhat is the potential
scale reduction factor on split chains (at convergence, Rhat = 1).Convergence is OK; both models lead to the same conclusions and similar estimates
4.4 PPC, preTreat RotLat: all days
pla<-pp_check(model1,type='dens_overlay',ndraws = 50)
plb<-pp_check(model2,type='dens_overlay',ndraws = 50)
plot_grid(pla, plb, labels=c("A", "B"), ncol = 1, nrow = 2)
pla<-pp_check(model1,type='dens_overlay_grouped',ndraws = 50,group='group')
plb<-pp_check(model2,type='dens_overlay_grouped',ndraws = 50,group='group')
plot_grid(pla, plb, labels=c("A", "B"), ncol = 1, nrow = 2)
pla<-pp_check(model1,type='scatter_avg') Using all posterior draws for ppc type 'scatter_avg' by default.plb<-pp_check(model2,type='scatter_avg') Using all posterior draws for ppc type 'scatter_avg' by default.plot_grid(pla, plb, labels=c("A", "B"), ncol = 1, nrow = 2)
pla<-pp_check(model1,type='ecdf_overlay_grouped',ndraws = 50,group='group')
plb<-pp_check(model2,type='ecdf_overlay_grouped',ndraws = 50,group='group')
plot_grid(pla, plb, labels=c("A", "B"), ncol = 1, nrow = 2)
pla<-pp_check(model1,type="stat_2d", stat = c("max", "min"),ndraws=200)
plb<-pp_check(model2,type="stat_2d", stat = c("max", "min"),ndraws=200)
plot_grid(pla, plb, labels=c("A", "B"), ncol = 1, nrow = 2)
pla<-pp_check(model1,type="stat_2d", stat = c("mean", "sd"),ndraws=200)
plb<-pp_check(model2,type="stat_2d", stat = c("mean", "sd"),ndraws=200)
plot_grid(pla, plb, labels=c("A", "B"), ncol = 1, nrow = 2)
### Are there too influential observations?
par(mfrow=c(1,2))
plot(loo(model1));loo(model1)
Computed from 12000 by 960 log-likelihood matrix
Estimate SE
elpd_loo -4694.3 26.9
p_loo 73.1 4.1
looic 9388.6 53.8
------
Monte Carlo SE of elpd_loo is 0.1.
Pareto k diagnostic values:
Count Pct. Min. n_eff
(-Inf, 0.5] (good) 959 99.9% 2736
(0.5, 0.7] (ok) 1 0.1% 561
(0.7, 1] (bad) 0 0.0% <NA>
(1, Inf) (very bad) 0 0.0% <NA>
All Pareto k estimates are ok (k < 0.7).
See help('pareto-k-diagnostic') for details.plot(loo(model2));loo(model2)
Computed from 12000 by 960 log-likelihood matrix
Estimate SE
elpd_loo -4759.1 28.9
p_loo 62.9 3.4
looic 9518.1 57.7
------
Monte Carlo SE of elpd_loo is 0.1.
All Pareto k estimates are good (k < 0.5).
See help('pareto-k-diagnostic') for details.Posterior predictive checks (PPC) shows that Gaussian model fits data clearly better than Gamma model. All models are OK, no pareto value > 0.7, indicating no overly influential observation
Lets to fit two other Gaussian models with different deficnition of pre-treatment rotarod latency
4.5 Setting prior for other models
## preTreat RotLat: 3rd day only
prior03 <- c(set_prior("student_t(3, 127, 229)", class = "b", coef = "Intercept"),
set_prior("normal(0,92)", class = "b", coef = "groupctrl.wt"),
set_prior("normal(0,92)", class = "b", coef = "groupeda.wt"),
set_prior("normal(0,46)" , class = "b", coef = "groupeda.sca"),
set_prior("normal(10,20)", class = "b", coef = "day_cen"),
set_prior("normal(0,20)", class = "b", coef = "session_cen"),
set_prior("normal(0,2)", class = "b", coef = "session_cen:day_cen"),
set_prior("normal(0.5,2)", class = "b", coef = "mean_lat"))
## preTreat RotLat: 2nd and 3rd day
prior04 <- c(set_prior("student_t(3, 127, 229)", class = "b", coef = "Intercept"),
set_prior("normal(0,92)", class = "b", coef = "groupctrl.wt"),
set_prior("normal(0,92)", class = "b", coef = "groupeda.wt"),
set_prior("normal(0,46)" , class = "b", coef = "groupeda.sca"),
set_prior("normal(10,20)", class = "b", coef = "day_cen"),
set_prior("normal(0,20)", class = "b", coef = "session_cen"),
set_prior("normal(0,2)", class = "b", coef = "session_cen:day_cen"),
set_prior("normal(0.5,2)", class = "b", coef = "mean_lat.y"))4.6 Fitting and diagnsotics, preTreat RotLat: 3rd (+2nd) days
## preTreat RotLat: 3rd day only
model3 <- run_model(
brm(rot_latency ~ 0 + Intercept + mean_lat + session_cen +
day_cen + session_cen:day_cen + group + (1|id),
family=gaussian(),
data = rot_post_co,
prior= prior03,
save_pars = save_pars(all = TRUE),
iter=8000, warmup=2000,chains=2,cores=1,seed=17,
control = list(adapt_delta = 0.99)),
'output/brm/model3', reuse = TRUE)
## preTreat RotLat: 2nd + 3rd day
model4 <- run_model(
brm(rot_latency ~ 0 + Intercept + mean_lat.y + session_cen +
day_cen + session_cen:day_cen + group + (1|id),
family=gaussian(),
data = rot_post_co,
prior= prior04,
save_pars = save_pars(all = TRUE),
iter=8000, warmup=2000,chains=2,cores=1,seed=17,
control = list(adapt_delta = 0.99)),
'output/brm/model4', reuse = TRUE)
## models summary
summary(model3) Family: gaussian
Links: mu = identity; sigma = identity
Formula: rot_latency ~ 0 + Intercept + mean_lat + session_cen + day_cen + session_cen:day_cen + group + (1 | id)
Data: rot_post_co (Number of observations: 960)
Draws: 2 chains, each with iter = 8000; warmup = 2000; thin = 1;
total post-warmup draws = 12000
Group-Level Effects:
~id (Number of levels: 80)
Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
sd(Intercept) 25.41 2.35 21.25 30.35 1.00 2368 4856
Population-Level Effects:
Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
Intercept 57.53 12.68 32.65 82.28 1.00 1533 2830
mean_lat 0.38 0.08 0.21 0.54 1.00 1635 3175
session_cen 5.89 0.90 4.12 7.66 1.00 12401 8102
day_cen 8.97 1.22 6.56 11.37 1.00 14924 9215
groupctrl.wt 26.77 8.58 10.17 43.81 1.00 1583 3191
groupeda.wt 27.79 8.89 10.47 45.04 1.00 1713 3408
groupeda.sca 2.87 8.39 -13.56 19.45 1.00 1604 3085
session_cen:day_cen -0.27 0.94 -2.14 1.56 1.00 12862 9054
Family Specific Parameters:
Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
sigma 30.89 0.73 29.50 32.35 1.00 11835 9004
Draws were sampled using sampling(NUTS). For each parameter, Bulk_ESS
and Tail_ESS are effective sample size measures, and Rhat is the potential
scale reduction factor on split chains (at convergence, Rhat = 1).summary(model4) Family: gaussian
Links: mu = identity; sigma = identity
Formula: rot_latency ~ 0 + Intercept + mean_lat.y + session_cen + day_cen + session_cen:day_cen + group + (1 | id)
Data: rot_post_co (Number of observations: 960)
Draws: 2 chains, each with iter = 8000; warmup = 2000; thin = 1;
total post-warmup draws = 12000
Group-Level Effects:
~id (Number of levels: 80)
Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
sd(Intercept) 24.55 2.31 20.48 29.48 1.00 2335 4617
Population-Level Effects:
Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
Intercept 47.30 13.27 21.52 73.68 1.00 1657 2816
mean_lat.y 0.47 0.09 0.29 0.65 1.00 1464 3143
session_cen 5.89 0.90 4.13 7.64 1.00 13934 8304
day_cen 8.96 1.23 6.53 11.38 1.00 12177 8586
groupctrl.wt 25.76 8.42 9.39 42.52 1.00 1587 3152
groupeda.wt 25.71 8.75 8.25 42.69 1.00 1504 2833
groupeda.sca 1.51 8.23 -15.36 17.35 1.00 1433 2399
session_cen:day_cen -0.26 0.97 -2.18 1.64 1.00 13103 8714
Family Specific Parameters:
Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
sigma 30.91 0.73 29.52 32.39 1.00 10656 8660
Draws were sampled using sampling(NUTS). For each parameter, Bulk_ESS
and Tail_ESS are effective sample size measures, and Rhat is the potential
scale reduction factor on split chains (at convergence, Rhat = 1).All models converged well (Rhat values = 1), ESSs are sufficient
4.7 PPC, preTreat RotLat: 3rd (+2nd) days
pla<-pp_check(model3,type='dens_overlay',ndraws = 50)
plb<-pp_check(model4,type='dens_overlay',ndraws = 50)
plot_grid(pla, plb, labels=c("A", "B"), ncol = 1, nrow = 2)
pla<-pp_check(model3,type='dens_overlay_grouped',ndraws = 50,group='group')
plb<-pp_check(model4,type='dens_overlay_grouped',ndraws = 50,group='group')
plot_grid(pla, plb, labels=c("A", "B"), ncol = 1, nrow = 2)
pla<-pp_check(model3,type='scatter_avg') Using all posterior draws for ppc type 'scatter_avg' by default.plb<-pp_check(model4,type='scatter_avg') Using all posterior draws for ppc type 'scatter_avg' by default.plot_grid(pla, plb, labels=c("A", "B"), ncol = 1, nrow = 2)
pla<-pp_check(model3,type="stat_2d", stat = c("max", "min"),ndraws=200)
plb<-pp_check(model4,type="stat_2d", stat = c("max", "min"),ndraws=200)
plot_grid(pla, plb, labels=c("A", "B"), ncol = 1, nrow = 2)
pla<-pp_check(model3,type="stat_2d", stat = c("mean", "sd"),ndraws=200)
plb<-pp_check(model4,type="stat_2d", stat = c("mean", "sd"),ndraws=200)
plot_grid(pla, plb, labels=c("A", "B"), ncol = 1, nrow = 2)
### Are there too influential observations?
par(mfrow=c(1,2))
plot(loo(model3));loo(model3)
Computed from 12000 by 960 log-likelihood matrix
Estimate SE
elpd_loo -4694.8 26.8
p_loo 74.0 4.1
looic 9389.6 53.6
------
Monte Carlo SE of elpd_loo is 0.1.
All Pareto k estimates are good (k < 0.5).
See help('pareto-k-diagnostic') for details.plot(loo(model4));loo(model4)
Computed from 12000 by 960 log-likelihood matrix
Estimate SE
elpd_loo -4694.7 26.8
p_loo 73.5 4.1
looic 9389.5 53.6
------
Monte Carlo SE of elpd_loo is 0.1.
All Pareto k estimates are good (k < 0.5).
See help('pareto-k-diagnostic') for details.Everything seems fine.
4.8 Models comparison
Next, we will compare all three Gaussian models using LOO-CV
model1 <- add_criterion(model1, criterion = "loo")
model3 <- add_criterion(model3, criterion = "loo")
model4 <- add_criterion(model4, criterion = "loo")
loo_compare(model1, model3, model4) elpd_diff se_diff
model1 0.0 0.0
model4 -0.5 0.6
model3 -0.5 1.1 Using baseline rotarod latency, averaged from all pre-treatment rotarod sessions, is comparable to using covariate from specifically 3rd or 2nd + 3rd day.
The model we currently have is, however, relatively simple. Likely the most importantly, it does not allow modelling the situation when group differ in day-to-day gain (implying motor learning).
Thus, let to create another model, which looks like model1 but also include interaction between group and effect of time
4.9 Model including group : day interaction
## fitting model
model5 <- run_model(
brm(rot_latency ~ 0 + Intercept + mean_lat.x + session_cen +
day_cen + session_cen:day_cen + group + day_cen:group + (1|id),
family=gaussian(),
data = rot_post_co,
prior= prior01,
save_pars = save_pars(all = TRUE),
iter=8000, warmup=2000,chains=2,cores=1,seed=17,
control = list(adapt_delta = 0.99, max_treedepth = 15)),
'output/brm/model5', reuse = TRUE)
## model summary
summary(model5) Family: gaussian
Links: mu = identity; sigma = identity
Formula: rot_latency ~ 0 + Intercept + mean_lat.x + session_cen + day_cen + session_cen:day_cen + group + day_cen:group + (1 | id)
Data: rot_post_co (Number of observations: 960)
Draws: 2 chains, each with iter = 8000; warmup = 2000; thin = 1;
total post-warmup draws = 12000
Group-Level Effects:
~id (Number of levels: 80)
Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
sd(Intercept) 24.24 2.30 20.01 29.09 1.00 2994 4963
Population-Level Effects:
Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS
Intercept 39.08 13.91 11.61 66.22 1.00 1919
mean_lat.x 0.56 0.11 0.36 0.77 1.00 1954
session_cen 5.90 0.90 4.16 7.63 1.00 13428
day_cen 3.66 2.40 -1.02 8.35 1.00 4552
groupctrl.wt 29.77 8.20 13.69 45.72 1.00 1876
groupeda.wt 27.86 8.40 11.12 44.26 1.00 1603
groupeda.sca 3.42 8.10 -12.55 19.03 1.00 1815
session_cen:day_cen -0.27 0.95 -2.17 1.59 1.00 15380
day_cen:groupctrl.wt 13.32 3.40 6.60 19.90 1.00 5779
day_cen:groupeda.wt 4.96 3.42 -1.73 11.68 1.00 5807
day_cen:groupeda.sca 3.01 3.38 -3.59 9.48 1.00 6298
Tail_ESS
Intercept 3810
mean_lat.x 3990
session_cen 9003
day_cen 6931
groupctrl.wt 3583
groupeda.wt 3157
groupeda.sca 3156
session_cen:day_cen 8538
day_cen:groupctrl.wt 7919
day_cen:groupeda.wt 7749
day_cen:groupeda.sca 8316
Family Specific Parameters:
Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
sigma 30.65 0.74 29.26 32.15 1.00 11851 8434
Draws were sampled using sampling(NUTS). For each parameter, Bulk_ESS
and Tail_ESS are effective sample size measures, and Rhat is the potential
scale reduction factor on split chains (at convergence, Rhat = 1).## prior summary
prior_summary(model5) prior class coef group resp dpar nlpar lb ub
(flat) b
normal(10,20) b day_cen
(flat) b day_cen:groupctrl.wt
(flat) b day_cen:groupeda.sca
(flat) b day_cen:groupeda.wt
normal(0,92) b groupctrl.wt
normal(0,46) b groupeda.sca
normal(0,92) b groupeda.wt
student_t(3, 127, 229) b Intercept
normal(0.5,2) b mean_lat.x
normal(0,20) b session_cen
normal(0,2) b session_cen:day_cen
student_t(3, 0, 45.3) sd 0
student_t(3, 0, 45.3) sd id 0
student_t(3, 0, 45.3) sd Intercept id 0
student_t(3, 0, 45.3) sigma 0
source
default
user
(vectorized)
(vectorized)
(vectorized)
user
user
user
user
user
user
user
default
(vectorized)
(vectorized)
default## PPC
pp_check(model5,type='dens_overlay',ndraws = 50)
pp_check(model5,type='dens_overlay_grouped',ndraws = 50,group='group') 
pp_check(model5,type='scatter_avg') Using all posterior draws for ppc type 'scatter_avg' by default.
pp_check(model5,type="stat_2d", stat = c("max", "min"),ndraws=200)
pp_check(model5,type="stat_2d", stat = c("mean", "sd"),ndraws=200)
## Are there too influential observations?
plot(loo(model5));loo(model5)
Computed from 12000 by 960 log-likelihood matrix
Estimate SE
elpd_loo -4688.7 26.3
p_loo 76.3 4.1
looic 9377.3 52.6
------
Monte Carlo SE of elpd_loo is 0.1.
All Pareto k estimates are good (k < 0.5).
See help('pareto-k-diagnostic') for details.## comparing with original model via LOO-CV
model5 <- add_criterion(model5, criterion = "loo")
loo_compare(model1, model5) elpd_diff se_diff
model5 0.0 0.0
model1 -5.6 4.8 Models comparison indicate that more complex model, allowing group-specific change in rotarod performance over day, has better predictive accuracy. The last complex model (Gaussian distribution, ignoring sex, including day : group interaction) thus will be used as the final model.
The model coverged well, ESSs are sufficient, PPC seems fine.
4.10 Posterior draws extraction
post_fix<-as.data.frame(model5, variable = c("b_Intercept","b_groupctrl.wt",
"b_groupeda.wt","b_groupeda.sca"))
names(post_fix)[1]<-"sca_ctrl"
post_fix$wt_ctrl<-post_fix$sca_ctrl+post_fix$b_groupctrl.wt
post_fix$wt_eda<-post_fix$sca_ctrl+post_fix$b_groupeda.wt
post_fix$sca_eda<-post_fix$sca_ctrl+post_fix$b_groupeda.sca
post_fix$wt_contrast<-post_fix$wt_eda-post_fix$wt_ctrl
summary(post_fix) sca_ctrl b_groupctrl.wt b_groupeda.wt b_groupeda.sca
Min. :-16.00 Min. : 0.5698 Min. :-5.929 Min. :-28.281
1st Qu.: 29.81 1st Qu.:24.2420 1st Qu.:22.191 1st Qu.: -2.051
Median : 39.13 Median :29.8007 Median :27.892 Median : 3.351
Mean : 39.08 Mean :29.7691 Mean :27.858 Mean : 3.416
3rd Qu.: 48.53 3rd Qu.:35.2160 3rd Qu.:33.536 3rd Qu.: 9.024
Max. : 91.93 Max. :66.4481 Max. :62.362 Max. : 33.938
wt_ctrl wt_eda sca_eda wt_contrast
Min. : 9.15 Min. : 2.729 Min. :-14.06 Min. :-40.123
1st Qu.: 58.85 1st Qu.: 56.104 1st Qu.: 32.93 1st Qu.: -7.248
Median : 69.06 Median : 67.066 Median : 42.76 Median : -1.859
Mean : 68.85 Mean : 66.939 Mean : 42.50 Mean : -1.911
3rd Qu.: 78.94 3rd Qu.: 77.698 3rd Qu.: 52.27 3rd Qu.: 3.473
Max. :126.46 Max. :124.232 Max. :103.39 Max. : 30.281 groupquant<-sapply(post_fix[,c(5,6,1,7)],function(p) quantile(p, probs = c(0.025,0.975,0.5)))
groupquant <- groupquant5 Data visualisation
rot_post_co$group <- factor(rot_post_co$group,levels = c('ctrl.wt','eda.wt','ctrl.sca','eda.sca') )
m= matrix(c(1,2), nrow = 1, ncol=2, byrow = TRUE)
layout(mat = m, heights = c(1),
widths = c(0.65, 0.4))
par(mar=c(2,3,0,0))
par(mgp=c(2,0.6,0))
range<-c(0, 180);scal<-range[2]-range[1]
xrange<-c(0.5,4.5);xscal=xrange[2]-xrange[1]
plot(NULL,xlim=c(xrange[1],xrange[2]),ylim=c(range[1],range[2]),xlab="",
ylab="Rotarod latency (s)",las=1, axes=FALSE,cex.lab=1.2)
rect(xrange[1],range[2],xrange[2],range[1],col="grey92",border=NA)
x<-range[1]
repeat{
lines(c(xrange[1],xrange[2]),c(x,x),col="white",lwd=0.7)
x=x+30;if(x>range[2]){break}}
x<-0.5
repeat{
lines(c(x,x),c(range[1],range[2]),col="white",lwd=0.7)
x=x+0.5;if(x>4.5){break}}
subdat <- subset(rotarod, pre_post == 'pre' & day == 1)
preD1 <- tapply(subdat$rot_latency, list(subdat$group, subdat$session), mean)
subdat <- subset(rotarod, pre_post == 'pre' & day == 2)
preD2 <- tapply(subdat$rot_latency, list(subdat$group, subdat$session), mean)
subdat <- subset(rotarod, pre_post == 'pre' & day == 3)
preD3 <- tapply(subdat$rot_latency, list(subdat$group, subdat$session), mean)
subdat <- subset(rotarod, pre_post == 'post' & day == 1)
postD1 <- tapply(subdat$rot_latency, list(subdat$group, subdat$session), mean)
subdat <- subset(rotarod, pre_post == 'post' & day == 2)
postD2 <- tapply(subdat$rot_latency, list(subdat$group, subdat$session), mean)
subdat <- subset(rotarod, pre_post == 'post' & day == 3)
postD3 <- tapply(subdat$rot_latency, list(subdat$group, subdat$session), mean)
xl <- c(0.85, 0.95, 1.05, 1.15)
for (x in 1:4){
lines(preD1[x,]~c(xl),col=cola[x], lwd = 1.6, type = 'o', pch= 16)
lines(preD2[x,]~c(xl+0.5),col= cola[x], lwd = 1.6, type = 'o', pch= 16)
lines(preD3[x,]~c(xl+1),col= cola[x], lwd = 1.6,, type = 'o', pch= 16)
lines(postD1[x,]~c(xl+2),col= cola[x], lwd = 1.6, type = 'o', pch= 16)
lines(postD2[x,]~c(xl+2.5),col= cola[x], lwd = 1.6, type = 'o', pch= 16)
lines(postD3[x,]~c(xl+3),col= cola[x], lwd = 1.6, type = 'o', pch= 16)
}
tckk=-0.02
axis(2,las=2,cex.axis=1.1,at=seq(range[1],range[2],by=30),pos=xrange[1],tck=tckk)
axis(side=1,las=1,cex.axis=1.1,at=c(1,1.5,2,3,3.5,4),
labels=c(rep("",6)),pos=range[1],tck=tckk)
lines(c(0.5,4.5), c(0,0))
text(c(1, 1.5, 2, 3, 3.5, 4) ,c(rep(range[1]-0.05*scal)),
c('D1', 'D2', 'D3', 'D1', 'D2', 'D3'), xpd= TRUE)
text(c(1.5, 3.5),c(range[1]-0.117*scal, range[1]-0.117*scal),
c('PRE', 'POST'), font = 2, xpd = TRUE)
lines(c(1,2),c(range[1]-0.085*scal, range[1]-0.085*scal), xpd = TRUE)
lines(c(3,4),c(range[1]-0.085*scal, range[1]-0.085*scal), xpd = TRUE)
x = 1
y = 58;yde <- 10
lines(c(1.2,1.6),c(y,y), col = cola[x], pch = 16, lwd = 1.6)
text(2.05, y, 'WT 0', col = cola[x])
points(1.4, y, pch = 16, col = cola[x]); x = x+1; y = y-yde
lines(c(1.2,1.6),c(y,y), col = cola[x], pch = 16, lwd = 1.6)
text(2.05, y, 'WT edv', col = cola[x])
points(1.4, y, pch = 16, col = cola[x]); x = x+1; y = y-yde
lines(c(1.2,1.6),c(y,y), col = cola[x], pch = 16, lwd = 1.6)
text(2.05, y, 'SCA1 0', col = cola[x])
points(1.4, y, pch = 16, col = cola[x]); x = x+1; y = y-yde
lines(c(1.2,1.6),c(y,y), col = cola[x], pch = 16, lwd = 1.6)
text(2.05, y, 'SCA1 edv', col = cola[x])
points(1.4, y, pch = 16, col = cola[x]); x = x+1; y = y-yde
## probability distributions
par(mar=c(2,0,0,0))
dif<-data.frame(post_fix$wt_contrast,post_fix$b_groupeda.sca,
post_fix$b_groupctrl.wt)
dif<-(dif)
yla<-"Difference in rotarod latency (s)"
tckk=-0.018
ste<-seq(-40,60,by=20)
mons_poste(dif, 0)
zpos=seq(0,1,1/3)
xx=1;ind=0.4;xpol<- -32
text(xpol,zpos[xx]+zpos[2]*ind,
"Wt: ed x 0 ",cex=1)
xx=xx+1
text(xpol,zpos[xx]+zpos[2]*ind,
"SCA1: ed x 0 ",cex=1)
xx=xx+1
text(xpol,zpos[xx]+zpos[2]*ind,
"0: Wt x SCA1",cex=1)