# installing and loading packages install.packages("MCMCglmm"); library(MCMCglmm) install.packages("BayesFactor"); library(BayesFactor) # generating scores for 3 variables for the Correlation and Regression analyses set.seed(123) # fix the random number generator seed in order to produce the same data on every run N = 200 X = rnorm(N, mean = 0, sd = 1) Y = rnorm(N, mean = 0, sd = 1) Z = rnorm(N, mean = 0, sd = 1) pc.cr = princomp( cbind(X,Y,Z) ) fscores = pc.cr$scores # imposing correlations between the factor scores hyp = ' 1 .4 .2 .4 1 .3 .2 .3 1' rtarget = data.matrix( read.table(text=hyp)) dataset1 = fscores %*% chol(rtarget) dataset1 = data.frame(scale(dataset1)) # standardizing colnames(dataset1) = c('varA','varB','varDV') round(cor(dataset1),2) # conventional Pearson correlations # Correlation # conventional Pearson correlation & NHST for varB & varDV cor.test(dataset1$varB,dataset1$varDV) # correlation through MCMCglmm # the variables are already standardized model1 = MCMCglmm(varDV ~ varB, dat=dataset1, nitt=13000, burnin=3000, thin=10) summary(model1) autocorr(model1$VCV) plot(model1$Sol) # REGRESSION # conventional multiple linear regression with NHST fit = lm(varDV ~ varA + varB, data=dataset1) summary(fit) confint(fit, level=0.95) anova(fit) # Bayesian REGRESSION using MCMCglmm model2 = MCMCglmm(varDV ~ varA + varB, data=dataset1, family='gaussian', nitt=13000, burnin=3000, thin=10) summary(model2) autocorr(model2$VCV) plot(model2$Sol) # Bayesian REGRESSION using MCMCglmm -- with different priors # default priors model2 = MCMCglmm(varDV ~ varA + varB, data=dataset1, family='gaussian', nitt=13000, burnin=3000, thin=10); summary(model2) plotdat <- model2$Sol[,3] # for MCMCglmm priors, B is for fixed effects. B is a list containing the expected value (mu) # and a (co)variance matrix (V) representing the strength of belief: the defaults are B$mu=0 # and B$V=I*1e+10, where where I is an identity matrix of appropriate dimension # mu values of 0 for the intercept, and 0 for the regression weights prior.list <- list(B = list(mu = c(0, 0, 0), V = diag(c(1e+10, 1e+10, 1e+1)))) model2 = MCMCglmm(varDV ~ varA + varB, data=dataset1, prior = prior.list, family='gaussian', nitt=13000, burnin=3000, thin=10, verbose = F); summary(model2) plotdat <- cbind( plotdat, model2$Sol[,3]) # mu values of 0 for the intercept, and -.5 for the regression weights prior.list <- list(B = list(mu = c(0, -.5, -.5), V = diag(c(1e+10, 1e+10, 1e+1)))) model2 = MCMCglmm(varDV ~ varA + varB, data=dataset1, prior = prior.list, family='gaussian', nitt=13000, burnin=3000, thin=10, verbose = F); summary(model2) plotdat <- cbind( plotdat, model2$Sol[,3]) # mu values of 0 for the intercept, and .5 for the regression weights prior.list <- list(B = list(mu = c(0, .5, .5), V = diag(c(1e+10, 1e+10, 1e+1)))) model2 = MCMCglmm(varDV ~ varA + varB, data=dataset1, prior = prior.list, family='gaussian', nitt=13000, burnin=3000, thin=10, verbose = F); summary(model2) plotdat <- cbind( plotdat, model2$Sol[,3]) # mu values of 0 for the intercept, and -1 for the regression weights prior.list <- list(B = list(mu = c(0, -1, -1), V = diag(c(1e+10, 1e+10, 1e+1)))) model2 = MCMCglmm(varDV ~ varA + varB, data=dataset1, prior = prior.list, family='gaussian', nitt=13000, burnin=3000, thin=10, verbose = F); summary(model2) plotdat <- cbind( plotdat, model2$Sol[,3]) # mu values of 0 for the intercept, and 1 for the regression weights prior.list <- list(B = list(mu = c(0, 1, 1), V = diag(c(1e+10, 1e+10, 1e+1)))) model2 = MCMCglmm(varDV ~ varA + varB, data=dataset1, prior = prior.list, family='gaussian', nitt=13000, burnin=3000, thin=10, verbose = F); summary(model2) plotdat <- cbind( plotdat, model2$Sol[,3]) plot (density(plotdat[,1]), yaxt="n", main="", ylim=c(0,6), ylab='Density', font.lab=1.8, cex.lab=1.2 ) lines (density(plotdat[,2])) lines (density(plotdat[,3])) lines (density(plotdat[,4])) lines (density(plotdat[,5])) lines (density(plotdat[,6])) title(main='Density Plots of Posteriors for varB Produced Using Differing Priors') axis(2, las=2, at=seq(0, 6, 1), tick = TRUE) # ANOVA # generating data for oneway ANOVA, 3 groups set.seed(123) # fix the random number generator seed in order to produce the same data on every run X = rnorm(20,40,7) Y = rnorm(20,60,7) Z = rnorm(20,43,7) dataset2 = data.frame(c(X,Y,Z)); colnames(dataset2) = c('DV') dataset2$Group = gl(3,20, labels = c("Group1", "Group2", "Group3")) # descriptives & plot tapply(dataset2$DV,dataset2$Group, summary) # display the group means tapply(dataset2$DV,dataset2$Group, sd) # display the group SDs plot(DV ~ Group, data=dataset2) # conventional ANOVA & multiple comparisons results = aov(DV ~ Group, data=dataset2) summary(results) pairwise.t.test(dataset2$DV,dataset2$Group, p.adjust='bonferroni') TukeyHSD(results,conf.level=.95) # Bayesian ANOVA using MCMCglmm model3 = MCMCglmm(DV ~ Group, data=dataset2, verbose=F, nitt=13000, thin=10, burnin=3000) summary(model3) autocorr(model3$VCV) plot(model3$Sol) # Bayes factors anovaBF(DV ~ Group, iterations = 10000, data=dataset2) # Bayes factor for the model anovaBF(DV ~ Group, iterations = 10000, data=subset(dataset2, Group != 'Group3') ) # Bayes factor for Groups 1 vs 2 anovaBF(DV ~ Group, iterations = 10000, data=subset(dataset2, Group != 'Group2') ) # Bayes factor for Groups 1 vs 3 anovaBF(DV ~ Group, iterations = 10000, data=subset(dataset2, Group != 'Group1') ) # Bayes factor for Groups 2 vs 3 > > # generating scores for 3 variables for the Correlation and Regression analyses > set.seed(123) # fix the random number generator seed in order to produce the same data on every run > N = 200 > X = rnorm(N, mean = 0, sd = 1) > Y = rnorm(N, mean = 0, sd = 1) > Z = rnorm(N, mean = 0, sd = 1) > pc.cr = princomp( cbind(X,Y,Z) ) > fscores = pc.cr$scores > # imposing correlations between the factor scores > hyp = ' + 1 .4 .2 + .4 1 .3 + .2 .3 1' > rtarget = data.matrix( read.table(text=hyp)) > dataset1 = fscores %*% chol(rtarget) > dataset1 = data.frame(scale(dataset1)) # standardizing > colnames(dataset1) = c('varA','varB','varDV') > > round(cor(dataset1),2) # conventional Pearson correlations varA varB varDV varA 1.00 0.42 0.22 varB 0.42 1.00 0.32 varDV 0.22 0.32 1.00 > > > > > # Correlation > > # conventional Pearson correlation & NHST for varB & varDV > cor.test(dataset1$varB,dataset1$varDV) Pearson's product-moment correlation data: dataset1$varB and dataset1$varDV t = 4.6897, df = 198, p-value = 5.091e-06 alternative hypothesis: true correlation is not equal to 0 95 percent confidence interval: 0.1855820 0.4358053 sample estimates: cor 0.316182 > > # correlation through MCMCglmm > # the variables are already standardized > model1 = MCMCglmm(varDV ~ varB, dat=dataset1, nitt=13000, burnin=3000, thin=10) MCMC iteration = 0 MCMC iteration = 1000 MCMC iteration = 2000 MCMC iteration = 3000 MCMC iteration = 4000 MCMC iteration = 5000 MCMC iteration = 6000 MCMC iteration = 7000 MCMC iteration = 8000 MCMC iteration = 9000 MCMC iteration = 10000 MCMC iteration = 11000 MCMC iteration = 12000 MCMC iteration = 13000 > summary(model1) Iterations = 3001:12991 Thinning interval = 10 Sample size = 1000 DIC: 551.5148 R-structure: ~units post.mean l-95% CI u-95% CI eff.samp units 0.9134 0.7483 1.112 1000 Location effects: varDV ~ varB post.mean l-95% CI u-95% CI eff.samp pMCMC (Intercept) 0.0007881 -0.1459576 0.1214560 720.4 0.942 varB 0.3126429 0.1746920 0.4318655 1127.6 <0.001 *** --- Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1 > autocorr(model1$VCV) , , units units Lag 0 1.00000000 Lag 10 0.02728845 Lag 50 0.01037955 Lag 100 0.06381656 Lag 500 -0.03396188 > plot(model1$Sol) > > > > > > > > # REGRESSION > > # conventional multiple linear regression with NHST > fit = lm(varDV ~ varA + varB, data=dataset1) > summary(fit) Call: lm(formula = varDV ~ varA + varB, data = dataset1) Residuals: Min 1Q Median 3Q Max -2.3847 -0.6331 -0.0599 0.5968 3.6897 Coefficients: Estimate Std. Error t value Pr(>|t|) (Intercept) 3.783e-18 6.709e-02 0.000 1.000000 varA 1.032e-01 7.393e-02 1.395 0.164487 varB 2.734e-01 7.393e-02 3.698 0.000282 *** --- Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1 Residual standard error: 0.9488 on 197 degrees of freedom Multiple R-squared: 0.1088, Adjusted R-squared: 0.09973 F-statistic: 12.02 on 2 and 197 DF, p-value: 1.185e-05 > confint(fit, level=0.95) 2.5 % 97.5 % (Intercept) -0.13231085 0.1323109 varA -0.04263971 0.2489561 varB 0.12756225 0.4191581 > anova(fit) Analysis of Variance Table Response: varDV Df Sum Sq Mean Sq F value Pr(>F) varA 1 9.339 9.3390 10.373 0.0014955 ** varB 1 12.308 12.3080 13.671 0.0002822 *** Residuals 197 177.353 0.9003 --- Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1 > > > # Bayesian REGRESSION using MCMCglmm > model2 = MCMCglmm(varDV ~ varA + varB, data=dataset1, family='gaussian', nitt=13000, burnin=3000, thin=10) MCMC iteration = 0 MCMC iteration = 1000 MCMC iteration = 2000 MCMC iteration = 3000 MCMC iteration = 4000 MCMC iteration = 5000 MCMC iteration = 6000 MCMC iteration = 7000 MCMC iteration = 8000 MCMC iteration = 9000 MCMC iteration = 10000 MCMC iteration = 11000 MCMC iteration = 12000 MCMC iteration = 13000 > summary(model2) Iterations = 3001:12991 Thinning interval = 10 Sample size = 1000 DIC: 551.5167 R-structure: ~units post.mean l-95% CI u-95% CI eff.samp units 0.9105 0.735 1.09 1000 Location effects: varDV ~ varA + varB post.mean l-95% CI u-95% CI eff.samp pMCMC (Intercept) -0.001926 -0.124576 0.128453 1107 0.994 varA 0.098120 -0.044695 0.247153 1038 0.198 varB 0.277865 0.127949 0.416467 1000 <0.001 *** --- Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1 > autocorr(model2$VCV) , , units units Lag 0 1.00000000 Lag 10 0.01371913 Lag 50 -0.06528034 Lag 100 0.01686856 Lag 500 0.02007962 > plot(model2$Sol) > > > > > > # Bayesian REGRESSION using MCMCglmm -- with different priors > > # default priors > model2 = MCMCglmm(varDV ~ varA + varB, data=dataset1, family='gaussian', nitt=13000, burnin=3000, thin=10); summary(model2) MCMC iteration = 0 MCMC iteration = 1000 MCMC iteration = 2000 MCMC iteration = 3000 MCMC iteration = 4000 MCMC iteration = 5000 MCMC iteration = 6000 MCMC iteration = 7000 MCMC iteration = 8000 MCMC iteration = 9000 MCMC iteration = 10000 MCMC iteration = 11000 MCMC iteration = 12000 MCMC iteration = 13000 Iterations = 3001:12991 Thinning interval = 10 Sample size = 1000 DIC: 551.5837 R-structure: ~units post.mean l-95% CI u-95% CI eff.samp units 0.9106 0.7409 1.093 1105 Location effects: varDV ~ varA + varB post.mean l-95% CI u-95% CI eff.samp pMCMC (Intercept) -0.001176 -0.136675 0.131468 1000 0.992 varA 0.107716 -0.027386 0.261374 1038 0.156 varB 0.270507 0.115525 0.408103 1000 0.002 ** --- Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1 > > > plotdat <- model2$Sol[,3] > > > # for MCMCglmm priors, B is for fixed effects. B is a list containing the expected value (mu) > # and a (co)variance matrix (V) representing the strength of belief: the defaults are B$mu=0 > # and B$V=I*1e+10, where where I is an identity matrix of appropriate dimension > > > # mu values of 0 for the intercept, and 0 for the regression weights > prior.list <- list(B = list(mu = c(0, 0, 0), V = diag(c(1e+10, 1e+10, 1e+1)))) > model2 = MCMCglmm(varDV ~ varA + varB, data=dataset1, prior = prior.list, family='gaussian', nitt=13000, burnin=3000, thin=10, verbose = F); summary(model2) Iterations = 3001:12991 Thinning interval = 10 Sample size = 1000 DIC: 551.5158 R-structure: ~units post.mean l-95% CI u-95% CI eff.samp units 0.9088 0.7496 1.129 1000 Location effects: varDV ~ varA + varB post.mean l-95% CI u-95% CI eff.samp pMCMC (Intercept) 0.0004871 -0.1224492 0.1356135 1018.0 0.998 varA 0.1009939 -0.0422584 0.2444908 1000.0 0.198 varB 0.2779698 0.1317281 0.4157104 901.3 <0.001 *** --- Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1 > plotdat <- cbind( plotdat, model2$Sol[,3]) > > # mu values of 0 for the intercept, and -.5 for the regression weights > prior.list <- list(B = list(mu = c(0, -.5, -.5), V = diag(c(1e+10, 1e+10, 1e+1)))) > model2 = MCMCglmm(varDV ~ varA + varB, data=dataset1, prior = prior.list, family='gaussian', nitt=13000, burnin=3000, thin=10, verbose = F); summary(model2) Iterations = 3001:12991 Thinning interval = 10 Sample size = 1000 DIC: 551.5169 R-structure: ~units post.mean l-95% CI u-95% CI eff.samp units 0.9099 0.732 1.086 855.4 Location effects: varDV ~ varA + varB post.mean l-95% CI u-95% CI eff.samp pMCMC (Intercept) -0.0008892 -0.1273231 0.1228500 1000.0 0.998 varA 0.1064190 -0.0406139 0.2458389 836.5 0.148 varB 0.2681316 0.1331627 0.4163012 1000.0 <0.001 *** --- Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1 > plotdat <- cbind( plotdat, model2$Sol[,3]) > > # mu values of 0 for the intercept, and .5 for the regression weights > prior.list <- list(B = list(mu = c(0, .5, .5), V = diag(c(1e+10, 1e+10, 1e+1)))) > model2 = MCMCglmm(varDV ~ varA + varB, data=dataset1, prior = prior.list, family='gaussian', nitt=13000, burnin=3000, thin=10, verbose = F); summary(model2) Iterations = 3001:12991 Thinning interval = 10 Sample size = 1000 DIC: 551.5282 R-structure: ~units post.mean l-95% CI u-95% CI eff.samp units 0.9094 0.7459 1.09 1107 Location effects: varDV ~ varA + varB post.mean l-95% CI u-95% CI eff.samp pMCMC (Intercept) 0.0002615 -0.1281780 0.1305094 1369 0.994 varA 0.1040694 -0.0464371 0.2361490 1000 0.144 varB 0.2728918 0.1394882 0.4133581 1102 <0.001 *** --- Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1 > plotdat <- cbind( plotdat, model2$Sol[,3]) > > # mu values of 0 for the intercept, and -1 for the regression weights > prior.list <- list(B = list(mu = c(0, -1, -1), V = diag(c(1e+10, 1e+10, 1e+1)))) > model2 = MCMCglmm(varDV ~ varA + varB, data=dataset1, prior = prior.list, family='gaussian', nitt=13000, burnin=3000, thin=10, verbose = F); summary(model2) Iterations = 3001:12991 Thinning interval = 10 Sample size = 1000 DIC: 551.4709 R-structure: ~units post.mean l-95% CI u-95% CI eff.samp units 0.9119 0.7412 1.09 1000 Location effects: varDV ~ varA + varB post.mean l-95% CI u-95% CI eff.samp pMCMC (Intercept) -0.0004289 -0.1298651 0.1347542 1000 0.966 varA 0.0998849 -0.0375466 0.2545398 1000 0.190 varB 0.2712774 0.1341460 0.4098331 1000 <0.001 *** --- Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1 > plotdat <- cbind( plotdat, model2$Sol[,3]) > > # mu values of 0 for the intercept, and 1 for the regression weights > prior.list <- list(B = list(mu = c(0, 1, 1), V = diag(c(1e+10, 1e+10, 1e+1)))) > model2 = MCMCglmm(varDV ~ varA + varB, data=dataset1, prior = prior.list, family='gaussian', nitt=13000, burnin=3000, thin=10, verbose = F); summary(model2) Iterations = 3001:12991 Thinning interval = 10 Sample size = 1000 DIC: 551.6378 R-structure: ~units post.mean l-95% CI u-95% CI eff.samp units 0.9093 0.7461 1.096 811.3 Location effects: varDV ~ varA + varB post.mean l-95% CI u-95% CI eff.samp pMCMC (Intercept) 0.003203 -0.134391 0.124581 1114 0.968 varA 0.104216 -0.048183 0.252562 1000 0.198 varB 0.269379 0.123461 0.413749 1000 <0.001 *** --- Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1 > plotdat <- cbind( plotdat, model2$Sol[,3]) > > > plot (density(plotdat[,1]), yaxt="n", main="", ylim=c(0,6), ylab='Density', font.lab=1.8, cex.lab=1.2 ) > lines (density(plotdat[,2])) > lines (density(plotdat[,3])) > lines (density(plotdat[,4])) > lines (density(plotdat[,5])) > lines (density(plotdat[,6])) > > title(main='Density Plots of Posteriors for varB Produced Using Differing Priors') > axis(2, las=2, at=seq(0, 6, 1), tick = TRUE) > > > > > > > > > # ANOVA > > # generating data for oneway ANOVA, 3 groups > set.seed(123) # fix the random number generator seed in order to produce the same data on every run > X = rnorm(20,40,7) > Y = rnorm(20,60,7) > Z = rnorm(20,43,7) > dataset2 = data.frame(c(X,Y,Z)); colnames(dataset2) = c('DV') > dataset2$Group = gl(3,20, labels = c("Group1", "Group2", "Group3")) > > > # descriptives & plot > tapply(dataset2$DV,dataset2$Group, summary) # display the group means $Group1 Min. 1st Qu. Median Mean 3rd Qu. Max. 26.23 36.55 40.84 40.99 43.84 52.51 $Group2 Min. 1st Qu. Median Mean 3rd Qu. Max. 48.19 55.44 59.02 59.64 65.05 68.78 $Group3 Min. 1st Qu. Median Mean 3rd Qu. Max. 32.16 40.07 42.75 43.75 47.43 58.18 > tapply(dataset2$DV,dataset2$Group, sd) # display the group SDs Group1 Group2 Group3 6.808657 5.809571 6.701384 > plot(DV ~ Group, data=dataset2) > > > # conventional ANOVA & multiple comparisons > results = aov(DV ~ Group, data=dataset2) > summary(results) Df Sum Sq Mean Sq F value Pr(>F) Group 2 4054 2026.9 48.64 4.74e-13 *** Residuals 57 2375 41.7 --- Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1 > pairwise.t.test(dataset2$DV,dataset2$Group, p.adjust='bonferroni') Pairwise comparisons using t tests with pooled SD data: dataset2$DV and dataset2$Group Group1 Group2 Group2 2.8e-12 - Group3 0.55 4.7e-10 P value adjustment method: bonferroni > TukeyHSD(results,conf.level=.95) Tukey multiple comparisons of means 95% family-wise confidence level Fit: aov(formula = DV ~ Group, data = dataset2) $Group diff lwr upr p adj Group2-Group1 18.64983 13.737410 23.562257 0.0000000 Group3-Group1 2.75403 -2.158394 7.666454 0.3743536 Group3-Group2 -15.89580 -20.808227 -10.983380 0.0000000 > > > # Bayesian ANOVA using MCMCglmm > model3 = MCMCglmm(DV ~ Group, data=dataset2, verbose=F, nitt=13000, thin=10, burnin=3000) > summary(model3) Iterations = 3001:12991 Thinning interval = 10 Sample size = 1000 DIC: 399.021 R-structure: ~units post.mean l-95% CI u-95% CI eff.samp units 43.26 29.75 60.02 1000 Location effects: DV ~ Group post.mean l-95% CI u-95% CI eff.samp pMCMC (Intercept) 41.010 38.213 44.114 1096 <0.001 *** GroupGroup2 18.564 14.230 22.670 1000 <0.001 *** GroupGroup3 2.750 -1.235 6.548 1000 0.17 --- Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1 > autocorr(model3$VCV) , , units units Lag 0 1.00000000 Lag 10 -0.03502126 Lag 50 -0.04229835 Lag 100 -0.04766084 Lag 500 -0.05137525 > plot(model3$Sol) > > > # Bayes factors > anovaBF(DV ~ Group, iterations = 10000, data=dataset2) # Bayes factor for the model |==================================================================================================================================================================================| 100% Bayes factor analysis -------------- [1] Group : 13017907395 ±0% Against denominator: Intercept only --- Bayes factor type: BFlinearModel, JZS > anovaBF(DV ~ Group, iterations = 10000, data=subset(dataset2, Group != 'Group3') ) # Bayes factor for Groups 1 vs 2 |==================================================================================================================================================================================| 100% Bayes factor analysis -------------- [1] Group : 242879306 ±0% Against denominator: Intercept only --- Bayes factor type: BFlinearModel, JZS > anovaBF(DV ~ Group, iterations = 10000, data=subset(dataset2, Group != 'Group2') ) # Bayes factor for Groups 1 vs 3 |==================================================================================================================================================================================| 100% Bayes factor analysis -------------- [1] Group : 0.5945621 ±0.01% Against denominator: Intercept only --- Bayes factor type: BFlinearModel, JZS > anovaBF(DV ~ Group, iterations = 10000, data=subset(dataset2, Group != 'Group1') ) # Bayes factor for Groups 2 vs 3 |==================================================================================================================================================================================| 100% Bayes factor analysis -------------- [1] Group : 6623975 ±0% Against denominator: Intercept only --- Bayes factor type: BFlinearModel, JZS > > >