jgx65 / hierfstat

the hierfstat package
24 stars 14 forks source link

problem performing sex.biased.test() function #59

Closed danielsangarci closed 2 years ago

danielsangarci commented 2 years ago

Hello,

I am trying to perform the test with the function sex.biased.test(), but I am getting an error I don't know how to solve.

When I perform the test with my whole dataset (2 metapopulations with 3 patches per metapopulation) it works.

microsatelites_dataset_example.xlsx

sexbias.test(dat=microsat_df[,c(2,4:25)],sex=microsat_df$sex)

But it actually doesn't make sense for me due to I need to analyse both population separately.

When I subset the dataset it doesn't work anymore and I get an error:

microsat_A_subset <- microsat_df$population == "A" sexbias.test(dat=microsat_df[microsat_A_subset,c(2,4:25)],sex=microsat_df[microsat_A_subset,]$sex)

Error in if (stderr < 10 * .Machine$double.eps * max(abs(mx), abs(my))) stop("data are essentially constant") : missing value where TRUE/FALSE needed

Do you know what's wrong here? I have rewrite the code in different ways and I cannot find the mistake.

Thank you very much in advance, Dani

jgx65 commented 2 years ago

Hi Dani,

The issue is probably how you've loaded the data set. You have separate columns for each allele, whereas hierfstat expects the two alleles at a locus to be collated together. You can generate the correct format using

dat<-read.csv("microsatelites_dataset_example.csv",sep=",",header=TRUE)
genot<-dat[,4:25]
genotb<-genot[,(0:10)*2+1]*1000+genot[,(0:10)*2+2]
genotb[genotb==0]<-NA

and then

library(hierfstat)
sexbias.test(data.frame(dat$meadow,genotb),dat$sex)
x<-which(dat$population=="A")
sexbias.test(data.frame(dat$meadow,genotb)[x,],dat$sex[x])

should give you meaningful results:

> sexbias.test(data.frame(dat$meadow,genotb)[x,],dat$sex[x])
$call
sexbias.test(dat = data.frame(dat$meadow, genotb)[x, ], sex = dat$sex[x])

$statistic
        t 
0.3829125 

$p.value
[1] 0.702486
danielsangarci commented 2 years ago

Thank you very much for your help.

All the best, Dani

El mié, 25 may 2022 a las 17:40, Jerome Goudet @.***>) escribió:

Hi Dani,

The issue is probably how you've loaded the data set. You have separate column for each allele, whereas hierfstat expects the two alleles at a locus to be collated together. You can generate the correct format using

dat<-read.csv("microsatelites_dataset_example.csv",sep=",",header=TRUE) genot<-dat[,4:25] genotb<-genot[,(0:10)2+1]1000+genot[,(0:10)*2+2] genotb[genotb==0]<-NA

and then

library(hierfstat) sexbias.test(data.frame(dat$meadow,genotb),dat$sex) x<-which(dat$population=="A") sexbias.test(data.frame(dat$meadow,genotb)[x,],dat$sex[x])

should give you meaningful results:

sexbias.test(data.frame(dat$meadow,genotb)[x,],dat$sex[x]) $call sexbias.test(dat = data.frame(dat$meadow, genotb)[x, ], sex = dat$sex[x])

$statistic t 0.3829125

$p.value [1] 0.702486

— Reply to this email directly, view it on GitHub https://github.com/jgx65/hierfstat/issues/59#issuecomment-1137453198, or unsubscribe https://github.com/notifications/unsubscribe-auth/AZKV25YDUQSE5PSPCWEQS6DVLZCXPANCNFSM5W5P3EFA . You are receiving this because you authored the thread.Message ID: @.***>

danielsangarci commented 2 years ago

Thank you very much for your help.

Dani