betas.Rd

[yasfoster]

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

[jgx65]

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

[yasfoster]

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

[jgx65]

weird... could you paste here the results of typing betas (without parentheses)? Thanks

[yasfoster]

{
    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

[jgx65]

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.

[yasfoster]

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 ```

[jgx65]

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.

[yasfoster]

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