{{ message }}

Instantly share code, notes, and snippets.

# zkamvar/statphi.R

Last active Dec 18, 2017
Get Phi statistics from amova object from either pegas or ade4
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters. Learn more about bidirectional Unicode characters
 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

### zkamvar commented Dec 14, 2017

```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
}
Statphi.default <- function(sig, method = "ade4") {
}

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")
#>
#>
#>    > 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")
#>
#> Attaching package: 'pegas'
#> The following object is masked from 'package:ape':
#>
#>     mst
#>
#>     amova
library("hierfstat")
#>
#> Attaching package: 'hierfstat'
#> The following objects are masked from 'package:ape':
#>
#>     pcoa, varcomp
#>
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).