non-numeric argument error in boot.ppfis

[jaflury]

Dear Jerome,

I had the same issue with the population name in basic.stats and it was resolved with the updated function, however I get the same error in boot.ppfis (error message "data[dim(data)[1], 1] + 1 : non-numeric argument to binary operator"). Unfortunately I could not figure out what needs to be changed, could you help me there?

Best, Jana

Originally posted by @jaflury in https://github.com/jgx65/hierfstat/issues/65#issuecomment-1743122646

[jgx65]

Dear Jana,

Could you post a reproducible example please?

A quick fix is to duplicate your unique population:

#load gtrunchier dataset
data(gtrunchier)
#extract individuals from 1st pop
dum<-gtrunchier[gtrunchier[,1]==1,-2]
#produces an error
boot.ppfis(dum)
#duplicates dum
dum1<-dum
dum1[,1]<-2
dumt<-rbind(dum,dum1)
#works, repeating twice the results
boot.ppfis(dumt)

[jgx65]

@jaflury , have you managed to solve the problem?

[jaflury]

Dear Jerome, sorry for the late reply. Here's a subset of the data I used. I tried your solution:

subset <- read.table("dat.txt")

dum<-subset[subset[,1]==1,-2]

produces an error

boot.ppfis(dum)

duplicates dum

dum1<-dum dum1[,1]<-2 dumt<-rbind(dum,dum1)

works, repeating twice the results

boot.ppfis(dumt)

but it didn't work. I know it's coming from the fact that there is only one population, but I don't find a way to fix it. It was fixed in this basic.stat function:

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 }

and this function works with the dat.txt. basic.stats(subset)

Thank you very much for your help!

Best, Jana

[jaflury]

and here's the file: dat.txt

[jgx65]

Dear @jaflury , sorry for my delayed answer. the following works for me:

library(hierfstat)
dat<-read.table("dat.txt",header=T)
tmp<-dat
tmp[,1]<-rep("pop2",dim(tmp)[1])
dat2<-rbind(dat,tmp)
boot.ppfis(dat2)

and I obtain:

$fis.ci
      ll     hl
1 0.0863 0.3339
2 0.0863 0.3339

Let me know if this solves the issue. I will look into a more permanent solution when I find time. Many thanks!