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
Hippocampal and cerebellar volume
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
This script analysis of hippocampal and cerebellar (molecular layer) volumes. For hippocampus, volume were analysed separately for both hippocampus with subject (mouse, id) as random-effect factor.
For hippocampus, volume was already estimated to be smaller in SCA mice (Tichanek et al.) and this will be reflected in prior specification for the genotype effect
1 Upload of packages
2 Data upload and wrangling
## Data upload
load("source_data/myEnvironment.RData")
## new and modified variable for hippocampal dataset
hp_vol$id <- as.factor(hp_vol$id)
hp_vol$group <- factor(hp_vol$group, levels=c('ctrl.wt','ctrl.sca','eda.sca'))
## new variable for cerebellar dataset
cbml$group<-relevel(cbml$group,ref='ctrl.sca')
## data summary
summary(hp_vol) id treatment genotype sex side volume group
9 : 2 ctrl:16 sca:16 f:12 l:12 Min. :1.184 ctrl.wt :8
11 : 2 eda : 8 wt : 8 m:12 p:12 1st Qu.:1.688 ctrl.sca:8
13 : 2 Median :1.899 eda.sca :8
15 : 2 Mean :1.827
18 : 2 3rd Qu.:2.011
19 : 2 Max. :2.303
(Other):12
sex_num group_sex
Min. :-0.5 ctrl.wt.f :4
1st Qu.:-0.5 ctrl.wt.m :4
Median : 0.0 ctrl.sca.f:4
Mean : 0.0 ctrl.sca.m:4
3rd Qu.: 0.5 eda.sca.f :4
Max. : 0.5 eda.sca.m :4
summary(cbml) id treatment genotype sex cb_ml group
9 :1 ctrl:8 sca:8 f:6 Min. :16.56 ctrl.sca:4
11 :1 eda :4 wt :4 m:6 1st Qu.:18.22 ctrl.wt :4
13 :1 Median :18.56 eda.sca :4
15 :1 Mean :18.89
18 :1 3rd Qu.:19.47
19 :1 Max. :21.64
(Other):6
sex_num group_sex volume
Min. :-0.5 ctrl.wt.f :2 Min. :16.56
1st Qu.:-0.5 ctrl.wt.m :2 1st Qu.:18.22
Median : 0.0 ctrl.sca.f:2 Median :18.56
Mean : 0.0 ctrl.sca.m:2 Mean :18.89
3rd Qu.: 0.5 eda.sca.f :2 3rd Qu.:19.47
Max. : 0.5 eda.sca.m :2 Max. :21.64
3 Hippocampal volume
3.1 Data exploration
## setting spacing
par(mfrow=c(1,3))
par(mar=c(4,3.5,1,1))
par(mgp=c(4,0.6,0))
## both sex together
plot(hp_vol$volume~hp_vol$group,
col=cola, cex.axis=0.8)
## by both group and sex
plot(hp_vol$volume~hp_vol$group_sex,
col=cola[c(1,1,2,2,3,3,4,4)],cex.axis=0.8, las=2)
## data distribution
hist(hp_vol$volume, 12)
Some values seem to be weirdly smaller and distribution seems to include outliers on low values, or heavy left-tail. This suggest that volume may be better modelled with robust regression.
3.2 Modelling
3.2.1 Prirs specifications
Summary statistics to get priors
hp_vol$group<-relevel(hp_vol$group,ref='ctrl.sca')
hp_scirep <- glmmTMB(DG_ML_volume~genotype,data=young_scirep)
summary(hp_scirep) Family: gaussian ( identity )
Formula: DG_ML_volume ~ genotype
Data: young_scirep
AIC BIC logLik deviance df.resid
2.9 7.6 1.6 -3.1 33
Dispersion estimate for gaussian family (sigma^2): 0.0537
Conditional model:
Estimate Std. Error z value Pr(>|z|)
(Intercept) 2.66767 0.05182 51.48 < 2e-16 ***
genotypew 0.52591 0.07773 6.77 1.32e-11 ***
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1sd(young_scirep$DG_ML_volume, na.rm=TRUE)[1] 0.354233sd(hp_vol$volume)*c(5, 2, 1.2, 1)[1] 1.4549249 0.5819699 0.3491820 0.2909850mean(hp_vol$volume)[1] 1.827024prios(hp_scirep) Estimate
0.3735612 0.5603419 In previous publication, the volume was measured in molecular layer of dentate gyrus, so the methodology was slightly different and SD differ. The estimated prior will be thus adjusted for different variability
0.37 *(0.29/0.354)[1] 0.30310730.56 *(0.29/0.354)[1] 0.4587571Setting priors
prior01 <- c(set_prior("student_t(3, 1.83, 1.45)", class = "b", coef = "Intercept"),
set_prior("normal(0.3, 0.46)", class = "b", coef = "groupctrl.wt"),
set_prior("normal(0, 0.34)" , class = "b", coef = "groupeda.sca"))3.2.2 Model fitting
Robust will be preffered, but compare it also with simpler Gaussian
m_hp_vol <- run_model(
brm(volume~0+Intercept+group+ (1|id),
data=hp_vol, prior = prior01,
family=student(),
save_pars = save_pars(all = TRUE),
iter=8000, warmup=2000,chains=2,cores=2,seed=17,
control = list(adapt_delta = 0.99)),
'output/brm/m_hp_vol', reuse = TRUE)
m_hp_vol_g <- run_model(
brm(volume~0+Intercept+group+ (1|id),
data=hp_vol, prior = prior01,
family= gaussian(),
save_pars = save_pars(all = TRUE),
iter=8000, warmup=2000,chains=2,cores=2,seed=17,
control = list(adapt_delta = 0.99)),
'output/brm/m_hp_vol_g', reuse = TRUE)
## chains for crucial fixed-effect predictors
mcmc_trace(m_hp_vol, pars = c('b_Intercept',
'b_groupctrl.wt',
'b_groupeda.sca'))
mcmc_trace(m_hp_vol_g, pars = c('b_Intercept',
'b_groupctrl.wt',
'b_groupeda.sca'))
## summary of models
summary(m_hp_vol) Family: student
Links: mu = identity; sigma = identity; nu = identity
Formula: volume ~ 0 + Intercept + group + (1 | id)
Data: hp_vol (Number of observations: 24)
Draws: 2 chains, each with iter = 8000; warmup = 2000; thin = 1;
total post-warmup draws = 12000
Group-Level Effects:
~id (Number of levels: 12)
Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
sd(Intercept) 0.09 0.07 0.00 0.26 1.00 6548 6181
Population-Level Effects:
Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
Intercept 1.69 0.11 1.47 1.91 1.00 8355 8230
groupctrl.wt 0.26 0.15 -0.04 0.56 1.00 9048 7727
groupeda.sca 0.16 0.15 -0.14 0.46 1.00 9868 8821
Family Specific Parameters:
Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
sigma 0.28 0.05 0.19 0.40 1.00 13725 8363
nu 21.80 14.22 4.03 57.72 1.00 14526 6822
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(m_hp_vol_g) Family: gaussian
Links: mu = identity; sigma = identity
Formula: volume ~ 0 + Intercept + group + (1 | id)
Data: hp_vol (Number of observations: 24)
Draws: 2 chains, each with iter = 8000; warmup = 2000; thin = 1;
total post-warmup draws = 12000
Group-Level Effects:
~id (Number of levels: 12)
Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
sd(Intercept) 0.09 0.07 0.00 0.26 1.00 5322 4899
Population-Level Effects:
Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
Intercept 1.69 0.11 1.48 1.92 1.00 7954 7076
groupctrl.wt 0.25 0.16 -0.05 0.56 1.00 8493 8242
groupeda.sca 0.16 0.15 -0.14 0.45 1.00 8773 7281
Family Specific Parameters:
Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
sigma 0.30 0.05 0.22 0.42 1.00 11829 7875
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).Both shows good convergence, sufficient ESS and practically the same results
3.2.3 PPC
pla<-pp_check(m_hp_vol,type='dens_overlay',ndraws = 50)
plb<-pp_check(m_hp_vol_g,type='dens_overlay',ndraws = 50)
plot_grid(pla, plb, labels=c("A", "B"), ncol = 1, nrow = 2)
pla<-pp_check(m_hp_vol,type='dens_overlay_grouped',ndraws = 50,group='group')
plb<-pp_check(m_hp_vol_g,type='dens_overlay_grouped',ndraws = 50,group='group')
plot_grid(pla, plb, labels=c("A", "B"), ncol = 1, nrow = 2)
pla<-pp_check(m_hp_vol,type='scatter_avg') Using all posterior draws for ppc type 'scatter_avg' by default.plb<-pp_check(m_hp_vol_g,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(m_hp_vol,type="stat_2d", stat = c("max", "min"),ndraws=200)
plb<-pp_check(m_hp_vol_g,type="stat_2d", stat = c("max", "min"),ndraws=200)
plot_grid(pla, plb, labels=c("A", "B"), ncol = 1, nrow = 2)
pla<-pp_check(m_hp_vol,type="stat_2d", stat = c("mean", "sd"),ndraws=200)
plb<-pp_check(m_hp_vol_g,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(m_hp_vol));loo(m_hp_vol)
Computed from 12000 by 24 log-likelihood matrix
Estimate SE
elpd_loo -7.3 3.2
p_loo 5.1 1.0
looic 14.6 6.3
------
Monte Carlo SE of elpd_loo is 0.0.
Pareto k diagnostic values:
Count Pct. Min. n_eff
(-Inf, 0.5] (good) 20 83.3% 3609
(0.5, 0.7] (ok) 4 16.7% 1882
(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(m_hp_vol_g));loo(m_hp_vol_g)Warning: Found 1 observations with a pareto_k > 0.7 in model 'm_hp_vol_g'. 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 'm_hp_vol_g'. It
is recommended to set 'moment_match = TRUE' in order to perform moment matching
for problematic observations.
Computed from 12000 by 24 log-likelihood matrix
Estimate SE
elpd_loo -7.2 3.1
p_loo 4.9 1.1
looic 14.3 6.2
------
Monte Carlo SE of elpd_loo is NA.
Pareto k diagnostic values:
Count Pct. Min. n_eff
(-Inf, 0.5] (good) 17 70.8% 1473
(0.5, 0.7] (ok) 6 25.0% 1648
(0.7, 1] (bad) 1 4.2% 1964
(1, Inf) (very bad) 0 0.0% <NA>
See help('pareto-k-diagnostic') for details.Both show good PPC. However, as expected, Gaussian model show overly influential observations (low-values outliers). Robust regression is thus natural choice.
3.2.4 Posterior draws extraction
post_fix_hp_vol <- as.data.frame(m_hp_vol, variable = c("b_Intercept","b_groupctrl.wt",
"b_groupeda.sca"), )
names(post_fix_hp_vol)[1]<-"sca_ctrl"
post_fix_hp_vol$wt_ctrl<-post_fix_hp_vol$sca_ctrl+post_fix_hp_vol$b_groupctrl.wt
post_fix_hp_vol$sca_eda<-post_fix_hp_vol$sca_ctrl+post_fix_hp_vol$b_groupeda.sca
summary(post_fix_hp_vol) sca_ctrl b_groupctrl.wt b_groupeda.sca wt_ctrl
Min. :1.253 Min. :-0.3184 Min. :-0.44631 Min. :1.409
1st Qu.:1.621 1st Qu.: 0.1555 1st Qu.: 0.06417 1st Qu.:1.876
Median :1.692 Median : 0.2562 Median : 0.16171 Median :1.951
Mean :1.693 Mean : 0.2568 Mean : 0.16213 Mean :1.950
3rd Qu.:1.765 3rd Qu.: 0.3559 3rd Qu.: 0.26127 3rd Qu.:2.024
Max. :2.181 Max. : 0.9485 Max. : 0.77453 Max. :2.523
sca_eda
Min. :1.193
1st Qu.:1.783
Median :1.856
Mean :1.855
3rd Qu.:1.928
Max. :2.511 groupquant_hp_vol<-sapply(post_fix_hp_vol[,c(4,1,5)],
function(p) quantile(p, probs = c(0.025,0.975,0.5)))3.3 Data and model visualization
In data visualization, individual points represent average per subject. Violin plot (area indicating data distribution) was, in contrast, constructed on the basis of individual data points (non-averaged, 1 animals thus generated two data-points)
suppressWarnings(suppressMessages( {
range<-c(0,2.5);scal<-range[2]-range[1]
par(mfrow=c(1,2))
par(mar=c(2,3,0,0))
par(mgp=c(2,0.6,0))
hp_vol$group<-factor(hp_vol$group,levels=c('ctrl.wt','ctrl.sca','eda.sca'))
xrange<-c(0.5,3.5);xscal=xrange[2]-xrange[1]
plot(NULL,xlim=c(xrange[1],xrange[2]),ylim=c(range[1],range[2]),xlab="",
ylab="Hippocampal (DG-ML) volume (mm^3)",las=1, axes=FALSE,cex.lab=1.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.5;if(x>range[2]){break}}
x<-0
repeat{
lines(c(x,x),c(range[1],range[2]),col="white",lwd=0.7)
x=x+1;if(x>3.5){break}}
vioplot(hp_vol$vol~hp_vol$group,col=colb[-2],at=c(1:4),add=T,
border=cola[-2],axes=FALSE,drawRect=F,lwd=1,wex=0.8)
hp_vol_sum <- hp_vol %>% group_by(id, group, sex_num) %>% summarize(hpvoll = mean(volume))
beeswarm(hp_vol_sum$hpvoll~hp_vol_sum$group,col=cola[-2],at=c(1:3),
add=T,pwpch=(16.5+hp_vol_sum$sex_num), cex=1.2 )
tp<-(groupquant_hp_vol[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_hp_vol[1,xx],groupquant_hp_vol[2,xx]),lwd=1.7,col="black")
xx<-xx+1
if(xx>3){break}}
tckk=-0.02
axis(2,las=2,cex.axis=1.1,at=seq(range[1],range[2],by=0.5),
labels=c(rep("",length(seq(range[1],range[2],by=0.5)))),pos=xrange[1],tck=tckk)
axis(2,las=2,cex.axis=1.1,at=seq(range[1],range[2],by=0.5),pos=xrange[1],tck=tckk)
axis(side=1,las=1,cex.axis=1.1,at=c(seq(1,3,by=1)),
labels=c(rep("",length(seq(1,3,by=1)))),pos=range[1],tck=tckk)
lines(c(xrange[1],xrange[2]),c(range[1],range[1]))
text(c(1:3),c(rep(range[1]-0.046*scal,3)),xpd=T,cex=1.1,col=cola[-2],
c("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[-2],
c("0", "0", "edv"))
ypos<-0.6;xpos=1
points(xpos,ypos,pch=16,cex=1.4,col=rgb(0.4,0.4,0.4,alpha=0.6))
text(xpos+0.24*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.24*xscal,ypos-0.08*scal,"Males",col="grey40",cex=1.2)
### probability distributions --------
par(mar=c(2.2,1,0,0))
dif<-data.frame(post_fix_hp_vol$b_groupeda.sca,
post_fix_hp_vol$b_groupctrl.wt)
dif<-(dif)
yla<-"Difference in HP-ML volume"
tckk=-0.018
ste<-seq(-0.8, 1.2,by=0.4)
mons_poste(dif,0)
zpos=seq(-0.05,1,1/2)
xx=1;ind=0.6;xpol<- -0.32
text(xpol,zpos[xx]+zpos[2]*ind,
"SCA1: ed x 0 ",cex=1, srt = 0)
xx=xx+1
text(xpol,zpos[xx]+zpos[2]*ind,
"0: Wt x SCA1",cex=1, srt = 0)
} ) )
4 Cerebellar volumes
The analysis will be done in the same way as for hippocampus
5 Data exploration
## setting spacing
par(mfrow=c(1,3))
par(mar=c(4,3.5,1,1))
par(mgp=c(4,0.6,0))
## both sex together
plot(cbml$volume~cbml$group,
col=cola, cex.axis=0.8)
## by both group and sex
plot(cbml$volume~cbml$group_sex,
col=cola[c(1,1,2,2,3,3,4,4)],cex.axis=0.8, las=2)
## data distribution
hist(cbml$volume, 12)
Similar patter as for the hippocampus
5.1 Modelling
5.1.1 Setting priors
Summary statistics for prior specification
sd(cbml$volume)*c(5, 2, 1.2, 1)[1] 7.871047 3.148419 1.889051 1.574209 mean(cbml$volume)[1] 18.89216Prior specification
cbml$group<-relevel(cbml$group,ref='ctrl.sca')
prior01 <- c(set_prior("student_t(3, 19, 8)", class = "b", coef = "Intercept"),
set_prior("normal(0, 3.15)", class = "b", coef = "groupctrl.wt"),
set_prior("normal(0, 1.9)" , class = "b", coef = "groupeda.sca"))5.1.2 Fitting the model
m_cbml <- run_model(brm(volume~0+Intercept+group,
data=cbml, prior = prior01,
family=student(),
save_pars = save_pars(all = TRUE),
iter=8000, warmup=2000,chains=2,cores=2,seed=17,
control = list(adapt_delta = 0.999)),
'output/brm/m_cbml', reuse = TRUE)
## chains mixing and convergence
mcmc_trace(m_cbml, pars = c('b_Intercept',
'b_groupctrl.wt',
'b_groupeda.sca'))
## model summary
summary(m_cbml) Family: student
Links: mu = identity; sigma = identity; nu = identity
Formula: volume ~ 0 + Intercept + group + (1 | id)
Data: cbml (Number of observations: 12)
Draws: 2 chains, each with iter = 8000; warmup = 2000; thin = 1;
total post-warmup draws = 12000
Group-Level Effects:
~id (Number of levels: 12)
Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
sd(Intercept) 0.88 0.51 0.05 1.96 1.00 1030 2780
Population-Level Effects:
Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
Intercept 18.28 0.64 17.02 19.57 1.00 4384 6571
groupctrl.wt 2.12 0.94 0.15 3.92 1.00 4765 3678
groupeda.sca -0.23 0.87 -1.98 1.49 1.00 5277 5789
Family Specific Parameters:
Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
sigma 0.84 0.48 0.09 1.87 1.01 568 649
nu 21.01 14.12 3.25 55.83 1.00 4931 2463
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 are fine
5.1.3 PPC
pp_check(m_cbml,type='dens_overlay',ndraws = 50)
pp_check(m_cbml,type='dens_overlay_grouped',ndraws = 50,group='group')
pp_check(m_cbml,type='scatter_avg', ndraws = 50)
pp_check(m_cbml,type="stat_2d", stat = c("max", "min"),ndraws=100)
pp_check(m_cbml,type="stat_2d", stat = c("mean", "sd"),ndraws=100)
5.1.4 Extraction of posterior samples
post_fix_cbml <- as.data.frame(m_cbml, variable = c("b_Intercept","b_groupctrl.wt",
"b_groupeda.sca"), )
names(post_fix_cbml)[1]<-"sca_ctrl"
post_fix_cbml$wt_ctrl<-post_fix_cbml$sca_ctrl+post_fix_cbml$b_groupctrl.wt
post_fix_cbml$sca_eda<-post_fix_cbml$sca_ctrl+post_fix_cbml$b_groupeda.sca
summary(post_fix_cbml) sca_ctrl b_groupctrl.wt b_groupeda.sca wt_ctrl
Min. :14.98 Min. :-2.906 Min. :-4.3060 Min. :16.49
1st Qu.:17.86 1st Qu.: 1.537 1st Qu.:-0.7981 1st Qu.:19.97
Median :18.27 Median : 2.134 Median :-0.2444 Median :20.42
Mean :18.28 Mean : 2.118 Mean :-0.2324 Mean :20.40
3rd Qu.:18.69 3rd Qu.: 2.726 3rd Qu.: 0.3276 3rd Qu.:20.85
Max. :21.09 Max. : 7.393 Max. : 3.3551 Max. :23.57
sca_eda
Min. :14.65
1st Qu.:17.63
Median :18.04
Mean :18.05
3rd Qu.:18.46
Max. :21.31 groupquant_cbml<-sapply(post_fix_cbml[,c(4,1,5)],
function(p) quantile(p, probs = c(0.025,0.975,0.5)))5.2 Data and model visualization
suppressWarnings(suppressMessages( {
range<-c(12,24);scal<-range[2]-range[1]
par(mfrow=c(1,2))
par(mar=c(2,3,0.2,0.2))
par(mgp=c(2,0.6,0))
cbml$group<-factor(cbml$group,levels=c('ctrl.wt','ctrl.sca','eda.sca'))
xrange<-c(0.5,3.5);xscal=xrange[2]-xrange[1]
plot(NULL,xlim=c(xrange[1],xrange[2]),ylim=c(range[1],range[2]),xlab="",
ylab="Cerebellar (ML) volume (mm^3)",las=1, axes=FALSE,cex.lab=1.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+3;if(x>range[2]){break}}
x<-0
repeat{
lines(c(x,x),c(range[1],range[2]),col="white",lwd=0.7)
x=x+1;if(x>3.5){break}}
vioplot(cbml$vol~cbml$group,col=colb[-2],at=c(1:4),add=T,
border=cola[-2],axes=FALSE,drawRect=F,lwd=1,wex=0.8)
cbml_sum <- cbml %>% group_by(id, group, sex_num) %>% summarize(cbvoll = mean(volume))
beeswarm(cbml_sum$cbvoll~cbml_sum$group,col=cola[-2],at=c(1:3),
add=T,pwpch=(16.5+cbml_sum$sex_num), cex=1.2 )
tp<-(groupquant_cbml[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_cbml[1,xx],groupquant_cbml[2,xx]),lwd=1.7,col="black")
xx<-xx+1
if(xx>3){break}}
tckk=-0.02
axis(2,las=2,cex.axis=1.1,at=c(seq(range[1],14,by=3),14),
labels=c('0','2'),pos=xrange[1],tck=tckk)
axis(2,las=2,cex.axis=1.1,at=c(seq(15,range[2],by=3)),,pos=xrange[1],tck=tckk)
axis(side=1,las=1,cex.axis=1.1,at=c(seq(1,3,by=1)),
labels=c(rep("",length(seq(1,3,by=1)))),pos=range[1],tck=tckk)
lines(c(xrange[1],xrange[2]),c(range[1],range[1]))
text(c(1:3),c(rep(range[1]-0.046*scal,3)),xpd=T,cex=1.1,col=cola[-2],
c("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[-2],
c("0", "0", "edv"))
ypos<-14;xpos=1
points(xpos,ypos,pch=16,cex=1.4,col=rgb(0.4,0.4,0.4,alpha=0.6))
text(xpos+0.24*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.24*xscal,ypos-0.08*scal,"Males",col="grey40",cex=1.2)
### probability distributions --------
ypos<-0.6;xpos=1
range=c(0, 24)
scal = 12
par(mar=c(2.2,1,0,0))
dif<-data.frame(post_fix_cbml$b_groupeda.sca,
post_fix_cbml$b_groupctrl.wt)
dif<-(dif)
yla<- "Difference in CB-ML volume"
tckk= -0.018
ste<- seq(-4, 6, by=2)
mons_poste(dif,0)
zpos=seq(-0.05,1,1/2)
xx=1;ind=0.8;xpol<- -3.05
text(xpol,zpos[xx]+zpos[2]*ind,
"SCA1: ed x 0 ",cex=1, srt = 0)
xx=xx+1
text(xpol,zpos[xx]+zpos[2]*ind,
"0: Wt x SCA1",cex=1, srt = 0)
} ) )