21 April 2016
Bayesian ANOVA
\[ \begin{aligned} y_i &= \mu + \pmb{\theta}^{\top} \pmb{x}_i + \varepsilon_i \\ \varepsilon_i &\sim N(0, \sigma) \end{aligned} \] Same as BayesFactor package, I used Zellner and Siow (1980) inspired g-priors: \[ \begin{aligned} \mu & \sim N(0, 100,000) \\ 1/\sigma^2 & \sim \Gamma(0.001, 0.001) \\ \theta_j &\sim N(0, g_j \sigma^2) \\ 1/g_j &\sim \chi^2(1) \end{aligned} \]
JAGS model (1/3)
# General Priors mu ~ dnorm(0,0.00001) tau ~ dgamma(0.001, 0.001) sigma2 <- 1/tau sigma <- sqrt(sigma2)
JAGS model (2/3)
# Theta coefficients:
for (i in 1:Nx){
theta[i] ~ dnorm(0, g_inv[i] * tau)
g[i] <- 1/g_inv[i]
g_inv[i] ~ dchisqr(1)
# Sample from priors:
theta_prior[i] ~ dnorm(0, g_inv[i] * tau)
}
JAGS model (3/3)
### Likelihood ###
for (p in 1:N){
mean_y[p] <- mu + X[p,] %*% theta
Y[p] ~ dnorm(mean_y[p], tau)
}
source("http://sachaepskamp.com/files/BayesianAnova.R")
Two functions:
BayesianAnova(formula, data,...)- Run the Bayesian ANOVA
- dots sent to JAGS
BF(results, ..., d=0.1)- Compute inequality and equality hypotheses
- Dots is any number of hyptoheses
- d is distance around equality hypotheses to test
Example 1: Tooth Growth
Data1 <- read.csv("http://sachaepskamp.com/files/ToothGrowth.csv")
Data1$dose <- as.factor(Data1$dose)
head(Data1)
## len supp dose ## 1 4.2 VC 500 ## 2 11.5 VC 500 ## 3 7.3 VC 500 ## 4 5.8 VC 500 ## 5 6.4 VC 500 ## 6 10.0 VC 500
library("BayesFactor")
anovaBF(len ~ supp * dose, data = Data1)
## Bayes factor analysis ## -------------- ## [1] supp : 1.198757 ±0.01% ## [2] dose : 4.983636e+12 ±0% ## [3] supp + dose : 2.76512e+14 ±1.09% ## [4] supp + dose + supp:dose : 8.090124e+14 ±3.53% ## ## Against denominator: ## Intercept only ## --- ## Bayes factor type: BFlinearModel, JZS
Res1 <- BayesianAnova(len ~ supp, data = Data1, n.iter = 10000,
n.chain = 3)
Res1
## parameter description mean sd ## 1 mu mean 20.457281 1.348873 ## 2 sigma2 variance 57.068776 10.913069 ## 3 theta suppVC -3.309619 1.875508 ## 4 mean[1] supp_OJ 20.457281 1.348873 ## 5 mean[2] supp_VC 17.147662 1.339423
Test if first mean is larger than second mean:
BF(Res1,mean[1] > mean[2])
## ## Bayes Factors computed on following hypotheses: ## mean[,1] > mean[,2] ## ## BF hypotheses vs unconstrained: 1.872 (1/BF: 0.534) ## BF hypotheses vs complement: 106.183 (1/BF: 0.009)
Test if means are equal:
BF(Res1,mean[1] == mean[2], d = 0.1)
## ## Bayes Factors computed on following hypotheses: ## mean[,1] > mean[,2] - d & ## mean[,1] < mean[,2] + d ## ## d = 0.1 ## ## BF hypotheses vs unconstrained: 0.694 (1/BF: 1.44) ## BF hypotheses vs complement: 0.709 (1/BF: 1.411)
BF(Res1,theta == 0, d = 0.1)
## ## Bayes Factors computed on following hypotheses: ## theta > 0 - d & ## theta < 0 + d ## ## d = 0.1 ## ## BF hypotheses vs unconstrained: 0.694 (1/BF: 1.44) ## BF hypotheses vs complement: 0.709 (1/BF: 1.411)
Res2 <- BayesianAnova(len ~ dose, data = Data1, n.iter = 10000,
n.chain = 3)
Res2
## parameter description mean sd ## 1 mu mean 10.907664 0.9894679 ## 2 sigma2 variance 19.359119 3.8733407 ## 3 theta[1] dose1000 8.628202 1.4129838 ## 4 theta[2] dose2000 15.072330 1.4001984 ## 5 mean[1] dose_500 10.907664 0.9894679 ## 6 mean[2] dose_1000 19.535866 0.9981120 ## 7 mean[3] dose_2000 25.979994 1.0021827
BF(Res2,theta[1] == 0, theta[2] == 0, d = 1)
## ## Bayes Factors computed on following hypotheses: ## theta[,1] > 0 - d & ## theta[,1] < 0 + d & ## theta[,2] > 0 - d & ## theta[,2] < 0 + d ## ## d = 1 ## ## BF hypotheses vs unconstrained: 0 (1/BF: Inf) ## BF hypotheses vs complement: 0 (1/BF: Inf)
Res3 <- BayesianAnova(len ~ supp*dose, data = Data1, n.iter = 10000,
n.chain = 3)
Res3
## parameter description mean sd ## 1 mu mean 13.1295467 1.1421550 ## 2 sigma2 variance 14.0681708 2.8114890 ## 3 theta[1] suppVC -4.5689279 1.5365381 ## 4 theta[2] dose1000 9.2232468 1.5549255 ## 5 theta[3] dose2000 13.0419921 1.6769859 ## 6 theta[4] suppVC:dose1000 -0.8963823 1.9510461 ## 7 theta[5] suppVC:dose2000 4.2131745 2.2205569 ## 8 mean[1] supp_OJ 20.5512930 0.6857014 ## 9 mean[2] supp_VC 17.0879625 0.6825468 ## 10 mean[3] dose_500 10.8450828 0.8435177 ## 11 mean[4] dose_1000 19.6201384 0.8407224 ## 12 mean[5] dose_2000 25.9936621 0.8316541 ## 13 mean[6] supp_OJ + dose_500 13.1295467 1.1421550 ## 14 mean[7] supp_OJ + dose_1000 22.3527935 1.1541648 ## 15 mean[8] supp_OJ + dose_2000 26.1715388 1.2120353 ## 16 mean[9] supp_VC + dose_500 8.5606188 1.1397372 ## 17 mean[10] supp_VC + dose_1000 16.8874833 1.1242269 ## 18 mean[11] supp_VC + dose_2000 25.8157854 1.1396110
BF(Res3, d = 2)
## ## Bayes Factors computed on following hypotheses: ## theta[,1] > 0 - d & ## theta[,1] < 0 + d & ## theta[,2] > 0 - d & ## theta[,2] < 0 + d & ## theta[,3] > 0 - d & ## theta[,3] < 0 + d & ## theta[,4] > 0 - d & ## theta[,4] < 0 + d & ## theta[,5] > 0 - d & ## theta[,5] < 0 + d ## ## d = 2 ## ## BF hypotheses vs unconstrained: 0 (1/BF: Inf) ## BF hypotheses vs complement: 0 (1/BF: Inf)
Example 2: Social Judgement
Data2 <- read.csv("http://sachaepskamp.com/files/OSF.csv")
head(Data2)
## aggressivenessScore mapCondition puzzleCondition ## 1 5.0 Local Neutral ## 2 8.0 Global Neutral ## 3 7.5 Global Neutral ## 4 6.5 Global Neutral ## 5 8.0 Local Aggression ## 6 8.0 Local Neutral
anovaBF(aggressivenessScore ~ mapCondition*puzzleCondition, data = Data2)
## Bayes factor analysis ## -------------- ## [1] mapCondition : 4.044688 ±0.02% ## [2] puzzleCondition : 0.4929054 ±0% ## [3] mapCondition + puzzleCondition : 2.198317 ±0.95% ## [4] mapCondition + puzzleCondition + mapCondition:puzzleCondition : 0.7951444 ±1.46% ## ## Against denominator: ## Intercept only ## --- ## Bayes factor type: BFlinearModel, JZS
Res4 <- BayesianAnova(aggressivenessScore ~ mapCondition, data = Data2,
n.iter = 10000, n.chain = 3)
Res4
## parameter description mean sd ## 1 mu mean 6.6680734 0.2983566 ## 2 sigma2 variance 1.9523843 0.3298221 ## 3 theta[1] mapConditionGlobal 1.0108898 0.4000795 ## 4 theta[2] mapConditionLocal 0.0957045 0.3908386 ## 5 mean[1] mapCondition_Control 6.6680734 0.2983566 ## 6 mean[2] mapCondition_Global 7.6789633 0.2658056 ## 7 mean[3] mapCondition_Local 6.7637779 0.2713693
BF(Res4, mean[2] > mean[1], mean[2] > mean[3])
## ## Bayes Factors computed on following hypotheses: ## mean[,2] > mean[,1] & ## mean[,2] > mean[,3] ## ## BF hypotheses vs unconstrained: 2.598 (1/BF: 0.385) ## BF hypotheses vs complement: 306.432 (1/BF: 0.003)
BF(Res4, mean[1] == mean[3], d = 0.1)
## ## Bayes Factors computed on following hypotheses: ## mean[,1] > mean[,3] - d & ## mean[,1] < mean[,3] + d ## ## d = 0.1 ## ## BF hypotheses vs unconstrained: 2.693 (1/BF: 0.371) ## BF hypotheses vs complement: 3.68 (1/BF: 0.272)
BF(Res4, d = 0.05)
## ## Bayes Factors computed on following hypotheses: ## theta[,1] > 0 - d & ## theta[,1] < 0 + d & ## theta[,2] > 0 - d & ## theta[,2] < 0 + d ## ## d = 0.05 ## ## BF hypotheses vs unconstrained: 0.333 (1/BF: 3) ## BF hypotheses vs complement: 0.334 (1/BF: 2.996)
Res5 <- BayesianAnova(aggressivenessScore ~ puzzleCondition, data = Data2,
n.iter = 10000, n.chain = 3)
Res5
## parameter description mean sd ## 1 mu mean 6.902633 0.2286833 ## 2 sigma2 variance 2.155039 0.3824632 ## 3 theta puzzleConditionNeutral 0.404074 0.3317627 ## 4 mean[1] puzzleCondition_Aggression 6.902633 0.2286833 ## 5 mean[2] puzzleCondition_Neutral 7.306707 0.2487358
BF(Res5, d = 0.1)
## ## Bayes Factors computed on following hypotheses: ## theta > 0 - d & ## theta < 0 + d ## ## d = 0.1 ## ## BF hypotheses vs unconstrained: 1.749 (1/BF: 0.572) ## BF hypotheses vs complement: 2.119 (1/BF: 0.472)
Conclusion
- My Bayesian ANOVA seems to work well when the evidence is overwelming
- Doesn't give results near BayesFactor when evidence is small
Codes: http://sachaepskamp.com/files/BayesianAnova.R