## -----------------------------------------------------------------------------------
## DÉBUT DE DÉFINITION DES FONCTIONS
## -----------------------------------------------------------------------------------

## PARALLEL ANALYSIS N(0,1)
parallel <- function(sujets=100,var=10,rep=100,cent=0.05)
 {
  r       = sujets
  c       = var
  y       = matrix(c(1:r*c), nrow=r, ncol=c)
  ycor    = matrix(c(1:c*c), nrow=c, ncol=c)
  evpea   = NULL
  leg.txt = "Pearson"
  # Simulation of k samples to obtain k random eigenvalues vectors
  # for Pearson correlation coefficients
  for (k in c(1:rep)) {
   y     = rnorm(y)
   y     = matrix(y, nrow=r, ncol=c)
   evpea = rbind(evpea, eigen(cor(y))[[1]]) 
  }
  # Summary statistics
  sprob = c(cent)
  mevpea = sapply(as.data.frame(evpea),mean)
  quant  <- function(x,sprobs=sprobs) {as.vector(quantile(x, probs = sprob)) }
  qevpea =  sapply(as.data.frame(evpea),quant)
  result <- list(eigen=data.frame(mevpea,qevpea),
                 sujets=r, variables=c, centile=cent)
  class(result) <- 'parallel'
  return(result)
 }

## USUAL SCREE TEST
uscree <- function(eig, main="", xlab="Component", ylab="Eigenvalue") {
 par(mfrow=c(1,1))
 nk  = length(eig)
 k   = 1:nk
 plot(k,eig, type = 'b',col="black", pch=1,
      ylab = ylab,
      xlab = xlab,
      main = main) }


## TWO-POINTS REGRESSION
tworeg <- function(eig) {
 nk  = length(eig)
 k   = 1:nk
 tag=0;flag = rep(1,nk); rep = rep(NA,nk); rd = as.character(rep("",nk))
 plot(k,vp, type = 'b',col="black", 
       ylab = "Eigenvalue",
       xlab = "Component",
       main = "TWO-POINTS REGRESSION")
 for (i in 2:(nk-1)) {
  ind = k[c(i,nk)]
  vp.p = lm(vp[c(i,nk)] ~ ind)
  vp.préc = rep[i-1] = sum(c(1,i-1)* coef(vp.p))  
  if (eig[i-1] < rep[i-1]) flag[i-1] = 0
  if ((flag[i-1] == 0) & (flag[i] == 1) & (tag == 0)) {rd[i-2] = "***";tag=1;nc=(i-2)}
  par(col=10+i)
  x =  sum(c(1,1)* coef(vp.p))
  y =  sum(c(1,nk)* coef(vp.p))
  lines(k[c(1,nk)],c(x,y)) }
 # Output Component, eigenvalue, predicted eigenvalue, Star for the value choosen
 list(summary=data.frame(Component       = k, 
                         Eigenvalue      = eig, 
                         Predicted.Eigen = rep, 
                         Flag            = flag, 
                         Criteria        = rd),
                         ncretained      = nc)
}
# tworeg(vp)

tworegW <- function(sujets=100,eig, graph=F,aparallel=rep(1,length(eig))) {
 nk  = length(eig)
 k   = 1:nk
 # if (sum(aparallel) != nk) parallel(sujets,var = nk, rep = 10)
 cond1 = T; cond2 = T; i = 0; pred.eig = af =rep(NA,nk)
 while ((cond1 == T) && (cond2 == T) && (i < nk)) {
  i       = i + 1
  ind     = k[c(i+1,nk)]
  # Optimal coordinate based on the next eigenvalue regression
  vp.p    = lm(eig[c(i+1,nk)] ~ ind)
  vp.préc = pred.eig[i] = sum(c(1,i)* coef(vp.p)) 
  cond1   = (eig[i] >= vp.préc)
  cond2   = (eig[i] >= aparallel[i])
  nc = i-1
  # print(c(i,eig[i], vp.préc, aparallel[i], af[i]))
  } 
  # Second derivative at the i eigenvalue (acceleration factor)
  # See Yakowitz and Szidarovszky (1986, p. 84)
  tag = 1
  for (j in 2:(nk-1)) {
    if (eig[j-1] >= aparallel[j-1]) af[j] = (eig[j+1] -2* eig[j]) + eig[j-1]
   }
  
  par(col=1)
  p.vec = which(eig >= aparallel,T)
  npar = sum(p.vec == (1:length(p.vec)))
  nkeyser = sum(eig >= rep(1,nk))
  # p.af = which(af > 0,T)
  # m.af = min(p.af)
  naf = which(af == max(af,na.rm=T),T) - 1

if (graph == T){
   par(mfrow = c(1,1))
   uscree(eig)
   vp.p = lm(eig[c(i,nk)] ~ k[c(i,nk)])
   x    = sum(c(1,1) * coef(vp.p))
   y    = sum(c(1,nk)* coef(vp.p))
   par(col=10)
   lines(k[c(1,nk)],c(x,y))
   par(col=11, pch=2)
   lines(1:nk, aparallel, type = "b")
   leg.txt <- c(paste("Eigenvalues ........(nkeyser =",nkeyser,")"), 
              c(paste("Parallel Analysis ..(n =",npar,")")),
              c(paste("Optimal Coordinates (n =",nc,")")),        
              c(paste("Acceleration factor (n =",naf,")")) )
   legend("topright",
          legend = leg.txt, 
          pch = c(1,2,16,16), 
          text.col= "black", col="black")
   }

  return(list(Components = data.frame(noc           = nc, 
                               naf           = naf,
                               npar.analysis = npar, 
                               nkeyser       = nkeyser), 
       Analysis   = data.frame(Eigenvalues   = eig, 
                               Par.Analysis  = aparallel,
                               Pred.eig      = pred.eig,
                               Acc.factor    = af)))
}


## -----------------------------
## INITIALISATION
sujets = 30
var    = 11
rep    = 100  
cent   = 0.95         # Centile value of the parallel analysis
vp1 = c(3.12,1.46,1.22,1.06,0.88,0.76,0.70,0.59,0.45,0.40,0.35)
vp2 = c(3.12,1.46,1.22,1.16,0.88,0.76,0.70,0.59,0.45,0.40,0.35)
vp3 = c(3.12,1.23,1.22,1.16,0.88,0.76,0.70,0.59,0.45,0.40,0.35)
vp4 = c(3.12,2.7,1.22,1.16,0.88,0.76,0.70,0.59,0.45,0.40,0.35)
vp  = eig = vp4

aparallel = parallel(var=length(eig),sujets=sujets, rep=rep, cent=cent)$eigen[,1]
tworegW(eig=eig,graph=T,aparallel=aparallel)
# tworegW(eig=eig,graph=T)
# tworeg(eig)

# par(mfrow=c(1,1))
# uscree(eig)


# par(mfrow=c(1,1))
# print(nkreg(vp))
# print(tworeg(vp))
# par(mfrow=c(1,1))