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 ELISA (brain IL6 and BDNF levels)
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
ELISA was used to measure:
- hippocampal levels of interleukin 6 (IL6, marker of inflammation)
- cerebellar level of IL6
- cerebellar level of Brain-derived neurotrophic factor (BDNF, neuroplasticity-related protein)
1 Upload of packages
2 Data upload and wrangling
## Data upload
load("source_data/myEnvironment.RData")
elisa<-read.csv("source_data/elisa_data.csv",stringsAsFactors = T)
summary(elisa) id treatment genotype sex tissue il6
Min. : 1.00 ctrl:32 sca:32 f:32 crbl:32 Min. :24.09
1st Qu.: 42.75 eda :32 wt :32 m:32 hip :32 1st Qu.:34.57
Median : 62.50 Median :47.29
Mean : 58.34 Mean :43.97
3rd Qu.: 76.25 3rd Qu.:51.47
Max. :121.00 Max. :58.45
bdnf X
Min. : 9.907 Mode:logical
1st Qu.:13.477 NA's:64
Median :15.211
Mean :15.702
3rd Qu.:17.591
Max. :25.690
NA's :32 ## Defining new variables
elisa$group<-interaction(elisa$treatment,elisa$genotype)
elisa$group<-factor(elisa$group,levels=c('ctrl.wt','eda.wt','ctrl.sca','eda.sca'))
elisa$sex_num<-ifelse(elisa$sex=="f",-0.5,0.5)
elisa$group_sex<-interaction(elisa$treatment,elisa$genotype,elisa$sex)
elisa$group_sex<-factor(elisa$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'))
## subseting datasets
elisa_hp<-subset(elisa,elisa$tissue=='hip')
elisa_cb<-subset(elisa,elisa$tissue=='crbl')
elisa_bdnf <- subset(elisa,elisa$bdnf != 'NA')3 Hippocampal IL6
3.1 Data exploration
par(mar=c(5,3.5,1,1))
par(mgp=c(1.9,0.7,0))
par(mfrow=c(3,1))
plot(elisa_hp$il6~elisa_hp$group,las=2,xlab="",
col=cola)
plot(elisa_hp$il6~elisa_hp$group_sex,las=2,xlab="",
col=cola[c(1,1,2,2,3,3,4,4)])
hist(elisa_hp$il6, 12)
There is no consistent sex-related difference. Groups with larger value show larger spread and the distribution seems right-tailed. Thus gamma model seems as natural choice. It will be used along Gaussian models and they will be compared with posterior predictive check (PPC)
3.2 Modelling
3.2.1 Setting priors
Summary statistics and other information used for prior specification
elisa_hp$group<-relevel(elisa_hp$group,ref='ctrl.sca')
## prior for gamma models
log(mean(elisa_hp$il6))[1] 3.626034sd(log(elisa_hp$il6))*5[1] 1.233038## prior for Gaussian model
mean(elisa_hp$il6)[1] 37.56354sd(elisa_hp[elisa_hp$genotype=="sca",]$il6) * c(5, 2, 1.2)[1] 41.88630 16.75452 10.05271sd(elisa_hp[elisa_hp$genotype=="wt",]$il6) * c(5, 2, 1.2)[1] 54.89251 21.95700 13.17420sd(elisa_hp$il6)* c(5, 2, 1.2)[1] 48.54912 19.41965 11.65179Defining priors
## Priors for gamma model
prior01 <- c(set_prior("student_t(3, 3.62, 1.23)", 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"))
## Priors for Gaussian model
prior02 <- c(set_prior("student_t(3, 37, 45)", class = "b", coef = "Intercept"),
set_prior("normal(0,22)", class = "b", coef = "groupctrl.wt"),
set_prior("normal(0,22)", class = "b", coef = "groupeda.wt"),
set_prior("normal(0,11)" , class = "b", coef = "groupeda.sca"))3.2.2 Fitting models and convergence check
## gamma model
elisa_hp_model_01 <- run_model(
brm(il6 ~0+Intercept+group,
data=elisa_hp, 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/elisa_hp_model_01', reuse = TRUE)
## Gaussian model
elisa_hp_model_02 <- run_model(
brm(il6~0+Intercept+group,
data=elisa_hp, prior = prior02,
family=gaussian(),
save_pars = save_pars(all = TRUE),
iter=8000, warmup=2000,chains=2,cores=1,seed=17,
control = list(adapt_delta = 0.99)),
'output/brm/elisa_hp_model_02', reuse = TRUE)
mcmc_trace(elisa_hp_model_01)
mcmc_trace(elisa_hp_model_02)
summary(elisa_hp_model_01) Family: gamma
Links: mu = log; shape = identity
Formula: il6 ~ 0 + Intercept + group
Data: elisa_hp (Number of observations: 32)
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 3.60 0.10 3.42 3.79 1.00 3408 4858
groupctrl.wt 0.04 0.13 -0.23 0.30 1.00 3865 5309
groupeda.wt 0.09 0.14 -0.18 0.36 1.00 3777 5656
groupeda.sca -0.02 0.13 -0.28 0.24 1.00 4223 5287
Family Specific Parameters:
Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
shape 14.87 3.92 8.20 23.64 1.00 5957 6776
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(elisa_hp_model_02) Family: gaussian
Links: mu = identity; sigma = identity
Formula: il6 ~ 0 + Intercept + group
Data: elisa_hp (Number of observations: 32)
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 36.48 3.35 29.99 43.13 1.00 3551 5467
groupctrl.wt 1.50 4.84 -8.09 10.81 1.00 4549 6723
groupeda.wt 3.39 4.80 -6.16 12.72 1.00 4494 6085
groupeda.sca -0.59 4.60 -9.66 8.32 1.00 4388 6173
Family Specific Parameters:
Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
sigma 10.17 1.35 7.91 13.17 1.00 5930 6789
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 of both models is fine
3.2.3 PPC
pla<-pp_check(elisa_hp_model_01,type='dens_overlay',ndraws = 50)
plb<-pp_check(elisa_hp_model_02,type='dens_overlay',ndraws = 50)
plot_grid(pla, plb, labels=c("A", "B"), ncol = 1, nrow = 2)
pla<-pp_check(elisa_hp_model_01,type='dens_overlay_grouped',ndraws = 50,group='group')
plb<-pp_check(elisa_hp_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(elisa_hp_model_01,type='scatter_avg') Using all posterior draws for ppc type 'scatter_avg' by default.plb<-pp_check(elisa_hp_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(elisa_hp_model_01,type='ecdf_overlay',ndraws = 200) 
pp_check(elisa_hp_model_02,type='ecdf_overlay',ndraws = 200) 
pla<-pp_check(elisa_hp_model_01,type='ecdf_overlay_grouped',ndraws = 50,group='group')
plb<-pp_check(elisa_hp_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(elisa_hp_model_01,type="stat_2d", stat = c("max", "min"),ndraws=200)
plb<-pp_check(elisa_hp_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(elisa_hp_model_01,type="stat_2d", stat = c("mean", "sd"),ndraws=200)
plb<-pp_check(elisa_hp_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(elisa_hp_model_01));loo(elisa_hp_model_01)
Computed from 12000 by 32 log-likelihood matrix
Estimate SE
elpd_loo -120.5 3.3
p_loo 4.3 0.7
looic 241.1 6.5
------
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(elisa_hp_model_02));loo(elisa_hp_model_02)
Computed from 12000 by 32 log-likelihood matrix
Estimate SE
elpd_loo -122.0 2.9
p_loo 4.1 0.7
looic 244.0 5.7
------
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 fit the data relatively well. For simplicity, we will use Gaussian models
3.2.4 Extraction of posterior samples
post_fix<-as.data.frame(elisa_hp_model_02, 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. :23.88 Min. :-19.364 Min. :-18.8974 Min. :-18.5055
1st Qu.:34.24 1st Qu.: -1.665 1st Qu.: 0.2454 1st Qu.: -3.6530
Median :36.45 Median : 1.461 Median : 3.3817 Median : -0.6670
Mean :36.48 Mean : 1.497 Mean : 3.3911 Mean : -0.5931
3rd Qu.:38.73 3rd Qu.: 4.718 3rd Qu.: 6.5097 3rd Qu.: 2.5098
Max. :48.88 Max. : 20.800 Max. : 25.5243 Max. : 16.4948
wt_ctrl wt_eda sca_eda wt_contrast
Min. :22.59 Min. :22.99 Min. :23.53 Min. :-17.600
1st Qu.:35.66 1st Qu.:37.53 1st Qu.:33.56 1st Qu.: -1.412
Median :37.94 Median :39.87 Median :35.84 Median : 1.891
Mean :37.97 Mean :39.87 Mean :35.88 Mean : 1.894
3rd Qu.:40.34 3rd Qu.:42.22 3rd Qu.:38.18 3rd Qu.: 5.264
Max. :52.70 Max. :57.23 Max. :50.85 Max. : 22.082 groupquant<-sapply(post_fix[,c(5,6,1,7)],function(p) quantile(p, probs = c(0.025,0.975,0.5)))
groupquant<-(groupquant)3.3 Vizualisation
## showing data
elisa_hp$group<-factor(elisa_hp$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,60);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="IL6 (pg/mg of proteins)",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+10;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(elisa_hp$il6~elisa_hp$group,col=colb,at=c(1:4),add=T,
border=cola,axes=FALSE,drawRect=F,lwd=1,wex=0.8)
beeswarm(elisa_hp$il6~elisa_hp$group,col=colc,at=c(1:4),
add=T,pwpch=(16.5+elisa_hp$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=20),
labels=c(rep("",length(seq(range[1],range[2],by=20)))),pos=xrange[1],tck=tckk)
axis(2,las=2,cex.axis=1.1,at=seq(range[1],range[2],by=20),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<-17;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)
## posterior probability distributions for between-groups difference
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 hippocampal IL6"
tckk=-0.018
ste<-seq(-20,15,by=5)
mons_poste(dif,0)
zpos=seq(0,1,1/3)
xx=1;ind=0.4;xpol<--15
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)
4 Cerebellar IL6
4.1 Data exploration
par(mar=c(5,3.5,1,1))
par(mgp=c(1.9,0.7,0))
par(mfrow=c(3,1))
plot(elisa_cb$il6~elisa_cb$group,las=2,xlab="",
col=cola)
plot(elisa_cb$il6~elisa_cb$group_sex,las=2,xlab="",
col=cola[c(1,1,2,2,3,3,4,4)])
hist(elisa_cb$il6, 12)
No consistent sex effect. Distribution seems relatively symetrical. We will again fit both gamma and Gaussian models
4.2 Modelling
4.2.1 Setting priors
Summary statistics and other information used for prior specification
elisa_cb$group<-relevel(elisa_cb$group,ref='ctrl.sca')
## prior for gamma models
log(mean(elisa_cb$il6))[1] 3.919395sd(log(elisa_cb$il6))*5[1] 0.5294394## prior for Gaussian model
mean(elisa_cb$il6)[1] 50.36994sd(elisa_cb[elisa_cb$genotype=="sca",]$il6) * c(5, 2, 1.2)[1] 22.569138 9.027655 5.416593sd(elisa_cb[elisa_cb$genotype=="wt",]$il6) * c(5, 2, 1.2)[1] 28.908948 11.563579 6.938147sd(elisa_cb$il6)* c(5, 2, 1.2)[1] 25.699076 10.279630 6.167778Defining priors
## Priors for gamma model
prior01 <- c(set_prior("student_t(3, 3.9, 0.53)", 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"))
## Priors for Gaussian model
prior02 <- c(set_prior("student_t(3, 50, 26)", class = "b", coef = "Intercept"),
set_prior("normal(0, 11.6)", class = "b", coef = "groupctrl.wt"),
set_prior("normal(0, 11.6)", class = "b", coef = "groupeda.wt"),
set_prior("normal(0, 5.4)" , class = "b", coef = "groupeda.sca"))4.2.2 Fitting models and convergence check
## gamma model
elisa_cb_model_01 <- run_model(
brm(il6 ~0+Intercept+group,
data=elisa_cb, 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/elisa_cb_model_01', reuse = TRUE)
## Gaussian model
elisa_cb_model_02 <- run_model(
brm(il6~0+Intercept+group,
data=elisa_cb, prior = prior02,
family=gaussian(),
save_pars = save_pars(all = TRUE),
iter=8000, warmup=2000,chains=2,cores=1,seed=17,
control = list(adapt_delta = 0.99)),
'output/brm/elisa_cb_model_02', reuse = TRUE)
mcmc_trace(elisa_cb_model_01)
mcmc_trace(elisa_cb_model_02)
summary(elisa_cb_model_01) Family: gamma
Links: mu = log; shape = identity
Formula: il6 ~ 0 + Intercept + group
Data: elisa_cb (Number of observations: 32)
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 3.96 0.04 3.88 4.04 1.00 3200 4824
groupctrl.wt -0.05 0.06 -0.16 0.06 1.00 3825 5585
groupeda.wt -0.06 0.06 -0.17 0.06 1.00 4043 5766
groupeda.sca -0.06 0.06 -0.18 0.05 1.00 3974 5394
Family Specific Parameters:
Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
shape 83.47 22.45 44.92 131.80 1.00 5794 6455
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(elisa_cb_model_02) Family: gaussian
Links: mu = identity; sigma = identity
Formula: il6 ~ 0 + Intercept + group
Data: elisa_cb (Number of observations: 32)
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 52.17 1.74 48.76 55.60 1.00 4169 6061
groupctrl.wt -2.20 2.50 -7.11 2.74 1.00 5204 6390
groupeda.wt -2.50 2.51 -7.48 2.43 1.00 4934 6637
groupeda.sca -2.55 2.36 -7.18 2.12 1.00 4943 6309
Family Specific Parameters:
Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
sigma 5.32 0.72 4.14 6.95 1.00 7125 7285
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 of both models is fine, ESSs are sufficient
4.2.3 PPC
pla<-pp_check(elisa_cb_model_01,type='dens_overlay',ndraws = 50)
plb<-pp_check(elisa_cb_model_02,type='dens_overlay',ndraws = 50)
plot_grid(pla, plb, labels=c("A", "B"), ncol = 1, nrow = 2)
pla<-pp_check(elisa_cb_model_01,type='dens_overlay_grouped',ndraws = 50,group='group')
plb<-pp_check(elisa_cb_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(elisa_cb_model_01,type='scatter_avg') Using all posterior draws for ppc type 'scatter_avg' by default.plb<-pp_check(elisa_cb_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(elisa_cb_model_01,type='ecdf_overlay',ndraws = 200) 
pp_check(elisa_cb_model_02,type='ecdf_overlay',ndraws = 200) 
pla<-pp_check(elisa_cb_model_01,type='ecdf_overlay_grouped',ndraws = 50,group='group')
plb<-pp_check(elisa_cb_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(elisa_cb_model_01,type="stat_2d", stat = c("max", "min"),ndraws=200)
plb<-pp_check(elisa_cb_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(elisa_cb_model_01,type="stat_2d", stat = c("mean", "sd"),ndraws=200)
plb<-pp_check(elisa_cb_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(elisa_cb_model_01));loo(elisa_cb_model_01)Warning: Found 1 observations with a pareto_k > 0.7 in model
'elisa_cb_model_01'. It is recommended to set 'moment_match = TRUE' in order to
perform moment matching for problematic observations.
Warning: Found 1 observations with a pareto_k > 0.7 in model
'elisa_cb_model_01'. It is recommended to set 'moment_match = TRUE' in order to
perform moment matching for problematic observations.
Computed from 12000 by 32 log-likelihood matrix
Estimate SE
elpd_loo -102.8 5.5
p_loo 5.1 2.0
looic 205.5 11.1
------
Monte Carlo SE of elpd_loo is NA.
Pareto k diagnostic values:
Count Pct. Min. n_eff
(-Inf, 0.5] (good) 31 96.9% 1417
(0.5, 0.7] (ok) 0 0.0% <NA>
(0.7, 1] (bad) 1 3.1% 216
(1, Inf) (very bad) 0 0.0% <NA>
See help('pareto-k-diagnostic') for details.plot(loo(elisa_cb_model_02));loo(elisa_cb_model_02)
Computed from 12000 by 32 log-likelihood matrix
Estimate SE
elpd_loo -101.6 5.4
p_loo 4.9 1.8
looic 203.1 10.8
------
Monte Carlo SE of elpd_loo is 0.1.
Pareto k diagnostic values:
Count Pct. Min. n_eff
(-Inf, 0.5] (good) 31 96.9% 1650
(0.5, 0.7] (ok) 1 3.1% 330
(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.Both models seem to fit the data relatively well. One observation is apparently outlier, with overly large impact on estimation. The large influence is seen particularly slightly more strongly in gamma model. Again, we will prefer Gaussian model, although both draw the same conlcusions.
4.2.4 Extraction of posterior samples
post_fix <- as.data.frame(elisa_cb_model_02, 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. :44.42 Min. :-12.4358 Min. :-12.8186 Min. :-12.1552
1st Qu.:51.02 1st Qu.: -3.8251 1st Qu.: -4.1481 1st Qu.: -4.1027
Median :52.17 Median : -2.1950 Median : -2.4996 Median : -2.5641
Mean :52.17 Mean : -2.1976 Mean : -2.4955 Mean : -2.5513
3rd Qu.:53.31 3rd Qu.: -0.5862 3rd Qu.: -0.8297 3rd Qu.: -0.9773
Max. :58.56 Max. : 8.1742 Max. : 8.6824 Max. : 6.8822
wt_ctrl wt_eda sca_eda wt_contrast
Min. :42.80 Min. :42.24 Min. :41.30 Min. :-10.4207
1st Qu.:48.75 1st Qu.:48.42 1st Qu.:48.45 1st Qu.: -2.0354
Median :49.97 Median :49.69 Median :49.62 Median : -0.3192
Mean :49.97 Mean :49.67 Mean :49.62 Mean : -0.2980
3rd Qu.:51.18 3rd Qu.:50.92 3rd Qu.:50.79 3rd Qu.: 1.4553
Max. :58.05 Max. :57.08 Max. :56.67 Max. : 10.1149 groupquant<-sapply(post_fix[,c(5,6,1,7)],function(p) quantile(p, probs = c(0.025,0.975,0.5)))
groupquant<-(groupquant)4.3 Vizualisation
## showing data
elisa_cb$group<-factor(elisa_cb$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,60);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="IL6 (pg/mg of proteins)",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+10;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(elisa_cb$il6~elisa_cb$group,col=colb,at=c(1:4),add=T,
border=cola,axes=FALSE,drawRect=F,lwd=1,wex=0.8)
beeswarm(elisa_cb$il6~elisa_cb$group,col=colc,at=c(1:4),
add=T,pwpch=(16.5+elisa_cb$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=20),
labels=c(rep("",length(seq(range[1],range[2],by=20)))),pos=xrange[1],tck=tckk)
axis(2,las=2,cex.axis=1.1,at=seq(range[1],range[2],by=20),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<-17;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)
## posterior probability distributions for between-groups difference
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 cerebellar IL6"
tckk=-0.018
ste<-seq(-20,15,by=5)
mons_poste(dif,0)
zpos=seq(0,1,1/3)
xx=1;ind=0.4;xpol<--15
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)
5 Cerebellar BDNF
5.1 Data exploration
par(mar=c(5,3.5,1,1))
par(mgp=c(1.9,0.7,0))
par(mfrow=c(3,1))
plot(elisa_bdnf$bdnf~elisa_bdnf$group,las=2,xlab="",
col=cola)
plot(elisa_bdnf$bdnf~elisa_bdnf$group_sex,las=2,xlab="",
col=cola[c(1,1,2,2,3,3,4,4)])
hist(elisa_bdnf$bdnf, 12)
There is no consistent sex-related difference, although specifically in WT, males show trend of higher BDNF levels. Distribution is symmetrical and spread is not larger in groups with larger values. However, we will be consistent with previous ELISA analyses and will also fit both gamma and Gaussian regression
5.2 Modelling
5.2.1 Setting priors
Summary statistics and other information used for prior specification
elisa_bdnf$group<-relevel(elisa_bdnf$group,ref='ctrl.sca')
## prior for gamma models
log(mean(elisa_bdnf$bdnf))[1] 2.753806sd(log(elisa_bdnf$bdnf))*5[1] 1.003424## prior for Gaussian model
mean(elisa_bdnf$bdnf)[1] 15.70229sd(elisa_bdnf[elisa_bdnf$genotype=="sca",]$bdnf) * c(5, 2, 1.2)[1] 17.708605 7.083442 4.250065sd(elisa_bdnf[elisa_bdnf$genotype=="wt",]$bdnf) * c(5, 2, 1.2)[1] 13.823583 5.529433 3.317660sd(elisa_bdnf$bdnf)* c(5, 2, 1.2)[1] 15.963888 6.385555 3.831333Defining priors
## Priors for gamma model
prior01 <- c(set_prior("student_t(3, 2.8, 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"))
## Priors for Gaussian model
prior02 <- c(set_prior("student_t(3, 16, 16)", class = "b", coef = "Intercept"),
set_prior("normal(0, 5.5)", class = "b", coef = "groupctrl.wt"),
set_prior("normal(0, 5.5)", class = "b", coef = "groupeda.wt"),
set_prior("normal(0, 4.3)" , class = "b", coef = "groupeda.sca"))5.2.2 Fitting models and convergence check
## gamma model
elisa_bdnf_model_01 <- run_model(
brm(bdnf ~0+Intercept+group,
data=elisa_bdnf, 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/elisa_bdnf_model_01', reuse = TRUE)
## Gaussian model
elisa_bdnf_model_02 <- run_model(
brm(bdnf~0+Intercept+group,
data=elisa_bdnf, prior = prior02,
family=gaussian(),
save_pars = save_pars(all = TRUE),
iter=8000, warmup=2000,chains=2,cores=1,seed=17,
control = list(adapt_delta = 0.99)),
'output/brm/elisa_bdnf_model_02', reuse = TRUE)
mcmc_trace(elisa_bdnf_model_01)
mcmc_trace(elisa_bdnf_model_02)
summary(elisa_bdnf_model_01) Family: gamma
Links: mu = log; shape = identity
Formula: bdnf ~ 0 + Intercept + group
Data: elisa_bdnf (Number of observations: 32)
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 2.77 0.08 2.62 2.92 1.00 3768 5459
groupctrl.wt -0.03 0.11 -0.24 0.18 1.00 4584 5944
groupeda.wt -0.08 0.11 -0.29 0.13 1.00 4329 5334
groupeda.sca 0.05 0.11 -0.16 0.26 1.00 4631 6046
Family Specific Parameters:
Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
shape 23.55 6.21 13.08 37.23 1.00 5260 5868
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(elisa_bdnf_model_02) Family: gaussian
Links: mu = identity; sigma = identity
Formula: bdnf ~ 0 + Intercept + group
Data: elisa_bdnf (Number of observations: 32)
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 15.90 1.10 13.74 18.10 1.00 3528 5061
groupctrl.wt -0.40 1.57 -3.56 2.68 1.00 4177 5464
groupeda.wt -1.20 1.55 -4.30 1.81 1.00 4350 5923
groupeda.sca 0.81 1.54 -2.24 3.86 1.00 4462 5476
Family Specific Parameters:
Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
sigma 3.32 0.45 2.57 4.36 1.00 6257 7185
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 of both models is fine, ESSs are sufficient
5.2.3 PPC
pla<-pp_check(elisa_bdnf_model_01,type='dens_overlay',ndraws = 50)
plb<-pp_check(elisa_bdnf_model_02,type='dens_overlay',ndraws = 50)
plot_grid(pla, plb, labels=c("A", "B"), ncol = 1, nrow = 2)
pla<-pp_check(elisa_bdnf_model_01,type='dens_overlay_grouped',ndraws = 50,group='group')
plb<-pp_check(elisa_bdnf_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(elisa_bdnf_model_01,type='scatter_avg') Using all posterior draws for ppc type 'scatter_avg' by default.plb<-pp_check(elisa_bdnf_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(elisa_bdnf_model_01,type='ecdf_overlay',ndraws = 200) 
pp_check(elisa_bdnf_model_02,type='ecdf_overlay',ndraws = 200) 
pla<-pp_check(elisa_bdnf_model_01,type='ecdf_overlay_grouped',ndraws = 50,group='group')
plb<-pp_check(elisa_bdnf_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(elisa_bdnf_model_01,type="stat_2d", stat = c("max", "min"),ndraws=200)
plb<-pp_check(elisa_bdnf_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(elisa_bdnf_model_01,type="stat_2d", stat = c("mean", "sd"),ndraws=200)
plb<-pp_check(elisa_bdnf_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(elisa_bdnf_model_01));loo(elisa_bdnf_model_01)
Computed from 12000 by 32 log-likelihood matrix
Estimate SE
elpd_loo -85.7 4.9
p_loo 4.8 1.7
looic 171.4 9.7
------
Monte Carlo SE of elpd_loo is 0.1.
Pareto k diagnostic values:
Count Pct. Min. n_eff
(-Inf, 0.5] (good) 31 96.9% 2619
(0.5, 0.7] (ok) 1 3.1% 243
(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(elisa_bdnf_model_02));loo(elisa_bdnf_model_02)Warning: Found 1 observations with a pareto_k > 0.7 in model
'elisa_bdnf_model_02'. It is recommended to set 'moment_match = TRUE' in order
to perform moment matching for problematic observations.
Warning: Found 1 observations with a pareto_k > 0.7 in model
'elisa_bdnf_model_02'. It is recommended to set 'moment_match = TRUE' in order
to perform moment matching for problematic observations.
Computed from 12000 by 32 log-likelihood matrix
Estimate SE
elpd_loo -86.8 6.2
p_loo 5.3 2.4
looic 173.6 12.4
------
Monte Carlo SE of elpd_loo is NA.
Pareto k diagnostic values:
Count Pct. Min. n_eff
(-Inf, 0.5] (good) 31 96.9% 2646
(0.5, 0.7] (ok) 0 0.0% <NA>
(0.7, 1] (bad) 1 3.1% 46
(1, Inf) (very bad) 0 0.0% <NA>
See help('pareto-k-diagnostic') for details.Both models seem to fit the data relatively well. For simplicity, we will use Gaussian models, both models show the same - no evidence for the edaravone effect
5.2.4 Extraction of posterior samples
post_fix <- as.data.frame(elisa_bdnf_model_02,
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. :11.18 Min. :-7.3248 Min. :-7.1888 Min. :-5.4116
1st Qu.:15.18 1st Qu.:-1.4531 1st Qu.:-2.2105 1st Qu.:-0.1994
Median :15.90 Median :-0.3961 Median :-1.2011 Median : 0.8149
Mean :15.90 Mean :-0.4038 Mean :-1.2003 Mean : 0.8101
3rd Qu.:16.62 3rd Qu.: 0.6532 3rd Qu.:-0.1657 3rd Qu.: 1.8247
Max. :20.13 Max. : 6.0358 Max. : 5.3296 Max. : 7.7811
wt_ctrl wt_eda sca_eda wt_contrast
Min. :10.68 Min. :10.46 Min. :11.64 Min. :-8.3769
1st Qu.:14.75 1st Qu.:13.93 1st Qu.:15.97 1st Qu.:-1.8619
Median :15.49 Median :14.70 Median :16.71 Median :-0.7779
Mean :15.50 Mean :14.70 Mean :16.71 Mean :-0.7965
3rd Qu.:16.25 3rd Qu.:15.46 3rd Qu.:17.46 3rd Qu.: 0.2754
Max. :20.69 Max. :19.98 Max. :21.45 Max. : 5.5642 groupquant <- sapply(post_fix[,c(5,6,1,7)],function(p) quantile(p, probs = c(0.025,0.975,0.5)))5.3 Vizualisation
## showing data
## Final plot ------------
elisa_bdnf$group<-factor(elisa_bdnf$group,levels=
c('ctrl.wt','eda.wt','ctrl.sca','eda.sca'))
summary(elisa) id treatment genotype sex tissue il6
Min. : 1.00 ctrl:32 sca:32 f:32 crbl:32 Min. :24.09
1st Qu.: 42.75 eda :32 wt :32 m:32 hip :32 1st Qu.:34.57
Median : 62.50 Median :47.29
Mean : 58.34 Mean :43.97
3rd Qu.: 76.25 3rd Qu.:51.47
Max. :121.00 Max. :58.45
bdnf X group sex_num group_sex
Min. : 9.907 Mode:logical ctrl.wt :16 Min. :-0.5 ctrl.wt.f : 8
1st Qu.:13.477 NA's:64 eda.wt :16 1st Qu.:-0.5 ctrl.wt.m : 8
Median :15.211 ctrl.sca:16 Median : 0.0 eda.wt.f : 8
Mean :15.702 eda.sca :16 Mean : 0.0 eda.wt.m : 8
3rd Qu.:17.591 3rd Qu.: 0.5 ctrl.sca.f: 8
Max. :25.690 Max. : 0.5 ctrl.sca.m: 8
NA's :32 (Other) :16 par(mfrow=c(1,2))
par(mar=c(2,3,0,0))
par(mgp=c(2,0.6,0))
range<-c(0,30);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="BDNF in cerebellum (pg/mg of proteins)",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+5;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(elisa_bdnf$bdnf~elisa_bdnf$group,col=colb,at=c(1:4),add=T,
border=cola,axes=FALSE,drawRect=F,lwd=1,wex=0.8)
beeswarm(elisa_bdnf$bdnf~elisa_bdnf$group,col=colc,at=c(1:4),add=T,pwpch=(16.5+elisa_bdnf$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=5),
labels=c(rep("",length(seq(range[1],range[2],by=5)))),pos=xrange[1],tck=tckk)
axis(2,las=2,cex.axis=1.1,at=seq(range[1],range[2],by=5),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<-27;xpos=1
points(xpos,ypos,pch=16,cex=1.4,col=rgb(0.4,0.4,0.4,alpha=0.6))
text(xpos+0.23*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.23*xscal,ypos-0.08*scal,"Males",col="grey40",cex=1.2)
## posterior probability distributions for between-groups difference
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 cerebellar bdnf"
tckk=-0.018
ste<-seq(-10,10,by=2)
mons_poste(dif,0)
zpos=seq(0,1,1/3)
xx=1;ind=0.4;xpol<- -6
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)