betas.Rd
yasfoster opened this issue · 8 comments
Hi Jerome,
I get an unused argument error for (betaij = TRUE) when trying to run the below script. My input data works for basic.stats(dat) and wc(dat). It also works for betas() when I don't include the betaijT = TRUE, but I would like coancestries/kinships and inbreeding coefficients in matrix form as per Weir and Goudet 2017.
ind.coan<-betas(cbind(1:120,dat[,-1]),betaij=TRUE)
Error in betas(cbind(1:120, dat[, -1]), betaij = TRUE) :
unused argument (betaij = TRUE)
According to R, the package is up-to-date. Am I missing a step or is there something else I should try?
Many thanks,
Yasmin
Hi Jasmin,
can you run the example of the betas function? e.g.:
dat<-sim.genot(size=20,nbloc=100,nbal=20,mig=0.01,f=c(0,0.3,0.7))
ind.coan<-betas(cbind(1:60,dat[,-1]),betaijT=TRUE)
what version of hierfstat do you have installed (current is 0.04-30, you can get it from sessionInfo() )
Cheers
Hi Jerome,
I'm using version 0.04-30. This is the output after running the betas function:
dat<-sim.genot(size=20,nbloc=100,nbal=20,mig=0.01,f=c(0,0.3,0.7))
ind.coan<-betas(cbind(1:60,dat[,-1]),betaijT=TRUE)
Error in betas(cbind(1:60, dat[, -1]), betaijT = TRUE) :
unused argument (betaijT = TRUE)`
Cheers
Yasmin
weird... could you paste here the results of typing betas (without parentheses)? Thanks
{
ratio.Hi.Hb <- function(x) {
dum <- which(!is.na(x))
sum(x[dum])/sum(Hb[dum])
}
if (is.genind(dat))
dat <- genind2hierfstat(dat)
pfr <- pop.freq(dat, diploid)
pfr2 <- lapply(pfr, function(x) t(x) %*% x)
nl <- dim(dat)[2] - 1
np <- length(table(dat[, 1]))
ns <- ind.count.n(dat)
if (diploid)
ns <- ns * 2
ns[is.na(ns)] <- 0
Hi <- matrix(numeric(np * nl), ncol = nl)
Hb <- numeric(nl)
for (il in 1:nl) {
Hi[, il] <- ns[, il]/(ns[, il] - 1) * (1 - diag(pfr2[[il]]))
npl <- sum(ns[, il] > 0, na.rm = TRUE)
Hb[il] <- 1 - 1/npl/(npl - 1) * sum((pfr2[[il]] - diag(diag(pfr2[[il]]))),
na.rm = TRUE)
}
betai <- 1 - apply(Hi, 1, ratio.Hi.Hb)
if (nboot < 100) {
return(list(betaiovl = betai, Hi = Hi, Hb = Hb))
}
else {
if (nl < 10) {
warning("Less than 10 loci, can't estimate Conf. Int.")
return(list(betaiovl = betai, Hi = Hi, Hb = Hb))
}
boot.bi <- matrix(numeric(nboot * np), nrow = nboot)
nls <- apply(ns, 1, function(x) which(x > 0))
if (is.matrix(nls)) {
for (ib in 1:nboot) {
for (ip in 1:np) {
dum <- sample(nls[, ip], replace = TRUE)
boot.bi[ib, ip] <- 1 - sum(Hi[ip, dum])/sum(Hb[dum])
}
}
}
else {
for (ib in 1:nboot) {
for (ip in 1:np) {
dum <- sample(nls[[ip]], replace = TRUE)
boot.bi[ib, ip] <- 1 - sum(Hi[ip, dum])/sum(Hb[dum])
}
}
}
bi.ci <- apply(boot.bi, 2, quantile, lim, na.rm = TRUE)
return(list(betaiovl = betai, ci = bi.ci, Hi = Hi, Hb = Hb))
}
}
<bytecode: 0x7f95040da4b8>
<environment: namespace:hierfstat>
Cheers
Yasmin
Hmmm, it is strange, as this version of betas is an old one. The one you should have installed is this
https://github.com/jgx65/hierfstat/blob/master/R/betas.R
can you show the results of sessionInfo() please? I am puzzled that you have this old version of betas.
Platform: x86_64-apple-darwin15.6.0 (64-bit)
Running under: macOS Sierra 10.12.6
Matrix products: default
BLAS: /System/Library/Frameworks/Accelerate.framework/Versions/A/Frameworks/vecLib.framework/Versions/A/libBLAS.dylib
LAPACK: /Library/Frameworks/R.framework/Versions/3.6/Resources/lib/libRlapack.dylib
Random number generation:
RNG: Mersenne-Twister
Normal: Inversion
Sample: Rounding
locale:
[1] en_NZ.UTF-8/en_NZ.UTF-8/en_NZ.UTF-8/C/en_NZ.UTF-8/en_NZ.UTF-8
attached base packages:
[1] stats graphics grDevices utils datasets methods base
other attached packages:
[1] devtools_2.3.0 usethis_1.6.1 ggpubr_0.4.0 ggplot2_3.3.2 readr_1.3.1
[6] hierfstat_0.04-30 pegas_0.13 adegenet_2.1.3 ade4_1.7-15 ape_5.4
loaded via a namespace (and not attached):
[1] colorspace_1.4-1 ggsignif_0.6.0 seqinr_3.6-1 deldir_0.1-25 ellipsis_0.3.1
[6] class_7.3-15 rio_0.5.16 rprojroot_1.3-2 fs_1.4.1 rstudioapi_0.11
[11] farver_2.0.3 remotes_2.1.1 fansi_0.4.1 codetools_0.2-16 splines_3.6.3
[16] pkgload_1.1.0 broom_0.5.6 cluster_2.1.0 shiny_1.5.0 compiler_3.6.3
[21] backports_1.1.8 assertthat_0.2.1 Matrix_1.2-18 fastmap_1.0.1 cli_2.0.2
[26] later_1.1.0.1 prettyunits_1.1.1 htmltools_0.5.0 tools_3.6.3 igraph_1.2.5
[31] coda_0.19-3 gtable_0.3.0 glue_1.4.1 reshape2_1.4.4 dplyr_1.0.0
[36] gmodels_2.18.1 Rcpp_1.0.4 carData_3.0-4 cellranger_1.1.0 raster_3.3-7
[41] vctrs_0.3.1 spdep_1.1-3 gdata_2.18.0 nlme_3.1-144 stringr_1.4.0
[46] ps_1.3.3 testthat_2.3.2 openxlsx_4.1.5 mime_0.9 lifecycle_0.2.0
[51] gtools_3.8.2 rstatix_0.6.0 LearnBayes_2.15.1 MASS_7.3-51.5 scales_1.1.1
[56] hms_0.5.3 promises_1.1.1 parallel_3.6.3 expm_0.999-4 curl_4.3
[61] memoise_1.1.0 stringi_1.4.6 desc_1.2.0 e1071_1.7-3 permute_0.9-5
[66] boot_1.3-24 pkgbuild_1.0.8 zip_2.0.4 spData_0.3.5 rlang_0.4.6
[71] pkgconfig_2.0.3 lattice_0.20-38 purrr_0.3.4 sf_0.9-4 labeling_0.3
[76] tidyselect_1.1.0 processx_3.4.2 plyr_1.8.6 magrittr_1.5 R6_2.4.1
[81] generics_0.0.2 DBI_1.1.0 pillar_1.4.4 haven_2.3.1 foreign_0.8-75
[86] withr_2.2.0 mgcv_1.8-31 units_0.6-7 abind_1.4-5 sp_1.4-2
[91] tibble_3.0.1 crayon_1.3.4 car_3.0-8 KernSmooth_2.23-16 grid_3.6.3
[96] readxl_1.3.1 data.table_1.12.8 callr_3.4.3 vegan_2.5-6 forcats_0.5.0
[101] digest_0.6.25 classInt_0.4-3 xtable_1.8-4 tidyr_1.1.0 httpuv_1.5.4
[106] munsell_0.5.0 sessioninfo_1.1.1 ```Could you try reinstalling hierfstat please? If this fails, you can copy and paste the code from the URL I sent previously, but I must admit I don't understand what is going on.
Hi Jerome
I redefined the betas function using the URL you sent. It then had an ind.count.n error so I also redefined that function with:
#should replace ind.count for ill behaved loci, e.g. many weird alleles
dum<-function(x){
a<-!is.na(x)
tapply(a,data[,1],sum)
}
if (is.genind(data)) data<-genind2hierfstat(data)
data[,1]<-factor(data[,1])
apply(data[,-1],2,dum)
}
It seems to work as long as I redefine betas and ind.count.n first.
Thanks for your help!
Yasmin