plot v

plot v

#You will need to load the R package "hypergeo" to use the vstat function

library(hypergeo)

#Below, I'm vectorizing the function so that I can plot curves. 
#The rest is unchanged from the vstat function by Stober-Davis & Dana.
#If you want to use R unbiased, remove the # before the Rsq adjustment calculation below
vstat <- Vectorize(function(n,p,Rsq)
{
  #Rsq = Re(1-((n-2)/(n-p))*(1-Rsq)*hypergeo(1,1,(n-p+2)*.5,1-Rsq))
  if (Rsq<=0) {Rsq = .0001}
  r = ((p-1)*(1-Rsq))/((n-p)*Rsq)
  g = min(r,1)
  if (g<.5001 && g>.4999) {g = .5001}
  z = (g - sqrt(g-g^2))/(2*g - 1)
  alpha = acos((1-z)/sqrt(1-2*z*(1-z)))
  v = Re((((2*cos(alpha)*gamma((p+2)/2))/(sqrt(pi)*gamma((p+1)/2)))*(hypergeo(.5,(1-p)/2, 3/2, cos(alpha)^2) - sin(alpha)^(p-1))))
  return(v)
}
)

#Plot all curves (there's probably a cleaner way to do this, if so, let me know)
curve(vstat(Rsq=x, n=300, p=2), 0.01, 0.25, type="l", col="orange", ylim=c(0, 1), xlab="R-squared when Estimating 2 Parameters", ylab="v-stat")
par(new=TRUE)
curve(vstat(Rsq=x, n=200, p=2), 0.01, 0.25, type="l", col="red", ylim=c(0, 1), xaxt = "n", yaxt = "n", xlab="", ylab="")
par(new=TRUE)
curve(vstat(Rsq=x, n=100, p=2), 0.01, 0.25, type="l", col="green", ylim=c(0, 1), xaxt = "n", yaxt = "n", xlab="", ylab="")
par(new=TRUE)
curve(vstat(Rsq=x, n=50, p=2), 0.01, 0.25, type="l", col="purple", ylim=c(0, 1), xaxt = "n", yaxt = "n", xlab="", ylab="")
par(new=TRUE)
curve(vstat(Rsq=x, n=30, p=2), 0.01, 0.25, type="l", col="black", ylim=c(0, 1), xaxt = "n", yaxt = "n", xlab="", ylab="")

#draw horizontal line at 0.5 cut-off
abline(h=0.5, col="azure4", lty=5)
#add legend
legend(0.2,0.44,c("n=300","n=200","n=100","n=50","n=30"), lty=c(1,1,1,1,1), lwd=c(2.5,2.5,2.5,2.5,2.5), col=c("orange","red","green","purple","black")) 

curve(vstat(Rsq=x, n=300, p=3), 0.01, 0.25, type="l", col="orange", ylim=c(0, 1), xlab="R-squared when Estimating 3 Parameters", ylab="v-stat")
par(new=TRUE)
curve(vstat(Rsq=x, n=200, p=3), 0.01, 0.25, type="l", col="red", ylim=c(0, 1), xaxt = "n", yaxt = "n", xlab="", ylab="")
par(new=TRUE)
curve(vstat(Rsq=x, n=100, p=3), 0.01, 0.25, type="l", col="green", ylim=c(0, 1), xaxt = "n", yaxt = "n", xlab="", ylab="")
par(new=TRUE)
curve(vstat(Rsq=x, n=50, p=3), 0.01, 0.25, type="l", col="purple", ylim=c(0, 1), xaxt = "n", yaxt = "n", xlab="", ylab="")
par(new=TRUE)
curve(vstat(Rsq=x, n=30, p=3), 0.01, 0.25, type="l", col="black", ylim=c(0, 1), xaxt = "n", yaxt = "n", xlab="", ylab="")

#draw horizontal line at 0.5 cut-off
abline(h=0.5, col="azure4", lty=5)
#add legend
legend(0.2,0.44,c("n=300","n=200","n=100","n=50","n=30"), lty=c(1,1,1,1,1), lwd=c(2.5,2.5,2.5,2.5,2.5), col=c("orange","red","green","purple","black")) 


curve(vstat(Rsq=x, n=300, p=4), 0.01, 0.25, type="l", col="orange", ylim=c(0, 1), xlab="R-squared when Estimating 4 Parameters", ylab="v-stat")
par(new=TRUE)
curve(vstat(Rsq=x, n=200, p=4), 0.01, 0.25, type="l", col="red", ylim=c(0, 1), xaxt = "n", yaxt = "n", xlab="", ylab="")
par(new=TRUE)
curve(vstat(Rsq=x, n=100, p=4), 0.01, 0.25, type="l", col="green", ylim=c(0, 1), xaxt = "n", yaxt = "n", xlab="", ylab="")
par(new=TRUE)
curve(vstat(Rsq=x, n=50, p=4), 0.01, 0.25, type="l", col="purple", ylim=c(0, 1), xaxt = "n", yaxt = "n", xlab="", ylab="")
par(new=TRUE)
curve(vstat(Rsq=x, n=30, p=4), 0.01, 0.25, type="l", col="black", ylim=c(0, 1), xaxt = "n", yaxt = "n", xlab="", ylab="")

#draw horizontal line at 0.5 cut-off
abline(h=0.5, col="azure4", lty=5)
#add legend
legend(0.2,0.44,c("n=300","n=200","n=100","n=50","n=30"), lty=c(1,1,1,1,1), lwd=c(2.5,2.5,2.5,2.5,2.5), col=c("orange","red","green","purple","black")) 


#Plot all curves (there's probably a cleaner way to overlay them)
curve(vstat(Rsq=x, n=300, p=6), 0.01, 0.25, type="l", col="orange", ylim=c(0, 1), xlab="R-squared when Estimating 6 Parameters", ylab="v-stat")
par(new=TRUE)
curve(vstat(Rsq=x, n=200, p=6), 0.01, 0.25, type="l", col="red", ylim=c(0, 1), xaxt = "n", yaxt = "n", xlab="", ylab="")
par(new=TRUE)
curve(vstat(Rsq=x, n=100, p=6), 0.01, 0.25, type="l", col="green", ylim=c(0, 1), xaxt = "n", yaxt = "n", xlab="", ylab="")
par(new=TRUE)
curve(vstat(Rsq=x, n=50, p=6), 0.01, 0.25, type="l", col="purple", ylim=c(0, 1), xaxt = "n", yaxt = "n", xlab="", ylab="")
par(new=TRUE)
curve(vstat(Rsq=x, n=30, p=6), 0.01, 0.25, type="l", col="black", ylim=c(0, 1), xaxt = "n", yaxt = "n", xlab="", ylab="")

#draw horizontal line at 0.5 cut-off
abline(h=0.5, col="azure4", lty=5)
#add legend
legend(0.2,0.44,c("n=300","n=200","n=100","n=50","n=30"), lty=c(1,1,1,1,1), lwd=c(2.5,2.5,2.5,2.5,2.5), col=c("orange","red","green","purple","black"))