# R/fourth.corner.r In zdealveindy/weimea: Weighted Mean Analysis

```#' calculates the fourth-corner statistic (according to Legendre et al. 1997, compiled according to Peres-Neto et al. 2012)
#' @param R,env Matrix of sample attributes.
#' @param L,sitspe Matrix of species composition.
#' @param Q,speatt Matrix of species attributes
#' @param chessel Logical; should the Chessel's variant of the fourth corner statistic be returned? Default is \code{FALSE}
#' @param fc.test Test to be chosen for fourthcorner analysis.
#' @param perm Number of permutations.
#' @param x,object Object of 'cwm' class.
#' @param ... Other arguments passing into \code{print} or \code{summary} functions (not implemented yet).
#' @details
#' Chessel's version of the fourth corner statistic has been introduced by Peres-Neto et al. (2016), to rescale the fourth corner r statistic into the range of -1 and 1. This is achieved by dividing the observed fourth corner r by square root of the first eigenvalue for correspondance analysis of L (sitspe) matrix.
#' @export
fourth.corner <- function (R, L, Q, chessel = FALSE)
{
R <- as.matrix (R)
L <- as.matrix (L)
Q <- as.matrix (Q)

missing.Q <- ifelse (any (is.na (Q)), TRUE, FALSE)
missing.R <- ifelse (any (is.na (R)), TRUE, FALSE)

if (!missing.Q & !missing.R) r.h <- fourth.corner.0 (R, L, Q, chessel = chessel)
if (missing.Q & !missing.R)
r.h <- apply (Q, 2, FUN = function (q) {L <- L[,!is.na (q)]; q <- q[!is.na (q)]; fourth.corner.0 (R, L, q, chessel = chessel)})
if (!missing.Q & missing.R)
r.h <- t(apply (R, 2, FUN = function (r) {L <- L[!is.na (r),]; r <- r[!is.na (r)]; fourth.corner.0 (r, L, Q, chessel = chessel)}))
if (missing.Q & missing.R)
r.h <- apply (Q, 2, FUN = function (q) apply (R, 2, FUN = function (r) {L <- L[!is.na (r), !is.na (q)]; q <- q[!is.na (q)]; r <- r[!is.na (r)]; fourth.corner.0 (r, L, q, chessel = chessel)}))
if (is.null (colnames (r.h))) colnames (r.h) <- colnames (Q)
if (is.null (rownames (r.h))) rownames (r.h) <- colnames (R)
class (r.h) <- 'fc'
return (r.h)
}

fourth.corner.0 <- function (R, L, Q, chessel)  # this can be used only in case of no missing values in Q
{
W.k <- diag (colSums (L))
Q <- as.matrix (Q)
R <- as.matrix (R)
L <- as.matrix (L)
W.n <- diag (rowSums (L))
I.k <- matrix (rep (1, ncol (L)))
I.n <- matrix (rep (1, nrow (L)))
I.g <- matrix (rep (1, ncol (Q)))
trace.m <- function (M) sum (diag (M))

Q.std.1 <- (Q - I.k %*% t(I.k) %*% W.k %*% Q * (1/trace.m (W.k)))
Q.std.2 <- (I.k %*% ((t(I.k) %*% W.k %*% (Q * Q) - ((t(I.k) %*% W.k %*% Q) * (t(I.k) %*% W.k %*% Q))/trace.m (W.k)) / trace.m (W.k)) ^ 0.5)
Q.std <- Q.std.1 / Q.std.2

R.std.1 <- (R - I.n %*% t(I.n) %*% W.n %*% R * (1/trace.m (W.n)))
R.std.2 <- (I.n %*% ((t(I.n) %*% W.n %*% (R * R) - ((t(I.n) %*% W.n %*% R) * (t(I.n) %*% W.n %*% R))/trace.m (W.n)) / trace.m (W.n)) ^ 0.5)
R.std <- R.std.1 / R.std.2

#Q.L <- apply (Q.std, 2, FUN = function (Q.std.1) I.n %*% t (Q.std.1) * L)
X.P <- L %*% Q.std / W.n %*% I.n %*% t(I.g)

r.h <- t(apply (R.std, 2, FUN = function (R.h.std) (t(R.h.std) %*% W.n %*% R.h.std) ^-1 %*% t(R.h.std) %*% W.n %*% X.P))
colnames (r.h) <- colnames (Q)
if (chessel) r.h <- r.h/sqrt (vegan::cca (L)\$CA\$eig[1])
return (r.h)
}

#' @export
#' @rdname fourth.corner
print.fc <- function (x, ...) print.default (as.matrix (unclass (x)))

#' @export
#' @rdname fourth.corner
summary.fc <- function (object, ...) print (object)

#' @export
#' @rdname fourth.corner
fourth.corner.ade <- function (sitspe, speatt, env, fc.test, perm, chessel = FALSE)
{
R <- as.data.frame (env)
L <- as.data.frame (sitspe)
Q <- as.data.frame (speatt)

missing.Q <- ifelse (any (is.na (Q)), TRUE, FALSE)
missing.R <- ifelse (any (is.na (R)), TRUE, FALSE)

if (!missing.Q & !missing.R) r.h <- fourth.corner.ade0 (R = R, L = L, Q = Q, fc.test = fc.test, perm = perm, chessel = chessel)
if (missing.Q & !missing.R){
r.h <- sapply (names (Q), FUN = function (q) {L1 <- L[,!is.na (Q[,q]), drop = FALSE]; Q1 <- Q[!is.na (Q[,q]), q, drop = FALSE]; fourth.corner.ade0 (R = R, L = L1, Q = Q1, fc.test = fc.test, perm = perm, chessel = chessel)}, simplify = FALSE)
r.h <- do.call (rbind.data.frame, r.h)
rownames (r.h) <- r.h\$names
r.h <- r.h[,-1]
}

if (!missing.Q & missing.R){
r.h <- sapply (names (R), FUN = function (r) {L1 <- L[!is.na (R[,r]),, drop = FALSE]; R1 <- R[!is.na (R[,r]), r, drop = FALSE]; fourth.corner.ade0 (R = R1, L = L1, Q = Q, fc.test = fc.test, perm = perm, chessel = chessel)}, simplify = FALSE)
r.h <- do.call (rbind.data.frame, r.h)
rownames (r.h) <- r.h\$names
r.h <- r.h[,-1]
}

if (missing.Q & missing.R){
r.h <- sapply (names (Q), FUN = function (q) sapply (names (R), FUN = function (r) {L1 <- L[!is.na (R[,r]), !is.na (Q[,q])]; Q1 <- Q[!is.na (Q[,q]), q, drop = FALSE]; R1 <- R[!is.na (R[,r]), r, drop = FALSE]; fourth.corner.ade0 (R = R1, L = L1, Q = Q1, fc.test = fc.test, perm = perm, chessel = chessel)}, simplify = FALSE), simplify = FALSE)
r.h <- do.call (rbind.data.frame, lapply (r.h, FUN = function (y) do.call (rbind.data.frame, y)))
rownames (r.h) <- r.h\$names
r.h <- r.h[,-1]
}
return (r.h)
}

fourth.corner.ade0 <- function (R, L, Q, fc.test, perm, chessel)
{
temp <- ade4::fourthcorner (tabR = R, tabL = L, tabQ = Q, modeltype = fc.test, nrepet = perm)
temp <- do.call (cbind.data.frame, temp\$tabD[c('names', 'obs', 'pvalue')])
names (temp) <-  c('names', 'fourthcorner', 'P.value')
if (chessel) temp[,'fourthcorner'] <- temp[,'fourthcorner']/sqrt (vegan::cca (L)\$CA\$eig[1])
return (temp)
}
```
zdealveindy/weimea documentation built on June 14, 2018, 3:45 a.m.