From ebb48fe1f2ad600c16cee9958a353a3b71e6e687 Mon Sep 17 00:00:00 2001 From: Jerome Goudet Date: Mon, 13 Feb 2023 18:16:17 +0100 Subject: [PATCH] Update misc.R in basic.stats, renamed the dummy pop when only one sample as DumPop. Cleaned up function output so that Dumpop is not showing --- R/misc.R | 165 ++++++++++++++++++++++++++++++------------------------- 1 file changed, 91 insertions(+), 74 deletions(-) diff --git a/R/misc.R b/R/misc.R index fa9109e..c7b7685 100644 --- a/R/misc.R +++ b/R/misc.R @@ -228,80 +228,97 @@ allelic.richness<-function (data, min.n = NULL, diploid = TRUE) ################ #' @export ################ -basic.stats<-function(data,diploid=TRUE,digits=4){ -# TODO : define plot functions for basic.stats -if (is.genind(data)) data<-genind2hierfstat(data) -if (length(table(data[,1]))<2) data[dim(data)[1]+1,1]<-data[dim(data)[1],1]+1 -if(dim(data)[2]==2) data<-data.frame(data,dummy.loc=data[,2]) -loc.names<-names(data)[-1] -p<-pop.freq(data,diploid) -n<-t(ind.count(data)) -if (diploid){ -dum<-getal.b(data[,-1]) -Ho<-dum[,,1]==dum[,,2] -#sHo<-(n-t(apply(Ho,2,fun<-function(x) tapply(x,data[,1],sum,na.rm=TRUE))))/n -sHo<-(1-t(apply(Ho,2,fun<-function(x) tapply(x,data[,1],mean,na.rm=TRUE)))) - -mHo<-apply(sHo,1,mean,na.rm=TRUE) -} -else { -sHo<-NA -mHo<-NA -} -sp2<-lapply(p,fun<-function(x) apply(x,2,fun2<-function(x) sum(x^2))) -sp2<-matrix(unlist(sp2),nrow=dim(data[,-1])[2],byrow=TRUE) -if (diploid) { -Hs<-(1-sp2-sHo/2/n) -Hs<-n/(n-1)*Hs -Fis=1-sHo/Hs -} -else { -Hs<-n/(n-1)*(1-sp2) -Fis<-NA -} -np<-apply(n,1,fun<-function(x) sum(!is.na(x))) -mn<-apply(n,1,fun<-function(x) {np<-sum(!is.na(x)); np/sum(1/x[!is.na(x)])}) -msp2<-apply(sp2,1,mean,na.rm=TRUE) -mp<-lapply(p,fun<-function(x) apply(x,1,mean,na.rm=TRUE)) -mp2<-unlist(lapply(mp,fun1<-function(x) sum(x^2))) -if (diploid){ -mHs<-mn/(mn-1)*(1-msp2-mHo/2/mn) -Ht<-1-mp2+mHs/mn/np-mHo/2/mn/np -mFis=1-mHo/mHs -} -else{ -mHs<-mn/(mn-1)*(1-msp2) -Ht<-1-mp2+mHs/mn/np -mFis<-NA -} -Dst<-Ht-mHs -Dstp<-np/(np-1)*Dst -Htp=mHs+Dstp -Fst=Dst/Ht -Fstp=Dstp/Htp -Dest<-Dstp/(1-mHs)#*np/(np-1) -res<-data.frame(cbind(mHo,mHs,Ht,Dst,Htp,Dstp,Fst,Fstp,mFis,Dest)) -names(res)<-c("Ho","Hs","Ht","Dst","Htp","Dstp","Fst","Fstp","Fis","Dest") -if (diploid){ -rownames(sHo)<-loc.names -rownames(Fis)<-loc.names -} -is.na(res)<-do.call(cbind,lapply(res,is.infinite)) -overall<-apply(res,2,mean,na.rm=TRUE) -overall[7]<-overall[4]/overall[3] -overall[8]<-overall[6]/overall[5] -overall[9]<-1-overall[1]/overall[2] -overall[10]<-overall[6]/(1-overall[2]) -names(overall)<-names(res) -if(!diploid){ -# res[,-2]<-NA - overall[c(1,9)]<-NA -} -all.res<-list(n.ind.samp=n,pop.freq=lapply(p,round,digits), -Ho=round(sHo,digits),Hs=round(Hs,digits),Fis=round(Fis,digits), -perloc=round(res,digits),overall=round(overall,digits)) -class(all.res)<-"basic.stats" -all.res +basic.stats<-function (data, diploid = TRUE, digits = 4) +{ + if (adegenet::is.genind(data)) + data <- genind2hierfstat(data) + if (length(table(data[, 1])) < 2){ + data[dim(data)[1] + 1, 1] <- "DumPop" + dum.pop<-TRUE + } + if (dim(data)[2] == 2) + data <- data.frame(data, dummy.loc = data[, 2]) + loc.names <- names(data)[-1] + p <- pop.freq(data, diploid) + n <- t(ind.count(data)) + if (diploid) { + dum <- getal.b(data[, -1]) + Ho <- dum[, , 1] == dum[, , 2] + sHo <- (1 - t(apply(Ho, 2, fun <- function(x) tapply(x, + data[, 1], mean, na.rm = TRUE)))) + mHo <- apply(sHo, 1, mean, na.rm = TRUE) + } + else { + sHo <- NA + mHo <- NA + } + sp2 <- lapply(p, fun <- function(x) apply(x, 2, fun2 <- function(x) sum(x^2))) + sp2 <- matrix(unlist(sp2), nrow = dim(data[, -1])[2], byrow = TRUE) + if (diploid) { + Hs <- (1 - sp2 - sHo/2/n) + Hs <- n/(n - 1) * Hs + Fis = 1 - sHo/Hs + } + else { + Hs <- n/(n - 1) * (1 - sp2) + Fis <- NA + } + np <- apply(n, 1, fun <- function(x) sum(!is.na(x))) + mn <- apply(n, 1, fun <- function(x) { + np <- sum(!is.na(x)) + np/sum(1/x[!is.na(x)]) + }) + msp2 <- apply(sp2, 1, mean, na.rm = TRUE) + mp <- lapply(p, fun <- function(x) apply(x, 1, mean, na.rm = TRUE)) + mp2 <- unlist(lapply(mp, fun1 <- function(x) sum(x^2))) + if (diploid) { + mHs <- mn/(mn - 1) * (1 - msp2 - mHo/2/mn) + Ht <- 1 - mp2 + mHs/mn/np - mHo/2/mn/np + mFis = 1 - mHo/mHs + } + else { + mHs <- mn/(mn - 1) * (1 - msp2) + Ht <- 1 - mp2 + mHs/mn/np + mFis <- NA + } + Dst <- Ht - mHs + Dstp <- np/(np - 1) * Dst + Htp = mHs + Dstp + Fst = Dst/Ht + Fstp = Dstp/Htp + Dest <- Dstp/(1 - mHs) + res <- data.frame(cbind(mHo, mHs, Ht, Dst, Htp, Dstp, Fst, + Fstp, mFis, Dest)) + names(res) <- c("Ho", "Hs", "Ht", "Dst", "Htp", "Dstp", + "Fst", "Fstp", "Fis", "Dest") + if (diploid) { + rownames(sHo) <- loc.names + rownames(Fis) <- loc.names + } + is.na(res) <- do.call(cbind, lapply(res, is.infinite)) + overall <- apply(res, 2, mean, na.rm = TRUE) + overall[7] <- overall[4]/overall[3] + overall[8] <- overall[6]/overall[5] + overall[9] <- 1 - overall[1]/overall[2] + overall[10] <- overall[6]/(1 - overall[2]) + names(overall) <- names(res) + if (!diploid) { + overall[c(1, 9)] <- NA + } + if(dum.pop){ + ToBeRemoved<-which(dimnames(Hs)[[2]]=="DumPop") + n<-n[,-ToBeRemoved,drop=FALSE] + Hs<-Hs[,-ToBeRemoved,drop=FALSE] + if (diploid) sHo<-sHo[,-ToBeRemoved,drop=FALSE] + Fis<-Fis[,-ToBeRemoved,drop=FALSE] + p<-lapply(p,function(x) x[,-ToBeRemoved,drop=FALSE]) + } + all.res <- list(n.ind.samp = n, pop.freq = lapply(p, round, + digits), Ho = round(sHo, digits), Hs = round(Hs, digits), + Fis = round(Fis, digits), perloc = round(res, digits), + overall = round(overall, digits)) + class(all.res) <- "basic.stats" + all.res } #' @method print basic.stats