Priors for regression

This document reproduces the plots from Rouder and Morey (2013) [pdf]. In the R code, the quantity \(\tau^2\) is refered to as r2, and the quantity \(\tau\) is referred to as r.

An interactive shiny app to demonstrate these priors can be found here.

The following functions define the density and cumulative density functions on \(\tau^2\), (analogous to \(R^2\)), given a prior scale.

dr2 = Vectorize(
  function(r2,scale) dcauchy(sqrt(r2/(1-r2)),scale=scale) * r2^(-.5)*(1-r2)^(-1.5),
  "r2"
  )

pr2 = Vectorize(
  function(r2,scale) integrate(dr2,.000001,r2,scale=scale)[[1]],
  "r2"
  )

The following functions define the density and cumulative density functions on \(\tau\), (analogous to \(\rho\), the true correlation), given a prior scale.

dr = Vectorize(
  function(r,scale) dr2(r^2,scale)*abs(r),
  "r"
  )
  

pr = Vectorize(
  function(r,scale) integrate(dr,-.999999,r,scale=scale)[[1]],
  "r"
  )

The variable otherscale gives the \(r\) scale to be compared to the default of 1. Change this to change the figures.

otherscale=.5

Sequences of \(x\) values.

r2r2 = seq(0,1,len=500)
rr = seq(-1,1,len=500)

Sometimes at the end points of the space, the density can be very high. These variables give the maximum values for the \(y\) axis in the plots for r2 and r respectively.

ymax1 = 7
ymax2 = 3.5

Plots from Rouder and Morey (2013)

par(mfrow=c(2,2),las=1,mgp=c(2,.7,0),mar=c(4,3,1.5,.1))

plot(r2r2,dr2(r2r2,otherscale),col="red",ty='l',lty=2,lwd=2,xlab=expression(paste("Prop. of nonerror variance, ",tau^2)),ylab="Prior density",ylim=c(0,ymax1),axes=FALSE)
lines(r2r2,dr2(r2r2,1),col="black",lty=1)
abline(h=0,col="gray")
abline(v=c(0,1),col="gray")
axis(1)
box()
mtext("A.",3,.1,adj=1,cex=1.2)
legend(.05,.98*ymax1,legend=c(1,otherscale),title="Cauchy prior scale",lty=c(1,2),
       col=c("black","red"),lwd=c(1,2))

plot(rr,dr(rr,otherscale),col="red",ty='l',lty=2,,lwd=2,xlab=expression(paste("Sqrt. prop. of nonerror variance, ",sqrt(tau^2))),ylab="",ylim=c(0,ymax2),axes=FALSE)
axis(1)
box()
lines(rr,dr(rr,1),col="black",lty=1)
abline(h=0,col="gray")
abline(v=c(-1,1),col="gray")
#abline(v=0,col="gray",lty=3)
mtext("C.",3,.1,adj=1,cex=1.2)


r2r2 = seq(0,1,len=100)
rr = seq(-1,1,len=100)


plot(r2r2,pr2(r2r2,otherscale),col="red",ty='l',lty=2,lwd=2,xlab=expression(paste("Prop. of nonerror variance, ",tau^2)),ylab="Cumulative prior probability",ylim=c(0,1),axes=FALSE)
lines(r2r2,pr2(r2r2,1),col="black",lty=1)
abline(h=c(0,1),col="gray")
abline(v=c(0,1),col="gray")
axis(1)
box()
axis(2,at=c(0,.5,1))
mtext("B.",3,.1,adj=1,cex=1.2)


plot(rr,pr(rr,otherscale),col="red",ty='l',lty=2,lwd=2,xlab=expression(paste("Sqrt. prop. of nonerror variance, ",sqrt(tau^2))),ylab="",ylim=c(0,1),axes=FALSE)
axis(1)
box()
lines(rr,pr(rr,1),col="black",lty=1)
abline(h=c(0,1),col="gray")
abline(v=c(-1,1),col="gray")
#abline(v=0,col="gray",lty=3)
mtext("D.",3,.1,adj=1,cex=1.2)
axis(2,at=c(0,.5,1))

Example showing the same thing, by sampling from the prior.

We can also show the prior distributions by sampling and transforming.

# How many simulation draws
M = 10000

# Draw from the prior
b = rcauchy(M,scale=otherscale)
# Transform to \tau^2
r2 = 1-1/(1+b^2) 
# Transform to \tau
r = sqrt(r2)*sample(c(-1,1),M,TRUE)

par(mfrow=c(1,2))

hist(r2, freq=FALSE,col="lightsteelblue1", ylab="Density", axes=FALSE, xlab=expression(tau^2))
axis(1)
lines(r2r2,dr2(r2r2,otherscale),col="red",lty=1,lwd=2)

abline(h=0,col="gray")
abline(v=c(0,1),col="gray")

hist(r, freq=FALSE,col="lightsteelblue1", ylab="Density", axes=FALSE, xlab=expression(sqrt(tau^2)))
axis(1)
lines(rr,dr(rr,otherscale),col="red",lty=1, lwd=2)

abline(h=0,col="gray")
abline(v=c(0,1),col="gray")