# R/Cperiodogram.R In adespatial: Multivariate Multiscale Spatial Analysis

#### Documented in Cperiodogram

#'Contingency periodogram
#'
#'Function to compute a contingency periodogram for a univariate series of
#'qualitative data
#'
#'The contingency periodogram of Legendre et al. (1981) identifies periodic
#'components in qualitative data vectors. The vector may contain classes of a
#'qualitative variable or the classes obtained by hierarchical clustering or
#'partitioning of a multivariate data table. The method is also described in
#'Legendre & Legendre (2012). The optional graph produced by the function shows
#'the following information:
#'
#'\itemize{ \item In red: the B statistics (information in common). \item In
#'blue: Confidence limits for B without correction. \item In green:
#'Bonferroni-corrected limits of the confidence intervals. \item In black:
#'Confidence limits with progressive Bonferroni correction. }
#'
#'@param x a qualitative variable (\code{factor})
#'@param T1 first period included in the calculations (default: T1 = 2)
#'@param T2 last period included in the calculations (default: T2 = n/2)
#'@param nperm Number of permutations for the chi-square test. For chi-square
#'  tests using the chi- square distribution, use the default \code{nperm=NULL}
#'@param alpha significance level for computation of the confidence limits
#'@param graph a logical indicating if a graph is requested, by default
#'  \code{TRUE}.
#'
#'@author Pierre Legendre \email{[email protected]@umontreal.ca}
#'
#'@return A table with the statistics for the selected periods: \itemize{ \item
#'  Wilks’ chi-square statistic (Wilks.chisq) \item information in common (B),
#'  \item degrees of freedom (df), \item p-value (prob) } Confidence interval
#'  limits: \itemize{ \item critical value of B without correction (B.crit),
#'  \item critical value of B with Bonferroni correction based on the number of
#'  periods studied in the periodogram (B.crit.Bonf), \item critical value of B
#'  with progressive Bonferroni correction (B.prog.Bonf). }
#'
#'@references Legendre, L., M. Fréchette & P. Legendre. 1981. The contingency
#'periodogram: a method of identifying rhythms in series on nonmetric ecological
#'data. Journal of Ecology 69: 965-979.
#'
#'Legendre, P. and Legendre, L. 2012. Numerical Ecology. 3rd English ed.
#'Elsevier, Amsterdam
#'
#' @examples
#' # Data from the numerical example of Subsection 12.4.2 of Legendre and Legendre (2012).
#' test.vec <- c(1,1,2,3,3,2,1,2,3,2,1,1,2,3,3,1)
#' # Periodogram with tests using the chi-square distribution
#' res <- Cperiodogram(test.vec)
#' # Periodogram with permutation tests
#' res <- Cperiodogram(test.vec, nperm=2000, graph=FALSE)
#'
#'@importFrom graphics lines
#'@importFrom stats chisq.test qchisq
#'@export
Cperiodogram <- function(x, T1 = 2, T2 = NULL, nperm = NULL, alpha = 0.05, graph = TRUE)
# Contingency periodogram (Legendre et al. 1981).
#
# Author:: Pierre Legendre, November 2010
{
### Internal functions
###
Wilks.chisq <- function(Obs, Exp)
{
rem <- which(Obs==0)
if(length(rem)==0)
{
res <- 2*sum(Obs * log(Obs/Exp))
} else {
res <- 2*sum(Obs[-rem] * log(Obs[-rem]/Exp[-rem]))
}
}

plo <- function(out)
{
plot(out[,1], out[,3], col="red", xlab="Period",
ylab="Information in common (B)", ylim=c(0, max(out[,3], out[,7])))
lines(out[,1], out[,3], col="red")     # B statistics
# points(out[,1], out[,6], col="blue")
lines(out[,1], out[,6], col="blue")    # Ordinary confidence limits
# points(out[,1], out[,7], col="green")
lines(out[,1], out[,7], col="green")   # Bonferroni-corrected conf. limits
# points(out[,1], out[,8], col="black")
lines(out[,1], out[,8], col="black")   # Progressive Bonferroni conf. limits
}
###
### End internal functions

if(!is.factor(x)) if(length(which(x <= 0)) > 0)
stop("x contains 0 or negative levels")
x <- as.factor(x)
nlev.x <- nlevels(x)
n <- length(x)
T.sup <- n %/% 2
if(is.null(T2)) T2 <- T.sup
if(T2 > T.sup) T2 <- T.sup
nk <- T2-T1+1
if(is.null(nperm)) {
simulate <- FALSE
} else {
simulate <- TRUE
}
out <- matrix(NA,nk,8)
colnames(out) <-c("Period","Wilks.chisq","B","df","prob","B.crit","B.crit.Bonf","B.prog.Bonf")
for(k in T1:T2) {
# cat("Period =",k,"\n")
vec <- as.factor(rep(1:k, ceiling(n/k)))
tab <- table(x,vec[1:n])
chisq.out <- chisq.test(tab, simulate.p.value=simulate, B=nperm)
chisq <- Wilks.chisq(chisq.out$observed, chisq.out$expected)
B <- chisq/(2*n)
if(simulate) {df <- (nrow(tab)-1)*(ncol(tab)-1)}
else {df <- chisq.out$parameter} out[(k-T1+1),1] <- k out[(k-T1+1),2] <- chisq out[(k-T1+1),3] <- B out[(k-T1+1),4] <- df out[(k-T1+1),5] <- chisq.out$p.value
out[(k-T1+1),6] <- qchisq(alpha, df, lower.tail=FALSE)/(2*n)
out[(k-T1+1),7] <- qchisq(alpha/nk, df, lower.tail=FALSE)/(2*n)
out[(k-T1+1),8] <- qchisq(alpha/(k-T1+1), df, lower.tail=FALSE)/(2*n)
}
if(graph) plo(out)
out
}