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)
} ) )Experimental treatment with edaravone in a mouse model of spinocerebellar ataxia 1
Analysis of gait abnormalities (CatWalk)
This page shows R code for the study Sucha et al. (2023, International Journal of Molecular Science).
Citation:
Sucha M, Benediktova S, Tichanek F, Jedlicka J, Kapl S, Jelinkova D, Purkartova Z, Tuma J, Kuncova J, Cendelin J. Experimental Treatment with Edaravone in a Mouse Model of Spinocerebellar Ataxia 1. International Journal of Molecular Sciences. 2023; 24(13):10689. https://doi.org/10.3390/ijms241310689
GitHub page: https://github.com/filip-tichanek/edaravonSCA1
1 Upload of packages
2 Data import, wrangling and summarization
## import of data
load("source_data/myEnvironment.RData")
## selection of gait parameters
names(behav) [1] "id" "treatment"
[3] "genotype" "sex"
[5] "sugar_water_ratio" "immobile_min"
[7] "speed_95th" "distance_01"
[9] "distance_02" "distance_03"
[11] "distance_04" "distance_05"
[13] "distance_06" "distance_07"
[15] "distance_08" "distance_09"
[17] "distance_10" "inactivity_01"
[19] "inactivity_02" "inactivity_03"
[21] "inactivity_04" "inactivity_05"
[23] "inactivity_06" "inactivity_07"
[25] "inactivity_08" "inactivity_09"
[27] "inactivity_10" "thigmo_propdist"
[29] "center_propdist" "center_avedist"
[31] "thigmo_proptime" "center_proptime"
[33] "center_avetime" "kurva_index"
[35] "t_session1" "t_session2"
[37] "t_session3" "t_session4"
[39] "t_session5" "t_session6"
[41] "t_session7" "t_session8"
[43] "t_session9" "t_session10"
[45] "t_session11" "t_session12"
[47] "grip_strength1" "grip_strength2"
[49] "grip_strength3" "grip_strength4"
[51] "Run_Average_Speed_.cm.s._Mean" "Run_Average_Speed_.cm.s._SD"
[53] "RF_Stand_.s._Mean" "RF_Stand_.s._SD"
[55] "RF_Swing_.s._Mean" "RF_Swing_.s._SD"
[57] "RF_SwingSpeed_.cm.s._Mean" "RF_SwingSpeed_.cm.s._SD"
[59] "RF_StrideLength_.cm._Mean" "RF_StrideLength_.cm._SD"
[61] "RH_Stand_.s._Mean" "RH_Stand_.s._SD"
[63] "RH_Swing_.s._Mean" "RH_Swing_.s._SD"
[65] "RH_SwingSpeed_.cm.s._Mean" "RH_SwingSpeed_.cm.s._SD"
[67] "RH_StrideLength_.cm._Mean" "RH_StrideLength_.cm._SD"
[69] "LF_Stand_.s._Mean" "LF_Stand_.s._SD"
[71] "LF_Swing_.s._Mean" "LF_Swing_.s._SD"
[73] "LF_SwingSpeed_.cm.s._Mean" "LF_SwingSpeed_.cm.s._SD"
[75] "LF_StrideLength_.cm._Mean" "LF_StrideLength_.cm._SD"
[77] "LH_Stand_.s._Mean" "LH_Stand_.s._SD"
[79] "LH_Swing_.s._Mean" "LH_Swing_.s._SD"
[81] "LH_SwingSpeed_.cm.s._Mean" "LH_SwingSpeed_.cm.s._SD"
[83] "LH_StrideLength_.cm._Mean" "LH_StrideLength_.cm._SD"
[85] "StepSequence_CA_..." "StepSequence_CB_..."
[87] "StepSequence_AA_..." "StepSequence_AB_..."
[89] "StepSequence_RA_..." "StepSequence_RB_..."
[91] "StepSequence_RegularityIndex_..." "BOS_FrontPaws_Mean_.cm."
[93] "BOS_HindPaws_Mean_.cm." "OtherStatistics_Average_Speed_Mean"
[95] "OtherStatistics_Average_Speed_SD" "PrintPositions_RightPaws_Mean_.cm."
[97] "PrintPositions_LeftPaws_Mean_.cm." "Support_Zero_..."
[99] "Support_Single_..." "Support_Diagonal_..."
[101] "Support_Girdle_..." "Support_Lateral_..."
[103] "Support_Three_..." "Support_Four_..."
[105] "group" "sex_num"
[107] "group_sex" catwalk <- behav[,c(1,105:107,seq(52,84,by=2),seq(51,83,by=2),85:88,98:104,92:93,96:97)]
names(catwalk)[5:53] [1] "Run_Average_Speed_.cm.s._SD" "RF_Stand_.s._SD"
[3] "RF_Swing_.s._SD" "RF_SwingSpeed_.cm.s._SD"
[5] "RF_StrideLength_.cm._SD" "RH_Stand_.s._SD"
[7] "RH_Swing_.s._SD" "RH_SwingSpeed_.cm.s._SD"
[9] "RH_StrideLength_.cm._SD" "LF_Stand_.s._SD"
[11] "LF_Swing_.s._SD" "LF_SwingSpeed_.cm.s._SD"
[13] "LF_StrideLength_.cm._SD" "LH_Stand_.s._SD"
[15] "LH_Swing_.s._SD" "LH_SwingSpeed_.cm.s._SD"
[17] "LH_StrideLength_.cm._SD" "Run_Average_Speed_.cm.s._Mean"
[19] "RF_Stand_.s._Mean" "RF_Swing_.s._Mean"
[21] "RF_SwingSpeed_.cm.s._Mean" "RF_StrideLength_.cm._Mean"
[23] "RH_Stand_.s._Mean" "RH_Swing_.s._Mean"
[25] "RH_SwingSpeed_.cm.s._Mean" "RH_StrideLength_.cm._Mean"
[27] "LF_Stand_.s._Mean" "LF_Swing_.s._Mean"
[29] "LF_SwingSpeed_.cm.s._Mean" "LF_StrideLength_.cm._Mean"
[31] "LH_Stand_.s._Mean" "LH_Swing_.s._Mean"
[33] "LH_SwingSpeed_.cm.s._Mean" "LH_StrideLength_.cm._Mean"
[35] "StepSequence_CA_..." "StepSequence_CB_..."
[37] "StepSequence_AA_..." "StepSequence_AB_..."
[39] "Support_Zero_..." "Support_Single_..."
[41] "Support_Diagonal_..." "Support_Girdle_..."
[43] "Support_Lateral_..." "Support_Three_..."
[45] "Support_Four_..." "BOS_FrontPaws_Mean_.cm."
[47] "BOS_HindPaws_Mean_.cm." "PrintPositions_RightPaws_Mean_.cm."
[49] "PrintPositions_LeftPaws_Mean_.cm." ## summary of gait parameters
catwalk2 <- catwalk[,c(2, 5:53)]
catwalk2 %>% tbl_summary(by = 'group') %>% add_p() %>% add_q(method = "fdr")add_q: Adjusting p-values with
`stats::p.adjust(x$table_body$p.value, method = "fdr")`| Characteristic | ctrl.wt, N = 201 | eda.wt, N = 201 | ctrl.sca, N = 201 | eda.sca, N = 201 | p-value2 | q-value3 |
|---|---|---|---|---|---|---|
| Run_Average_Speed_.cm.s._SD | 7.4 (3.7, 10.9) | 8.0 (4.6, 10.2) | 4.9 (2.6, 7.4) | 5.8 (4.4, 8.0) | 0.4 | 0.5 |
| RF_Stand_.s._SD | 0.036 (0.027, 0.045) | 0.042 (0.027, 0.047) | 0.044 (0.034, 0.061) | 0.037 (0.032, 0.052) | 0.3 | 0.4 |
| RF_Swing_.s._SD | 0.026 (0.018, 0.033) | 0.028 (0.024, 0.036) | 0.028 (0.022, 0.035) | 0.026 (0.021, 0.030) | 0.8 | 0.8 |
| RF_SwingSpeed_.cm.s._SD | 16 (8, 22) | 18 (12, 29) | 20 (16, 23) | 22 (18, 28) | 0.13 | 0.2 |
| RF_StrideLength_.cm._SD | 0.76 (0.61, 1.09) | 0.87 (0.65, 1.14) | 0.88 (0.60, 1.37) | 0.99 (0.76, 1.18) | 0.7 | 0.8 |
| RH_Stand_.s._SD | 0.04 (0.03, 0.04) | 0.04 (0.03, 0.05) | 0.04 (0.03, 0.07) | 0.03 (0.03, 0.06) | 0.3 | 0.4 |
| RH_Swing_.s._SD | 0.029 (0.025, 0.042) | 0.036 (0.026, 0.043) | 0.042 (0.028, 0.054) | 0.039 (0.025, 0.059) | 0.4 | 0.5 |
| RH_SwingSpeed_.cm.s._SD | 13 (7, 18) | 14 (10, 21) | 14 (11, 18) | 16 (11, 20) | 0.6 | 0.7 |
| RH_StrideLength_.cm._SD | 0.80 (0.68, 1.19) | 1.02 (0.75, 1.11) | 0.94 (0.71, 1.31) | 1.00 (0.75, 1.32) | 0.6 | 0.7 |
| LF_Stand_.s._SD | 0.03 (0.02, 0.05) | 0.04 (0.03, 0.05) | 0.04 (0.03, 0.07) | 0.04 (0.03, 0.06) | 0.3 | 0.4 |
| LF_Swing_.s._SD | 0.025 (0.019, 0.036) | 0.029 (0.022, 0.033) | 0.027 (0.022, 0.034) | 0.025 (0.021, 0.032) | 0.9 | 0.9 |
| LF_SwingSpeed_.cm.s._SD | 17 (8, 23) | 20 (14, 25) | 21 (16, 24) | 24 (16, 28) | 0.2 | 0.3 |
| LF_StrideLength_.cm._SD | 0.73 (0.61, 1.14) | 0.86 (0.73, 1.14) | 0.86 (0.67, 1.27) | 1.01 (0.71, 1.29) | 0.6 | 0.7 |
| LH_Stand_.s._SD | 0.032 (0.026, 0.042) | 0.039 (0.026, 0.046) | 0.038 (0.032, 0.072) | 0.040 (0.033, 0.063) | 0.2 | 0.3 |
| LH_Swing_.s._SD | 0.028 (0.024, 0.046) | 0.032 (0.026, 0.045) | 0.038 (0.029, 0.050) | 0.036 (0.028, 0.052) | 0.4 | 0.5 |
| LH_SwingSpeed_.cm.s._SD | 13 (7, 17) | 13 (10, 16) | 14 (10, 18) | 18 (11, 22) | 0.2 | 0.4 |
| LH_StrideLength_.cm._SD | 0.90 (0.70, 1.19) | 0.88 (0.71, 1.23) | 0.96 (0.65, 1.25) | 1.10 (0.72, 1.29) | 0.8 | 0.9 |
| Run_Average_Speed_.cm.s._Mean | 26 (23, 28) | 28 (22, 31) | 23 (18, 27) | 27 (23, 32) | 0.11 | 0.2 |
| RF_Stand_.s._Mean | 0.15 (0.14, 0.16) | 0.14 (0.13, 0.15) | 0.17 (0.14, 0.20) | 0.15 (0.12, 0.16) | 0.048 | 0.12 |
| RF_Swing_.s._Mean | 0.128 (0.118, 0.138) | 0.114 (0.109, 0.131) | 0.112 (0.094, 0.139) | 0.093 (0.082, 0.123) | 0.006 | 0.040 |
| RF_SwingSpeed_.cm.s._Mean | 55 (51, 62) | 65 (55, 71) | 59 (48, 73) | 70 (59, 82) | 0.040 | 0.12 |
| RF_StrideLength_.cm._Mean | 6.67 (6.34, 7.02) | 6.95 (6.31, 7.30) | 6.31 (5.91, 6.62) | 6.35 (6.18, 6.75) | 0.017 | 0.074 |
| RH_Stand_.s._Mean | 0.119 (0.108, 0.126) | 0.112 (0.103, 0.120) | 0.129 (0.111, 0.160) | 0.108 (0.086, 0.118) | 0.048 | 0.12 |
| RH_Swing_.s._Mean | 0.16 (0.15, 0.18) | 0.15 (0.13, 0.17) | 0.15 (0.14, 0.19) | 0.14 (0.11, 0.18) | 0.4 | 0.5 |
| RH_SwingSpeed_.cm.s._Mean | 44 (40, 51) | 51 (43, 58) | 43 (34, 52) | 51 (40, 62) | 0.3 | 0.4 |
| RH_StrideLength_.cm._Mean | 6.63 (6.25, 7.01) | 6.92 (6.32, 7.38) | 6.31 (5.96, 6.59) | 6.37 (6.14, 6.90) | 0.022 | 0.091 |
| LF_Stand_.s._Mean | 0.15 (0.13, 0.16) | 0.14 (0.13, 0.15) | 0.17 (0.14, 0.19) | 0.15 (0.12, 0.16) | 0.066 | 0.14 |
| LF_Swing_.s._Mean | 0.127 (0.118, 0.141) | 0.118 (0.109, 0.128) | 0.106 (0.092, 0.142) | 0.099 (0.084, 0.117) | 0.003 | 0.025 |
| LF_SwingSpeed_.cm.s._Mean | 55 (51, 64) | 65 (53, 70) | 59 (47, 71) | 70 (61, 85) | 0.040 | 0.12 |
| LF_StrideLength_.cm._Mean | 6.65 (6.25, 7.07) | 6.88 (6.38, 7.34) | 6.33 (5.91, 6.58) | 6.37 (6.14, 6.77) | 0.012 | 0.060 |
| LH_Stand_.s._Mean | 0.120 (0.107, 0.138) | 0.113 (0.095, 0.122) | 0.125 (0.110, 0.149) | 0.106 (0.094, 0.120) | 0.061 | 0.14 |
| LH_Swing_.s._Mean | 0.15 (0.15, 0.18) | 0.15 (0.14, 0.17) | 0.14 (0.13, 0.19) | 0.14 (0.11, 0.16) | 0.3 | 0.4 |
| LH_SwingSpeed_.cm.s._Mean | 45 (40, 50) | 48 (43, 57) | 44 (34, 54) | 52 (41, 62) | 0.3 | 0.4 |
| LH_StrideLength_.cm._Mean | 6.61 (6.30, 7.04) | 6.98 (6.36, 7.29) | 6.25 (5.99, 6.58) | 6.34 (6.10, 6.89) | 0.011 | 0.060 |
| StepSequence_CA_... | 23 (14, 35) | 23 (13, 35) | 10 (8, 19) | 13 (7, 19) | 0.025 | 0.10 |
| StepSequence_CB_... | 22 (12, 28) | 22 (13, 34) | 21 (7, 26) | 14 (8, 23) | 0.4 | 0.5 |
| StepSequence_AA_... | 12 (5, 22) | 11 (8, 14) | 1 (0, 5) | 1 (0, 6) | <0.001 | <0.001 |
| StepSequence_AB_... | 38 (20, 57) | 34 (25, 48) | 63 (53, 74) | 65 (54, 78) | <0.001 | <0.001 |
| Support_Zero_... | 0.06 (0.00, 0.64) | 0.19 (0.00, 0.49) | 0.11 (0.00, 0.17) | 0.00 (0.00, 0.12) | 0.13 | 0.2 |
| Support_Single_... | 13 (9, 18) | 16 (11, 18) | 9 (5, 15) | 8 (6, 14) | 0.055 | 0.13 |
| Support_Diagonal_... | 74 (71, 78) | 72 (69, 74) | 69 (62, 73) | 64 (61, 74) | 0.010 | 0.060 |
| Support_Girdle_... | 3.20 (2.43, 5.02) | 3.71 (3.07, 5.22) | 5.86 (4.49, 7.14) | 6.80 (4.64, 8.98) | <0.001 | 0.006 |
| Support_Lateral_... | 2.16 (1.47, 3.06) | 2.52 (1.61, 3.53) | 1.23 (0.76, 1.96) | 1.70 (1.33, 2.40) | 0.047 | 0.12 |
| Support_Three_... | 7 (5, 8) | 4 (4, 6) | 11 (7, 19) | 12 (9, 15) | <0.001 | <0.001 |
| Support_Four_... | 0.31 (0.10, 0.62) | 0.07 (0.00, 0.35) | 0.93 (0.18, 3.25) | 0.39 (0.11, 1.53) | 0.031 | 0.11 |
| BOS_FrontPaws_Mean_.cm. | 1.23 (1.15, 1.31) | 1.17 (1.12, 1.26) | 1.12 (1.07, 1.24) | 1.12 (1.04, 1.16) | 0.049 | 0.12 |
| BOS_HindPaws_Mean_.cm. | 2.22 (2.09, 2.36) | 2.21 (2.12, 2.28) | 2.21 (2.00, 2.46) | 2.16 (2.06, 2.28) | 0.7 | 0.8 |
| PrintPositions_RightPaws_Mean_.cm. | 0.72 (0.47, 0.95) | 0.52 (0.32, 0.63) | 0.58 (0.26, 0.85) | 0.43 (0.26, 0.57) | 0.12 | 0.2 |
| PrintPositions_LeftPaws_Mean_.cm. | 0.69 (0.56, 0.92) | 0.49 (0.37, 0.63) | 0.54 (0.40, 0.85) | 0.33 (0.18, 0.48) | 0.004 | 0.035 |
| 1 Median (IQR) | ||||||
| 2 Kruskal-Wallis rank sum test | ||||||
| 3 False discovery rate correction for multiple testing | ||||||
Transformation of gait parameters (proportions by logit, SD of gait parameters with log-transformation)
## tranformation
catwalk[,c(5:21)]<-log(catwalk[,c(5:21)])
catwalk[,c(39:49)]<- logit( (catwalk[,c(39:49)]/100) +0.01)
## scaling
catwalk[,5:53]<-scale(catwalk[,5:53])
## getting principal components and % explained variance
pca<-princomp(catwalk[,5:53])
plot(pca)
round( ((pca$sdev**2)/sum(pca$sdev**2)),3) Comp.1 Comp.2 Comp.3 Comp.4 Comp.5 Comp.6 Comp.7 Comp.8 Comp.9 Comp.10
0.297 0.268 0.117 0.064 0.039 0.032 0.027 0.023 0.018 0.016
Comp.11 Comp.12 Comp.13 Comp.14 Comp.15 Comp.16 Comp.17 Comp.18 Comp.19 Comp.20
0.015 0.013 0.011 0.008 0.007 0.005 0.005 0.004 0.004 0.004
Comp.21 Comp.22 Comp.23 Comp.24 Comp.25 Comp.26 Comp.27 Comp.28 Comp.29 Comp.30
0.003 0.003 0.002 0.002 0.002 0.002 0.001 0.001 0.001 0.001
Comp.31 Comp.32 Comp.33 Comp.34 Comp.35 Comp.36 Comp.37 Comp.38 Comp.39 Comp.40
0.001 0.001 0.001 0.001 0.000 0.000 0.000 0.000 0.000 0.000
Comp.41 Comp.42 Comp.43 Comp.44 Comp.45 Comp.46 Comp.47 Comp.48 Comp.49
0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.000 sum( round( ((pca$sdev**2)/sum(pca$sdev**2)),3) [1:3])[1] 0.682## extraction the PCs as another variables
catwalk$pca1<- c(scale(pca$scores[,1]))
catwalk$pca2<- c(scale(pca$scores[,2]))
catwalk$pca3<- c(scale(pca$scores[,3]))3 principal components represent majority of variability in data of gait
2.1 data exploration
## setting spacing
par(mar=c(7,3.5,1,1))
par(mgp=c(5,0.7,0))
par(mfrow = c(1,3))
## plotting values of the 3 main PCs across groups
plot(catwalk$pca1~catwalk$group_sex,
col=cola[c(1,1,2,2,3,3,4,4)], las=2)
plot(catwalk$pca2~catwalk$group_sex,
col=cola[c(1,1,2,2,3,3,4,4)], las=2)
plot(catwalk$pca3~catwalk$group_sex,
col=cola[c(1,1,2,2,3,3,4,4)], las=2)
par(mfrow = c(1, 3))
plot(catwalk$pca1~catwalk$group,col=cola, las=2)
plot(catwalk$pca2~catwalk$group,col=cola, las=2)
plot(catwalk$pca3~catwalk$group,col=cola, las=2)
Integration of PCs score to general table and assessing their correlations with individual gait parameters
## Extraction
behav$walk_pc1<- c(scale(pca$scores[,1]))
behav$walk_pc2<- c(scale(pca$scores[,2]))
behav$walk_pc3<- c(scale(pca$scores[,3]))
## Correlation of PCs with gait parameters
cor1<-c(cor(catwalk$pca1,catwalk[,5:53],method="pearson"))
cor2<-c(cor(catwalk$pca2,catwalk[,5:53],method="pearson"))
cor3<-c(cor(catwalk$pca3,catwalk[,5:53],method="pearson"))
pc_gait_cor<-data.frame(cbind(names(catwalk[,5:53]),
cor1,cor2,cor3))
pc_gait_cor$cor1<-as.numeric(pc_gait_cor$cor1)
pc_gait_cor$cor2<-as.numeric(pc_gait_cor$cor2)
pc_gait_cor$cor3<-as.numeric(pc_gait_cor$cor3)
pc_gait_cor<-data.frame(parameter = pc_gait_cor$V1,
PC1_cor = cor1,
PC2_cor = cor2,
PC3_cor = cor3)
## Correlations of all axes with all parameters
knitr::kable(pc_gait_cor)| parameter | PC1_cor | PC2_cor | PC3_cor |
|---|---|---|---|
| Run_Average_Speed_.cm.s._SD | 0.2672053 | 0.6095708 | 0.3561662 |
| RF_Stand_.s._SD | -0.6625579 | 0.6548039 | 0.1116333 |
| RF_Swing_.s._SD | -0.4970953 | 0.5323204 | 0.4379428 |
| RF_SwingSpeed_.cm.s._SD | 0.2303728 | 0.8915512 | 0.1911648 |
| RF_StrideLength_.cm._SD | -0.1581137 | 0.8213267 | 0.2811808 |
| RH_Stand_.s._SD | -0.5495105 | 0.6637357 | 0.0013374 |
| RH_Swing_.s._SD | -0.5164171 | 0.6159595 | 0.2649329 |
| RH_SwingSpeed_.cm.s._SD | 0.2379274 | 0.8582441 | 0.2234228 |
| RH_StrideLength_.cm._SD | -0.1944708 | 0.7945607 | 0.3271262 |
| LF_Stand_.s._SD | -0.6342372 | 0.6517018 | 0.0938350 |
| LF_Swing_.s._SD | -0.4188181 | 0.5059832 | 0.5451368 |
| LF_SwingSpeed_.cm.s._SD | 0.1613490 | 0.8980204 | 0.1137190 |
| LF_StrideLength_.cm._SD | -0.2103648 | 0.7978682 | 0.3177006 |
| LH_Stand_.s._SD | -0.6582108 | 0.6496784 | -0.0457254 |
| LH_Swing_.s._SD | -0.3323126 | 0.6148714 | 0.4479048 |
| LH_SwingSpeed_.cm.s._SD | 0.2507938 | 0.8393571 | 0.0126718 |
| LH_StrideLength_.cm._SD | -0.0651910 | 0.7609980 | 0.4168295 |
| Run_Average_Speed_.cm.s._Mean | 0.9004667 | 0.3414745 | 0.1116079 |
| RF_Stand_.s._Mean | -0.9110190 | -0.1901959 | 0.0645795 |
| RF_Swing_.s._Mean | -0.5214467 | -0.6390176 | 0.4954172 |
| RF_SwingSpeed_.cm.s._Mean | 0.7708241 | 0.5397599 | -0.1303897 |
| RF_StrideLength_.cm._Mean | 0.7666779 | -0.2154543 | 0.3955947 |
| RH_Stand_.s._Mean | -0.7881435 | -0.1053852 | -0.0662461 |
| RH_Swing_.s._Mean | -0.6218096 | -0.4758645 | 0.4438371 |
| RH_SwingSpeed_.cm.s._Mean | 0.8066497 | 0.4620823 | -0.1070256 |
| RH_StrideLength_.cm._Mean | 0.7576182 | -0.1996631 | 0.3866716 |
| LF_Stand_.s._Mean | -0.9249427 | -0.1788519 | 0.0182499 |
| LF_Swing_.s._Mean | -0.4384682 | -0.6459233 | 0.5276477 |
| LF_SwingSpeed_.cm.s._Mean | 0.7461173 | 0.5609859 | -0.2107481 |
| LF_StrideLength_.cm._Mean | 0.7653601 | -0.2310794 | 0.3892210 |
| LH_Stand_.s._Mean | -0.8352683 | -0.0606679 | -0.1319910 |
| LH_Swing_.s._Mean | -0.5268990 | -0.5088357 | 0.5057004 |
| LH_SwingSpeed_.cm.s._Mean | 0.7775097 | 0.4774909 | -0.2267086 |
| LH_StrideLength_.cm._Mean | 0.7533525 | -0.1833291 | 0.4148233 |
| StepSequence_CA_… | 0.0991969 | -0.1318993 | 0.2078292 |
| StepSequence_CB_… | 0.3363822 | -0.1149716 | 0.3563170 |
| StepSequence_AA_… | 0.2104824 | -0.2899867 | 0.4769519 |
| StepSequence_AB_… | -0.3590790 | 0.2682774 | -0.5431625 |
| Support_Zero_… | 0.2829771 | 0.0144798 | 0.5997906 |
| Support_Single_… | 0.2402412 | -0.2521296 | 0.7595718 |
| Support_Diagonal_… | 0.4158972 | -0.5117275 | -0.2713623 |
| Support_Girdle_… | -0.2942092 | 0.1913824 | -0.1087525 |
| Support_Lateral_… | 0.0790337 | 0.1132000 | 0.4384725 |
| Support_Three_… | -0.5224694 | 0.4856310 | -0.5824357 |
| Support_Four_… | -0.4867327 | 0.5440529 | -0.3823720 |
| BOS_FrontPaws_Mean_.cm. | -0.0516410 | 0.0871550 | -0.0738046 |
| BOS_HindPaws_Mean_.cm. | -0.0957587 | -0.5028307 | 0.1940434 |
| PrintPositions_RightPaws_Mean_.cm. | -0.6804820 | -0.2246993 | 0.1745176 |
| PrintPositions_LeftPaws_Mean_.cm. | -0.6513525 | -0.2344030 | 0.1784472 |
## correlates of PC1
knitr::kable(pc_gait_cor[abs(pc_gait_cor$PC1_cor)>0.75,])| parameter | PC1_cor | PC2_cor | PC3_cor | |
|---|---|---|---|---|
| 18 | Run_Average_Speed_.cm.s._Mean | 0.9004667 | 0.3414745 | 0.1116079 |
| 19 | RF_Stand_.s._Mean | -0.9110190 | -0.1901959 | 0.0645795 |
| 21 | RF_SwingSpeed_.cm.s._Mean | 0.7708241 | 0.5397599 | -0.1303897 |
| 22 | RF_StrideLength_.cm._Mean | 0.7666779 | -0.2154543 | 0.3955947 |
| 23 | RH_Stand_.s._Mean | -0.7881435 | -0.1053852 | -0.0662461 |
| 25 | RH_SwingSpeed_.cm.s._Mean | 0.8066497 | 0.4620823 | -0.1070256 |
| 26 | RH_StrideLength_.cm._Mean | 0.7576182 | -0.1996631 | 0.3866716 |
| 27 | LF_Stand_.s._Mean | -0.9249427 | -0.1788519 | 0.0182499 |
| 30 | LF_StrideLength_.cm._Mean | 0.7653601 | -0.2310794 | 0.3892210 |
| 31 | LH_Stand_.s._Mean | -0.8352683 | -0.0606679 | -0.1319910 |
| 33 | LH_SwingSpeed_.cm.s._Mean | 0.7775097 | 0.4774909 | -0.2267086 |
| 34 | LH_StrideLength_.cm._Mean | 0.7533525 | -0.1833291 | 0.4148233 |
## Correlates of PC2
knitr::kable(pc_gait_cor[abs(pc_gait_cor$PC2_cor)>0.75,])| parameter | PC1_cor | PC2_cor | PC3_cor | |
|---|---|---|---|---|
| 4 | RF_SwingSpeed_.cm.s._SD | 0.2303728 | 0.8915512 | 0.1911648 |
| 5 | RF_StrideLength_.cm._SD | -0.1581137 | 0.8213267 | 0.2811808 |
| 8 | RH_SwingSpeed_.cm.s._SD | 0.2379274 | 0.8582441 | 0.2234228 |
| 9 | RH_StrideLength_.cm._SD | -0.1944708 | 0.7945607 | 0.3271262 |
| 12 | LF_SwingSpeed_.cm.s._SD | 0.1613490 | 0.8980204 | 0.1137190 |
| 13 | LF_StrideLength_.cm._SD | -0.2103648 | 0.7978682 | 0.3177006 |
| 16 | LH_SwingSpeed_.cm.s._SD | 0.2507938 | 0.8393571 | 0.0126718 |
| 17 | LH_StrideLength_.cm._SD | -0.0651910 | 0.7609980 | 0.4168295 |
## correlates of PC3
knitr::kable(pc_gait_cor[abs(pc_gait_cor$PC3_cor)>0.55,])| parameter | PC1_cor | PC2_cor | PC3_cor | |
|---|---|---|---|---|
| 39 | Support_Zero_… | 0.2829771 | 0.0144798 | 0.5997906 |
| 40 | Support_Single_… | 0.2402412 | -0.2521296 | 0.7595718 |
| 44 | Support_Three_… | -0.5224694 | 0.4856310 | -0.5824357 |
3 Modelling
Statistical analysis will consist from 2 parts:
- Multivariate analysis in frequentist framework, using PERMANOVA
- Bayesian multivariate regression
1st part will assees which groups differed ‘statistically significanlty’ in terms of conglomerate of all parameters
2nd part will employ Bayesian regression and will estimate how much the groups differ in terms of the three main principal components (PCs explaining together majority of variance in the gait data)
As the overall structure of gait is affected by walk speed (which differ between WT and SCA1 mice, but also differ substatially from animal to animal), the 2nd analysis will be also validated with fitting same model but adjusted for speed (walk speed will be included as a covariate)
3.1 PERMANOVA
## defining outcome
outcome <- catwalk[,5:53]
## PERMANOVA, all groups
if(file.exists('perm_all') == FALSE){
perm_all <- adonis2(outcome ~ group, data=catwalk, perm=999,method="euclidian")}
perm_allPermutation test for adonis under reduced model
Terms added sequentially (first to last)
Permutation: free
Number of permutations: 999
adonis2(formula = outcome ~ group, data = catwalk, permutations = 999, method = "euclidian")
Df SumOfSqs R2 F Pr(>F)
group 3 376 0.09712 2.7251 0.002 **
Residual 76 3495 0.90288
Total 79 3871 1.00000
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1### edaravon in SCA1
catwalk_sca<-subset(catwalk,catwalk$group=="ctrl.sca"|catwalk$group=="eda.sca")
outcome<-catwalk_sca[,5:53]
if(file.exists('perm_sca1') == FALSE){
perm_sca1 <- adonis2(outcome ~ group, data=catwalk_sca, perm=99999,method="euclidian") }
perm_sca1Permutation test for adonis under reduced model
Terms added sequentially (first to last)
Permutation: free
Number of permutations: 99999
adonis2(formula = outcome ~ group, data = catwalk_sca, permutations = 99999, method = "euclidian")
Df SumOfSqs R2 F Pr(>F)
group 1 70.69 0.03463 1.3631 0.2203
Residual 38 1970.74 0.96537
Total 39 2041.43 1.00000 ### edaravon in WT -
catwalk_wt<-subset(catwalk,catwalk$group=="eda.wt"|catwalk$group=="ctrl.wt")
outcome<-catwalk_wt[,5:53]
if(file.exists('perm_wt') == FALSE){
perm_wt <- adonis2(outcome ~ group, data=catwalk_wt, perm=9999,method="euclidian") }
perm_wtPermutation test for adonis under reduced model
Terms added sequentially (first to last)
Permutation: free
Number of permutations: 9999
adonis2(formula = outcome ~ group, data = catwalk_wt, permutations = 9999, method = "euclidian")
Df SumOfSqs R2 F Pr(>F)
group 1 38.24 0.02447 0.9532 0.4241
Residual 38 1524.30 0.97553
Total 39 1562.54 1.00000 ### geenotype in ctrl
set.seed(17)
catwalk_ctrl<-subset(catwalk,catwalk$group=="ctrl.sca"|catwalk$group=="ctrl.wt")
outcome<-catwalk_ctrl[,5:53]
if(file.exists('perm_ctrl') == FALSE){
perm_ctrl <- adonis2(outcome ~ group, data=catwalk_sca, perm=9999,method="euclidian")}
perm_ctrlPermutation test for adonis under reduced model
Terms added sequentially (first to last)
Permutation: free
Number of permutations: 9999
adonis2(formula = outcome ~ group, data = catwalk_sca, permutations = 9999, method = "euclidian")
Df SumOfSqs R2 F Pr(>F)
group 1 123.51 0.06211 2.5167 0.0303 *
Residual 38 1864.97 0.93789
Total 39 1988.49 1.00000
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1As seen, groups generally differ, with ‘significant’ difference between WT and SCA1 mice on saline. Edaravone treatment does not seem to have effect.
3.2 Bayesian models - without speed adjustment
3.2.1 Priors and formulas
## re-setting reference levels and creating new variable
behav$group <- relevel(behav$group,ref='ctrl.sca')
behav$speed <- scale(behav$Run_Average_Speed_.cm.s._Mean)
## formula
bf_pc1<-bf(walk_pc1~0+Intercept+group)
bf_pc2<-bf(walk_pc2~0+Intercept+group)
bf_pc3<-bf(walk_pc3~0+Intercept+group)
## prior
prior_1 <- get_prior(bf_pc1+bf_pc2+bf_pc3+
set_rescor(TRUE), data = behav,family=gaussian())
prior_1 prior class coef group resp dpar nlpar lb ub
(flat) b
lkj(1) rescor
(flat) b walkpc1
(flat) b groupctrl.wt walkpc1
(flat) b groupeda.sca walkpc1
(flat) b groupeda.wt walkpc1
(flat) b Intercept walkpc1
student_t(3, 0, 2.5) sigma walkpc1 0
(flat) b walkpc2
(flat) b groupctrl.wt walkpc2
(flat) b groupeda.sca walkpc2
(flat) b groupeda.wt walkpc2
(flat) b Intercept walkpc2
student_t(3, 0, 2.5) sigma walkpc2 0
(flat) b walkpc3
(flat) b groupctrl.wt walkpc3
(flat) b groupeda.sca walkpc3
(flat) b groupeda.wt walkpc3
(flat) b Intercept walkpc3
student_t(3, 0, 2.5) sigma walkpc3 0
source
default
default
default
(vectorized)
(vectorized)
(vectorized)
(vectorized)
default
default
(vectorized)
(vectorized)
(vectorized)
(vectorized)
default
default
(vectorized)
(vectorized)
(vectorized)
(vectorized)
defaultprior_1$prior[c(4,6,10,12,16,18)]<-"normal(0,2)"
prior_1$prior[c(5,11,17)]<-"normal(0,1.2)"
prior_1$prior[c(7,13,19)]<-"student_t(3,0,5)"3.2.2 Model fit and diagnostics
## without adjustment for speed
set.seed(17)
pc_mult1 <- run_model(
brm(bf_pc1+bf_pc2+bf_pc3
+set_rescor(TRUE),
data=behav, prior = prior_1,
family=gaussian(),
save_pars = save_pars(all = TRUE),
iter=8000, warmup=2000,chains=2,cores=1,seed=17,
control = list(adapt_delta = 0.99)),
'output/brm/pc_mult1', reuse = TRUE)
mcmc_trace(pc_mult1) 
summary(pc_mult1) Family: MV(gaussian, gaussian, gaussian)
Links: mu = identity; sigma = identity
mu = identity; sigma = identity
mu = identity; sigma = identity
Formula: walk_pc1 ~ 0 + Intercept + group
walk_pc2 ~ 0 + Intercept + group
walk_pc3 ~ 0 + Intercept + group
Data: behav (Number of observations: 80)
Draws: 2 chains, each with iter = 8000; warmup = 2000; thin = 1;
total post-warmup draws = 12000
Population-Level Effects:
Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS
walkpc1_Intercept -0.47 0.21 -0.89 -0.04 1.00 6796
walkpc1_groupctrl.wt 0.52 0.31 -0.08 1.12 1.00 8242
walkpc1_groupeda.wt 0.82 0.31 0.21 1.43 1.00 8474
walkpc1_groupeda.sca 0.53 0.30 -0.07 1.11 1.00 8149
walkpc2_Intercept 0.15 0.21 -0.28 0.57 1.00 5957
walkpc2_groupctrl.wt -0.60 0.30 -1.20 0.00 1.00 7397
walkpc2_groupeda.wt -0.27 0.31 -0.87 0.34 1.00 7407
walkpc2_groupeda.sca 0.27 0.29 -0.31 0.85 1.00 7705
walkpc3_Intercept -0.39 0.21 -0.79 0.01 1.00 6349
walkpc3_groupctrl.wt 0.74 0.29 0.17 1.32 1.00 8287
walkpc3_groupeda.wt 0.85 0.29 0.27 1.42 1.00 8133
walkpc3_groupeda.sca -0.05 0.29 -0.61 0.52 1.00 8158
Tail_ESS
walkpc1_Intercept 7958
walkpc1_groupctrl.wt 8608
walkpc1_groupeda.wt 8339
walkpc1_groupeda.sca 8095
walkpc2_Intercept 6609
walkpc2_groupctrl.wt 8175
walkpc2_groupeda.wt 7616
walkpc2_groupeda.sca 8570
walkpc3_Intercept 7940
walkpc3_groupctrl.wt 8642
walkpc3_groupeda.wt 7921
walkpc3_groupeda.sca 8638
Family Specific Parameters:
Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
sigma_walkpc1 0.99 0.08 0.85 1.17 1.00 13278 8666
sigma_walkpc2 0.99 0.08 0.84 1.17 1.00 13803 7693
sigma_walkpc3 0.95 0.08 0.81 1.12 1.00 14664 8168
Residual Correlations:
Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS
rescor(walkpc1,walkpc2) 0.03 0.11 -0.19 0.26 1.00 12091
rescor(walkpc1,walkpc3) -0.10 0.11 -0.31 0.12 1.00 12416
rescor(walkpc2,walkpc3) 0.13 0.11 -0.09 0.34 1.00 11722
Tail_ESS
rescor(walkpc1,walkpc2) 8352
rescor(walkpc1,walkpc3) 8543
rescor(walkpc2,walkpc3) 8247
Draws were sampled using sampling(NUTS). For each parameter, Bulk_ESS
and Tail_ESS are effective sample size measures, and Rhat is the potential
scale reduction factor on split chains (at convergence, Rhat = 1).Convergence is fine
3.3 Bayesian model with speed adjustment
3.3.1 formula and prior
bf_adj_pc1<-bf(walk_pc1~0+Intercept+group+s(speed,k=3))
bf_adj_pc2<-bf(walk_pc2~0+Intercept+group+s(speed,k=3))
bf_adj_pc3<-bf(walk_pc3~0+Intercept+group+s(speed,k=3))
prior_2 <- get_prior(bf_adj_pc1+bf_adj_pc2+bf_adj_pc3+
set_rescor(TRUE), data = behav,family=gaussian())
prior_2 prior class coef group resp dpar nlpar lb ub
(flat) b
lkj(1) rescor
(flat) b walkpc1
(flat) b groupctrl.wt walkpc1
(flat) b groupeda.sca walkpc1
(flat) b groupeda.wt walkpc1
(flat) b Intercept walkpc1
(flat) b sspeed_1 walkpc1
student_t(3, 0, 2.5) sds walkpc1 0
student_t(3, 0, 2.5) sds s(speed, k = 3) walkpc1 0
student_t(3, 0, 2.5) sigma walkpc1 0
(flat) b walkpc2
(flat) b groupctrl.wt walkpc2
(flat) b groupeda.sca walkpc2
(flat) b groupeda.wt walkpc2
(flat) b Intercept walkpc2
(flat) b sspeed_1 walkpc2
student_t(3, 0, 2.5) sds walkpc2 0
student_t(3, 0, 2.5) sds s(speed, k = 3) walkpc2 0
student_t(3, 0, 2.5) sigma walkpc2 0
(flat) b walkpc3
(flat) b groupctrl.wt walkpc3
(flat) b groupeda.sca walkpc3
(flat) b groupeda.wt walkpc3
(flat) b Intercept walkpc3
(flat) b sspeed_1 walkpc3
student_t(3, 0, 2.5) sds walkpc3 0
student_t(3, 0, 2.5) sds s(speed, k = 3) walkpc3 0
student_t(3, 0, 2.5) sigma walkpc3 0
source
default
default
default
(vectorized)
(vectorized)
(vectorized)
(vectorized)
(vectorized)
default
(vectorized)
default
default
(vectorized)
(vectorized)
(vectorized)
(vectorized)
(vectorized)
default
(vectorized)
default
default
(vectorized)
(vectorized)
(vectorized)
(vectorized)
(vectorized)
default
(vectorized)
defaultprior_2$prior[c(4,6,13,15,22,24)]<-"normal(0,2)"
prior_2$prior[c(5,14,23)]<-"normal(0,1.2)"
prior_2$prior[c(7,16,25)]<-"student_t(3,0,5)"
prior_2$prior[c(8,17,26)]<-"normal(0,2)"3.3.2 Model fit and diagnostics
pc_mult2 <- run_model(
brm(bf_adj_pc1+bf_adj_pc2+bf_adj_pc3+
set_rescor(TRUE),
data=behav, prior = prior_2,
family=gaussian(),
save_pars = save_pars(all = TRUE),
iter=8000, warmup=2000,chains=2,cores=1,seed=17,
control = list(adapt_delta = 0.99)),
'output/brm/pc_mult2', reuse = TRUE)
mcmc_trace(pc_mult2)
summary(pc_mult2)Warning: There were 3 divergent transitions after warmup. Increasing
adapt_delta above 0.99 may help. See
http://mc-stan.org/misc/warnings.html#divergent-transitions-after-warmup Family: MV(gaussian, gaussian, gaussian)
Links: mu = identity; sigma = identity
mu = identity; sigma = identity
mu = identity; sigma = identity
Formula: walk_pc1 ~ 0 + Intercept + group + s(speed, k = 3)
walk_pc2 ~ 0 + Intercept + group + s(speed, k = 3)
walk_pc3 ~ 0 + Intercept + group + s(speed, k = 3)
Data: behav (Number of observations: 80)
Draws: 2 chains, each with iter = 8000; warmup = 2000; thin = 1;
total post-warmup draws = 12000
Smooth Terms:
Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS
sds(walkpc1_sspeed_1) 2.54 2.00 0.22 7.67 1.00 6258
sds(walkpc2_sspeed_1) 2.05 2.15 0.07 6.96 1.00 7763
sds(walkpc3_sspeed_1) 3.46 2.86 0.18 10.66 1.00 6619
Tail_ESS
sds(walkpc1_sspeed_1) 4523
sds(walkpc2_sspeed_1) 4351
sds(walkpc3_sspeed_1) 4665
Population-Level Effects:
Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS
walkpc1_Intercept -0.12 0.10 -0.32 0.07 1.00 3203
walkpc1_groupctrl.wt 0.19 0.14 -0.08 0.46 1.00 3848
walkpc1_groupeda.wt 0.21 0.14 -0.06 0.49 1.00 3895
walkpc1_groupeda.sca 0.09 0.14 -0.18 0.36 1.00 3961
walkpc2_Intercept 0.29 0.20 -0.11 0.70 1.00 3466
walkpc2_groupctrl.wt -0.74 0.29 -1.30 -0.17 1.00 4166
walkpc2_groupeda.wt -0.52 0.29 -1.09 0.06 1.00 4378
walkpc2_groupeda.sca 0.09 0.29 -0.48 0.64 1.00 4130
walkpc3_Intercept -0.40 0.21 -0.82 0.01 1.00 4993
walkpc3_groupctrl.wt 0.79 0.30 0.21 1.38 1.00 5664
walkpc3_groupeda.wt 0.85 0.30 0.27 1.44 1.00 6321
walkpc3_groupeda.sca -0.02 0.29 -0.61 0.56 1.00 6383
walkpc1_sspeed_1 -0.84 0.06 -0.94 -0.73 1.00 7357
walkpc2_sspeed_1 -0.36 0.11 -0.58 -0.14 1.00 8563
walkpc3_sspeed_1 -0.13 0.12 -0.37 0.12 1.00 10248
Tail_ESS
walkpc1_Intercept 5802
walkpc1_groupctrl.wt 6260
walkpc1_groupeda.wt 7120
walkpc1_groupeda.sca 5778
walkpc2_Intercept 5985
walkpc2_groupctrl.wt 7036
walkpc2_groupeda.wt 7231
walkpc2_groupeda.sca 6404
walkpc3_Intercept 7073
walkpc3_groupctrl.wt 8156
walkpc3_groupeda.wt 8349
walkpc3_groupeda.sca 8120
walkpc1_sspeed_1 8640
walkpc2_sspeed_1 8626
walkpc3_sspeed_1 8654
Family Specific Parameters:
Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
sigma_walkpc1 0.43 0.04 0.37 0.51 1.00 6726 7677
sigma_walkpc2 0.91 0.07 0.78 1.07 1.00 7220 8190
sigma_walkpc3 0.93 0.08 0.79 1.11 1.00 11531 8651
Residual Correlations:
Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS
rescor(walkpc1,walkpc2) -0.76 0.05 -0.84 -0.65 1.00 8240
rescor(walkpc1,walkpc3) -0.30 0.10 -0.49 -0.09 1.00 12339
rescor(walkpc2,walkpc3) 0.12 0.11 -0.11 0.33 1.00 12441
Tail_ESS
rescor(walkpc1,walkpc2) 8532
rescor(walkpc1,walkpc3) 8824
rescor(walkpc2,walkpc3) 8512
Draws were sampled using sampling(NUTS). For each parameter, Bulk_ESS
and Tail_ESS are effective sample size measures, and Rhat is the potential
scale reduction factor on split chains (at convergence, Rhat = 1).Convergence is fine
3.3.3 PPC and models compariron
### Bayesian checking - without speed
pp_check(pc_mult1,type='dens_overlay',ndraws = 50,resp="walkpc1")
pp_check(pc_mult1,type='dens_overlay_grouped',ndraws = 50,group='group',resp="walkpc1")
pp_check(pc_mult1,type="stat_2d", stat = c("max", "min"),ndraws=100,resp="walkpc1")
pp_check(pc_mult1,type="stat_2d", stat = c("mean", "sd"),ndraws=200,resp="walkpc1")
pp_check(pc_mult1,type='dens_overlay',ndraws = 50,resp="walkpc2")
pp_check(pc_mult1,type='dens_overlay_grouped',ndraws = 50,group='group',resp="walkpc2")
pp_check(pc_mult1,type="stat_2d", stat = c("max", "min"),ndraws=100,resp="walkpc2")
pp_check(pc_mult1,type="stat_2d", stat = c("mean", "sd"),ndraws=200,resp="walkpc2")
pp_check(pc_mult1,type='dens_overlay',ndraws = 50,resp="walkpc3")
pp_check(pc_mult1,type='dens_overlay_grouped',ndraws = 50,group='group',resp="walkpc3")
pp_check(pc_mult1,type="stat_2d", stat = c("max", "min"),ndraws=100,resp="walkpc3")
pp_check(pc_mult1,type="stat_2d", stat = c("mean", "sd"),ndraws=200,resp="walkpc3")
### Bayesian checking - speed included
pp_check(pc_mult2,type='dens_overlay',ndraws = 50,resp="walkpc1")
pp_check(pc_mult2,type='dens_overlay_grouped',ndraws = 200,group='group',resp="walkpc1")
pp_check(pc_mult2,type="stat_2d", stat = c("max", "min"),ndraws=100,resp="walkpc1")
pp_check(pc_mult2,type="stat_2d", stat = c("mean", "sd"),ndraws=200,resp="walkpc1")
pp_check(pc_mult2,type='dens_overlay',ndraws = 50,resp="walkpc2")
pp_check(pc_mult2,type='dens_overlay_grouped',ndraws = 50,group='group',resp="walkpc2")
pp_check(pc_mult2,type="stat_2d", stat = c("max", "min"),ndraws=100,resp="walkpc2")
pp_check(pc_mult2,type="stat_2d", stat = c("mean", "sd"),ndraws=200,resp="walkpc2")
pp_check(pc_mult2,type='dens_overlay',ndraws = 50,resp="walkpc3")
pp_check(pc_mult2,type='dens_overlay_grouped',ndraws = 50,group='group',resp="walkpc3")
pp_check(pc_mult2,type="stat_2d", stat = c("max", "min"),ndraws=100,resp="walkpc3")
pp_check(pc_mult2,type="stat_2d", stat = c("mean", "sd"),ndraws=200,resp="walkpc3")
plot(loo(pc_mult1));loo(pc_mult1)
Computed from 12000 by 80 log-likelihood matrix
Estimate SE
elpd_loo -339.9 12.8
p_loo 17.4 2.9
looic 679.8 25.5
------
Monte Carlo SE of elpd_loo is 0.1.
Pareto k diagnostic values:
Count Pct. Min. n_eff
(-Inf, 0.5] (good) 79 98.8% 2395
(0.5, 0.7] (ok) 1 1.2% 210
(0.7, 1] (bad) 0 0.0% <NA>
(1, Inf) (very bad) 0 0.0% <NA>
All Pareto k estimates are ok (k < 0.7).
See help('pareto-k-diagnostic') for details.plot(loo(pc_mult2));loo(pc_mult2)
Computed from 12000 by 80 log-likelihood matrix
Estimate SE
elpd_loo -230.0 11.2
p_loo 22.9 3.8
looic 459.9 22.5
------
Monte Carlo SE of elpd_loo is NA.
Pareto k diagnostic values:
Count Pct. Min. n_eff
(-Inf, 0.5] (good) 79 98.8% 1655
(0.5, 0.7] (ok) 0 0.0% <NA>
(0.7, 1] (bad) 1 1.2% 68
(1, Inf) (very bad) 0 0.0% <NA>
See help('pareto-k-diagnostic') for details.### comparison of predictive ability via leave-one-out CV
pc_mult1 <- add_criterion(pc_mult1, criterion = "loo")
pc_mult2 <- add_criterion(pc_mult2, criterion = "loo")
loo_compare(pc_mult1,pc_mult2) elpd_diff se_diff
pc_mult2 0.0 0.0
pc_mult1 -109.9 10.3 Models are generally OK. Speed adjustment highly improves predictive accuracy of model.
4 Visualisation
4.1 PCA plot
## plotting multidimensional-------------
behav$group<-factor(behav$group,levels=c('ctrl.wt','eda.wt','ctrl.sca','eda.sca'))
par(mfrow=c(1,1))
par(mar=c(2.5,2.5,0,0))
par(mgp=c(1.3,0.6,0))
range<-c(-2.2,1.8);scal<-range[2]-range[1]
xrange<-c(-2.8,4.2);xscal=xrange[2]-xrange[1]
plot(NULL,xlim=c(xrange[1],xrange[2]),ylim=c(range[1],range[2]),
xlab="PC3 (explaining 12% of var.)",
ylab="PC2 (explaining 27% of var.)",
las=1, axes=FALSE,cex.lab=1)
rect(xrange[1],range[2],xrange[2],range[1],col="grey92",border=NA)
x<-range[1]
repeat{
lines(c(xrange[1],xrange[2]),c(x,x),col="white",lwd=0.7)
x=x+0.5;if(x>range[2]){break}}
x<-xrange[1]
repeat{
lines(c(x,x),c(range[1],range[2]),col="white",lwd=0.7)
x=x+0.5;if(x>xrange[2]){break}}
x=1
points(behav[behav$group=="ctrl.wt",]$walk_pc3,behav[behav$group=="ctrl.wt",]$walk_pc2,
pch=16,col=colc[x],cex=2.3);x=x+1
points(behav[behav$group=="eda.wt",]$walk_pc3,behav[behav$group=="eda.wt",]$walk_pc2,
pch=16,col=colc[x],cex=2.3);x=x+1
points(behav[behav$group=="ctrl.sca",]$walk_pc3,behav[behav$group=="ctrl.wt",]$walk_pc2,
pch=16,col=colc[x],cex=2.3);x=x+1
points(behav[behav$group=="eda.sca",]$walk_pc3,behav[behav$group=="eda.wt",]$walk_pc2,
pch=16,col=colc[x],cex=2.3);x=x+1
tckk=-0.02
axis(2,las=2,cex.axis=1,at=seq(range[1],range[2],by=0.5),
labels=c(rep("",length(seq(range[1],range[2],by=0.5)))),pos=xrange[1],tck=tckk)
axis(2,las=2,cex.axis=1,at=seq(range[1],range[2],by=1),pos=xrange[1],tck=tckk)
axis(side=1,las=1,cex.axis=1.1,at=seq(xrange[1],xrange[2],by=0.5)
,pos=range[1],tck=tckk)
pos=range[1]+scal*0.27
xx=1;xp=1.7;lwd=1.6
pos1=xrange[1]+xscal*0.72
pos3=xrange[1]+xscal*0.97
rect(pos1,pos+0.03*scal,pos3,range[1]+0.035*scal,
col="white",border="grey50",lwd=0.8)
ypos=0.8;xpos=xrange[1]+1.01*xscal;cons=0.065
text(xpos,range[1]+scal*ypos,expression(P[Wt_0 : SCA1_0]~ '= 0.027'),pos=2);ypos=ypos-cons
text(xpos,range[1]+scal*ypos,expression(P[SCA1_0 : SCA1_eda]~ '= 0.2'),pos=2);ypos=ypos-cons
text(xpos,range[1]+scal*ypos,expression(P[Wt_0 : Wt_eda]~ '= 0.4'),pos=2)
## legend
pos=range[1]+scal*0.27
xx=1;xp=1.7;lwd=1.6
pos1=xrange[1]+xscal*0.72
pos3=xrange[1]+xscal*0.97
rect(pos1,pos+0.03*scal,pos3,range[1]+0.035*scal,
col="white",border="grey50",lwd=0.8)
points(2.4,-1.15,pch=16,col=cola[xx],cex=2)
text(3.3,-1.15,"Wt_0",pch=16,col=cola[xx],cex=1.1);xx=xx+1
points(2.4,-1.4,pch=16,col=cola[xx],cex=2)
text(3.3,-1.4,"Wt_eda",pch=16,col=cola[xx],cex=1.1);xx=xx+1
points(2.4,-1.65,pch=16,col=cola[xx],cex=2)
text(3.3,-1.65,"SCA1_0",pch=16,col=cola[xx],cex=1.1);xx=xx+1
points(2.4,-1.9,pch=16,col=cola[xx],cex=2)
text(3.3,-1.9,"SCA1_eda",pch=16,col=cola[xx],cex=1.1);xx=xx+1
4.2 Posterior draws extraction
## Posterior extraction PC1
post_fix_na_pc1<-as.data.frame(pc_mult1, variable = c("b_walkpc1_Intercept","b_walkpc1_groupctrl.wt",
"b_walkpc1_groupeda.wt","b_walkpc1_groupeda.sca"))
names(post_fix_na_pc1)[1]<-"sca_ctrl"
post_fix_na_pc1$wt_ctrl<-post_fix_na_pc1$sca_ctrl+post_fix_na_pc1$b_walkpc1_groupctrl.wt
post_fix_na_pc1$wt_eda<-post_fix_na_pc1$sca_ctrl+post_fix_na_pc1$b_walkpc1_groupeda.wt
post_fix_na_pc1$sca_eda<-post_fix_na_pc1$sca_ctrl+post_fix_na_pc1$b_walkpc1_groupeda.sca
post_fix_na_pc1$wt_contrast<-post_fix_na_pc1$wt_eda-post_fix_na_pc1$wt_ctrl
post_fix_adj_pc1<-as.data.frame(pc_mult2, variable = c("b_walkpc1_Intercept","b_walkpc1_groupctrl.wt",
"b_walkpc1_groupeda.wt","b_walkpc1_groupeda.sca"))
names(post_fix_adj_pc1)[1]<-"sca_ctrl"
post_fix_adj_pc1$wt_ctrl<-post_fix_adj_pc1$sca_ctrl+post_fix_adj_pc1$b_walkpc1_groupctrl.wt
post_fix_adj_pc1$wt_eda<-post_fix_adj_pc1$sca_ctrl+post_fix_adj_pc1$b_walkpc1_groupeda.wt
post_fix_adj_pc1$sca_eda<-post_fix_adj_pc1$sca_ctrl+post_fix_adj_pc1$b_walkpc1_groupeda.sca
post_fix_adj_pc1$wt_contrast<-post_fix_adj_pc1$wt_eda-post_fix_adj_pc1$wt_ctrl
## Posterior extraction PC2
post_fix_na_pc2<-as.data.frame(pc_mult1, variable = c("b_walkpc2_Intercept","b_walkpc2_groupctrl.wt",
"b_walkpc2_groupeda.wt","b_walkpc2_groupeda.sca"))
names(post_fix_na_pc2)[1]<-"sca_ctrl"
post_fix_na_pc2$wt_ctrl<-post_fix_na_pc2$sca_ctrl+post_fix_na_pc2$b_walkpc2_groupctrl.wt
post_fix_na_pc2$wt_eda<-post_fix_na_pc2$sca_ctrl+post_fix_na_pc2$b_walkpc2_groupeda.wt
post_fix_na_pc2$sca_eda<-post_fix_na_pc2$sca_ctrl+post_fix_na_pc2$b_walkpc2_groupeda.sca
post_fix_na_pc2$wt_contrast<-post_fix_na_pc2$wt_eda-post_fix_na_pc2$wt_ctrl
post_fix_adj_pc2<-as.data.frame(pc_mult2, variable = c("b_walkpc2_Intercept","b_walkpc2_groupctrl.wt",
"b_walkpc2_groupeda.wt","b_walkpc2_groupeda.sca"))
names(post_fix_adj_pc2)[1]<-"sca_ctrl"
post_fix_adj_pc2$wt_ctrl<-post_fix_adj_pc2$sca_ctrl+post_fix_adj_pc2$b_walkpc2_groupctrl.wt
post_fix_adj_pc2$wt_eda<-post_fix_adj_pc2$sca_ctrl+post_fix_adj_pc2$b_walkpc2_groupeda.wt
post_fix_adj_pc2$sca_eda<-post_fix_adj_pc2$sca_ctrl+post_fix_adj_pc2$b_walkpc2_groupeda.sca
post_fix_adj_pc2$wt_contrast<-post_fix_adj_pc2$wt_eda-post_fix_adj_pc2$wt_ctrl
## Posterior extraction PC3
post_fix_na_pc3<-as.data.frame(pc_mult1, variable = c("b_walkpc3_Intercept","b_walkpc3_groupctrl.wt",
"b_walkpc3_groupeda.wt","b_walkpc3_groupeda.sca"))
names(post_fix_na_pc3)[1]<-"sca_ctrl"
post_fix_na_pc3$wt_ctrl<-post_fix_na_pc3$sca_ctrl+post_fix_na_pc3$b_walkpc3_groupctrl.wt
post_fix_na_pc3$wt_eda<-post_fix_na_pc3$sca_ctrl+post_fix_na_pc3$b_walkpc3_groupeda.wt
post_fix_na_pc3$sca_eda<-post_fix_na_pc3$sca_ctrl+post_fix_na_pc3$b_walkpc3_groupeda.sca
post_fix_na_pc3$wt_contrast<-post_fix_na_pc3$wt_eda-post_fix_na_pc3$wt_ctrl
post_fix_adj_pc3<-as.data.frame(pc_mult2, variable = c("b_walkpc3_Intercept","b_walkpc3_groupctrl.wt",
"b_walkpc3_groupeda.wt","b_walkpc3_groupeda.sca"))
names(post_fix_adj_pc3)[1]<-"sca_ctrl"
post_fix_adj_pc3$wt_ctrl<-post_fix_adj_pc3$sca_ctrl+post_fix_adj_pc3$b_walkpc3_groupctrl.wt
post_fix_adj_pc3$wt_eda<-post_fix_adj_pc3$sca_ctrl+post_fix_adj_pc3$b_walkpc3_groupeda.wt
post_fix_adj_pc3$sca_eda<-post_fix_adj_pc3$sca_ctrl+post_fix_adj_pc3$b_walkpc3_groupeda.sca
post_fix_adj_pc3$wt_contrast<-post_fix_adj_pc3$wt_eda-post_fix_adj_pc3$wt_ctrl4.3 Showing posterior probabilities
## plot setting and labeling
behav$group<-factor(behav$group,
levels=c('ctrl.wt','eda.wt','ctrl.sca','eda.sca'))
m <- matrix(c(9,1,2,3,
4,6,7,8,
5,6,7,8), nrow = 3, ncol=4 ,byrow = TRUE)
layout(mat = m,heights = c(0.06,0.94/2,0.94/2),
widths = c(0.07,0.31,0.31,0.31))
par(mgp=c(1.6,0.58,0))
par(mar=c(0,0,0,0))
plot(NULL,xlim=c(-1,1),ylim=c(-1,1),xlab="",ylab="",axes=F)
text(0,0.5,"PC1",xpd=T,cex=1.8)
text(0,-0.65,"fast walk & swings, long strides",xpd=T,cex=1.6,font=3)
plot(NULL,xlim=c(-1,1),ylim=c(-1,1),xlab="",ylab="",axes=F)
text(0.0,0.5,"PC2",xpd=T,cex=1.8)
text(0.0,-0.65,"inconsistent walk",xpd=T,cex=1.5,font=3)
plot(NULL,xlim=c(-1,1),ylim=c(-1,1),xlab="",ylab="",axes=F)
text(0.0,0.5,"PC3",xpd=T,cex=1.8)
text(0.0,-0.65,"support with 0 or 1 legs",xpd=T,cex=1.6,font=3)
plot(NULL,xlim=c(-1,1),ylim=c(-1,1),xlab="",ylab="",axes=F)
text(-0.5,0,"Original PC values",xpd=T,cex=1.5,srt=90)
lines(c(0,0),c(-1,1))
lines(c(0,0.2),c(-1,-1))
lines(c(0,0.2),c(1,1))
plot(NULL,xlim=c(-1,1),ylim=c(-1,1),xlab="",ylab="",axes=F)
text(-0.5,0,"Adjusted for speed",xpd=T,cex=1.5,srt=90)
lines(c(0,0),c(-1,1))
lines(c(0,0.2),c(-1,-1))
lines(c(0,0.2),c(1,1))
## posterior probability distributions
par(mar=c(2,0,0,0))
dif<-data.frame(post_fix_adj_pc1$wt_contrast,post_fix_adj_pc1$b_walkpc1_groupeda.sca,
post_fix_adj_pc1$b_walkpc1_groupctrl.wt,
post_fix_na_pc1$wt_contrast,post_fix_na_pc1$b_walkpc1_groupeda.sca,
post_fix_na_pc1$b_walkpc1_groupctrl.wt)
tckk=-0.018
ste<-seq(-2,2,by=0.5)
yla = ''
mons_poste6(dif, 0)
text(0,-0.06,expression(beta[PC1]),xpd=T,cex=1.7)
zpos=seq(0,1,1/6)
xx=1;ind=0.3;xpol<- -0.6
text(xpol,zpos[xx]+zpos[2]*ind,
"Wt: ed x 0 ",cex=1.3,xpd=T,font=2)
xx=xx+1
text(xpol,zpos[xx]+zpos[2]*ind,
"SCA1: ed x 0 ",cex=1.3,font=2)
xx=xx+1
text(xpol,zpos[xx]+zpos[2]*ind,
"0: Wt x SCA1",cex=1.3,font=2)
xx=xx+1
text(xpol,zpos[xx]+zpos[2]*ind,
"Wt: ed x 0 ",cex=1.3,xpd=T,font=2)
xx=xx+1
text(xpol,zpos[xx]+zpos[2]*ind,
"SCA1: ed x 0 ",cex=1.3,font=2)
xx=xx+1
text(xpol,zpos[xx]+zpos[2]*ind,
"0: Wt x SCA1",cex=1.3,font=2)
dif<-data.frame(post_fix_adj_pc2$wt_contrast,post_fix_adj_pc2$b_walkpc2_groupeda.sca,
post_fix_adj_pc2$b_walkpc2_groupctrl.wt,
post_fix_na_pc2$wt_contrast,post_fix_na_pc2$b_walkpc2_groupeda.sca,
post_fix_na_pc2$b_walkpc2_groupctrl.wt)
tckk=-0.018
ste<-seq(-1.5,1.5,by=0.5)
mons_poste6(dif,0)
text(0,-0.06,expression(beta[PC2]),xpd=T,cex=1.7, xpd = TRUE)
dif<-data.frame(post_fix_adj_pc3$wt_contrast,post_fix_adj_pc3$b_walkpc3_groupeda.sca,
post_fix_adj_pc3$b_walkpc3_groupctrl.wt,
post_fix_na_pc3$wt_contrast,post_fix_na_pc3$b_walkpc3_groupeda.sca,
post_fix_na_pc3$b_walkpc3_groupctrl.wt)
tckk=-0.018
ste<-seq(-1.5,1.5,by=0.5)
mons_poste6(dif,0)
text(0,-0.06,expression(beta[PC3]),xpd=T,cex=1.7, xpd = TRUE)