Get Phi statistics from amova object from either pegas or ade4

statphi.R

ade4phi <- function(sig){
  n <- length(sig)
  f <- rep(0, n - 1)
  if (length(sig) == 3) {
    f <- rep(0, 1)
  }
  f[1] <- (sig[n] - sig[n - 1]) / sig[n]
  if (length(f) > 1) {
    s1 <- cumsum(sig[(n - 1):2])[-1]
    s2 <- sig[(n - 2):2]
    f[length(f)] <- sig[1] / sig[n]
    f[2:(length(f) - 1)] <- s2 / s1
  }
  return(f)
}

hfsphi <- function(sig){
  tot      <- sig
  nblevels <- length(tot)
  n.names  <- names(sig)
  n.cols   <- n.names[-nblevels]
  n.rows   <- c("Total", n.names[-(nblevels - 0:1)])
  f <- matrix(rep(0, (nblevels - 1)^2), 
              ncol = (nblevels - 1),
              dimnames = list(n.rows, n.cols))
  for (i in 1:(nblevels - 1)) {
    for (j in i:(nblevels - 1)) {
      f[i, j] <- sum(tot[i:j])/sum(tot[i:nblevels])
    }
  }
  f
}
# adapted from the ade4 amova code
Statphi.default <- function(sig, method = "ade4") {
  method <- match.arg(method, c("ade4", "hierfstat"))
  if (method == "ade4") ade4phi(sig) else hfsphi(sig[-length(sig)])
}

is_pegas_amova <- function(sig) identical(names(sig), c("tab", "varcoef", "varcomp", "call"))

Statphi.amova <- function(sig, ...){
  if (is_pegas_amova(sig)){
    if (is.data.frame(sig$varcomp)){
      siggies <- setNames(sig$varcomp$sigma2, rownames(sig$varcomp))
    } else {
      siggies <- sig$varcomp
    }
    Statphi(c(siggies, sum(siggies)), ...)
  } else {
    Statphi(sig$componentsofcovariance$Sigma, ...)
  }
}

Statphi <- function(sig, ...){
  UseMethod("Statphi")
}

library("poppr")
library("pegas")
library("hierfstat")
data(nancycats)
strata(nancycats) <- data.frame(pop = pop(nancycats))
a <- poppr.amova(nancycats, ~pop, method = "pegas", nperm = 0)
Statphi(a)
Statphi(a, "h")

data(microbov)
strata(microbov) <- data.frame(other(microbov))
b <- poppr.amova(microbov, ~coun/spe, method = "pegas", nperm = 0)
Statphi(b)
Statphi(b, "h")

# PhiST is always going to be in the top right corner of the matrix
# all other Phi values will be the diagonal. 
bc <- poppr.amova(microbov, ~coun/spe/breed, method = "pegas", nperm = 0)
pha <- Statphi(bc)
phh <- Statphi(bc, "h")
phh <- c(phh[1, ncol(phh)], rev(diag(phh)))
phh
pha

ade4phi <- function(sig){
  n <- length(sig)
  f <- rep(0, n - 1)
  if (length(sig) == 3) {
    f <- rep(0, 1)
  }
  f[1] <- (sig[n] - sig[n - 1]) / sig[n]
  if (length(f) > 1) {
    s1 <- cumsum(sig[(n - 1):2])[-1]
    s2 <- sig[(n - 2):2]
    f[length(f)] <- sig[1] / sig[n]
    f[2:(length(f) - 1)] <- s2 / s1
  }
  return(f)
}

hfsphi <- function(sig){
  tot      <- sig
  nblevels <- length(tot)
  f <- matrix(rep(0, (nblevels - 1)^2), 
              ncol = (nblevels - 1))
  for (i in 1:(nblevels - 1)) {
    for (j in i:(nblevels - 1)) {
      f[i, j] <- sum(tot[i:j])/sum(tot[i:nblevels])
    }
  }
  f
}
# adapted from the ade4 amova code
Statphi.default <- function(sig, method = "ade4") {
  method <- match.arg(method, c("ade4", "hierfstat"))
  if (method == "ade4") ade4phi(sig) else hfsphi(sig[-length(sig)])
}

is_pegas_amova <- function(sig) identical(names(sig), c("tab", "varcoef", "varcomp", "call"))

Statphi.amova <- function(sig, ...){
  if (is_pegas_amova(sig)){
    if (is.data.frame(sig$varcomp)){
      siggies <- sig$varcomp$sigma2
    } else {
      siggies <- sig$varcomp
    }
    Statphi(c(siggies, sum(siggies)), ...)
  } else {
    Statphi(sig$componentsofcovariance$Sigma, ...)
  }
}

Statphi <- function(sig, ...){
  UseMethod("Statphi")
}

library("poppr")
#> Loading required package: adegenet
#> Loading required package: ade4
#> 
#>    /// adegenet 2.1.0 is loaded ////////////
#> 
#>    > overview: '?adegenet'
#>    > tutorials/doc/questions: 'adegenetWeb()' 
#>    > bug reports/feature requests: adegenetIssues()
#> This is poppr version 2.5.0.99.64. To get started, type package?poppr
#> OMP parallel support: available
#> 
#> This version of poppr is under development.
#> If you find any bugs, please report them at https://github.com/grunwaldlab/poppr/issues
library("pegas")
#> Loading required package: ape
#> 
#> Attaching package: 'pegas'
#> The following object is masked from 'package:ape':
#> 
#>     mst
#> The following object is masked from 'package:ade4':
#> 
#>     amova
library("hierfstat")
#> 
#> Attaching package: 'hierfstat'
#> The following objects are masked from 'package:ape':
#> 
#>     pcoa, varcomp
#> The following object is masked from 'package:adegenet':
#> 
#>     read.fstat
data(nancycats)
strata(nancycats) <- data.frame(pop = pop(nancycats))
a <- poppr.amova(nancycats, ~pop, method = "pegas", nperm = 0)
#> 
#> Found 1234 missing values.
#> 
#> 2 loci contained missing values greater than 5%
#> 
#> Removing 2 loci: fca8, fca45
#> Warning in is.euclid(xdist): Zero distance(s)
#> Distance matrix is non-euclidean.
#> Utilizing quasieuclid correction method. See ?quasieuclid for details.
#> Warning in is.euclid(distmat): Zero distance(s)
Statphi(a)
#> [1] 0.20563714 0.13469696 0.08198305
Statphi(a, "h")
#>            [,1]      [,2]
#> [1,] 0.08198305 0.2056371
#> [2,] 0.00000000 0.1346970

data(microbov)
strata(microbov) <- data.frame(other(microbov))
b <- poppr.amova(microbov, ~coun/spe, method = "pegas", nperm = 0)
#> 
#> Found 12650 missing values.
#> 
#> 2 loci contained missing values greater than 5%
#> 
#> Removing 2 loci: ILSTS6, ETH185
#> Distance matrix is non-euclidean.
#> Utilizing quasieuclid correction method. See ?quasieuclid for details.
Statphi(b)
#> [1] 0.19187860 0.10167494 0.07099783 0.03166337
Statphi(b, "h")
#>            [,1]       [,2]      [,3]
#> [1,] 0.03166337 0.10041316 0.1918786
#> [2,] 0.00000000 0.07099783 0.1654541
#> [3,] 0.00000000 0.00000000 0.1016749

# PhiST is always going to be in the top right corner of the matrix
# all other Phi values will be the diagonal. 
bc <- poppr.amova(microbov, ~coun/spe/breed, method = "pegas", nperm = 0)
#> 
#> Found 12650 missing values.
#> 
#> 2 loci contained missing values greater than 5%
#> 
#> Removing 2 loci: ILSTS6, ETH185
#> Distance matrix is non-euclidean.
#> Utilizing quasieuclid correction method. See ?quasieuclid for details.
pha <- Statphi(bc)
phh <- Statphi(bc, "h")
phh <- c(phh[1, ncol(phh)], rev(diag(phh)))
phh
#> [1] 0.16650510 0.05090004 0.06588374 0.02812636 0.03265754
pha
#> [1] 0.16650510 0.05090004 0.06588374 0.02812636 0.03265754

Created on 2017-12-14 by the reprex
package (v0.1.1.9000).