rm(list = ls())
suppressWarnings(suppressMessages( {
library(brms)
library(beeswarm)
library(vioplot)
library(glmmTMB)
library(car)
library(cowplot)
library(ggplot2)
library(tidyverse)
library(bayesplot)
} ) )Experimental treatment with edaravone in a mouse model of spinocerebellar ataxia 1
Analysis of grip strength
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 import and wrangling
## import of data
load("source_data/myEnvironment.RData")
## averaging grip strength (4 measuremnts)
behav$grip_strength<-rowSums(behav[,c(47:50)])/43 Data exploration
par(mar=c(5,3.5,1,1))
par(mgp=c(1.9,0.7,0))
par(mfrow=c(3,1))
plot(behav$grip_strength~behav$group,las=2,xlab="",
col=cola)
plot(behav$grip_strength~behav$group_sex,las=2,xlab="",
col=cola[c(1,1,2,2,3,3,4,4)])
hist(behav$grip_strength, 20)
Males may have generally larger strength. Spread seems to be larger for larger values.
4 Modelling
We will fit two gamma models, one with sex as a covariate
4.1 Prior specification
behav$group<-relevel(behav$group,ref='ctrl.sca')
log(mean(behav$grip_strength))[1] -0.1267686sd(log(behav$grip_strength))*5[1] 1.582227Setting priors
prior01 <- c(set_prior("student_t(3,-0.13,1.6)", 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"))
prior02 <- c(set_prior("student_t(3,-0.13,1.6)", 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,1.2)" , class = "b", coef = "sex_num"))4.2 Fitting models and diagnostics
## fitting model without sex
grip_strength_model_01 <- run_model(
brm(grip_strength ~0+Intercept+group,
data=behav, prior = prior01,
family=Gamma(link="log"),
save_pars = save_pars(all = TRUE),
iter=8000, warmup=2000,chains=2,cores=1,seed=17,
control = list(adapt_delta = 0.99)),
'output/brm/grip_strength_model_01', reuse = TRUE)
## fitting model with 'sex' as a covariate
grip_strength_model_02 <- run_model(
brm(grip_strength~0+Intercept+group+sex_num,
data=behav, prior = prior02,
family=Gamma(link="log"),
save_pars = save_pars(all = TRUE),
iter=8000, warmup=2000,chains=2,cores=1,seed=17,
control = list(adapt_delta = 0.99)),
'output/brm/grip_strength_model_02', reuse = TRUE)
mcmc_trace(grip_strength_model_01)
mcmc_trace(grip_strength_model_02)
summary(grip_strength_model_01) Family: gamma
Links: mu = log; shape = identity
Formula: grip_strength ~ 0 + Intercept + group
Data: behav (Number of observations: 80)
Draws: 2 chains, each with iter = 8000; warmup = 2000; thin = 1;
total post-warmup draws = 12000
Population-Level Effects:
Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
Intercept -0.25 0.07 -0.38 -0.11 1.00 3857 5146
groupctrl.wt 0.20 0.10 0.01 0.39 1.00 4779 6147
groupeda.wt 0.27 0.10 0.08 0.46 1.00 4374 5682
groupeda.sca -0.01 0.10 -0.20 0.18 1.00 4623 5893
Family Specific Parameters:
Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
shape 11.40 1.81 8.10 15.17 1.00 6682 6444
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(grip_strength_model_02) Family: gamma
Links: mu = log; shape = identity
Formula: grip_strength ~ 0 + Intercept + group + sex_num
Data: behav (Number of observations: 80)
Draws: 2 chains, each with iter = 8000; warmup = 2000; thin = 1;
total post-warmup draws = 12000
Population-Level Effects:
Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
Intercept -0.25 0.07 -0.38 -0.12 1.00 4298 5795
groupctrl.wt 0.20 0.09 0.01 0.38 1.00 5134 6848
groupeda.wt 0.27 0.09 0.09 0.45 1.00 5119 6744
groupeda.sca -0.01 0.09 -0.19 0.18 1.00 5217 6836
sex_num 0.14 0.07 0.02 0.27 1.00 8843 8056
Family Specific Parameters:
Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
shape 11.94 1.93 8.43 15.98 1.00 7661 6844
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 and ESS of both models are OK
4.3 PPC
pla<-pp_check(grip_strength_model_01,type='dens_overlay',ndraws = 50)
plb<-pp_check(grip_strength_model_02,type='dens_overlay',ndraws = 50)
plot_grid(pla, plb, labels=c("A", "B"), ncol = 1, nrow = 2)
pla<-pp_check(grip_strength_model_01,type='dens_overlay_grouped',ndraws = 50,group='group')
plb<-pp_check(grip_strength_model_02,type='dens_overlay_grouped',ndraws = 50,group='group')
plot_grid(pla, plb, labels=c("A", "B"), ncol = 1, nrow = 2)
pla<-pp_check(grip_strength_model_01,type='scatter_avg') Using all posterior draws for ppc type 'scatter_avg' by default.plb<-pp_check(grip_strength_model_02,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)
pp_check(grip_strength_model_01,type='ecdf_overlay',ndraws = 200) 
pp_check(grip_strength_model_02,type='ecdf_overlay',ndraws = 200) 
pla<-pp_check(grip_strength_model_01,type='ecdf_overlay_grouped',ndraws = 50,group='group')
plb<-pp_check(grip_strength_model_02,type='ecdf_overlay_grouped',ndraws = 50,group='group')
plot_grid(pla, plb, labels=c("A", "B"), ncol = 1, nrow = 2)
pla<-pp_check(grip_strength_model_01,type="stat_2d", stat = c("max", "min"),ndraws=200)
plb<-pp_check(grip_strength_model_02,type="stat_2d", stat = c("max", "min"),ndraws=200)
plot_grid(pla, plb, labels=c("A", "B"), ncol = 1, nrow = 2)
pla<-pp_check(grip_strength_model_01,type="stat_2d", stat = c("mean", "sd"),ndraws=200)
plb<-pp_check(grip_strength_model_02,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(grip_strength_model_01));loo(grip_strength_model_01)
Computed from 12000 by 80 log-likelihood matrix
Estimate SE
elpd_loo -6.3 6.4
p_loo 4.8 0.8
looic 12.5 12.9
------
Monte Carlo SE of elpd_loo is 0.0.
All Pareto k estimates are good (k < 0.5).
See help('pareto-k-diagnostic') for details.plot(loo(grip_strength_model_02));loo(grip_strength_model_02)
Computed from 12000 by 80 log-likelihood matrix
Estimate SE
elpd_loo -4.7 6.0
p_loo 5.5 0.8
looic 9.4 11.9
------
Monte Carlo SE of elpd_loo is 0.0.
All Pareto k estimates are good (k < 0.5).
See help('pareto-k-diagnostic') for details.Both models seem to retrodict the data well
4.4 Leave-one out cross-validation to compare predicitve accuracy of both models
grip_strength_model_01 <- add_criterion(grip_strength_model_01, criterion = "loo")
grip_strength_model_02 <- add_criterion(grip_strength_model_02, criterion = "loo")
loo_compare(grip_strength_model_01,grip_strength_model_02) elpd_diff se_diff
grip_strength_model_02 0.0 0.0
grip_strength_model_01 -1.6 2.2 The difference between the models in terms predictive accuracy is minimal and both show basically the same estimates. We will use the 1st (simpler) model
4.5 Extraction of posterior samples
post_fix<-as.data.frame(grip_strength_model_01, 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
groupquant<-sapply(post_fix[,c(5,6,1,7)],function(p) quantile(p, probs = c(0.025,0.975,0.5)))
groupquant<-exp(groupquant)4.6 Data and model visualization
behav$group<-factor(behav$group,levels=c('ctrl.wt','eda.wt','ctrl.sca','eda.sca'))
par(mfrow=c(1,2))
par(mar=c(2,3,0,0))
par(mgp=c(2,0.6,0))
range<-c(0,2);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="Grip strength",las=1, axes=FALSE,cex.lab=1)
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+0.2;if(x>range[2]){break}}
x<-1
repeat{
lines(c(x,x),c(range[1],range[2]),col="white",lwd=0.7)
x=x+1;if(x>4.5){break}}
vioplot(behav$grip_strength~behav$group,col=colb,at=c(1:4),add=T,
border=cola,axes=FALSE,drawRect=F,lwd=1,wex=0.8)
beeswarm(behav$grip_strength~behav$group,col=colc,at=c(1:4),add=T,pwpch=(16.5+behav$sex_num) )
tp<-(groupquant[3,])
xx<-1;wid=0.25
repeat{
lines(c(xx-wid,xx+wid),c(tp[xx],tp[xx]),lwd=2,col="black");
lines(c(xx,xx),c(groupquant[1,xx],groupquant[2,xx]),lwd=1.7,col="black")
xx<-xx+1
if(xx>4){break}}
tckk=-0.02
axis(2,las=2,cex.axis=1.1,at=seq(range[1],range[2],by=0.4),
labels=c(rep("",length(seq(range[1],range[2],by=0.4)))),pos=xrange[1],tck=tckk)
axis(2,las=2,cex.axis=1.1,at=seq(range[1],range[2],by=0.4),pos=xrange[1],tck=tckk)
axis(side=1,las=1,cex.axis=1.1,at=c(seq(1,4,by=1)),
labels=c(rep("",length(seq(1,4,by=1)))),pos=range[1],tck=tckk)
lines(c(xrange[1],xrange[2]),c(range[1],range[1]))
text(c(1:4),c(rep(range[1]-0.046*scal,4)),xpd=T,cex=1.1,col=cola,
c("Wt","Wt", "SCA1", "SCA1"))
ypo<-c(range[1]-0.1*scal,range[1]-0.1*scal)
text(c(1:4),c(ypo[1],ypo[2],ypo[1],ypo[2]),xpd=T,cex=1.1,col=cola,
c("0","edv", "0", "edv"))
ypos<-0.3;xpos=1.4
points(xpos,ypos,pch=16,cex=1.4,col=rgb(0.4,0.4,0.4,alpha=0.6))
text(xpos+0.19*xscal,ypos,"Females",col="grey40",cex=1.2)
points(xpos,ypos-0.08*scal,pch=17,cex=1.4,col=rgb(0.4,0.4,0.4,alpha=0.6))
text(xpos+0.19*xscal,ypos-0.08*scal,"Males",col="grey40",cex=1.2)
### 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<-exp(dif)
yla<-"Fold-difference grip strength"
tckk=-0.018
ste<-seq(0.4,1.6,by=0.2)
mons_poste(dif,1)
zpos=seq(0,1,1/3)
xx=1;ind=0.3;xpol<-0.6
text(xpol+1,zpos[xx]+zpos[2]*ind,
"Wt: ed x 0 ",cex=1)
xx=xx+1
text(xpol+1,zpos[xx]+zpos[2]*ind,
"SCA1: ed x 0 ",cex=1)
xx=xx+1
text(xpol+0.2,zpos[xx]+zpos[2]*ind,
"0: Wt x SCA1",cex=1)