# dwinter/mmod

diff_stat functions use base R for all calculations -2x speed up

1 parent c4362ab commit fb2612a528f712e059e9e20080e6e400f6c49945 committed Feb 1, 2012
Showing with 19 additions and 12 deletions.
1. +2 −1 R/D_Jost.R
2. +1 −1 R/D_Jost.R~
3. +2 −1 R/Gst_Hedrick.R
4. +1 −1 R/Gst_Hedrick.R~
5. +3 −2 R/Gst_Nei.R
6. +2 −1 R/diff_stats.R
7. +5 −1 R/diff_stats.R~
8. +1 −1 R/summarise_bootstrap.R
9. +2 −3 R/summarise_bootstrap.R~
3 R/D_Jost.R
 @@ -24,8 +24,9 @@ D_Jost <- function(x){ n <- length(unique(pop(x))) harmN <- harmonic_mean(table(pop(x))) + pops <- pop(x) D.per.locus <- function(g) { - a <- makefreq(genind2genpop(g, quiet=T), quiet=T)[[1]] + a <- apply(g@tab,2,function(row) tapply(row, pops, mean, na.rm=TRUE)) HpS <- sum(1 - apply(a^2, 1, sum, na.rm=TRUE)) / n Hs_est <- (2*harmN/(2*harmN-1))*HpS HpT <- 1 - sum(apply(a,2,mean, na.rm=TRUE)^2)
2 R/D_Jost.R~
 @@ -38,7 +38,7 @@ D_Jost <- function(x){ global_Hs <- mean(loci[,1], na.rm=T) global_Ht <- mean(loci[,2], na.rm=T) global_D <- (global_Ht - global_Hs)/(1 - global_Hs ) * (n/(n-1)) - harm_D <- harmonic_mean(loci) + harm_D <- harmonic_mean(loci[,3]) return(list("per.locus"=loci[,3], "global.het"=global_D, "global.harm_mean" = harm_D
3 R/Gst_Hedrick.R
 @@ -23,9 +23,10 @@ Gst_Hedrick <- function(x){ n <- length(unique(pop(x))) harmN <- harmonic_mean(table(pop(x))) + pops <- pop(x) D.per.locus <- function(g) { #what we need to calculate these stats - a <- makefreq(genind2genpop(g, quiet=T), quiet=T)[[1]] + a <- apply(g@tab,2,function(row) tapply(row, pops, mean, na.rm=TRUE)) HpS <- sum(1 - apply(a^2, 1, sum, na.rm=TRUE)) / n Hs_est <- (2*harmN/(2*harmN-1))*HpS HpT <- 1 - sum(apply(a,2,mean, na.rm=TRUE)^2)
2 R/Gst_Hedrick.R~
 @@ -3,7 +3,7 @@ #' This function calculates Hedrick's G'st from a genind object #' #' Takes a genind object with population information and calculates Hedrick's -#' G'st. This Returns a list with values for each locus as well as a global estimates +#' G''st. This Returns a list with values for each locus as well as a global estimates #' #' Because estimators of Hs and Ht are used, it's possible to have negative #' estimates of Gst. You should treat such results as zeros (or estimating a
5 R/Gst_Nei.R
 @@ -19,9 +19,10 @@ Gst_Nei <- function(x){ n <- length(unique(pop(x))) harmN <- harmonic_mean(table(pop(x))) + pops <- pop(x) Gst.per.locus <- function(g) { - #what we need to calculate these stats - a <- makefreq(genind2genpop(g, quiet=T), quiet=T)[[1]] + #what we need to calculate these + a <- apply(g@tab,2,function(row) tapply(row, pops, mean, na.rm=TRUE)) HpS <- sum(1 - apply(a^2, 1, sum, na.rm=TRUE)) / n Hs_est <- (2*harmN/(2*harmN-1))*HpS HpT <- 1 - sum(apply(a,2,mean, na.rm=TRUE)^2)
3 R/diff_stats.R
 @@ -26,9 +26,10 @@ diff_stats <- function(x){ n <- length(unique(pop(x))) harmN <- harmonic_mean(table(pop(x))) + pops <- pop(x) per.locus <- function(g) { #what we need to calculate these stats - a <- makefreq(genind2genpop(g, quiet=T), quiet=T)[[1]] + a <- apply(g@tab,2,function(row) tapply(row, pops, mean, na.rm=TRUE)) HpS <- sum(1 - apply(a^2, 1, sum, na.rm=TRUE)) / n Hs_est <- (2*harmN/(2*harmN-1))*HpS HpT <- 1 - sum(apply(a,2,mean, na.rm=TRUE)^2)
6 R/diff_stats.R~
 @@ -14,18 +14,22 @@ #' diff_stats(nancycats) #' @references #' Hedrick, PW. (2005), A Standardized Genetic Differentiation Measure. Evolution 59: 1633-1638. +#' @references #' Jost, L. (2008), GST and its relatives do not measure differentiation. Molecular Ecology, 17: 4015-4026. +#' @references #' Nei M. (1973) Analysis of gene diversity in subdivided populations. PNAS: 3321-3323. +#' @references #' Nei M, Chesser RK. (1983). Estimation of fixation indices and gene diversities. Annals of Human Genetics. 47: 253-259. #' @family diffstat diff_stats <- function(x){ n <- length(unique(pop(x))) harmN <- harmonic_mean(table(pop(x))) + pops <- pop(x) per.locus <- function(g) { #what we need to calculate these stats - a <- makefreq(genind2genpop(g, quiet=T), quiet=T)[[1]] + a <- apply(g@tab,2,function(row) tapply(row, pops, mean, na.rm=TRUE)) HpS <- sum(1 - apply(a^2, 1, sum, na.rm=TRUE)) / n Hs_est <- (2*harmN/(2*harmN-1))*HpS HpT <- 1 - sum(apply(a,2,mean, na.rm=TRUE)^2)
2 R/summarise_bootstrap.R
 @@ -34,7 +34,7 @@ summarise_bootsrap <- function(bs, statistic){ res\$global.harm <- unlist(stats[3,]) } summarise <- function(x){ - return(c(mean=mean(x), quantile(x, c(0.025, 0.975)))) + return(c(mean=mean(x), quantile(x, c(0.025, 0.975), na.rm=TRUE))) } res\$summary.loci <- apply(loc_stats, 2, summarise) res\$summary.global.het <- summarise(res\$global.het)
5 R/summarise_bootstrap.R~
 @@ -24,9 +24,8 @@ summarise_bootsrap <- function(bs, statistic){ nreps <- length(bs) stats <- sapply(bs, statistic) - loc_stats <- matrix(unlist(stats[1,]), nrow=nreps, - dimnames=list(paste("rep", 1:nreps, sep=""), bs[[1]]@loc.names) - ) + loc_stats <- do.call(rbind, stats["per.locus",]) + res <-list("per.locus"= loc_stats, "global.het"=unlist(stats[2,]) )