#' Estimate expected heterozygosity
#'
#' Returns the general expected total heterozygosity parameter
#' @param x A \code{data.frame} object with \code{locus} objects
#' @param stratum The name of the column representing the stratum variable (default=Population)
#' @export
#' @author Rodney J. Dyer \email{rjdyer@@vcu.edu}
#' @examples
#' loci1 <- c( locus( c("A","A") ), locus( c("A","A") ), locus( c("A","B")))
#' loci2 <- c( locus( c("A","A") ), locus( c("A","B")), locus(c("B","B")))
#' df <- data.frame( Population=c("One","One","One","Two","Two","Two"), Locus=c(loci1,loci2) )
#' Ht( df )
Ht <- function( x, stratum="Population" ) {
if( missing(x) || !(stratum %in% names(x)))
stop("You need to pass both a data.frame and the name of the stratum to this function")
locus_names <- column_class(x,class="locus")
if( length(locus_names)==0)
stop("Cannot estimate expected total heterozygosity if there are no loci...")
ret <- data.frame( Locus=locus_names, Ht=0 )
ho <- Hos( x, stratum=stratum )
hs <- Hes( x, stratum=stratum )
freqs <- frequencies( x, stratum=stratum)
cts <- genotype_counts(x, stratum)
K <- nrow(cts)
for( locus in locus_names ) {
nbar <- harmonic_mean(cts[[locus]])
xki <- freqs[ freqs$Locus==locus, ]
xibar <- unlist(by( xki$Frequency, xki$Allele, function(x) sum(x/K )))
ht <- 1 - sum(xibar^2) + hs$Hes[hs$Locus==locus]/(K*nbar) - ho$Hos[ho$Locus==locus]/(K*nbar*2)
ret$Ht[ ret$Locus == locus ] <- ht
}
return( ret )
}
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.