####Load functions from NCAP.R
source(file.choose())

#################################Holdfast data################################
#Species data are 351 (species) by  80 (samples)
HoldfastSpecies.df=read.table("NZHoldfastSpecies.dat",sep="",skip=1)
EnvVars.df=read.table("Environmentals.dat",head=T); attach(EnvVars.df)
D=distance(t(as.matrix(HoldfastSpecies.df)),measure="BC",trans="fourthroot")
pcoD=pco(D)

gradient=gradient.choice("vonB")

X=model(~Volume)
#X=model(~Volume+I(log(Density)*Volume)) #To include volume*density interaction
npar=ncol(X)

#########Sequence of fits for different number of pco's#####################
b0=rep(0.01,npar)
par(mfrow=c(2,2),mar=c(4,5,2,3))
fitseq=NLCCorSeq(b0,X,pcoD,gradient,m=25,stat="Rsq")
fitseq2=NLCCorSeq(b0,X,pcoD,gradient,m=25,stat="RDA")

m=5
############################Holdfast analyses###############################
criterion="Rsq"

#Linear fit
CCr2=LinCCor(X,pcoD,m,stat=criterion)

#Fit nonlinear CAP
NLCC=nlm(NLCCor,b0,X=X,pcoD=pcoD,gradient=gradient,m=m,stat=criterion,
         print.level=1,iterlim=500,typsize=b0,blow=rep(0,npar),bhigh=rep(1,npar),pwgt=0)
cat("\nMaximized",criterion,"statistic  is",-NLCC$min,"\n")
bhat=NLCC$est

#Pseudo F-stat for nonlinear gradient vs linear gradient
cat("\nPseudo F-stat for nonlinear vs linear gradient is",
     (-NLCC$min-CCr2)/(1+NLCC$min),"\n")

par(mfrow=c(2,1),mar=c(4,5,2,3))
plot.NCAP(NLCC$est,X,pcoD,gradient=gradient,m=m,stat=criterion,
          xpts=0:300,ylim=c(0,1))


###############################Bootstraps CI's#################################
nsamps=100 #Takes about 10 seconds on a 2GHz desktop
boot=BootNLCor(b0,X=X,D=D,gradient,m=m,stat=criterion,nsamps=nsamps,typsize=bhat,
      blow=rep(0.1*bhat,npar),bhigh=rep(10*bhat,npar),pwgt=0.0)

###############################Permutations####################################
nsamps=99 #Takes about 1 minute on a 2GHz desktop
perm=PermNLCor(b0,X=X,D=D,gradient,m=m,stat=criterion,nsamps=nsamps,typsize=bhat,
      blow=rep(0.1*bhat,npar),bhigh=rep(10*bhat,npar),pwgt=0.0)