#HOMOGENEITY OF MULTIVARIATE DISPERSIONS MODIFIED FOR COMPOSITIONAL FIDELITY #ANALYSES (LIVING VERSUS DEATH ASSEMBLAGES) FOR UNEQUAL NUMBER OF LIVING AND DEATH #ASSEMBLAGES. 


#LIBRARIES
require(vegan)
require(ade4)
###################################################################################

#THIS IS THE FUNCTION
#This function allows different number of living and death assemblages
#grouping = factor with sample assignment to "Live"or "Dead" (or any other equivalent)
fidelity<-function(livedead, grouping, method=method) {

livedead<-livedead; 
method=method;
livegroups=which(grouping=="Live")
deadgroups=which(grouping=="Dead")

if (method=="horn") {distance<-vegdist(t(livedead),method="horn")}
if (method=="bray") {distance<-vegdist(t(livedead),method="bray")}
if (method=="hellinger") {distance<-vegdist(t(sqrt(livedead)),method="bray")}
if (method=="sorenson") {distance<-vegdist(t(livedead),method="bray", binary=T)}
if (method=="jaccard") {distance<-vegdist(t(livedead),method="jaccard", binary=T)}
if (method=="gower") {distance<-vegdist(decostand(t(livedead), "log"), "altGower")}
if (method=="euclidean") {distance<-vegdist(t(livedead),method="euclidean")}

#PRINCIPAL COORDINATE ANALYSIS USING ADE4 PACKAGE
tol = -1;
distmat <- as.matrix(distance);n <- ncol(distmat)
row.w <- rep(1, n)
row.w <- row.w/sum(row.w)
delta <- -0.5 * bicenter.wt(distmat * distmat, row.wt = row.w, col.wt = row.w)
wsqrt <- sqrt(row.w)
delta <- delta * wsqrt
delta <- t(t(delta) * wsqrt)
eig <- eigen(delta, symmetric = TRUE)
lambda <- eig$values
w0 <- lambda[n]/lambda[1]
r <- sum(lambda > (lambda[1] * tol))
nf<-r
res <- list()
res$eig <- lambda[1:r]                  #list of eigenvalues
res$rank <- r
res$cw <- rep(1, r)
w <- t(t(eig$vectors[, 1:r]) * sqrt(abs(lambda[1:r])))/wsqrt            
w <- data.frame(w)
res$tab <- w
res$li <- data.frame(w[, 1:nf])         #eigenvectors - principal coordinates 

#ESTIMATE THE NUMBER OF POSITIVE AND NEGATIVE EIGENVALUES
positive<-sum(res$eig > 0); negative<-sum(res$eig < 0);

#ESTIMATE DISTANCES BETWEEN DEATH ASSEMBLAGES AND CENTROID OF LIVING ASSEMBLAGES #(TOTAL LIVE-DEAD VARIATION)
#COMPUTE CENTROID LOCATION OF LIVING ASSEMBLAGES
#FOR POSITIVE EIGENVALUES
livecentroid1<-rep(NA, length(NA))
for (i in 1:positive) 
                {livecentroid1[i]<-mean(res$li[1:length(livegroups),i])
                }
#FOR NEGATIVE EIGENVALUES
if (negative > 0) 
                {livecentroid2<-rep(NA, length(NA))
                        for (i in (positive+1):(positive+negative)) 
                        {livecentroid2[i]<-mean(res$li[1:length(livegroups),i])
                        }
                }

#DISTANCES BETWEEN DEATH ASSEMBLAGES AND CENTROID OF LIVING ASSEMBLAGES FOR #POSITIVE EIGENVALUES
if (length(livecentroid1) == 1) 
                {n1<-append(livecentroid1, res$li[(length(livegroups)+1):ncol(livedead),1:positive])
                }
if (length(livecentroid1) > 1) 
                {n1<-rbind(livecentroid1, res$li[(length(livegroups)+1):ncol(livedead),1:positive])
                }
nn1<-designdist(n1,method="(A+B-2*J)", terms=c("quadratic"))    #distances between centroid and all assemblages (living and death) 
positiveLDdistances<-as.vector(nn1[1:length(deadgroups)])       #extract distances between centroid of LA and death assemblages only

#DISTANCES BETWEEN DEATH ASSEMBLAGES AND CENTROID OF LIVING ASSEMBLAGES FOR #NEGATIVE EIGENVALUES
if (negative == 0) 
                {LCdead<-sqrt(positiveLDdistances)
                }
if (negative > 1) 
                {
                n2<-rbind(livecentroid2[(positive+1):(positive+negative)], res$li[(length(livegroups)+1):ncol(livedead),(positive+1):(positive+negative)])
                nn2<-designdist(n2,method="(A+B-2*J)", terms=c("quadratic"))    #distances between centroid and all assemblages (living and death) 
                negativeLDdistances<-as.vector(nn2[1:length(deadgroups)])       #extract distances between centroid of LA and death assemblages only
                LCdead<-sqrt(abs(positiveLDdistances-negativeLDdistances))
                }

if (negative == 1)      
                {n2<-append(livecentroid2[(positive+1):(positive+negative)], res$li[(length(livegroups)+1):ncol(livedead),(positive+1):(positive+negative)])
                nn2<-designdist(n2,method="(A+B-2*J)", terms=c("quadratic"))
                negativeLDdistances<-as.vector(nn2[1:length(deadgroups)])
                LCdead<-sqrt(abs(positiveLDdistances-negativeLDdistances))
                }

#DISTANCES BETWEEN LIVING ASSEMBLAGES AND CENTROID OF LIVING ASSEMBLAGES FOR #POSITIVE EIGENVALUES
if (length(livecentroid1) == 1) 
                {n1<-append(livecentroid1, res$li[1:length(livegroups),1:positive])}
if (length(livecentroid1) > 1) 
                {n1<-rbind(livecentroid1, res$li[1:length(livegroups),1:positive])}
nn1<-designdist(n1,method="(A+B-2*J)", terms=c("quadratic"))    #all distances between centroid and among samples for positive eigenvalues
positiveLLdistances<-as.vector(nn1[1:length(livegroups)])       

#DISTANCES BETWEEN LIVING ASSEMBLAGES AND CENTROID OF LIVING ASSEMBLAGES FOR #NEGATIVE EIGENVALUES 
if (negative == 0) 
                {LClive<-sqrt(positiveLLdistances)
                }
if (negative > 1) 
                {n2<-rbind(livecentroid2[(positive+1):(positive+negative)], res$li[1:length(livegroups),(positive+1):(positive+negative)])
                nn2<-designdist(n2,method="(A+B-2*J)", terms=c("quadratic"))
                negativeLLdistances<-as.vector(nn2[1:length(livegroups)])       
                LClive<-sqrt(abs(positiveLLdistances-negativeLLdistances))
                }
if (negative == 1) 
                {n2<-append(livecentroid2[(positive+1):(positive+negative)], res$li[1:length(livegroups),(positive+1):(positive+negative)])
                nn2<-designdist(n2,method="(A+B-2*J)", terms=c("quadratic"))
                negativeLLdistances<-as.vector(nn2[1:length(livegroups)])
                LClive<-sqrt(abs(positiveLLdistances-negativeLLdistances))}

################################################
#COMPUTE CENTROID LOCATION OF DEATH ASSEMBLAGES
#FOR POSITIVE EIGENVALUES
deadcentroid1<-rep(NA, length(NA))
for (i in 1:positive) 
                {deadcentroid1[i]<-mean(res$li[(length(livegroups)+1):ncol(livedead),i])
                }
#FOR NEGATIVE EIGENVALUES
if (negative > 0) 
                {deadcentroid2<-rep(NA, length(NA))
                        for (i in (positive+1):(positive+negative)) 
                        {deadcentroid2[i]<-mean(res$li[1:ncol(livedead),i])
                        }
                }

#DISTANCES BETWEEN LIVING AND DEAD CENTROID FOR POSITIVE EIGENVALUES
if (length(livecentroid1) == 1) 
                {n1<-append(livecentroid1, deadcentroid1)}
if (length(livecentroid1) > 1) 
                {n1<-rbind(livecentroid1, deadcentroid1)}
nn1<-designdist(n1,method="(A+B-2*J)", terms=c("quadratic"))    #all distances between centroid and among samples for positive eigenvalues
positivedistance<-as.vector(nn1)

#DISTANCES BETWEEN LIVING AND DEAD CENTROID FOR POSITIVE EIGENVALUES
if (negative == 0) 
                {LCDC<-sqrt(positivedistance)
                }
if (negative > 1) 
                {n2<-rbind(livecentroid2[(positive+1):(positive+negative)], deadcentroid2[(positive+1):(positive+negative)])
                nn2<-designdist(n2,method="(A+B-2*J)", terms=c("quadratic"))
                negativedistance<-as.vector(nn2)       
                LCDC<-sqrt(abs(positivedistance-negativedistance))
                }
if (negative == 1) 
                {n2<-append(livecentroid2[(positive+1):(positive+negative)],deadcentroid2[(positive+1):(positive+negative)])
                nn2<-designdist(n2,method="(A+B-2*J)", terms=c("quadratic"))
                negativedistance<-as.vector(nn2)
                LCDC<-sqrt(abs(positivedistance-negativedistance))
                }


#############################################
#95% BOOTSTRAPPED CONFIDENCE INTERVALS FOR DISTANCES OF LIVING ASSEMBLAGES TO #CENTROID OF LIVING ASSEMBLAGES-LClive
#AND FOR DISTANCES OF DEATH ASSEMBLAGES TO 
#CENTROID OF LIVING ASSEMBLAGES - LCdead

nLive<-length(livegroups)                   
nDead<-length(deadgroups)
LCdead.resample<-rep(NA,nLive);LClive.resample<-rep(NA,nDead);          
num.rep<-100                    
LCdeadmean.resample<-rep(NA,num.rep); LClivemean.resample<-rep(NA,num.rep);
for (i in 1:num.rep) {
        for (j in 1:nLive) {
                LClive.resample[j]<-sample(LClive,1,replace=T)
                }
        for (j in 1:nDead) {
                LCdead.resample[j]<-sample(LCdead,1,replace=T)
                }
        LCdeadmean.resample[i]<-mean(LCdead.resample)
        LClivemean.resample[i]<-mean(LClive.resample)
        }

LCdeadci<-quantile(LCdeadmean.resample, c(0.025,0.975),na.rm=T)
LCliveci<-quantile(LClivemean.resample, c(0.025,0.975),na.rm=T)
                                                        
#TOTAL LIVE-DEAD VARIATION IN TERMS OF SUM OF SQUARED DISTANCES
SST = sum(LCdead^2)

#PREMORTEM VARIATION IN TERMS OF SUM OF SQUARED DISTANCES
SSL = sum(LClive^2)

#PERMUTEST USING VEGAN PACKAGE
control = permControl(nperm = 1000)
xdistances<-append(LCdead,LClive)
xgroup<-grouping
nobs <- length(xdistances)
mod <- lm(xdistances ~ xgroup)
mod.Q <- mod$qr
p <- mod.Q$rank
resids <- qr.resid(mod.Q, xdistances)
res <- numeric(length = control$nperm + 1)
res[1] <- summary(mod)$fstatistic[1]
  
for (i in seq(along = res[-1])) {
        perm <- shuffle(nobs, control = control)
        perm.resid <- resids[perm]
        f <- qr.fitted(mod.Q, perm.resid)
        mss <- sum((f - mean(f))^2)
        r <- qr.resid(mod.Q, perm.resid)
        rss <- sum(r^2)
        rdf <- nobs - p
        resvar <- rss/rdf
        res[i + 1] <- (mss/(p - 1))/resvar
        }
pvalue <- sum(res >= res[1])/length(res)

out <- list("Distances between DAs and centroid of LA" = LCdead, "Distances between LAs and centroid of LA" = LClive, 
"Magnitude of total live-dead variation (mean centroid of LA-DA distance)" = mean(LCdead), 
"Magnitude of premortem variation (mean centroid of LA-LA distance)" = mean(LClive), 
"Magnitude of postmortem variation (can be negative)" = mean(LCdead)-mean(LClive), "Between-centroid distance" = LCDC,
"95% CI for mean centroid of LA-DA distance" = LCdeadci, "95% CI for mean centroid of LA-LA distance" = LCliveci, "F value"=res[1],
"P value for homogeneity of dispersions" = pvalue)
class(out) <- "fidelity"
out
}