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 mitochondrial respiration
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
Mitochodnrial respiration was analysed in cerebellar (cb) and hippocampal (hp) tissue. Respiration was standardized by either tissue wet weight of (mg) or by citrate synthase activity (cs), across different phases experiment
1 Upload of packages
2 Data upload and wrangling
## Data upload
load("source_data/myEnvironment.RData")
## connecting datasets
mitochondria <- mitochondria %>% filter(
C_iv != 'NA'
) %>% mutate(id = as.factor(id))
## defining new variables
mitochondria$group <- interaction(mitochondria$treatment,mitochondria$genotype)
mitochondria$group <- factor(mitochondria$group,
levels=c('ctrl.wt','eda.wt','ctrl.sca','eda.sca'))
mitochondria$sex_num <- ifelse(mitochondria$sex == "f", -0.5, 0.5)
mitochondria$group_sex <- interaction(mitochondria$treatment,
mitochondria$genotype,
mitochondria$sex)
mitochondria$group_sex <- factor(mitochondria$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'))
mitochondria$id <- as.factor(mitochondria$id)
names(mitochondria) [1] "id.mito" "id" "sex" "genotype"
[5] "treatment" "normalisation" "region" "Sample"
[9] "M_plus_G" "D" "c" "P"
[13] "S" "Fccp" "Rot" "Ama"
[17] "As.Tm" "Azd" "C_iv" "b"
[21] "out" "out_susp" "group" "sex_num"
[25] "group_sex" ## subseting according to sex
mitochondria_m <- subset(mitochondria, mitochondria$sex == "m")
mitochondria_f <- subset(mitochondria, mitochondria$sex == "f")
## subsetting according to brain region and standardization
mg_cb <- subset(mitochondria,
mitochondria$normalisation == 'mg' &
mitochondria$region == 'cerebellum')
cs_cb <- subset(mitochondria,
mitochondria$normalisation == 'cs' &
mitochondria$region == 'cerebellum')
mg_hp <- subset(mitochondria,
mitochondria$normalisation == 'mg' &
mitochondria$region == 'hippocampus')
cs_hp <- subset(mitochondria,
mitochondria$normalisation == 'cs' &
mitochondria$region=='hippocampus')
## re-scaling respiration standardized by citrate synthase
cs_cb[,c(9:19)] <- cs_cb[,c(9:19)]/1000
cs_hp[,c(9:19)] <- cs_hp[,c(9:19)]/10003 Hippocampus, weight standardization
3.1 Data exploration
par(mfrow=c(2,5))
par(mar=c(5,5,1,0))
par(mgp=c(3.6, 0.6,0))
## both sexes together
plot(mg_hp$P~mg_hp$group,
col=cola, cex.axis=0.8, las=2)
plot(mg_hp$S~mg_hp$group,
col=cola, cex.axis=0.8, las=2)
plot(mg_hp$Fccp~mg_hp$group,
col=cola, cex.axis=0.8, las=2)
plot(mg_hp$Rot~mg_hp$group,
col=cola, cex.axis=0.8, las=2)
plot(mg_hp$C_iv~mg_hp$group,
col=cola, cex.axis=0.8, las=2)
## by both group and sex
plot(mg_hp$P~mg_hp$group_sex,
col=cola[c(1,1,2,2,3,3,4,4)],cex.axis=0.8, las=2)
plot(mg_hp$S~mg_hp$group_sex,
col=cola[c(1,1,2,2,3,3,4,4)],cex.axis=0.8, las=2)
plot(mg_hp$Fccp~mg_hp$group_sex,
col=cola[c(1,1,2,2,3,3,4,4)],cex.axis=0.8, las=2)
plot(mg_hp$Rot~mg_hp$group_sex,
col=cola[c(1,1,2,2,3,3,4,4)],cex.axis=0.8, las=2)
plot(mg_hp$C_iv~mg_hp$group_sex,
col=cola[c(1,1,2,2,3,3,4,4)],cex.axis=0.8, las=2)
No consistent sex effect. Se will be thus ignored
Data distribution seems to be asymmetric From this reason, and also because the data are continuous and non-negative, we will employ gamma regression
3.2 Modelling
To model all 5 parameters (reflecting mitochondrial respiration capacity in difference phases of mitochondrial respiratory cycle) in a single model, we will fit multivariate regression.
As there are technical replicates (two measurements per animal and brain region), we will apply hierarchical model with random intercept.
3.2.1 Setting priors
Although we measured mitochondrial respiration in our previous publication (Tichanek et al., 2020), context was different in the sense that previous publication used group-housed mice not undergoing any test or manipulation. In contrast, these data are from mice which have been tested in behavioral tests and which were handled regularly.
For simplicity, we will use default prior for intercept
## re-leveling groups
mg_hp$group<-relevel(mg_hp$group,ref='ctrl.sca')
## formulas
bf_P<-bf(P~group+(1|id))
bf_S<-bf(S~group+(1|id))
bf_Fccp<-bf(Fccp~group+(1|id))
bf_Rot<-bf(Rot~group+(1|id))
bf_C_iv<-bf(C_iv~group+(1|id))
## prior setting
prior_1 <- get_prior(bf_P + bf_S + bf_Fccp + bf_Rot + bf_C_iv,
data = mg_hp,
family=Gamma(link="log"))Setting 'rescor' to FALSE by default for this modelprior_1 prior class coef group resp dpar nlpar lb ub
(flat) b
(flat) Intercept
(flat) b Civ
(flat) b groupctrl.wt Civ
(flat) b groupeda.sca Civ
(flat) b groupeda.wt Civ
student_t(3, 5, 2.5) Intercept Civ
student_t(3, 0, 2.5) sd Civ 0
student_t(3, 0, 2.5) sd id Civ 0
student_t(3, 0, 2.5) sd Intercept id Civ 0
gamma(0.01, 0.01) shape Civ 0
(flat) b Fccp
(flat) b groupctrl.wt Fccp
(flat) b groupeda.sca Fccp
(flat) b groupeda.wt Fccp
student_t(3, 4.4, 2.5) Intercept Fccp
student_t(3, 0, 2.5) sd Fccp 0
student_t(3, 0, 2.5) sd id Fccp 0
student_t(3, 0, 2.5) sd Intercept id Fccp 0
gamma(0.01, 0.01) shape Fccp 0
(flat) b P
(flat) b groupctrl.wt P
(flat) b groupeda.sca P
(flat) b groupeda.wt P
student_t(3, 3.1, 2.5) Intercept P
student_t(3, 0, 2.5) sd P 0
student_t(3, 0, 2.5) sd id P 0
student_t(3, 0, 2.5) sd Intercept id P 0
gamma(0.01, 0.01) shape P 0
(flat) b Rot
(flat) b groupctrl.wt Rot
(flat) b groupeda.sca Rot
(flat) b groupeda.wt Rot
student_t(3, 3.4, 2.5) Intercept Rot
student_t(3, 0, 2.5) sd Rot 0
student_t(3, 0, 2.5) sd id Rot 0
student_t(3, 0, 2.5) sd Intercept id Rot 0
gamma(0.01, 0.01) shape Rot 0
(flat) b S
(flat) b groupctrl.wt S
(flat) b groupeda.sca S
(flat) b groupeda.wt S
student_t(3, 4, 2.5) Intercept S
student_t(3, 0, 2.5) sd S 0
student_t(3, 0, 2.5) sd id S 0
student_t(3, 0, 2.5) sd Intercept id S 0
gamma(0.01, 0.01) shape S 0
source
default
default
default
(vectorized)
(vectorized)
(vectorized)
default
default
(vectorized)
(vectorized)
default
default
(vectorized)
(vectorized)
(vectorized)
default
default
(vectorized)
(vectorized)
default
default
(vectorized)
(vectorized)
(vectorized)
default
default
(vectorized)
(vectorized)
default
default
(vectorized)
(vectorized)
(vectorized)
default
default
(vectorized)
(vectorized)
default
default
(vectorized)
(vectorized)
(vectorized)
default
default
(vectorized)
(vectorized)
defaultprior_1$prior[c(4,6,13,15,22,24,31,33,40,42)] <- "normal(0,2)"
prior_1$prior[c(5,14,23,32,41)] <- "normal(0,1.2)"3.2.2 Fitting the model
mg_hp_model <- run_model(
brm(bf_P+bf_S+bf_Fccp+bf_Rot+bf_C_iv,
data=mg_hp, prior = prior_1,
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/mg_hp_model', reuse = TRUE)
summary(mg_hp_model) Family: MV(gamma, gamma, gamma, gamma, gamma)
Links: mu = log; shape = identity
mu = log; shape = identity
mu = log; shape = identity
mu = log; shape = identity
mu = log; shape = identity
Formula: P ~ group + (1 | id)
S ~ group + (1 | id)
Fccp ~ group + (1 | id)
Rot ~ group + (1 | id)
C_iv ~ group + (1 | id)
Data: mg_hp (Number of observations: 84)
Draws: 2 chains, each with iter = 8000; warmup = 2000; thin = 1;
total post-warmup draws = 12000
Group-Level Effects:
~id (Number of levels: 42)
Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
sd(P_Intercept) 0.40 0.05 0.31 0.52 1.00 2579 4599
sd(S_Intercept) 0.22 0.03 0.16 0.29 1.00 3106 5276
sd(Fccp_Intercept) 0.20 0.03 0.14 0.27 1.00 2885 4803
sd(Rot_Intercept) 0.17 0.03 0.12 0.23 1.00 3394 5334
sd(Civ_Intercept) 0.23 0.03 0.17 0.30 1.00 3309 5642
Population-Level Effects:
Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
P_Intercept 3.04 0.14 2.76 3.31 1.00 1682 3226
S_Intercept 3.95 0.08 3.79 4.11 1.00 2463 4866
Fccp_Intercept 4.36 0.07 4.22 4.51 1.00 2147 4278
Rot_Intercept 3.36 0.06 3.24 3.48 1.00 2395 4718
Civ_Intercept 4.94 0.08 4.77 5.10 1.00 2039 4149
P_groupctrl.wt 0.02 0.19 -0.36 0.39 1.00 1711 3089
P_groupeda.wt -0.03 0.19 -0.41 0.34 1.00 1727 3172
P_groupeda.sca 0.20 0.19 -0.18 0.59 1.00 1529 2976
S_groupctrl.wt 0.02 0.11 -0.20 0.23 1.00 2304 4444
S_groupeda.wt -0.03 0.11 -0.24 0.18 1.00 2542 4994
S_groupeda.sca 0.10 0.11 -0.11 0.33 1.00 2375 3959
Fccp_groupctrl.wt 0.01 0.10 -0.19 0.21 1.00 2299 4321
Fccp_groupeda.wt -0.01 0.10 -0.20 0.18 1.00 2307 4466
Fccp_groupeda.sca 0.14 0.10 -0.06 0.35 1.00 2188 4496
Rot_groupctrl.wt -0.04 0.09 -0.20 0.13 1.00 2655 4810
Rot_groupeda.wt -0.04 0.08 -0.20 0.12 1.00 2852 4423
Rot_groupeda.sca 0.08 0.09 -0.09 0.25 1.00 2694 4363
Civ_groupctrl.wt 0.25 0.11 0.03 0.47 1.00 2124 4123
Civ_groupeda.wt 0.11 0.11 -0.11 0.33 1.00 2207 4439
Civ_groupeda.sca 0.18 0.12 -0.05 0.40 1.00 2189 4146
Family Specific Parameters:
Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
shape_P 32.19 7.04 19.93 47.55 1.00 5060 7330
shape_S 46.98 10.35 28.66 69.44 1.00 4828 6254
shape_Fccp 51.46 11.36 31.74 76.17 1.00 4322 6382
shape_Rot 61.32 13.57 37.51 90.22 1.00 4132 5926
shape_Civ 45.46 9.77 27.90 66.34 1.00 4553 5892
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).Rhat values and ESS are fine, model converged well.
3.2.3 PPC
pp_check(mg_hp_model,type='dens_overlay',ndraws = 50,resp="P")
pp_check(mg_hp_model,type='dens_overlay_grouped',ndraws = 50,group='group',resp="P")
pp_check(mg_hp_model,type="stat_2d", stat = c("max", "min"),ndraws=100,resp="P")
pp_check(mg_hp_model,type="stat_2d", stat = c("mean", "sd"),ndraws=200,resp="P")
pp_check(mg_hp_model,type='dens_overlay',ndraws = 50,resp="S")
pp_check(mg_hp_model,type='dens_overlay_grouped',ndraws = 50,group='group',resp="S")
pp_check(mg_hp_model,type="stat_2d", stat = c("max", "min"),ndraws=100,resp="S")
pp_check(mg_hp_model,type="stat_2d", stat = c("mean", "sd"),ndraws=200,resp="S")
pp_check(mg_hp_model,type='dens_overlay',ndraws = 50,resp="Fccp")
pp_check(mg_hp_model,type='dens_overlay_grouped',ndraws = 50,group='group',resp="Fccp")
pp_check(mg_hp_model,type="stat_2d", stat = c("max", "min"),ndraws=100,resp="Fccp")
pp_check(mg_hp_model,type="stat_2d", stat = c("mean", "sd"),ndraws=200,resp="Fccp")
pp_check(mg_hp_model,type='dens_overlay',ndraws = 50,resp="Rot")
pp_check(mg_hp_model,type='dens_overlay_grouped',ndraws = 50,group='group',resp="Rot")
pp_check(mg_hp_model,type="stat_2d", stat = c("max", "min"),ndraws=100,resp="Rot")
pp_check(mg_hp_model,type="stat_2d", stat = c("mean", "sd"),ndraws=200,resp="Rot")
pp_check(mg_hp_model,type='dens_overlay',ndraws = 50,resp="Civ")
pp_check(mg_hp_model,type='dens_overlay_grouped',ndraws = 50,group='group',resp="Civ")
pp_check(mg_hp_model,type="stat_2d", stat = c("max", "min"),ndraws=100,resp="Civ")
pp_check(mg_hp_model,type="stat_2d", stat = c("mean", "sd"),ndraws=200,resp="Civ")
PPC looks fine, no serious issue detected.
3.2.4 Extraction of posterior samples
post_fix_mg_hp_P <- as.data.frame(mg_hp_model, variable = c("b_P_Intercept","b_P_groupctrl.wt",
"b_P_groupeda.wt","b_P_groupeda.sca"))
names(post_fix_mg_hp_P)[1] <- "sca_ctrl"
post_fix_mg_hp_P$wt_ctrl <- post_fix_mg_hp_P$sca_ctrl+post_fix_mg_hp_P$b_P_groupctrl.wt
post_fix_mg_hp_P$wt_eda <- post_fix_mg_hp_P$sca_ctrl+post_fix_mg_hp_P$b_P_groupeda.wt
post_fix_mg_hp_P$sca_eda <- post_fix_mg_hp_P$sca_ctrl+post_fix_mg_hp_P$b_P_groupeda.sca
post_fix_mg_hp_P$wt_contrast <- post_fix_mg_hp_P$wt_eda-post_fix_mg_hp_P$wt_ctrl
post_fix_mg_hp_S <- as.data.frame(mg_hp_model, variable = c("b_S_Intercept","b_S_groupctrl.wt",
"b_S_groupeda.wt","b_S_groupeda.sca"))
names(post_fix_mg_hp_S)[1] <- "sca_ctrl"
post_fix_mg_hp_S$wt_ctrl <- post_fix_mg_hp_S$sca_ctrl+post_fix_mg_hp_S$b_S_groupctrl.wt
post_fix_mg_hp_S$wt_eda <- post_fix_mg_hp_S$sca_ctrl+post_fix_mg_hp_S$b_S_groupeda.wt
post_fix_mg_hp_S$sca_eda <- post_fix_mg_hp_S$sca_ctrl+post_fix_mg_hp_S$b_S_groupeda.sca
post_fix_mg_hp_S$wt_contrast <- post_fix_mg_hp_S$wt_eda-post_fix_mg_hp_S$wt_ctrl
post_fix_mg_hp_Fccp <- as.data.frame(mg_hp_model, variable =
c("b_Fccp_Intercept","b_Fccp_groupctrl.wt",
"b_Fccp_groupeda.wt","b_Fccp_groupeda.sca"))
names(post_fix_mg_hp_Fccp)[1] <- "sca_ctrl"
post_fix_mg_hp_Fccp$wt_ctrl <- post_fix_mg_hp_Fccp$sca_ctrl+post_fix_mg_hp_Fccp$b_Fccp_groupctrl.wt
post_fix_mg_hp_Fccp$wt_eda <- post_fix_mg_hp_Fccp$sca_ctrl+post_fix_mg_hp_Fccp$b_Fccp_groupeda.wt
post_fix_mg_hp_Fccp$sca_eda <- post_fix_mg_hp_Fccp$sca_ctrl+post_fix_mg_hp_Fccp$b_Fccp_groupeda.sca
post_fix_mg_hp_Fccp$wt_contrast <- post_fix_mg_hp_Fccp$wt_eda-post_fix_mg_hp_Fccp$wt_ctrl
post_fix_mg_hp_Rot <- as.data.frame(mg_hp_model, variable = c("b_Rot_Intercept","b_Rot_groupctrl.wt",
"b_Rot_groupeda.wt","b_Rot_groupeda.sca"))
names(post_fix_mg_hp_Rot)[1] <- "sca_ctrl"
post_fix_mg_hp_Rot$wt_ctrl <- post_fix_mg_hp_Rot$sca_ctrl+post_fix_mg_hp_Rot$b_Rot_groupctrl.wt
post_fix_mg_hp_Rot$wt_eda <- post_fix_mg_hp_Rot$sca_ctrl+post_fix_mg_hp_Rot$b_Rot_groupeda.wt
post_fix_mg_hp_Rot$sca_eda <- post_fix_mg_hp_Rot$sca_ctrl+post_fix_mg_hp_Rot$b_Rot_groupeda.sca
post_fix_mg_hp_Rot$wt_contrast <- post_fix_mg_hp_Rot$wt_eda-post_fix_mg_hp_Rot$wt_ctrl
post_fix_mg_hp_Civ <- as.data.frame(mg_hp_model, variable = c("b_Civ_Intercept","b_Civ_groupctrl.wt",
"b_Civ_groupeda.wt","b_Civ_groupeda.sca"))
names(post_fix_mg_hp_Civ)[1] <- "sca_ctrl"
post_fix_mg_hp_Civ$wt_ctrl <- post_fix_mg_hp_Civ$sca_ctrl+post_fix_mg_hp_Civ$b_Civ_groupctrl.wt
post_fix_mg_hp_Civ$wt_eda <- post_fix_mg_hp_Civ$sca_ctrl+post_fix_mg_hp_Civ$b_Civ_groupeda.wt
post_fix_mg_hp_Civ$sca_eda <- post_fix_mg_hp_Civ$sca_ctrl+post_fix_mg_hp_Civ$b_Civ_groupeda.sca
post_fix_mg_hp_Civ$wt_contrast <- post_fix_mg_hp_Civ$wt_eda-post_fix_mg_hp_Civ$wt_ctrl3.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( {
### Plotting P ----------------
groupquant<-sapply(post_fix_mg_hp_P[,c(5,6,1,7)],function(p) quantile(p, probs = c(0.025,0.975,0.5)))
groupquant<-exp(groupquant)
mg_hp$group<-factor(mg_hp$group,levels=c('ctrl.wt','eda.wt','ctrl.sca','eda.sca'))
par(mfrow=c(2,5))
par(mar=c(2,3,1,0))
par(mgp=c(1.9,0.6,0))
range<-c(0,100);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="Respiration (pmol/s/mg) in hp",las=1, axes=FALSE,cex.lab=1.2)
rect(xrange[1],range[2],xrange[2],range[1],col="grey90",border=NA)
x<-range[1]
repeat{
lines(c(xrange[1],xrange[2]),c(x,x),col="white",lwd=0.7)
x=x+20;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(mg_hp$P~mg_hp$group,col=colb,at=c(1:4),add=T,
border=cola,axes=FALSE,drawRect=F,lwd=1,wex=0.8)
mg_hp_sum <- mg_hp %>% group_by(id, group, sex_num) %>% summarize(P = mean(P))
beeswarm(mg_hp_sum$P~mg_hp_sum$group,col=colc,at=c(1:4),
add=T,pwpch=(16.5+mg_hp_sum$sex_num),
corral = "wrap")
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<-140;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.25*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.25*xscal,ypos-0.08*scal,"Males",col="grey40",cex=1.2)
text(mean(xrange),range[2]+0.05*scal,"P I",cex=1.6,font=3,xpd=T)
### Plotting S ----------------
groupquant<-sapply(post_fix_mg_hp_S[,c(5,6,1,7)],function(p) quantile(p, probs = c(0.025,0.975,0.5)))
groupquant<-exp(groupquant)
range<-c(0,100);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="",las=1, axes=FALSE,cex.lab=1.2)
rect(xrange[1],range[2],xrange[2],range[1],col="grey90",border=NA)
x<-range[1]
repeat{
lines(c(xrange[1],xrange[2]),c(x,x),col="white",lwd=0.7)
x=x+20;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(mg_hp$S~mg_hp$group,col=colb,at=c(1:4),add=T,
border=cola,axes=FALSE,drawRect=F,lwd=1,wex=0.8)
mg_hp_sum <- mg_hp %>% group_by(id, group, sex_num) %>% summarize(S = mean(S))
beeswarm(mg_hp_sum$S~mg_hp_sum$group,col=colc,at=c(1:4),
add=T,pwpch=(16.5+mg_hp_sum$sex_num),
corral = "wrap")
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"))
text(mean(xrange),range[2]+0.05*scal,"P I+II",cex=1.6,font=3,xpd=T)
### Plotting Fccp ----------------
groupquant<-sapply(post_fix_mg_hp_Fccp[,c(5,6,1,7)],function(p) quantile(p, probs = c(0.025,0.975,0.5)))
groupquant<-exp(groupquant)
range<-c(0,200);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="",las=1, axes=FALSE,cex.lab=1.2)
rect(xrange[1],range[2],xrange[2],range[1],col="grey90",border=NA)
x<-range[1]
repeat{
lines(c(xrange[1],xrange[2]),c(x,x),col="white",lwd=0.7)
x=x+20;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(mg_hp$Fccp~mg_hp$group,col=colb,at=c(1:4),add=T,
border=cola,axes=FALSE,drawRect=F,lwd=1,wex=0.8)
mg_hp_sum <- mg_hp %>% group_by(id, group, sex_num) %>% summarize(P = mean(Fccp))
beeswarm(mg_hp_sum$P~mg_hp_sum$group,col=colc,at=c(1:4),
add=T,pwpch=(16.5+mg_hp_sum$sex_num),
corral = "wrap")
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=40),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"))
text(mean(xrange),range[2]+0.05*scal,"E I+II",cex=1.6,font=3,xpd=T)
### Plotting Rot ----------------
groupquant<-sapply(post_fix_mg_hp_Rot[,c(5,6,1,7)],function(p) quantile(p, probs = c(0.025,0.975,0.5)))
groupquant<-exp(groupquant)
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="",las=1, axes=FALSE,cex.lab=1.2)
rect(xrange[1],range[2],xrange[2],range[1],col="grey90",border=NA)
x<-range[1]
repeat{
lines(c(xrange[1],xrange[2]),c(x,x),col="white",lwd=0.7)
x=x+20;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(mg_hp$Rot~mg_hp$group,col=colb,at=c(1:4),add=T,
border=cola,axes=FALSE,drawRect=F,lwd=1,wex=0.8)
mg_hp_sum <- mg_hp %>% group_by(id, group, sex_num) %>% summarize(P = mean(Rot))
beeswarm(mg_hp_sum$P~mg_hp_sum$group,col=colc,at=c(1:4),
add=T,pwpch=(16.5+mg_hp_sum$sex_num),
corral = "wrap")
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"))
text(mean(xrange),range[2]+0.05*scal,"E II",cex=1.6,font=3,xpd=T)
### Plotting Civ ----------------
groupquant<-sapply(post_fix_mg_hp_Civ[,c(5,6,1,7)],function(p) quantile(p, probs = c(0.025,0.975,0.5)))
groupquant<-exp(groupquant)
range<-c(0,400);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="",las=1, axes=FALSE,cex.lab=1.2)
rect(xrange[1],range[2],xrange[2],range[1],col="grey90",border=NA)
x<-range[1]
repeat{
lines(c(xrange[1],xrange[2]),c(x,x),col="white",lwd=0.7)
x=x+20;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(mg_hp$C_iv~mg_hp$group,col=colb,at=c(1:4),add=T,
border=cola,axes=FALSE,drawRect=F,lwd=1,wex=0.8)
mg_hp_sum <- mg_hp %>% group_by(id, group, sex_num) %>% summarize(P = mean(C_iv))
beeswarm(mg_hp_sum$P~mg_hp_sum$group,col=colc,at=c(1:4),
add=T,pwpch=(16.5+mg_hp_sum$sex_num),
corral = "wrap")
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=80),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"))
text(mean(xrange),range[2]+0.05*scal,"Complex IV",cex=1.6,font=3,xpd=T)
## probability distributions -------------
### PD - P ----------
par(mar=c(2,1,0,0))
dif<-data.frame(post_fix_mg_hp_P$wt_contrast,post_fix_mg_hp_P$b_P_groupeda.sca,
post_fix_mg_hp_P$b_P_groupctrl.wt)
dif<-exp(dif)
yla<-"Fold-difference in respiration"
tckk=-0.018
ste<-seq(0.2,2.2,by=0.2)
mons_poste(dif,1)
zpos=seq(0,1,1/3)
xx=1;ind=0.5;xpol<-0.4;srt=55
text(xpol,zpos[xx]+zpos[2]*ind,
"Wt: ed x 0 ",cex=1,srt=srt,xpd=T)
xx=xx+1
text(xpol,zpos[xx]+zpos[2]*ind,
"SCA1: ed x 0 ",cex=1,srt=srt,xpd=T)
xx=xx+1
text(xpol,zpos[xx]+zpos[2]*ind,
"0: Wt x SCA1",cex=1,srt=srt,xpd=T)
### PD - _S ----------
dif<-data.frame(post_fix_mg_hp_S$wt_contrast,post_fix_mg_hp_S$b_S_groupeda.sca,
post_fix_mg_hp_S$b_S_groupctrl.wt)
dif<-exp(dif)
yla<-"Fold-difference in respiration"
tckk=-0.018
ste<-seq(0.2,2.2,by=0.2)
mons_poste(dif,1)
### PD - _Fccp ----------
dif<-data.frame(post_fix_mg_hp_Fccp$wt_contrast,post_fix_mg_hp_Fccp$b_Fccp_groupeda.sca,
post_fix_mg_hp_Fccp$b_Fccp_groupctrl.wt)
dif<-exp(dif)
yla<-"Fold-difference in respiration"
tckk=-0.018
ste<-seq(0.2,2.2,by=0.2)
mons_poste(dif,1)
### PD - _Rot ----------
dif<-data.frame(post_fix_mg_hp_Rot$wt_contrast,post_fix_mg_hp_Rot$b_Rot_groupeda.sca,
post_fix_mg_hp_Rot$b_Rot_groupctrl.wt)
dif<-exp(dif)
yla<-"Fold-difference in respiration"
tckk=-0.018
ste<-seq(0.2,1.7,by=0.2)
mons_poste(dif,1)
### PD - _Civ ----------
dif<-data.frame(post_fix_mg_hp_Civ$wt_contrast,post_fix_mg_hp_Civ$b_Civ_groupeda.sca,
post_fix_mg_hp_Civ$b_Civ_groupctrl.wt)
dif<-exp(dif)
yla<-"Fold-difference in respiration"
tckk=-0.018
ste<-seq(0.6,2,by=0.2)
mons_poste(dif,1)
} ) )
4 Cerebellum, weight standardization
4.1 Data exploration
par(mfrow=c(2,5))
par(mar=c(5,5,1,0))
par(mgp=c(3.6, 0.6,0))
## both sexes together
plot(mg_cb$P~mg_cb$group,
col=cola, cex.axis=0.8, las=2)
plot(mg_cb$S~mg_cb$group,
col=cola, cex.axis=0.8, las=2)
plot(mg_cb$Fccp~mg_cb$group,
col=cola, cex.axis=0.8, las=2)
plot(mg_cb$Rot~mg_cb$group,
col=cola, cex.axis=0.8, las=2)
plot(mg_cb$C_iv~mg_cb$group,
col=cola, cex.axis=0.8, las=2)
## by both group and sex
plot(mg_cb$P~mg_cb$group_sex,
col=cola[c(1,1,2,2,3,3,4,4)],cex.axis=0.8, las=2)
plot(mg_cb$S~mg_cb$group_sex,
col=cola[c(1,1,2,2,3,3,4,4)],cex.axis=0.8, las=2)
plot(mg_cb$Fccp~mg_cb$group_sex,
col=cola[c(1,1,2,2,3,3,4,4)],cex.axis=0.8, las=2)
plot(mg_cb$Rot~mg_cb$group_sex,
col=cola[c(1,1,2,2,3,3,4,4)],cex.axis=0.8, las=2)
plot(mg_cb$C_iv~mg_cb$group_sex,
col=cola[c(1,1,2,2,3,3,4,4)],cex.axis=0.8, las=2)
There seems to be consistent sex effect for complex iv respiration, but not in other parameters. As multivariate model will be fitted, and because it affects only one or two parameters, sex will be ignored.
Spread seems to be higher in groups and parameters with larger median values, implying heteroscedasticity. Moreover, data distribution seems to be asymmetric, with most outliers being above median. From this reason, and also because the data are continuous and non-negative, we will employ gamma regression
4.2 Modelling
To model all 5 parameters (reflecting mitochondrial respiration capacity in difference phases of mitochondrial respiratory cycle) in a single model, we will fit multivariate regression.
As there are technical replicates (two measurements per animal and brain region), we will apply hierarchical model with random intercept.
4.2.1 Setting priors
For simplicity, we will use default prior for intercept
## re-leveling groups
mg_cb$group<-relevel(mg_cb$group,ref='ctrl.sca')
## formulas
bf_P<-bf(P~group+(1|id))
bf_S<-bf(S~group+(1|id))
bf_Fccp<-bf(Fccp~group+(1|id))
bf_Rot<-bf(Rot~group+(1|id))
bf_C_iv<-bf(C_iv~group+(1|id))
## prior setting
prior_1 <- get_prior(bf_P + bf_S + bf_Fccp + bf_Rot + bf_C_iv,
data = mg_cb,
family=Gamma(link="log"))Setting 'rescor' to FALSE by default for this modelprior_1 prior class coef group resp dpar nlpar lb ub
(flat) b
(flat) Intercept
(flat) b Civ
(flat) b groupctrl.wt Civ
(flat) b groupeda.sca Civ
(flat) b groupeda.wt Civ
student_t(3, 5.3, 2.5) Intercept Civ
student_t(3, 0, 2.5) sd Civ 0
student_t(3, 0, 2.5) sd id Civ 0
student_t(3, 0, 2.5) sd Intercept id Civ 0
gamma(0.01, 0.01) shape Civ 0
(flat) b Fccp
(flat) b groupctrl.wt Fccp
(flat) b groupeda.sca Fccp
(flat) b groupeda.wt Fccp
student_t(3, 4.6, 2.5) Intercept Fccp
student_t(3, 0, 2.5) sd Fccp 0
student_t(3, 0, 2.5) sd id Fccp 0
student_t(3, 0, 2.5) sd Intercept id Fccp 0
gamma(0.01, 0.01) shape Fccp 0
(flat) b P
(flat) b groupctrl.wt P
(flat) b groupeda.sca P
(flat) b groupeda.wt P
student_t(3, 3.5, 2.5) Intercept P
student_t(3, 0, 2.5) sd P 0
student_t(3, 0, 2.5) sd id P 0
student_t(3, 0, 2.5) sd Intercept id P 0
gamma(0.01, 0.01) shape P 0
(flat) b Rot
(flat) b groupctrl.wt Rot
(flat) b groupeda.sca Rot
(flat) b groupeda.wt Rot
student_t(3, 3.7, 2.5) Intercept Rot
student_t(3, 0, 2.5) sd Rot 0
student_t(3, 0, 2.5) sd id Rot 0
student_t(3, 0, 2.5) sd Intercept id Rot 0
gamma(0.01, 0.01) shape Rot 0
(flat) b S
(flat) b groupctrl.wt S
(flat) b groupeda.sca S
(flat) b groupeda.wt S
student_t(3, 4.3, 2.5) Intercept S
student_t(3, 0, 2.5) sd S 0
student_t(3, 0, 2.5) sd id S 0
student_t(3, 0, 2.5) sd Intercept id S 0
gamma(0.01, 0.01) shape S 0
source
default
default
default
(vectorized)
(vectorized)
(vectorized)
default
default
(vectorized)
(vectorized)
default
default
(vectorized)
(vectorized)
(vectorized)
default
default
(vectorized)
(vectorized)
default
default
(vectorized)
(vectorized)
(vectorized)
default
default
(vectorized)
(vectorized)
default
default
(vectorized)
(vectorized)
(vectorized)
default
default
(vectorized)
(vectorized)
default
default
(vectorized)
(vectorized)
(vectorized)
default
default
(vectorized)
(vectorized)
defaultprior_1$prior[c(4,6,13,15,22,24,31,33,40,42)] <- "normal(0,2)"
prior_1$prior[c(5,14,23,32,41)] <- "normal(0,1.2)"4.2.2 Fitting the model
mg_cb_model <- run_model(
brm(bf_P+bf_S+bf_Fccp+bf_Rot+bf_C_iv,
data=mg_cb, prior = prior_1,
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/mg_cb_model', reuse = TRUE)
summary(mg_cb_model) Family: MV(gamma, gamma, gamma, gamma, gamma)
Links: mu = log; shape = identity
mu = log; shape = identity
mu = log; shape = identity
mu = log; shape = identity
mu = log; shape = identity
Formula: P ~ group + (1 | id)
S ~ group + (1 | id)
Fccp ~ group + (1 | id)
Rot ~ group + (1 | id)
C_iv ~ group + (1 | id)
Data: mg_cb (Number of observations: 84)
Draws: 2 chains, each with iter = 8000; warmup = 2000; thin = 1;
total post-warmup draws = 12000
Group-Level Effects:
~id (Number of levels: 42)
Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
sd(P_Intercept) 0.34 0.05 0.26 0.44 1.00 3789 6580
sd(S_Intercept) 0.23 0.04 0.17 0.31 1.00 3181 5752
sd(Fccp_Intercept) 0.22 0.03 0.16 0.29 1.00 3707 6410
sd(Rot_Intercept) 0.16 0.03 0.10 0.23 1.00 2690 3780
sd(Civ_Intercept) 0.26 0.03 0.20 0.33 1.00 3231 5918
Population-Level Effects:
Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
P_Intercept 3.69 0.11 3.47 3.91 1.00 4497 5765
S_Intercept 4.40 0.08 4.24 4.57 1.00 6003 6696
Fccp_Intercept 4.63 0.08 4.48 4.78 1.00 5961 6373
Rot_Intercept 3.72 0.06 3.60 3.85 1.00 6576 7633
Civ_Intercept 5.31 0.09 5.13 5.48 1.00 4424 5408
P_groupctrl.wt -0.09 0.16 -0.41 0.23 1.00 4710 6090
P_groupeda.wt -0.15 0.16 -0.45 0.16 1.00 4747 6052
P_groupeda.sca -0.22 0.16 -0.53 0.10 1.00 4868 5867
S_groupctrl.wt -0.01 0.12 -0.24 0.23 1.00 6425 6900
S_groupeda.wt -0.12 0.11 -0.34 0.10 1.00 6457 6733
S_groupeda.sca -0.22 0.12 -0.45 0.01 1.00 6161 6669
Fccp_groupctrl.wt -0.02 0.11 -0.23 0.20 1.00 6239 6671
Fccp_groupeda.wt -0.11 0.10 -0.32 0.09 1.00 6215 7251
Fccp_groupeda.sca -0.15 0.11 -0.36 0.07 1.00 5939 6904
Rot_groupctrl.wt 0.01 0.09 -0.17 0.19 1.00 7476 6927
Rot_groupeda.wt -0.08 0.09 -0.26 0.09 1.00 6870 8005
Rot_groupeda.sca -0.10 0.09 -0.28 0.08 1.00 7263 7472
Civ_groupctrl.wt 0.23 0.12 -0.01 0.47 1.00 4488 5500
Civ_groupeda.wt 0.06 0.12 -0.17 0.29 1.00 4315 5235
Civ_groupeda.sca 0.02 0.12 -0.23 0.27 1.00 4366 5184
Family Specific Parameters:
Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
shape_P 34.02 7.42 20.93 50.14 1.00 5028 6615
shape_S 34.15 7.54 20.98 50.23 1.00 4249 6929
shape_Fccp 41.43 8.93 25.92 61.00 1.00 4495 7269
shape_Rot 40.25 8.96 24.56 59.53 1.00 4034 4906
shape_Civ 75.39 16.38 47.01 111.07 1.00 4858 7584
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).Rhat values and ESS are fine, model converged well.
4.2.3 PPC
pp_check(mg_cb_model,type='dens_overlay',ndraws = 50,resp="P")
pp_check(mg_cb_model,type='dens_overlay_grouped',ndraws = 50,group='group',resp="P")
pp_check(mg_cb_model,type="stat_2d", stat = c("max", "min"),ndraws=100,resp="P")
pp_check(mg_cb_model,type="stat_2d", stat = c("mean", "sd"),ndraws=200,resp="P")
pp_check(mg_cb_model,type='dens_overlay',ndraws = 50,resp="S")
pp_check(mg_cb_model,type='dens_overlay_grouped',ndraws = 50,group='group',resp="S")
pp_check(mg_cb_model,type="stat_2d", stat = c("max", "min"),ndraws=100,resp="S")
pp_check(mg_cb_model,type="stat_2d", stat = c("mean", "sd"),ndraws=200,resp="S")
pp_check(mg_cb_model,type='dens_overlay',ndraws = 50,resp="Fccp")
pp_check(mg_cb_model,type='dens_overlay_grouped',ndraws = 50,group='group',resp="Fccp")
pp_check(mg_cb_model,type="stat_2d", stat = c("max", "min"),ndraws=100,resp="Fccp")
pp_check(mg_cb_model,type="stat_2d", stat = c("mean", "sd"),ndraws=200,resp="Fccp")
pp_check(mg_cb_model,type='dens_overlay',ndraws = 50,resp="Rot")
pp_check(mg_cb_model,type='dens_overlay_grouped',ndraws = 50,group='group',resp="Rot")
pp_check(mg_cb_model,type="stat_2d", stat = c("max", "min"),ndraws=100,resp="Rot")
pp_check(mg_cb_model,type="stat_2d", stat = c("mean", "sd"),ndraws=200,resp="Rot")
pp_check(mg_cb_model,type='dens_overlay',ndraws = 50,resp="Civ")
pp_check(mg_cb_model,type='dens_overlay_grouped',ndraws = 50,group='group',resp="Civ")
pp_check(mg_cb_model,type="stat_2d", stat = c("max", "min"),ndraws=100,resp="Civ")
pp_check(mg_cb_model,type="stat_2d", stat = c("mean", "sd"),ndraws=200,resp="Civ")
PPC looks fine, no serious issue detected.
4.2.4 Extraction of posterior samples
post_fix_mg_cb_P <- as.data.frame(mg_cb_model, variable = c("b_P_Intercept","b_P_groupctrl.wt",
"b_P_groupeda.wt","b_P_groupeda.sca"))
names(post_fix_mg_cb_P)[1] <- "sca_ctrl"
post_fix_mg_cb_P$wt_ctrl <- post_fix_mg_cb_P$sca_ctrl+post_fix_mg_cb_P$b_P_groupctrl.wt
post_fix_mg_cb_P$wt_eda <- post_fix_mg_cb_P$sca_ctrl+post_fix_mg_cb_P$b_P_groupeda.wt
post_fix_mg_cb_P$sca_eda <- post_fix_mg_cb_P$sca_ctrl+post_fix_mg_cb_P$b_P_groupeda.sca
post_fix_mg_cb_P$wt_contrast <- post_fix_mg_cb_P$wt_eda-post_fix_mg_cb_P$wt_ctrl
post_fix_mg_cb_S <- as.data.frame(mg_cb_model, variable = c("b_S_Intercept","b_S_groupctrl.wt",
"b_S_groupeda.wt","b_S_groupeda.sca"))
names(post_fix_mg_cb_S)[1] <- "sca_ctrl"
post_fix_mg_cb_S$wt_ctrl <- post_fix_mg_cb_S$sca_ctrl+post_fix_mg_cb_S$b_S_groupctrl.wt
post_fix_mg_cb_S$wt_eda <- post_fix_mg_cb_S$sca_ctrl+post_fix_mg_cb_S$b_S_groupeda.wt
post_fix_mg_cb_S$sca_eda <- post_fix_mg_cb_S$sca_ctrl+post_fix_mg_cb_S$b_S_groupeda.sca
post_fix_mg_cb_S$wt_contrast <- post_fix_mg_cb_S$wt_eda-post_fix_mg_cb_S$wt_ctrl
post_fix_mg_cb_Fccp <- as.data.frame(mg_cb_model, variable =
c("b_Fccp_Intercept","b_Fccp_groupctrl.wt",
"b_Fccp_groupeda.wt","b_Fccp_groupeda.sca"))
names(post_fix_mg_cb_Fccp)[1] <- "sca_ctrl"
post_fix_mg_cb_Fccp$wt_ctrl <- post_fix_mg_cb_Fccp$sca_ctrl+post_fix_mg_cb_Fccp$b_Fccp_groupctrl.wt
post_fix_mg_cb_Fccp$wt_eda <- post_fix_mg_cb_Fccp$sca_ctrl+post_fix_mg_cb_Fccp$b_Fccp_groupeda.wt
post_fix_mg_cb_Fccp$sca_eda <- post_fix_mg_cb_Fccp$sca_ctrl+post_fix_mg_cb_Fccp$b_Fccp_groupeda.sca
post_fix_mg_cb_Fccp$wt_contrast <- post_fix_mg_cb_Fccp$wt_eda-post_fix_mg_cb_Fccp$wt_ctrl
post_fix_mg_cb_Rot <- as.data.frame(mg_cb_model, variable = c("b_Rot_Intercept","b_Rot_groupctrl.wt",
"b_Rot_groupeda.wt","b_Rot_groupeda.sca"))
names(post_fix_mg_cb_Rot)[1] <- "sca_ctrl"
post_fix_mg_cb_Rot$wt_ctrl <- post_fix_mg_cb_Rot$sca_ctrl+post_fix_mg_cb_Rot$b_Rot_groupctrl.wt
post_fix_mg_cb_Rot$wt_eda <- post_fix_mg_cb_Rot$sca_ctrl+post_fix_mg_cb_Rot$b_Rot_groupeda.wt
post_fix_mg_cb_Rot$sca_eda <- post_fix_mg_cb_Rot$sca_ctrl+post_fix_mg_cb_Rot$b_Rot_groupeda.sca
post_fix_mg_cb_Rot$wt_contrast <- post_fix_mg_cb_Rot$wt_eda-post_fix_mg_cb_Rot$wt_ctrl
post_fix_mg_cb_Civ <- as.data.frame(mg_cb_model, variable = c("b_Civ_Intercept","b_Civ_groupctrl.wt",
"b_Civ_groupeda.wt","b_Civ_groupeda.sca"))
names(post_fix_mg_cb_Civ)[1] <- "sca_ctrl"
post_fix_mg_cb_Civ$wt_ctrl <- post_fix_mg_cb_Civ$sca_ctrl+post_fix_mg_cb_Civ$b_Civ_groupctrl.wt
post_fix_mg_cb_Civ$wt_eda <- post_fix_mg_cb_Civ$sca_ctrl+post_fix_mg_cb_Civ$b_Civ_groupeda.wt
post_fix_mg_cb_Civ$sca_eda <- post_fix_mg_cb_Civ$sca_ctrl+post_fix_mg_cb_Civ$b_Civ_groupeda.sca
post_fix_mg_cb_Civ$wt_contrast <- post_fix_mg_cb_Civ$wt_eda-post_fix_mg_cb_Civ$wt_ctrl4.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( {
### Plotting P ----------------
groupquant<-sapply(post_fix_mg_cb_P[,c(5,6,1,7)],
function(p) quantile(p, probs = c(0.025,0.975,0.5)))
groupquant<-exp(groupquant)
mg_cb$group<-factor(mg_cb$group,levels=c('ctrl.wt','eda.wt','ctrl.sca','eda.sca'))
#mat <- matrix(c(1:6,7:12),nrow = )
par(mfrow=c(2,5))
par(mar=c(2,3,1,0))
par(mgp=c(1.9,0.6,0))
range<-c(0,120);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="Respiration (pmol/s/mg) in cb",las=1, axes=FALSE,cex.lab=1.2)
rect(xrange[1],range[2],xrange[2],range[1],col="grey90",border=NA)
x<-range[1]
repeat{
lines(c(xrange[1],xrange[2]),c(x,x),col="white",lwd=0.7)
x=x+20;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(mg_cb$P~mg_cb$group,col=colb,at=c(1:4),add=T,
border=cola,axes=FALSE,drawRect=F,lwd=1,wex=0.8)
mg_cb_sum <- mg_cb %>% group_by(id, group, sex_num) %>% summarize(P = mean(P))
beeswarm(mg_cb_sum$P~mg_cb_sum$group,col=colc,at=c(1:4),
add=T,pwpch=(16.5+mg_cb_sum$sex_num),
corral = "wrap")
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<- 110 ;xpos=1.1
points(xpos,ypos,pch=16,cex=1.4,col=rgb(0.4,0.4,0.4,alpha=0.6))
text(xpos+0.25*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.25*xscal,ypos-0.08*scal,"Males",col="grey40",cex=1.2)
text(mean(xrange),range[2]+0.05*scal,"P I",cex=1.6,font=3,xpd=T)
### Plotting S ----------------
groupquant<-sapply(post_fix_mg_cb_S[,c(5,6,1,7)],function(p) quantile(p, probs = c(0.025,0.975,0.5)))
groupquant<-exp(groupquant)
range<-c(0,200);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="",las=1, axes=FALSE,cex.lab=1.2)
rect(xrange[1],range[2],xrange[2],range[1],col="grey90",border=NA)
x<-range[1]
repeat{
lines(c(xrange[1],xrange[2]),c(x,x),col="white",lwd=0.7)
x=x+20;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(mg_cb$S~mg_cb$group,col=colb,at=c(1:4),add=T,
border=cola,axes=FALSE,drawRect=F,lwd=1,wex=0.8)
mg_cb_sum <- mg_cb %>% group_by(id, group, sex_num) %>% summarize(P = mean(S))
beeswarm(mg_cb_sum$P~mg_cb_sum$group,col=colc,at=c(1:4),
add=T,pwpch=(16.5+mg_cb_sum$sex_num),
corral = "wrap")
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=40),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"))
text(mean(xrange),range[2]+0.05*scal,"P I+II",cex=1.6,font=3,xpd=T)
### Plotting Fccp ----------------
groupquant<-sapply(post_fix_mg_cb_Fccp[,c(5,6,1,7)],function(p) quantile(p, probs = c(0.025,0.975,0.5)))
groupquant<-exp(groupquant)
range<-c(0,240);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="",las=1, axes=FALSE,cex.lab=1.2)
rect(xrange[1],range[2],xrange[2],range[1],col="grey90",border=NA)
x<-range[1]
repeat{
lines(c(xrange[1],xrange[2]),c(x,x),col="white",lwd=0.7)
x=x+20;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(mg_cb$Fccp~mg_cb$group,col=colb,at=c(1:4),add=T,
border=cola,axes=FALSE,drawRect=F,lwd=1,wex=0.8)
mg_cb_sum <- mg_cb %>% group_by(id, group, sex_num) %>% summarize(P = mean(Fccp))
beeswarm(mg_cb_sum$P~mg_cb_sum$group,col=colc,at=c(1:4),
add=T,pwpch=(16.5+mg_cb_sum$sex_num),
corral = "wrap")
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=40),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"))
text(mean(xrange),range[2]+0.05*scal,"E I+II",cex=1.6,font=3,xpd=T)
### Plotting Rot ----------------
groupquant<-sapply(post_fix_mg_cb_Rot[,c(5,6,1,7)],function(p) quantile(p, probs = c(0.025,0.975,0.5)))
groupquant<-exp(groupquant)
range<-c(0,80);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="",las=1, axes=FALSE,cex.lab=1.2)
rect(xrange[1],range[2],xrange[2],range[1],col="grey90",border=NA)
x<-range[1]
repeat{
lines(c(xrange[1],xrange[2]),c(x,x),col="white",lwd=0.7)
x=x+20;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(mg_cb$Rot~mg_cb$group,col=colb,at=c(1:4),add=T,
border=cola,axes=FALSE,drawRect=F,lwd=1,wex=0.8)
mg_cb_sum <- mg_cb %>% group_by(id, group, sex_num) %>% summarize(P = mean(Rot))
beeswarm(mg_cb_sum$P~mg_cb_sum$group,col=colc,at=c(1:4),
add=T,pwpch=(16.5+mg_cb_sum$sex_num),
corral = "wrap")
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"))
text(mean(xrange),range[2]+0.05*scal,"E II",cex=1.6,font=3,xpd=T)
### Plotting Civ ----------------
groupquant<-sapply(post_fix_mg_cb_Civ[,c(5,6,1,7)],function(p) quantile(p, probs = c(0.025,0.975,0.5)))
groupquant<-exp(groupquant)
range<-c(0,500);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="",las=1, axes=FALSE,cex.lab=1.2)
rect(xrange[1],range[2],xrange[2],range[1],col="grey90",border=NA)
x<-range[1]
repeat{
lines(c(xrange[1],xrange[2]),c(x,x),col="white",lwd=0.7)
x=x+20;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(mg_cb$C_iv~mg_cb$group,col=colb,at=c(1:4),add=T,
border=cola,axes=FALSE,drawRect=F,lwd=1,wex=0.8)
mg_cb_sum <- mg_cb %>% group_by(id, group, sex_num) %>% summarize(P = mean(C_iv))
beeswarm(mg_cb_sum$P~mg_cb_sum$group,col=colc,at=c(1:4),
add=T,pwpch=(16.5+mg_cb_sum$sex_num),
corral = "wrap")
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=100),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"))
text(mean(xrange),range[2]+0.05*scal,"Complex IV",cex=1.6,font=3,xpd=T)
## probability distributions -------------
### PD - P ----------
par(mar=c(2,1,0,0))
dif<-data.frame(post_fix_mg_cb_P$wt_contrast,post_fix_mg_cb_P$b_P_groupeda.sca,
post_fix_mg_cb_P$b_P_groupctrl.wt)
dif<-exp(dif)
yla<-"Fold-difference in respiration"
tckk=-0.018
ste<-seq(0.2,2.2,by=0.2)
mons_poste(dif,1)
zpos=seq(0,1,1/3)
xx=1;ind=0.5;xpol<-0.4;srt=55
text(xpol,zpos[xx]+zpos[2]*ind,
"Wt: ed x 0 ",cex=1,srt=srt,xpd=T)
xx=xx+1
text(xpol,zpos[xx]+zpos[2]*ind,
"SCA1: ed x 0 ",cex=1,srt=srt,xpd=T)
xx=xx+1
text(xpol,zpos[xx]+zpos[2]*ind,
"0: Wt x SCA1",cex=1,srt=srt,xpd=T)
### PD - _S ----------
dif<-data.frame(post_fix_mg_cb_S$wt_contrast,post_fix_mg_cb_S$b_S_groupeda.sca,
post_fix_mg_cb_S$b_S_groupctrl.wt)
dif<-exp(dif)
yla<-"Fold-difference in respiration"
tckk=-0.018
ste<-seq(0.2,2.2,by=0.2)
mons_poste(dif,1)
### PD - _Fccp ----------
dif<-data.frame(post_fix_mg_cb_Fccp$wt_contrast,post_fix_mg_cb_Fccp$b_Fccp_groupeda.sca,
post_fix_mg_cb_Fccp$b_Fccp_groupctrl.wt)
dif<-exp(dif)
yla<-"Fold-difference in respiration"
tckk=-0.018
ste<-seq(0.2,2.2,by=0.2)
mons_poste(dif,1)
### PD - _Rot ----------
dif<-data.frame(post_fix_mg_cb_Rot$wt_contrast,post_fix_mg_cb_Rot$b_Rot_groupeda.sca,
post_fix_mg_cb_Rot$b_Rot_groupctrl.wt)
dif<-exp(dif)
yla<-"Fold-difference in respiration"
tckk=-0.018
ste<-seq(0.2,1.7,by=0.2)
mons_poste(dif,1)
### PD - _Civ ----------
dif<-data.frame(post_fix_mg_cb_Civ$wt_contrast,post_fix_mg_cb_Civ$b_Civ_groupeda.sca,
post_fix_mg_cb_Civ$b_Civ_groupctrl.wt)
dif<-exp(dif)
yla<-"Fold-difference in respiration"
tckk=-0.018
ste<-seq(0.6,2,by=0.2)
mons_poste(dif,1)
} ) )
5 Hippocampus, citrate standardization
5.1 Data exploration
par(mfrow=c(2,5))
par(mar=c(5,5,1,0))
par(mgp=c(3.6, 0.6,0))
## both sexes together
plot(cs_hp$P~cs_hp$group,
col=cola, cex.axis=0.8, las=2)
plot(cs_hp$S~cs_hp$group,
col=cola, cex.axis=0.8, las=2)
plot(cs_hp$Fccp~cs_hp$group,
col=cola, cex.axis=0.8, las=2)
plot(cs_hp$Rot~cs_hp$group,
col=cola, cex.axis=0.8, las=2)
plot(cs_hp$C_iv~cs_hp$group,
col=cola, cex.axis=0.8, las=2)
## by both group and sex
plot(cs_hp$P~cs_hp$group_sex,
col=cola[c(1,1,2,2,3,3,4,4)],cex.axis=0.8, las=2)
plot(cs_hp$S~cs_hp$group_sex,
col=cola[c(1,1,2,2,3,3,4,4)],cex.axis=0.8, las=2)
plot(cs_hp$Fccp~cs_hp$group_sex,
col=cola[c(1,1,2,2,3,3,4,4)],cex.axis=0.8, las=2)
plot(cs_hp$Rot~cs_hp$group_sex,
col=cola[c(1,1,2,2,3,3,4,4)],cex.axis=0.8, las=2)
plot(cs_hp$C_iv~cs_hp$group_sex,
col=cola[c(1,1,2,2,3,3,4,4)],cex.axis=0.8, las=2)
No consistent sex effect. Data distribution seems to be more symmetrical than for weight-standardized respirtation, but we will keep gamma regression for consistency.
5.2 Modelling
5.2.1 Setting priors
Although we measured mitochondrial respiration in our previous publication (Tichanek et al., 2020), context was different in the sense that previous publication used group-housed mice not undergoing any test or manipulation. In contrast, these data are from mice which have been tested in behavioral tests and which were handled regularly.
For simplicity, we will use default prior for intercept
## re-leveling groups
cs_hp$group<-relevel(cs_hp$group,ref='ctrl.sca')
## formulas
bf_P<-bf(P~group+(1|id))
bf_S<-bf(S~group+(1|id))
bf_Fccp<-bf(Fccp~group+(1|id))
bf_Rot<-bf(Rot~group+(1|id))
bf_C_iv<-bf(C_iv~group+(1|id))
## prior setting
prior_1 <- get_prior(bf_P + bf_S + bf_Fccp + bf_Rot + bf_C_iv,
data = cs_hp,
family=Gamma(link="log"))Setting 'rescor' to FALSE by default for this modelprior_1 prior class coef group resp dpar nlpar lb ub
(flat) b
(flat) Intercept
(flat) b Civ
(flat) b groupctrl.wt Civ
(flat) b groupeda.sca Civ
(flat) b groupeda.wt Civ
student_t(3, 2.8, 2.5) Intercept Civ
student_t(3, 0, 2.5) sd Civ 0
student_t(3, 0, 2.5) sd id Civ 0
student_t(3, 0, 2.5) sd Intercept id Civ 0
gamma(0.01, 0.01) shape Civ 0
(flat) b Fccp
(flat) b groupctrl.wt Fccp
(flat) b groupeda.sca Fccp
(flat) b groupeda.wt Fccp
student_t(3, 2.2, 2.5) Intercept Fccp
student_t(3, 0, 2.5) sd Fccp 0
student_t(3, 0, 2.5) sd id Fccp 0
student_t(3, 0, 2.5) sd Intercept id Fccp 0
gamma(0.01, 0.01) shape Fccp 0
(flat) b P
(flat) b groupctrl.wt P
(flat) b groupeda.sca P
(flat) b groupeda.wt P
student_t(3, 0.9, 2.5) Intercept P
student_t(3, 0, 2.5) sd P 0
student_t(3, 0, 2.5) sd id P 0
student_t(3, 0, 2.5) sd Intercept id P 0
gamma(0.01, 0.01) shape P 0
(flat) b Rot
(flat) b groupctrl.wt Rot
(flat) b groupeda.sca Rot
(flat) b groupeda.wt Rot
student_t(3, 1.1, 2.5) Intercept Rot
student_t(3, 0, 2.5) sd Rot 0
student_t(3, 0, 2.5) sd id Rot 0
student_t(3, 0, 2.5) sd Intercept id Rot 0
gamma(0.01, 0.01) shape Rot 0
(flat) b S
(flat) b groupctrl.wt S
(flat) b groupeda.sca S
(flat) b groupeda.wt S
student_t(3, 1.8, 2.5) Intercept S
student_t(3, 0, 2.5) sd S 0
student_t(3, 0, 2.5) sd id S 0
student_t(3, 0, 2.5) sd Intercept id S 0
gamma(0.01, 0.01) shape S 0
source
default
default
default
(vectorized)
(vectorized)
(vectorized)
default
default
(vectorized)
(vectorized)
default
default
(vectorized)
(vectorized)
(vectorized)
default
default
(vectorized)
(vectorized)
default
default
(vectorized)
(vectorized)
(vectorized)
default
default
(vectorized)
(vectorized)
default
default
(vectorized)
(vectorized)
(vectorized)
default
default
(vectorized)
(vectorized)
default
default
(vectorized)
(vectorized)
(vectorized)
default
default
(vectorized)
(vectorized)
defaultprior_1$prior[c(4,6,13,15,22,24,31,33,40,42)] <- "normal(0,2)"
prior_1$prior[c(5,14,23,32,41)] <- "normal(0,1.2)"5.2.2 Fitting the model
cs_hp_model <- run_model(
brm(bf_P+bf_S+bf_Fccp+bf_Rot+bf_C_iv,
data=cs_hp, prior = prior_1,
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/cs_hp_model', reuse = TRUE)
summary(cs_hp_model) Family: MV(gamma, gamma, gamma, gamma, gamma)
Links: mu = log; shape = identity
mu = log; shape = identity
mu = log; shape = identity
mu = log; shape = identity
mu = log; shape = identity
Formula: P ~ group + (1 | id)
S ~ group + (1 | id)
Fccp ~ group + (1 | id)
Rot ~ group + (1 | id)
C_iv ~ group + (1 | id)
Data: cs_hp (Number of observations: 77)
Draws: 2 chains, each with iter = 8000; warmup = 2000; thin = 1;
total post-warmup draws = 12000
Group-Level Effects:
~id (Number of levels: 39)
Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
sd(P_Intercept) 0.32 0.04 0.24 0.41 1.00 3126 6335
sd(S_Intercept) 0.13 0.03 0.08 0.18 1.00 2820 2524
sd(Fccp_Intercept) 0.10 0.03 0.05 0.16 1.00 1917 1187
sd(Rot_Intercept) 0.13 0.02 0.08 0.18 1.00 3291 5305
sd(Civ_Intercept) 0.21 0.03 0.15 0.28 1.00 3299 5838
Population-Level Effects:
Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
P_Intercept 0.83 0.11 0.61 1.04 1.00 2276 3515
S_Intercept 1.75 0.05 1.64 1.85 1.00 3512 5564
Fccp_Intercept 2.16 0.05 2.07 2.25 1.00 4017 6242
Rot_Intercept 1.16 0.05 1.06 1.26 1.00 3100 5165
Civ_Intercept 2.73 0.08 2.58 2.88 1.00 2640 4994
P_groupctrl.wt 0.05 0.15 -0.25 0.34 1.00 2308 3907
P_groupeda.wt -0.06 0.16 -0.37 0.25 1.00 2438 4320
P_groupeda.sca 0.11 0.15 -0.18 0.41 1.00 2252 3784
S_groupctrl.wt 0.04 0.07 -0.10 0.18 1.00 4114 6209
S_groupeda.wt 0.00 0.08 -0.15 0.15 1.00 3871 5574
S_groupeda.sca 0.02 0.07 -0.12 0.16 1.00 3948 5530
Fccp_groupctrl.wt 0.04 0.06 -0.09 0.16 1.00 4374 6575
Fccp_groupeda.wt 0.04 0.07 -0.09 0.17 1.00 4311 6400
Fccp_groupeda.sca 0.05 0.06 -0.07 0.18 1.00 4181 6988
Rot_groupctrl.wt -0.01 0.07 -0.14 0.13 1.00 3482 5450
Rot_groupeda.wt 0.02 0.07 -0.12 0.16 1.00 3458 5813
Rot_groupeda.sca -0.00 0.07 -0.14 0.13 1.00 3303 4998
Civ_groupctrl.wt 0.28 0.10 0.08 0.49 1.00 2883 4585
Civ_groupeda.wt 0.24 0.11 0.02 0.46 1.00 2951 5472
Civ_groupeda.sca 0.09 0.11 -0.11 0.31 1.00 2992 4797
Family Specific Parameters:
Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
shape_P 42.09 9.52 25.51 62.76 1.00 6103 7811
shape_S 63.08 14.62 37.34 95.36 1.00 3746 3571
shape_Fccp 71.62 16.97 42.39 107.85 1.00 2766 2335
shape_Rot 83.22 19.42 49.49 125.06 1.00 4263 6442
shape_Civ 48.45 10.96 29.35 72.53 1.00 4617 6040
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).Rhat values and ESS are fine, model converged well.
5.2.3 PPC
pp_check(cs_hp_model,type='dens_overlay',ndraws = 50,resp="P")
pp_check(cs_hp_model,type='dens_overlay_grouped',ndraws = 50,group='group',resp="P")
pp_check(cs_hp_model,type="stat_2d", stat = c("max", "min"),ndraws=100,resp="P")
pp_check(cs_hp_model,type="stat_2d", stat = c("mean", "sd"),ndraws=200,resp="P")
pp_check(cs_hp_model,type='dens_overlay',ndraws = 50,resp="S")
pp_check(cs_hp_model,type='dens_overlay_grouped',ndraws = 50,group='group',resp="S")
pp_check(cs_hp_model,type="stat_2d", stat = c("max", "min"),ndraws=100,resp="S")
pp_check(cs_hp_model,type="stat_2d", stat = c("mean", "sd"),ndraws=200,resp="S")
pp_check(cs_hp_model,type='dens_overlay',ndraws = 50,resp="Fccp")
pp_check(cs_hp_model,type='dens_overlay_grouped',ndraws = 50,group='group',resp="Fccp")
pp_check(cs_hp_model,type="stat_2d", stat = c("max", "min"),ndraws=100,resp="Fccp")
pp_check(cs_hp_model,type="stat_2d", stat = c("mean", "sd"),ndraws=200,resp="Fccp")
pp_check(cs_hp_model,type='dens_overlay',ndraws = 50,resp="Rot")
pp_check(cs_hp_model,type='dens_overlay_grouped',ndraws = 50,group='group',resp="Rot")
pp_check(cs_hp_model,type="stat_2d", stat = c("max", "min"),ndraws=100,resp="Rot")
pp_check(cs_hp_model,type="stat_2d", stat = c("mean", "sd"),ndraws=200,resp="Rot")
pp_check(cs_hp_model,type='dens_overlay',ndraws = 50,resp="Civ")
pp_check(cs_hp_model,type='dens_overlay_grouped',ndraws = 50,group='group',resp="Civ")
pp_check(cs_hp_model,type="stat_2d", stat = c("max", "min"),ndraws=100,resp="Civ")
pp_check(cs_hp_model,type="stat_2d", stat = c("mean", "sd"),ndraws=200,resp="Civ")
PPC looks generally fine, no serious issue detected. As in previous parts, minimal values are smaller than is predicted, suggesting below-outliers occurrence (technical issues with some replicates or wrong isolation of grey tissue from white nerves??)
5.2.4 Extraction of posterior samples
post_fix_cs_hp_P <- as.data.frame(
cs_hp_model, variable = c("b_P_Intercept","b_P_groupctrl.wt",
"b_P_groupeda.wt","b_P_groupeda.sca"))
names(post_fix_cs_hp_P)[1] <- "sca_ctrl"
post_fix_cs_hp_P$wt_ctrl <- post_fix_cs_hp_P$sca_ctrl+post_fix_cs_hp_P$b_P_groupctrl.wt
post_fix_cs_hp_P$wt_eda <- post_fix_cs_hp_P$sca_ctrl+post_fix_cs_hp_P$b_P_groupeda.wt
post_fix_cs_hp_P$sca_eda <- post_fix_cs_hp_P$sca_ctrl+post_fix_cs_hp_P$b_P_groupeda.sca
post_fix_cs_hp_P$wt_contrast <- post_fix_cs_hp_P$wt_eda-post_fix_cs_hp_P$wt_ctrl
post_fix_cs_hp_S <- as.data.frame(cs_hp_model, variable =
c("b_S_Intercept","b_S_groupctrl.wt",
"b_S_groupeda.wt","b_S_groupeda.sca"))
names(post_fix_cs_hp_S)[1] <- "sca_ctrl"
post_fix_cs_hp_S$wt_ctrl <- post_fix_cs_hp_S$sca_ctrl+post_fix_cs_hp_S$b_S_groupctrl.wt
post_fix_cs_hp_S$wt_eda <- post_fix_cs_hp_S$sca_ctrl+post_fix_cs_hp_S$b_S_groupeda.wt
post_fix_cs_hp_S$sca_eda <- post_fix_cs_hp_S$sca_ctrl+post_fix_cs_hp_S$b_S_groupeda.sca
post_fix_cs_hp_S$wt_contrast <- post_fix_cs_hp_S$wt_eda-post_fix_cs_hp_S$wt_ctrl
post_fix_cs_hp_Fccp <- as.data.frame(cs_hp_model, variable =
c("b_Fccp_Intercept","b_Fccp_groupctrl.wt",
"b_Fccp_groupeda.wt","b_Fccp_groupeda.sca"))
names(post_fix_cs_hp_Fccp)[1] <- "sca_ctrl"
post_fix_cs_hp_Fccp$wt_ctrl <- post_fix_cs_hp_Fccp$sca_ctrl+post_fix_cs_hp_Fccp$b_Fccp_groupctrl.wt
post_fix_cs_hp_Fccp$wt_eda <- post_fix_cs_hp_Fccp$sca_ctrl+post_fix_cs_hp_Fccp$b_Fccp_groupeda.wt
post_fix_cs_hp_Fccp$sca_eda <- post_fix_cs_hp_Fccp$sca_ctrl+post_fix_cs_hp_Fccp$b_Fccp_groupeda.sca
post_fix_cs_hp_Fccp$wt_contrast <- post_fix_cs_hp_Fccp$wt_eda-post_fix_cs_hp_Fccp$wt_ctrl
post_fix_cs_hp_Rot <- as.data.frame(cs_hp_model,
variable = c("b_Rot_Intercept","b_Rot_groupctrl.wt",
"b_Rot_groupeda.wt","b_Rot_groupeda.sca"))
names(post_fix_cs_hp_Rot)[1] <- "sca_ctrl"
post_fix_cs_hp_Rot$wt_ctrl <- post_fix_cs_hp_Rot$sca_ctrl+
post_fix_cs_hp_Rot$b_Rot_groupctrl.wt
post_fix_cs_hp_Rot$wt_eda <- post_fix_cs_hp_Rot$sca_ctrl+
post_fix_cs_hp_Rot$b_Rot_groupeda.wt
post_fix_cs_hp_Rot$sca_eda <- post_fix_cs_hp_Rot$sca_ctrl + post_fix_cs_hp_Rot$b_Rot_groupeda.sca
post_fix_cs_hp_Rot$wt_contrast <- post_fix_cs_hp_Rot$wt_eda-post_fix_cs_hp_Rot$wt_ctrl
post_fix_cs_hp_Civ <- as.data.frame(cs_hp_model,
variable = c("b_Civ_Intercept","b_Civ_groupctrl.wt",
"b_Civ_groupeda.wt","b_Civ_groupeda.sca"))
names(post_fix_cs_hp_Civ)[1] <- "sca_ctrl"
post_fix_cs_hp_Civ$wt_ctrl <- post_fix_cs_hp_Civ$sca_ctrl+post_fix_cs_hp_Civ$b_Civ_groupctrl.wt
post_fix_cs_hp_Civ$wt_eda <- post_fix_cs_hp_Civ$sca_ctrl+post_fix_cs_hp_Civ$b_Civ_groupeda.wt
post_fix_cs_hp_Civ$sca_eda <- post_fix_cs_hp_Civ$sca_ctrl +
post_fix_cs_hp_Civ$b_Civ_groupeda.sca
post_fix_cs_hp_Civ$wt_contrast <- post_fix_cs_hp_Civ$wt_eda-post_fix_cs_hp_Civ$wt_ctrl5.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( {
### Plotting P ----------------
groupquant<-sapply(post_fix_cs_hp_P[,c(5,6,1,7)],function(p) quantile(p, probs = c(0.025,0.975,0.5)))
groupquant<-exp(groupquant)
cs_hp$group<-factor(cs_hp$group,levels=c('ctrl.wt','eda.wt','ctrl.sca','eda.sca'))
par(mfrow=c(2,5))
par(mar=c(2,2.6,1,0))
par(mgp=c(1.5,0.6,0))
range<-c(0,6);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="Respiration per CS (umol/s/IU) in hp",las=1, axes=FALSE,cex.lab=1.2)
rect(xrange[1],range[2],xrange[2],range[1],col="grey90",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<-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(cs_hp$P~cs_hp$group,col=colb,at=c(1:4),add=T,
border=cola,axes=FALSE,drawRect=F,lwd=1,wex=0.8)
mg_cb_sum <- cs_hp %>% group_by(id, group, sex_num) %>% summarize(P = mean(P))
beeswarm(mg_cb_sum$P~mg_cb_sum$group,col=colc,at=c(1:4),
add=T,pwpch=(16.5+mg_cb_sum$sex_num),
corral = "wrap")
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=3),
labels=c(rep("",length(seq(range[1],range[2],by=3)))),pos=xrange[1],tck=tckk)
axis(2,las=2,cex.axis=1.1,at=seq(range[1],range[2],by=3),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<- 5.5 ;xpos=1.1
points(xpos,ypos,pch=16,cex=1.4,col=rgb(0.4,0.4,0.4,alpha=0.6))
text(xpos+0.25*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.25*xscal,ypos-0.08*scal,"Males",col="grey40",cex=1.2)
text(mean(xrange),range[2]+0.05*scal,"P I",cex=1.6,font=3,xpd=T)
### Plotting S ----------------
groupquant<-sapply(post_fix_cs_hp_S[,c(5,6,1,7)],function(p) quantile(p, probs = c(0.025,0.975,0.5)))
groupquant<-exp(groupquant)
cs_hp$group<-factor(cs_hp$group,levels=c('ctrl.wt','eda.wt','ctrl.sca','eda.sca'))
range<-c(0,9);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="",las=1, axes=FALSE,cex.lab=1.2)
rect(xrange[1],range[2],xrange[2],range[1],col="grey90",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<-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(cs_hp$S~cs_hp$group,col=colb,at=c(1:4),add=T,
border=cola,axes=FALSE,drawRect=F,lwd=1,wex=0.8)
mg_cb_sum <- cs_hp %>% group_by(id, group, sex_num) %>% summarize(P = mean(S))
beeswarm(mg_cb_sum$P~mg_cb_sum$group,col=colc,at=c(1:4),
add=T,pwpch=(16.5+mg_cb_sum$sex_num),
corral = "wrap")
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=3),
labels=c(rep("",length(seq(range[1],range[2],by=3)))),pos=xrange[1],tck=tckk)
axis(2,las=2,cex.axis=1.1,at=seq(range[1],range[2],by=3),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"))
text(mean(xrange),range[2]+0.05*scal,"P I+II",cex=1.6,font=3,xpd=T)
### Plotting Fccp ----------------
groupquant<-sapply(post_fix_cs_hp_Fccp[,c(5,6,1,7)],function(p) quantile(p, probs = c(0.025,0.975,0.5)))
groupquant<-exp(groupquant)
cs_hp$group<-factor(cs_hp$group,levels=c('ctrl.wt','eda.wt','ctrl.sca','eda.sca'))
range<-c(0,12);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="",las=1, axes=FALSE,cex.lab=1.2)
rect(xrange[1],range[2],xrange[2],range[1],col="grey90",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<-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(cs_hp$Fccp~cs_hp$group,col=colb,at=c(1:4),add=T,
border=cola,axes=FALSE,drawRect=F,lwd=1,wex=0.8)
mg_cb_sum <- cs_hp %>% group_by(id, group, sex_num) %>% summarize(P = mean(Fccp))
beeswarm(mg_cb_sum$P~mg_cb_sum$group,col=colc,at=c(1:4),
add=T,pwpch=(16.5+mg_cb_sum$sex_num),
corral = "wrap")
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=3),
labels=c(rep("",length(seq(range[1],range[2],by=3)))),pos=xrange[1],tck=tckk)
axis(2,las=2,cex.axis=1.1,at=seq(range[1],range[2],by=3),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"))
text(mean(xrange),range[2]+0.05*scal,"E I+II",cex=1.6,font=3,xpd=T)
### Plotting Rot ----------------
groupquant<-sapply(post_fix_cs_hp_Rot[,c(5,6,1,7)],function(p) quantile(p, probs = c(0.025,0.975,0.5)))
groupquant<-exp(groupquant)
cs_hp$group<-factor(cs_hp$group,levels=c('ctrl.wt','eda.wt','ctrl.sca','eda.sca'))
range<-c(0,6);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="",las=1, axes=FALSE,cex.lab=1.2)
rect(xrange[1],range[2],xrange[2],range[1],col="grey90",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<-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(cs_hp$Rot~cs_hp$group,col=colb,at=c(1:4),add=T,
border=cola,axes=FALSE,drawRect=F,lwd=1,wex=0.8)
mg_cb_sum <- cs_hp %>% group_by(id, group, sex_num) %>% summarize(P = mean(Rot))
beeswarm(mg_cb_sum$P~mg_cb_sum$group,col=colc,at=c(1:4),
add=T,pwpch=(16.5+mg_cb_sum$sex_num),
corral = "wrap")
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=3),
labels=c(rep("",length(seq(range[1],range[2],by=3)))),pos=xrange[1],tck=tckk)
axis(2,las=2,cex.axis=1.1,at=seq(range[1],range[2],by=3),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"))
text(mean(xrange),range[2]+0.05*scal,"E II",cex=1.6,font=3,xpd=T)
### Plotting Civ ----------------
groupquant<-sapply(post_fix_cs_hp_Civ[,c(5,6,1,7)],function(p) quantile(p, probs = c(0.025,0.975,0.5)))
groupquant<-exp(groupquant)
cs_hp$group<-factor(cs_hp$group,levels=c('ctrl.wt','eda.wt','ctrl.sca','eda.sca'))
range<-c(0,36);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="",las=1, axes=FALSE,cex.lab=1.2)
rect(xrange[1],range[2],xrange[2],range[1],col="grey90",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<-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(cs_hp$C_iv~cs_hp$group,col=colb,at=c(1:4),add=T,
border=cola,axes=FALSE,drawRect=F,lwd=1,wex=0.8)
mg_cb_sum <- cs_hp %>% group_by(id, group, sex_num) %>% summarize(P = mean(C_iv))
beeswarm(mg_cb_sum$P~mg_cb_sum$group,col=colc,at=c(1:4),
add=T,pwpch=(16.5+mg_cb_sum$sex_num),
corral = "wrap")
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=3),
labels=c(rep("",length(seq(range[1],range[2],by=3)))),pos=xrange[1],tck=tckk)
axis(2,las=2,cex.axis=1.1,at=seq(range[1],range[2],by=6),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"))
text(mean(xrange),range[2]+0.05*scal,"Complex IV",cex=1.6,font=3,xpd=T)
## probability distributions -------------
### PD - P ----------
par(mar=c(2,1,0,0))
dif<-data.frame(post_fix_cs_hp_P$wt_contrast,post_fix_cs_hp_P$b_P_groupeda.sca,
post_fix_cs_hp_P$b_P_groupctrl.wt)
dif<-exp(dif)
yla<-"Fold-difference in respiration"
tckk=-0.018
ste<-seq(0.2,2.2,by=0.2)
mons_poste(dif,1)
zpos=seq(0,1,1/3)
xx=1;ind=0.5;xpol<-0.4;srt=55
text(xpol,zpos[xx]+zpos[2]*ind,
"Wt: ed x 0 ",cex=1,srt=srt,xpd=T)
xx=xx+1
text(xpol,zpos[xx]+zpos[2]*ind,
"SCA1: ed x 0 ",cex=1,srt=srt,xpd=T)
xx=xx+1
text(xpol,zpos[xx]+zpos[2]*ind,
"0: Wt x SCA1",cex=1,srt=srt,xpd=T)
### PD - _S ----------
dif<-data.frame(post_fix_cs_hp_S$wt_contrast,post_fix_cs_hp_S$b_S_groupeda.sca,
post_fix_cs_hp_S$b_S_groupctrl.wt)
dif<-exp(dif)
yla<-"Fold-difference in respiration"
tckk=-0.018
ste<-seq(0.2,2.2,by=0.2)
mons_poste(dif,1)
### PD - _Fccp ----------
dif<-data.frame(post_fix_cs_hp_Fccp$wt_contrast,post_fix_cs_hp_Fccp$b_Fccp_groupeda.sca,
post_fix_cs_hp_Fccp$b_Fccp_groupctrl.wt)
dif<-exp(dif)
yla<-"Fold-difference in respiration"
tckk=-0.018
ste<-seq(0.2,2.2,by=0.2)
mons_poste(dif,1)
### PD - _Rot ----------
dif<-data.frame(post_fix_cs_hp_Rot$wt_contrast,post_fix_cs_hp_Rot$b_Rot_groupeda.sca,
post_fix_cs_hp_Rot$b_Rot_groupctrl.wt)
dif<-exp(dif)
yla<-"Fold-difference in respiration"
tckk=-0.018
ste<-seq(0.2,2.2,by=0.2)
mons_poste(dif,1)
### PD - _Civ ----------
dif<-data.frame(post_fix_cs_hp_Civ$wt_contrast,post_fix_cs_hp_Civ$b_Civ_groupeda.sca,
post_fix_cs_hp_Civ$b_Civ_groupctrl.wt)
dif<-exp(dif)
yla<-"Fold-difference in respiration"
tckk=-0.018
ste<-seq(0.2,2.1,by=0.2)
mons_poste(dif,1)
} ) )
6 Cerebellum, citrate standardization
6.1 Data exploration
par(mfrow=c(2,5))
par(mar=c(5,5,1,0))
par(mgp=c(3.6, 0.6,0))
## both sexes together
plot(cs_cb$P~cs_cb$group,
col=cola, cex.axis=0.8, las=2)
plot(cs_cb$S~cs_cb$group,
col=cola, cex.axis=0.8, las=2)
plot(cs_cb$Fccp~cs_cb$group,
col=cola, cex.axis=0.8, las=2)
plot(cs_cb$Rot~cs_cb$group,
col=cola, cex.axis=0.8, las=2)
plot(cs_cb$C_iv~cs_cb$group,
col=cola, cex.axis=0.8, las=2)
## by both group and sex
plot(cs_cb$P~cs_cb$group_sex,
col=cola[c(1,1,2,2,3,3,4,4)],cex.axis=0.8, las=2)
plot(cs_cb$S~cs_cb$group_sex,
col=cola[c(1,1,2,2,3,3,4,4)],cex.axis=0.8, las=2)
plot(cs_cb$Fccp~cs_cb$group_sex,
col=cola[c(1,1,2,2,3,3,4,4)],cex.axis=0.8, las=2)
plot(cs_cb$Rot~cs_cb$group_sex,
col=cola[c(1,1,2,2,3,3,4,4)],cex.axis=0.8, las=2)
plot(cs_cb$C_iv~cs_cb$group_sex,
col=cola[c(1,1,2,2,3,3,4,4)],cex.axis=0.8, las=2)
No consistent sex effect. Data distribution seems again more asymmetrical compared to hp data. Gamma seems more appropriate and will be used. Sex will be ignored
6.2 Modelling
6.2.1 Setting priors
Although we measured mitochondrial respiration in our previous publication (Tichanek et al., 2020), context was different in the sense that previous publication used group-housed mice not undergoing any test or manipulation. In contrast, these data are from mice which have been tested in behavioral tests and which were handled regularly.
For simplicity, we will use default prior for intercept
## re-leveling groups
cs_cb$group<-relevel(cs_cb$group,ref='ctrl.sca')
## formulas
bf_P<-bf(P~group+(1|id))
bf_S<-bf(S~group+(1|id))
bf_Fccp<-bf(Fccp~group+(1|id))
bf_Rot<-bf(Rot~group+(1|id))
bf_C_iv<-bf(C_iv~group+(1|id))
## prior setting
prior_1 <- get_prior(bf_P + bf_S + bf_Fccp + bf_Rot + bf_C_iv,
data = cs_cb,
family=Gamma(link="log"))Setting 'rescor' to FALSE by default for this modelprior_1 prior class coef group resp dpar nlpar lb ub
(flat) b
(flat) Intercept
(flat) b Civ
(flat) b groupctrl.wt Civ
(flat) b groupeda.sca Civ
(flat) b groupeda.wt Civ
student_t(3, 3.2, 2.5) Intercept Civ
student_t(3, 0, 2.5) sd Civ 0
student_t(3, 0, 2.5) sd id Civ 0
student_t(3, 0, 2.5) sd Intercept id Civ 0
gamma(0.01, 0.01) shape Civ 0
(flat) b Fccp
(flat) b groupctrl.wt Fccp
(flat) b groupeda.sca Fccp
(flat) b groupeda.wt Fccp
student_t(3, 2.4, 2.5) Intercept Fccp
student_t(3, 0, 2.5) sd Fccp 0
student_t(3, 0, 2.5) sd id Fccp 0
student_t(3, 0, 2.5) sd Intercept id Fccp 0
gamma(0.01, 0.01) shape Fccp 0
(flat) b P
(flat) b groupctrl.wt P
(flat) b groupeda.sca P
(flat) b groupeda.wt P
student_t(3, 1.4, 2.5) Intercept P
student_t(3, 0, 2.5) sd P 0
student_t(3, 0, 2.5) sd id P 0
student_t(3, 0, 2.5) sd Intercept id P 0
gamma(0.01, 0.01) shape P 0
(flat) b Rot
(flat) b groupctrl.wt Rot
(flat) b groupeda.sca Rot
(flat) b groupeda.wt Rot
student_t(3, 1.5, 2.5) Intercept Rot
student_t(3, 0, 2.5) sd Rot 0
student_t(3, 0, 2.5) sd id Rot 0
student_t(3, 0, 2.5) sd Intercept id Rot 0
gamma(0.01, 0.01) shape Rot 0
(flat) b S
(flat) b groupctrl.wt S
(flat) b groupeda.sca S
(flat) b groupeda.wt S
student_t(3, 2.2, 2.5) Intercept S
student_t(3, 0, 2.5) sd S 0
student_t(3, 0, 2.5) sd id S 0
student_t(3, 0, 2.5) sd Intercept id S 0
gamma(0.01, 0.01) shape S 0
source
default
default
default
(vectorized)
(vectorized)
(vectorized)
default
default
(vectorized)
(vectorized)
default
default
(vectorized)
(vectorized)
(vectorized)
default
default
(vectorized)
(vectorized)
default
default
(vectorized)
(vectorized)
(vectorized)
default
default
(vectorized)
(vectorized)
default
default
(vectorized)
(vectorized)
(vectorized)
default
default
(vectorized)
(vectorized)
default
default
(vectorized)
(vectorized)
(vectorized)
default
default
(vectorized)
(vectorized)
defaultprior_1$prior[c(4,6,13,15,22,24,31,33,40,42)] <- "normal(0,2)"
prior_1$prior[c(5,14,23,32,41)] <- "normal(0,1.2)"6.2.2 Fitting the model
cs_cb_model <- run_model(
brm(bf_P+bf_S+bf_Fccp+bf_Rot+bf_C_iv,
data=cs_cb, prior = prior_1,
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/cs_cb_model', reuse = TRUE)
summary(cs_cb_model) Family: MV(gamma, gamma, gamma, gamma, gamma)
Links: mu = log; shape = identity
mu = log; shape = identity
mu = log; shape = identity
mu = log; shape = identity
mu = log; shape = identity
Formula: P ~ group + (1 | id)
S ~ group + (1 | id)
Fccp ~ group + (1 | id)
Rot ~ group + (1 | id)
C_iv ~ group + (1 | id)
Data: cs_cb (Number of observations: 76)
Draws: 2 chains, each with iter = 8000; warmup = 2000; thin = 1;
total post-warmup draws = 12000
Group-Level Effects:
~id (Number of levels: 39)
Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
sd(P_Intercept) 0.29 0.04 0.22 0.38 1.00 2866 4709
sd(S_Intercept) 0.19 0.03 0.13 0.26 1.00 3027 4546
sd(Fccp_Intercept) 0.17 0.03 0.12 0.24 1.00 3439 5244
sd(Rot_Intercept) 0.14 0.03 0.07 0.20 1.00 2389 2365
sd(Civ_Intercept) 0.21 0.04 0.15 0.29 1.00 3035 5107
Population-Level Effects:
Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
P_Intercept 1.48 0.10 1.28 1.67 1.00 2077 3585
S_Intercept 2.20 0.07 2.06 2.33 1.00 2829 4918
Fccp_Intercept 2.42 0.06 2.30 2.55 1.00 2816 4777
Rot_Intercept 1.52 0.05 1.42 1.63 1.00 3754 6137
Civ_Intercept 3.10 0.08 2.95 3.26 1.00 3008 4555
P_groupctrl.wt -0.07 0.14 -0.35 0.21 1.00 1899 3034
P_groupeda.wt -0.11 0.14 -0.39 0.17 1.00 2340 4148
P_groupeda.sca -0.16 0.14 -0.43 0.11 1.00 2023 4074
S_groupctrl.wt 0.00 0.10 -0.19 0.19 1.00 2863 4991
S_groupeda.wt -0.04 0.10 -0.24 0.15 1.00 3310 5526
S_groupeda.sca -0.17 0.10 -0.36 0.02 1.00 3022 4513
Fccp_groupctrl.wt -0.00 0.09 -0.19 0.17 1.00 3089 5244
Fccp_groupeda.wt -0.04 0.09 -0.22 0.14 1.00 3106 5539
Fccp_groupeda.sca -0.10 0.09 -0.28 0.08 1.00 2816 4656
Rot_groupctrl.wt 0.02 0.08 -0.13 0.17 1.00 3598 5792
Rot_groupeda.wt 0.01 0.08 -0.15 0.17 1.00 3932 5274
Rot_groupeda.sca -0.06 0.08 -0.22 0.09 1.00 3889 5483
Civ_groupctrl.wt 0.25 0.11 0.04 0.46 1.00 3199 5086
Civ_groupeda.wt 0.23 0.11 0.01 0.46 1.00 3077 5350
Civ_groupeda.sca 0.08 0.11 -0.13 0.30 1.00 3171 5246
Family Specific Parameters:
Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
shape_P 49.65 11.45 29.85 74.86 1.00 4978 6146
shape_S 48.53 11.41 28.65 73.15 1.00 3977 5528
shape_Fccp 57.93 13.76 34.15 88.11 1.00 4361 5882
shape_Rot 52.43 12.73 30.43 80.31 1.00 2715 3120
shape_Civ 38.55 8.95 23.11 58.09 1.00 4295 6066
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).Rhat values and ESS are fine, model converged well.
6.2.3 PPC
pp_check(cs_cb_model,type='dens_overlay',ndraws = 50,resp="P")
pp_check(cs_cb_model,type='dens_overlay_grouped',ndraws = 50,group='group',resp="P")
pp_check(cs_cb_model,type="stat_2d", stat = c("max", "min"),ndraws=100,resp="P")
pp_check(cs_cb_model,type="stat_2d", stat = c("mean", "sd"),ndraws=200,resp="P")
pp_check(cs_cb_model,type='dens_overlay',ndraws = 50,resp="S")
pp_check(cs_cb_model,type='dens_overlay_grouped',ndraws = 50,group='group',resp="S")
pp_check(cs_cb_model,type="stat_2d", stat = c("max", "min"),ndraws=100,resp="S")
pp_check(cs_cb_model,type="stat_2d", stat = c("mean", "sd"),ndraws=200,resp="S")
pp_check(cs_cb_model,type='dens_overlay',ndraws = 50,resp="Fccp")
pp_check(cs_cb_model,type='dens_overlay_grouped',ndraws = 50,group='group',resp="Fccp")
pp_check(cs_cb_model,type="stat_2d", stat = c("max", "min"),ndraws=100,resp="Fccp")
pp_check(cs_cb_model,type="stat_2d", stat = c("mean", "sd"),ndraws=200,resp="Fccp")
pp_check(cs_cb_model,type='dens_overlay',ndraws = 50,resp="Rot")
pp_check(cs_cb_model,type='dens_overlay_grouped',ndraws = 50,group='group',resp="Rot")
pp_check(cs_cb_model,type="stat_2d", stat = c("max", "min"),ndraws=100,resp="Rot")
pp_check(cs_cb_model,type="stat_2d", stat = c("mean", "sd"),ndraws=200,resp="Rot")
pp_check(cs_cb_model,type='dens_overlay',ndraws = 50,resp="Civ")
pp_check(cs_cb_model,type='dens_overlay_grouped',ndraws = 50,group='group',resp="Civ")
pp_check(cs_cb_model,type="stat_2d", stat = c("max", "min"),ndraws=100,resp="Civ")
pp_check(cs_cb_model,type="stat_2d", stat = c("mean", "sd"),ndraws=200,resp="Civ")
PPC looks generally fine, no serious issue detected. Occurrence of unexpectedly low minimal values is OK here.
6.2.4 Extraction of posterior samples
post_fix_cs_cb_P <- as.data.frame(
cs_cb_model, variable = c("b_P_Intercept","b_P_groupctrl.wt",
"b_P_groupeda.wt","b_P_groupeda.sca"))
names(post_fix_cs_cb_P)[1] <- "sca_ctrl"
post_fix_cs_cb_P$wt_ctrl <- post_fix_cs_cb_P$sca_ctrl+post_fix_cs_cb_P$b_P_groupctrl.wt
post_fix_cs_cb_P$wt_eda <- post_fix_cs_cb_P$sca_ctrl+post_fix_cs_cb_P$b_P_groupeda.wt
post_fix_cs_cb_P$sca_eda <- post_fix_cs_cb_P$sca_ctrl+post_fix_cs_cb_P$b_P_groupeda.sca
post_fix_cs_cb_P$wt_contrast <- post_fix_cs_cb_P$wt_eda-post_fix_cs_cb_P$wt_ctrl
post_fix_cs_cb_S <- as.data.frame(cs_cb_model, variable =
c("b_S_Intercept","b_S_groupctrl.wt",
"b_S_groupeda.wt","b_S_groupeda.sca"))
names(post_fix_cs_cb_S)[1] <- "sca_ctrl"
post_fix_cs_cb_S$wt_ctrl <- post_fix_cs_cb_S$sca_ctrl+post_fix_cs_cb_S$b_S_groupctrl.wt
post_fix_cs_cb_S$wt_eda <- post_fix_cs_cb_S$sca_ctrl+post_fix_cs_cb_S$b_S_groupeda.wt
post_fix_cs_cb_S$sca_eda <- post_fix_cs_cb_S$sca_ctrl+post_fix_cs_cb_S$b_S_groupeda.sca
post_fix_cs_cb_S$wt_contrast <- post_fix_cs_cb_S$wt_eda-post_fix_cs_cb_S$wt_ctrl
post_fix_cs_cb_Fccp <- as.data.frame(cs_cb_model, variable =
c("b_Fccp_Intercept","b_Fccp_groupctrl.wt",
"b_Fccp_groupeda.wt","b_Fccp_groupeda.sca"))
names(post_fix_cs_cb_Fccp)[1] <- "sca_ctrl"
post_fix_cs_cb_Fccp$wt_ctrl <- post_fix_cs_cb_Fccp$sca_ctrl+post_fix_cs_cb_Fccp$b_Fccp_groupctrl.wt
post_fix_cs_cb_Fccp$wt_eda <- post_fix_cs_cb_Fccp$sca_ctrl+post_fix_cs_cb_Fccp$b_Fccp_groupeda.wt
post_fix_cs_cb_Fccp$sca_eda <- post_fix_cs_cb_Fccp$sca_ctrl+post_fix_cs_cb_Fccp$b_Fccp_groupeda.sca
post_fix_cs_cb_Fccp$wt_contrast <- post_fix_cs_cb_Fccp$wt_eda-post_fix_cs_cb_Fccp$wt_ctrl
post_fix_cs_cb_Rot <- as.data.frame(cs_cb_model,
variable = c("b_Rot_Intercept","b_Rot_groupctrl.wt",
"b_Rot_groupeda.wt","b_Rot_groupeda.sca"))
names(post_fix_cs_cb_Rot)[1] <- "sca_ctrl"
post_fix_cs_cb_Rot$wt_ctrl <- post_fix_cs_cb_Rot$sca_ctrl+
post_fix_cs_cb_Rot$b_Rot_groupctrl.wt
post_fix_cs_cb_Rot$wt_eda <- post_fix_cs_cb_Rot$sca_ctrl+
post_fix_cs_cb_Rot$b_Rot_groupeda.wt
post_fix_cs_cb_Rot$sca_eda <- post_fix_cs_cb_Rot$sca_ctrl + post_fix_cs_cb_Rot$b_Rot_groupeda.sca
post_fix_cs_cb_Rot$wt_contrast <- post_fix_cs_cb_Rot$wt_eda-post_fix_cs_cb_Rot$wt_ctrl
post_fix_cs_cb_Civ <- as.data.frame(cs_cb_model,
variable = c("b_Civ_Intercept","b_Civ_groupctrl.wt",
"b_Civ_groupeda.wt","b_Civ_groupeda.sca"))
names(post_fix_cs_cb_Civ)[1] <- "sca_ctrl"
post_fix_cs_cb_Civ$wt_ctrl <- post_fix_cs_cb_Civ$sca_ctrl+post_fix_cs_cb_Civ$b_Civ_groupctrl.wt
post_fix_cs_cb_Civ$wt_eda <- post_fix_cs_cb_Civ$sca_ctrl+post_fix_cs_cb_Civ$b_Civ_groupeda.wt
post_fix_cs_cb_Civ$sca_eda <- post_fix_cs_cb_Civ$sca_ctrl +
post_fix_cs_cb_Civ$b_Civ_groupeda.sca
post_fix_cs_cb_Civ$wt_contrast <- post_fix_cs_cb_Civ$wt_eda-post_fix_cs_cb_Civ$wt_ctrl6.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( {
### Plotting P ----------------
groupquant<-sapply(post_fix_cs_cb_P[,c(5,6,1,7)],function(p) quantile(p, probs = c(0.025,0.975,0.5)))
groupquant<-exp(groupquant)
cs_cb$group<-factor(cs_cb$group,levels=c('ctrl.wt','eda.wt','ctrl.sca','eda.sca'))
par(mfrow=c(2,5))
par(mar=c(2,2.6,1,0))
par(mgp=c(1.5,0.6,0))
range<-c(0,12);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="Respiration per CS (umol/s/IU) in cb",las=1, axes=FALSE,cex.lab=1.2)
rect(xrange[1],range[2],xrange[2],range[1],col="grey90",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<-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(cs_cb$P~cs_cb$group,col=colb,at=c(1:4),add=T,
border=cola,axes=FALSE,drawRect=F,lwd=1,wex=0.8)
mg_cb_sum <- cs_cb %>% group_by(id, group, sex_num) %>% summarize(P = mean(P))
beeswarm(mg_cb_sum$P~mg_cb_sum$group,col=colc,at=c(1:4),
add=T,pwpch=(16.5+mg_cb_sum$sex_num),
corral = "wrap")
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=3),
labels=c(rep("",length(seq(range[1],range[2],by=3)))),pos=xrange[1],tck=tckk)
axis(2,las=2,cex.axis=1.1,at=seq(range[1],range[2],by=3),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<- 11.5 ;xpos=1.1
points(xpos,ypos,pch=16,cex=1.4,col=rgb(0.4,0.4,0.4,alpha=0.6))
text(xpos+0.25*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.25*xscal,ypos-0.08*scal,"Males",col="grey40",cex=1.2)
text(mean(xrange),range[2]+0.05*scal,"P I",cex=1.6,font=3,xpd=T)
### Plotting S ----------------
groupquant<-sapply(post_fix_cs_cb_S[,c(5,6,1,7)],function(p) quantile(p, probs = c(0.025,0.975,0.5)))
groupquant<-exp(groupquant)
cs_cb$group<-factor(cs_cb$group,levels=c('ctrl.wt','eda.wt','ctrl.sca','eda.sca'))
range<-c(0,15);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="",las=1, axes=FALSE,cex.lab=1.2)
rect(xrange[1],range[2],xrange[2],range[1],col="grey90",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<-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(cs_cb$S~cs_cb$group,col=colb,at=c(1:4),add=T,
border=cola,axes=FALSE,drawRect=F,lwd=1,wex=0.8)
mg_cb_sum <- cs_cb %>% group_by(id, group, sex_num) %>% summarize(P = mean(S))
beeswarm(mg_cb_sum$P~mg_cb_sum$group,col=colc,at=c(1:4),
add=T,pwpch=(16.5+mg_cb_sum$sex_num),
corral = "wrap")
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=3),
labels=c(rep("",length(seq(range[1],range[2],by=3)))),pos=xrange[1],tck=tckk)
axis(2,las=2,cex.axis=1.1,at=seq(range[1],range[2],by=3),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"))
text(mean(xrange),range[2]+0.05*scal,"P I+II",cex=1.6,font=3,xpd=T)
### Plotting Fccp ----------------
groupquant<-sapply(post_fix_cs_cb_Fccp[,c(5,6,1,7)],function(p) quantile(p, probs = c(0.025,0.975,0.5)))
groupquant<-exp(groupquant)
cs_cb$group<-factor(cs_cb$group,levels=c('ctrl.wt','eda.wt','ctrl.sca','eda.sca'))
range<-c(0,21);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="",las=1, axes=FALSE,cex.lab=1.2)
rect(xrange[1],range[2],xrange[2],range[1],col="grey90",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<-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(cs_cb$Fccp~cs_cb$group,col=colb,at=c(1:4),add=T,
border=cola,axes=FALSE,drawRect=F,lwd=1,wex=0.8)
mg_cb_sum <- cs_cb %>% group_by(id, group, sex_num) %>% summarize(P = mean(Fccp))
beeswarm(mg_cb_sum$P~mg_cb_sum$group,col=colc,at=c(1:4),
add=T,pwpch=(16.5+mg_cb_sum$sex_num),
corral = "wrap")
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=3),
labels=c(rep("",length(seq(range[1],range[2],by=3)))),pos=xrange[1],tck=tckk)
axis(2,las=2,cex.axis=1.1,at=seq(range[1],range[2],by=3),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"))
text(mean(xrange),range[2]+0.05*scal,"E I+II",cex=1.6,font=3,xpd=T)
### Plotting Rot ----------------
groupquant<-sapply(post_fix_cs_cb_Rot[,c(5,6,1,7)],function(p) quantile(p, probs = c(0.025,0.975,0.5)))
groupquant<-exp(groupquant)
cs_cb$group<-factor(cs_cb$group,levels=c('ctrl.wt','eda.wt','ctrl.sca','eda.sca'))
range<-c(0,9);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="",las=1, axes=FALSE,cex.lab=1.2)
rect(xrange[1],range[2],xrange[2],range[1],col="grey90",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<-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(cs_cb$Rot~cs_cb$group,col=colb,at=c(1:4),add=T,
border=cola,axes=FALSE,drawRect=F,lwd=1,wex=0.8)
mg_cb_sum <- cs_cb %>% group_by(id, group, sex_num) %>% summarize(P = mean(Rot))
beeswarm(mg_cb_sum$P~mg_cb_sum$group,col=colc,at=c(1:4),
add=T,pwpch=(16.5+mg_cb_sum$sex_num),
corral = "wrap")
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=3),
labels=c(rep("",length(seq(range[1],range[2],by=3)))),pos=xrange[1],tck=tckk)
axis(2,las=2,cex.axis=1.1,at=seq(range[1],range[2],by=3),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"))
text(mean(xrange),range[2]+0.05*scal,"E II",cex=1.6,font=3,xpd=T)
### Plotting Civ ----------------
groupquant<-sapply(post_fix_cs_cb_Civ[,c(5,6,1,7)],function(p) quantile(p, probs = c(0.025,0.975,0.5)))
groupquant<-exp(groupquant)
cs_cb$group<-factor(cs_cb$group,levels=c('ctrl.wt','eda.wt','ctrl.sca','eda.sca'))
range<-c(0,54);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="",las=1, axes=FALSE,cex.lab=1.2)
rect(xrange[1],range[2],xrange[2],range[1],col="grey90",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<-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(cs_cb$C_iv~cs_cb$group,col=colb,at=c(1:4),add=T,
border=cola,axes=FALSE,drawRect=F,lwd=1,wex=0.8)
mg_cb_sum <- cs_cb %>% group_by(id, group, sex_num) %>% summarize(P = mean(C_iv))
beeswarm(mg_cb_sum$P~mg_cb_sum$group,col=colc,at=c(1:4),
add=T,pwpch=(16.5+mg_cb_sum$sex_num),
corral = "wrap")
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=3),
labels=c(rep("",length(seq(range[1],range[2],by=3)))),pos=xrange[1],tck=tckk)
axis(2,las=2,cex.axis=1.1,at=seq(range[1],range[2],by=9),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"))
text(mean(xrange),range[2]+0.05*scal,"Complex IV",cex=1.6,font=3,xpd=T)
## probability distributions -------------
### PD - P ----------
par(mar=c(2,1,0,0))
dif<-data.frame(post_fix_cs_cb_P$wt_contrast,post_fix_cs_cb_P$b_P_groupeda.sca,
post_fix_cs_cb_P$b_P_groupctrl.wt)
dif<-exp(dif)
yla<-"Fold-difference in respiration"
tckk=-0.018
ste<-seq(0.2,2.2,by=0.2)
mons_poste(dif,1)
zpos=seq(0,1,1/3)
xx=1;ind=0.5;xpol<-0.4;srt=55
text(xpol,zpos[xx]+zpos[2]*ind,
"Wt: ed x 0 ",cex=1,srt=srt,xpd=T)
xx=xx+1
text(xpol,zpos[xx]+zpos[2]*ind,
"SCA1: ed x 0 ",cex=1,srt=srt,xpd=T)
xx=xx+1
text(xpol,zpos[xx]+zpos[2]*ind,
"0: Wt x SCA1",cex=1,srt=srt,xpd=T)
### PD - _S ----------
dif<-data.frame(post_fix_cs_cb_S$wt_contrast,post_fix_cs_cb_S$b_S_groupeda.sca,
post_fix_cs_cb_S$b_S_groupctrl.wt)
dif<-exp(dif)
yla<-"Fold-difference in respiration"
tckk=-0.018
ste<-seq(0.2,2.2,by=0.2)
mons_poste(dif,1)
### PD - _Fccp ----------
dif<-data.frame(post_fix_cs_cb_Fccp$wt_contrast,post_fix_cs_cb_Fccp$b_Fccp_groupeda.sca,
post_fix_cs_cb_Fccp$b_Fccp_groupctrl.wt)
dif<-exp(dif)
yla<-"Fold-difference in respiration"
tckk=-0.018
ste<-seq(0.2,2.2,by=0.2)
mons_poste(dif,1)
### PD - _Rot ----------
dif<-data.frame(post_fix_cs_cb_Rot$wt_contrast,post_fix_cs_cb_Rot$b_Rot_groupeda.sca,
post_fix_cs_cb_Rot$b_Rot_groupctrl.wt)
dif<-exp(dif)
yla<-"Fold-difference in respiration"
tckk=-0.018
ste<-seq(0.2,2.2,by=0.2)
mons_poste(dif,1)
### PD - _Civ ----------
dif<-data.frame(post_fix_cs_cb_Civ$wt_contrast,post_fix_cs_cb_Civ$b_Civ_groupeda.sca,
post_fix_cs_cb_Civ$b_Civ_groupctrl.wt)
dif<-exp(dif)
yla<-"Fold-difference in respiration"
tckk=-0.018
ste<-seq(0.2,2.2,by=0.2)
mons_poste(dif,1)
} ) )