R/sexbiaseddisp.R

Defines functions sexbias.test AIc

Documented in AIc sexbias.test

################################################################################
subsamp.within<-function (lev,ni=10) {
#subsample ni observations per levels of lev
    y <- 1:length(lev)
    nlev <- nlevels(factor(lev))
    nl <- as.integer(factor(lev))
    x <- list()
    for (i in 1:nlev) if (length(y[nl == i]) > ni) 
        x[[i]] <- sample(y[nl == i],size=ni)
    else x[[i]] <- y[nl == i]
    return(unlist(x))
}
##################################################################################
#' @title Calculates corrected Assignment Index
#' @description Calculates corrected Assignment Index as described in 
#' \href{https://onlinelibrary.wiley.com/doi/abs/10.1046/j.1365-294X.2002.01496.x}{Goudet etal. (2002)}
#' @usage AIc(dat)
#' @param dat a data frane with nlocs+1 columns, 
#' @return aic  The corrected assignment index of each individual
#' @author Jerome Goudet \email{jerome.goudet@@unil.ch}
#' @references 
#' \href{https://onlinelibrary.wiley.com/doi/abs/10.1046/j.1365-294X.2002.01496.x}{Goudet J, Perrin N, Waser P (2002)} 
#' Tests for sex-biased dispersal 
#'  using bi-parentally inherited genetic markers 11, 1103:1114
#' @export
##################################################################################
AIc<-function(dat){
dat[,1]<-as.integer(as.factor(dat[,1]))
popfreq<-pop.freq(dat) #freq per pop
indfreq<-pop.freq(cbind(1:dim(dat)[1],dat[,-1])) # freq per ind
id.al<-lapply(indfreq,function(y) apply(y,2,fun2<-function(z) which(z>0))) #which allele in which ind
nloc<-dim(dat)[2]-1 
npop<-length(table(dat[,1]))
nind<-dim(dat)[1]
ail<-matrix(numeric(nind*nloc),ncol=nloc)
for (il in 1:nloc){
  for (ii in 1:nind){
    dum<-as.vector(id.al[[il]][[ii]])
    if (length(dum)!=0){
      if (length(dum)==1)
        ail[ii,il]<-2*log(popfreq[[il]][[dum,dat[ii,1]]])
      else
        ail[ii,il]<-log(2*popfreq[[il]][[dum[1],dat[ii,1]]]*popfreq[[il]][[dum[2],dat[ii,1]]])
    }
    else ail[ii,il]<-NA
  }
}
#mean assignment per locus per pop
mean.ail<-apply(ail,2,function(x) {dum<-tapply(x,dat[,1],mean,na.rm=T);return(rep(dum,table(dat[,1])))})

nas<-which(is.na(ail))
ail[nas]<-mean.ail[nas] #set missing values to population mean
aic<-apply(ail-mean.ail,1,sum)
return(aic)
}
##################################################################################


#' @title Test for sex biased dispersal
#' @description Test whether one
#' sex disperses more than the other using the method described in 
#' \href{https://onlinelibrary.wiley.com/doi/abs/10.1046/j.1365-294X.2002.01496.x}{Goudet etal. (2002)} 
#' @usage sexbias.test(dat,sex,nperm=NULL,test="mAIc",alternative="two.sided")
#' @param dat a data frame with n.locs+1 columns and n.inds rows
#' @param sex a vector containing the individual's sex
#' @param nperm the number of permutation to carry out
#' @param test one of "mAIc" (default), "vAIc","FIS" or "FST"
#' @param alternative one of "two.sided" (default),"less" or "greater"
#' @return call the function call
#' @return res the observation for each sex
#' @return statistic the observed statistic for the chosen test
#' @return p.value the p-value of the hypothesis
#' @author Jerome Goudet \email{jerome.goudet@@unil.ch}
#' @references \href{https://onlinelibrary.wiley.com/doi/abs/10.1046/j.1365-294X.2002.01496.x}{Goudet J, Perrin N, Waser P (2002)} Tests for sex-biased dispersal 
#'  using bi-parentally inherited genetic markers 11, 1103:1114
#' @examples 
#'   data(crocrussula)
#'   sexbias.test(crocrussula$genot,crocrussula$sex)
#'   dat<-qn2.read.fstat(system.file("extdata","qn2_sex.dat",package="hierfstat"))
#'   sexbias.test(dat[[1]],sex=dat[[2]])
#'   \dontrun{
#'   sexbias.test(crocrussula$genot,crocrussula$sex,nperm=1000)
#'   sexbias.test(dat[[1]],sex=dat[[2]],nperm=100,test="FST",alternative="greater")
#'   }
#' @export
###################################################################################
sexbias.test<-function(dat,sex,nperm=NULL,test="mAIc",
                        alternative="two.sided"){
#calculate the corrected assignment index for all individuals
#dat is a genotype data frame
#sex is a vector of sex labels
#nperm is the number of permutations to be carried out in the permutation test
#2 tailed test on mean AIc for the time being
#TODO: 
#1. options for one tail tests
#2. include test on the variance and Fis / Fst (call wc)

  call<-match.call()
  if(!is.na(pmatch(test,"mAIc")))
    test<-"mAIc"
  TESTS<-c("mAIc","vAIc","FIS","FST")
  test<-pmatch(test,TESTS)
  if (is.na(test)) stop("Invalid test")
  if (test==-1) stop("Ambiguous test")
  
  if(!is.na(pmatch(alternative,"two.sided")))
    alternative<-"two.sided"
  ALTS<-c("two.sided","less","greater")
  alt<-pmatch(alternative,ALTS)
  if(is.na(alt)) stop("Invalid alternative")
  if(alt==-1) stop("Ambiguous alternative")
  

sexes<-names(table(sex))
sex1<-which(sex==sexes[1])
sex2<-which(sex==sexes[2])

if (test==3 | test==4){
  a<-wc(dat[sex1,])
  b<-wc(dat[sex2,])
  FIS<-c(a$FIS,b$FIS)
  names(FIS)<-sexes
  FST<-c(a$FST,b$FST)
  names(FST)<-sexes
}  

if (test==1 | test==2){
aic<-AIc(dat)
mAic1<-mean(aic[sex1])
mAic2<-mean(aic[sex2])

vAic1<-stats::var(aic[sex1])
vAic2<-stats::var(aic[sex2])

mAics<-data.frame(mAic1,mAic2)
names(mAics)<-sexes
vAics<-data.frame(vAic1,vAic2)
names(vAics)<-sexes
}

if (is.null(nperm)){
  if (test==1){
  tmAic<-stats::t.test(aic[sex1],aic[sex2],alternative=alternative)
  res<-list(statistic=tmAic$statistic,p.value=tmAic$p.value)
  }
  if (test==2){
    vAic<-stats::var.test(aic[sex1],aic[sex2],alternative=alternative)
    res<-list(statistic=vAic$statistic,p.value=vAic$p.value)
  }
  if (test==3 | test==4) stop("You must specify nperm for this test. Exiting")
}
else{
res<-list()
if (test==1){
tstats<-numeric(nperm)
tstats[1]<-stats::t.test(aic[sex1],aic[sex2])$statistic

foo<-function(){
dum<-samp.within(dat[,1])
stats::t.test(aic[dum][sex1],aic[dum][sex2])$statistic
}
tstats[-1]<-replicate(nperm-1,foo())
res$statistic<-tstats[1]
if (alt==1){res$p.value<-sum(abs(tstats)>=abs(tstats[1]))/nperm}
if (alt==2){res$p.value<-sum(tstats<=tstats[1])/nperm}
if (alt==3){res$p.value<-sum(tstats>=tstats[1])/nperm}
}

if (test==2){
  fstats<-numeric(nperm)
  fstats[1]<-stats::var.test(aic[sex1],aic[sex2])$statistic
  
  foo<-function(){
    dum<-samp.within(dat[,1])
    return(stats::var.test(aic[dum][sex1],aic[dum][sex2])$statistic)
  }
  fstats[-1]<-replicate(nperm-1,foo())
  res$statistic<-fstats[1]
  if (alt==1){res$p.value<-sum(abs(fstats)>=abs(fstats[1]))/nperm}
  if (alt==2){res$p.value<-sum(fstats<=fstats[1])/nperm}
  if (alt==3){res$p.value<-sum(fstats>=fstats[1])/nperm}
  
}
if (test==3){
  fstats<-numeric(nperm)
  fstats[1]<-FIS[1]-FIS[2]
  foo<-function(){
    dum<-samp.within(dat[,1])
    d1<-dum[sex1]
    d2<-dum[sex2]
    fis1<-wc(dat[d1,])$FIS
    fis2<-wc(dat[d2,])$FIS
    return(fis1-fis2)
  }
  fstats[-1]<-replicate(nperm-1,foo())
  res$statistic<-fstats[1]
  if (alt==1){res$p.value<-sum(abs(fstats)>=abs(fstats[1]))/nperm}
  if (alt==2){res$p.value<-sum(fstats<=fstats[1])/nperm}
  if (alt==3){res$p.value<-sum(fstats>=fstats[1])/nperm}
  
}
if (test==4){
  fstats<-numeric(nperm)
  fstats[1]<-FST[1]-FST[2]
  foo<-function(){
    dum<-samp.within(dat[,1])
    d1<-dum[sex1]
    d2<-dum[sex2]
    fst1<-wc(dat[d1,])$FST
    fst2<-wc(dat[d2,])$FST
    return(fst1-fst2)
  }
  fstats[-1]<-replicate(nperm-1,foo())
  res$statistic<-fstats[1]
  if (alt==1){res$p.value<-sum(abs(fstats)>=abs(fstats[1]))/nperm}
  if (alt==2){res$p.value<-sum(fstats<=fstats[1])/nperm}
  if (alt==3){res$p.value<-sum(fstats>=fstats[1])/nperm}
} 
  
}
return(list(call=match.call(),statistic=res$statistic,p.value=res$p.value))

}

 
jgx65/hierfstat documentation built on April 20, 2023, 8:34 a.m.