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
jgx65 · Feb 13, 2023 · ebb48fe · ebb48fe
1 parent ad81c6a
commit ebb48fe
Showing 1 changed file with 91 additions and 74 deletions.
diff --git 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