# Empirical Critical Values for the Correlation Coefficient Test of Normality
#
# Reference: Looney and Gulledge (1985), "Use of the Correlation Coefficient with Normal
# Probability Plots," The American Statistician, 75-79.
#
N<-c(5,6,7,8,9,10,11,12,13,14,15,16,17,18,19,20,21,22,23,24,25,26,27,28,29,30,35,40,45,50,60,70,80,90,100)
p05<-c(.880,.888,.898,.906,.912,.918,.923,.928,.932,.935,.939,.941,.944,.946,.949,.951,.952,.954,.956,.957,.959,.960,.961,.962,.963,.964,.969,.972,.974,.977,.980,.983,.985,.986,.987)
p01<-c(.826,.838,.850,.861,.871,.879,.886,.892,.899,.905,.910,.913,.917,.920,.924,.926,.930,.933,.935,.937,.939,.941,.943,.944,.946,.947,.954,.959,.963,.966,.971,.975,.978,.980,.982)
#
# A second set from a technical report on the MINITAB website
#
nobs<-10*1:10
rcrit<-c(0.917,0.950,0.964,0.972,0.977,0.980,0.982,0.984,0.985,0.987)
#
# A routine to estimate the Filliben R
#
filliben <- function(x) {
  y<-sort(x)
  n<-length(y)
  p<-(1:n-0.3175)/(n+0.365)
  p[1]<-1-0.5^(1/n)
  p[n]<-0.5^(1/n)
  p<-qnorm(p)
  cor(y,p)
}
#
# Simulation based distribution of r_f
#
fillmoment <- function(n,reps) {
  moments <- function(data,r) sum((data-mean(data))^r)/length(data)
  skew <- function(data) moments(data,3)/(moments(data,2)*sqrt(moments(data,2)))
  kurt <- function(data) moments(data,4)/(moments(data,2)*moments(data,2))
  nr<-length(n)
  nc<-4
  result<-matrix(0,nrow=nr,ncol=nc)
  rownames(result)<-n
  colnames(result)<-c("Mean","SD","Skew","Kurt")

  x<-double(length=reps)  
  for(i in 1:length(n)) {
    for(j in 1:reps) { x[j]<-log(1-filliben(rnorm(n[i]))) }
    result[i,1]<-mean(x)
    result[i,2]<-sd(x)
    result[i,3]<-skew(x)
    result[i,4]<-kurt(x)
  }
  result
}
#
# Simulation based quantiles for range of Ns -- compare to those published
#
filltable <- function(n,reps,q=c(0.001,0.005,0.01,0.05,0.10,0.20,0.50)) {
  nr<-length(n)
  nc<-length(q)
  result<-matrix(0,nrow=nr,ncol=nc)
  rownames(result)<-n
  colnames(result)<-q

  x<-double(length=reps)  
  for(i in 1:length(n)) {
    for(j in 1:reps) { x[j]<-filliben(rnorm(n[i])) }
    result[i,]<-quantile(x,q)
  }
  result
}
#
# Critical r_f values based on fitting to simulation results 
#
fillcrit <- function(n) {
  a<- -3.7707
  b<-c(0.3044, 0.3805, 0.4673, 0.5211, 0.7210, 0.8656, 1.0861, 1.6813)
  p<-c(0.0001, 0.001,  0.005,  0.01,   0.05,   0.10,   0.20,   0.50)
  r<- 1- 1/(a+b*(13+n))
  names(r)<-p
  r
} 
#
# P-values based on similar approximations to Royston 1993 for W' ie log(1-r)
#
fillp <- function(r,n) {
  u<-log(n)
  v<-log(u)
  mu<- -1.99196+1.0402*(v-u)
  sd<- 0.31293 + 0.788392/u
  pnorm((log(1-r)-mu)/sd, lower=F)
}
#
# P-values based on approximation using 1/(1-r)
#
fillp2 <- function(r,n) {
   b<- (1/(1-r)+3.7707)/(13+n)
   z<- 8.0382-20.9143*b+26.6971*b^2-16.6688*b^3+3.8645*b^4
   pnorm(z, lower=F)
}
#
# test fillp
#
filltestp <- function(n,reps,alpha=0.05,fun=fillp) {
  p<-double(length(n))
  for(i in 1:reps) {
    x<-filliben(rnorm(n)) 
    p<-p+(fun(x,n)<alpha)
  }
  p/reps
}
       
          

#
# An example of use from German Rodriquez's notes on Stata
#
# country<-c("Bolivia","Brazil","Chile","Colombia","CostaRica","Cuba",
#            "DominicanRep","Ecuador","ElSalvador","Guatemala","Haiti",
#            "Honduras","Jamaica","Mexico","Nicaragua","Panama","Paraguay",
#            "Peru","TrinidadTobago","Venezuela")
# setting<-c(46,74,89,77,84,89,68,70,60,55,35,51,87,83,68,84,74,73,84,91)
# effort<-c(0,0,16,16,21,15,14,6,13,9,3,7,23,4,0,19,3,0,15,7)
# change<-c(1,10,29,25,29,40,21,0,13,4,0,7,21,9,7,22,6,2,29,11)
# peg<-cut(effort,c(-99,4,14,99),labels=c("0-4","5-14","15+"))
# 
# m1<-lm(change ~ setting + peg)
#
# Should be 0.9655
# 
# print(paste("Filliben r=",filliben(rstandard(m1))))