R/jzs_cor.R In MicheleNuijten/BayesMed: Default Bayesian Hypothesis Tests for Correlation, Partial Correlation, and Mediation

Documented in jzs_cor

```jzs_cor <-
function(V1,V2,
alternative=c("two.sided","less","greater"),
n.iter=10000,n.burnin=500,standardize=TRUE){

runif(1) # defines .Random.seed

# standardize variables
if(standardize==TRUE){
X <- (V1-mean(V1))/sd(V1)
Y <- (V2-mean(V2))/sd(V2)
}else {
X <- V1
Y <- V2
}

r <- cor(X,Y)
n <- length(X)

# main function to analytically calculate the BF for correlation
# see Wetzels, R. & Wagenmakers, E.-J. (2012). A default Bayesian hypothesis test for correlations and partial correlations. Psychonomic Bulletin & Review, 19, 1057-1064
# the jzs_corbf function is based on updated R code that can handle larger values of n

jzs_corbf=function(r,n){
int=function(r,n,g){
exp(
((n-2)/2)*log(1+g)+
(-(n-1)/2)*log(1+(1-r^2)*g)+
(-3/2)*log(g)+
-n/(2*g))
}

bf10 <- try(sqrt((n/2))/gamma(1/2)*
integrate(int,lower=0,upper=Inf,r=r,n=n)\$value)

# if the evidence is overwhelming, the BF will become infinite
# to avoid this resulting in an error, give back the max possible number
if(class(bf10)=="try-error" & grepl("non-finite function value",bf10[1])){
bf10 <- .Machine\$double.xmax
message("Note: the Error above was caused because the BF from the function jzs_corbf was approaching infinity. To avoid the function from crashing, the returned BF is the largest possible number in R.")
} else {
if(class(bf10)=="try-error" & !grepl("non-finite function value",bf10[1])){
bf10 <- NA
message("An error occurred. The BF could not be calculated")
}
}
return(bf10)
}

BF <- jzs_corbf(r,n)

### the next part is needed to impose an order restriction
### for the order restrictions we need to estimate the posterior samples

#==========================================================
#==========================================================

jagsmodelcorrelation <-

"####### Cauchy-prior on single regression coefficient #######
model

{

for (i in 1:n)

{
mu[i] <- intercept + alpha*x[i]
y[i]   ~ dnorm(mu[i],phi)

}

# uninformative prior on the intercept intercept,
# Jeffreys' prior on precision phi
intercept ~ dnorm(0,.0001)
phi   ~ dgamma(.0001,.0001)
#phi   ~ dgamma(0.0000001,0.0000001) #JAGS accepts even this
#phi   ~ dgamma(0.01,0.01)           #WinBUGS wants this

# inverse-gamma prior on g:
g       <- 1/invg
a.gamma <- 1/2
b.gamma <- n/2
invg     ~ dgamma(a.gamma,b.gamma)

# g-prior on beta:
vari <- (g/phi) * invSigma
prec <- 1/vari
alpha    ~ dnorm(0, prec)
}

# Explanation------------------------------------------------------------------
# Prior on g:
# We know that g ~ inverse_gamma(1/2, n/2), with 1/2 the shape
# parameter and n/2 the scale parameter.
# It follows that 1/g ~ gamma(1/2, 2/n).
# However, BUGS/JAGS uses the *rate parameterization* 1/theta instead of the
# scale parametrization theta. Hence we obtain, in de BUGS/JAGS rate notation:
# 1/g ~ dgamma(1/2, n/2)
#------------------------------------------------------------------------------
"

jags.model.file1 <- tempfile(fileext=".txt")
write(jagsmodelcorrelation,jags.model.file1)

#========================================================================
# Estimate Posterior Distribution for the Correlation Coefficient Alpha
#========================================================================

x <- X
y <- Y

invSigma <- solve(t(x)%*%x)

jags.data   <- list("n", "x", "y", "invSigma")
jags.params <- c("alpha", "g")
jags.inits  <-  list(
list(alpha = 0.0), #chain 1 starting value
list(alpha = -0.3), #chain 2 starting value
list(alpha = 0.3)) #chain 3 starting value

jagssamples <- jags(data=jags.data, inits=jags.inits, jags.params,
n.chains=3, n.iter=n.iter, DIC=T,
n.burnin=n.burnin, n.thin=1, model.file=jags.model.file1)

# estimate the posterior regression coefficient and scaling factor g
alpha <- jagssamples\$BUGSoutput\$sims.list\$alpha[,1]
g  <- jagssamples\$BUGSoutput\$sims.list\$g

#-------------------------------------------------------

# one-sided test?

# save BF for one-tailed test
# BF21 = 2*{proportion posterior samples of alpha < 0}

propposterior_less <- sum(alpha<0)/length(alpha)
propposterior_greater <- sum(alpha>0)/length(alpha)

# posterior proportion cannot be zero, because this renders a BF of zero
# none of the samples of the parameter follow the restriction
# ergo: the posterior proportion is smaller than 1/length(parameter)

if(propposterior_less==0){
propposterior_less <- 1/length(alpha)
}

if(propposterior_greater==0){
propposterior_greater <- 1/length(alpha)
}

BF21_less <- 2*propposterior_less
BF21_greater <- 2*propposterior_greater

if(alternative[1]=="less"){
# BF10 = p(D|a~cauchy(0,1))/p(D|a=0)
BF10 <- BF

# BF21 = p(D|a~cauchy-(0,1))/p(D|a~cauchy(0,1))
# BF21 = 2*{proportion posterior samples of alpha < 0}
BF21 <- BF21_less

BF <- BF10*BF21

} else if(alternative[1]=="greater"){
# BF10 = p(D|a~cauchy(0,1))/p(D|a=0)
BF10 <- BF

# BF21 = p(D|a~cauchy+(0,1))/p(D|a~cauchy(0,1))
# BF21 = 2*{proportion posterior samples of alpha > 0}
BF21 <- BF21_greater

BF <- BF10*BF21

}

#--------------------------------------------------------

# convert BFs to posterior probability
# prob cannot be exactly 1 or 0
prob_r <- BF/(BF+1)

if(prob_r == 1){
prob_r <- prob_r - .Machine\$double.eps
}
if(prob_r == 0){
prob_r <- prob_r + .Machine\$double.eps
}

#==================================================

# convert posterior samples for the regression coefficient x-y to correlation
cor_coef <- alpha*(sd(x)/sd(y))

#===================================================

res <- list(Correlation=mean(cor_coef),
BayesFactor=BF,
PosteriorProbability=prob_r,
cor_coef_samples=cor_coef,
jagssamples=jagssamples)

class(res) <- c("jzs_med","list")
class(res\$jagssamples) <- "rjags"
class(res\$cor_coef_samples) <- "CI"

return(res)
}
```
MicheleNuijten/BayesMed documentation built on May 7, 2019, 4:56 p.m.