TracePst {Pstat}R Documentation

Pst variations in function of c/h^2

Description

'TracePst' plots the curves of the functions that map c/h^2 onto Pst (for chosen quantitative measures). Indeed, Pst depends on the value of c/h^2, where c is the assumed additive genetic proportion of differences between populations and where h^2 is (narrow-sense heritability) the assumed additive genetic proportion of differences between individuals within populations.

Usage

TracePst(data,va=0,ci=1,boot=1000,pe=0.95,Fst=-1,Pw=0,Rp=0,Ri=0,xm=2,pts=30)

Arguments

data

a dataframe with as many rows as individuals. The first column contains the name of the population to which the individual belongs, the others contain quantitative variables.

va

a vector containing the selected variables names or numbers (i.e. those of the quantitative measures considered). If va=0 all the variables are selected.

ci

if ci=1 the confidence interval of Pst is plotted.

boot

the number of data frames generated to determine the confidence interval or to construct the dotted lines representing this confidence interval (using the bootstrap method).

pe

the confidence level of the calculated interval.

Fst

the value of Wright's Fst, if avalaible.

Pw

a vector containing the names of the two populations considered to obtain pairwise Pst.

Rp

a vector containing the names of the populations to be deleted.

Ri

a vector containing each number of individual to be deleted. The vector Ri must contain existent individuals, each of them once.

xm

the maximum on x-axis (values of c/h^2).

pts

number of points used to plot the curves.

Value

In any case, the sizes of each population considered. The expected curves.

Note

The time required to construct the dotted lines associated with the confidence intervals might be fairly long depending on the user choices.

Author(s)

Blondeau Da Silva Stephane - Da Silva Anne.

References

Brommer J., 2011. Whither Pst? The approximation of Qst by Pst in evolutionary and conservation biology. Journal of Evolutionary Biology, 24:1160-1168.

Lima M.R. et al., 2012. Genetic and Morphometric Divergence of an Invasive Bird: The Introduced House Sparrow (Passer domesticus) in Brazil. PloS One 7 (12).

On Fst : Wright S., 1951. The genetical structure of populations. Annals of Eugenics 15, 323-354.

Examples


data(test)
# TracePst(test)
# TracePst(test,boot=2000,va="QM7",Ri=18,pe=0.9,pts=40,xm=4)
TracePst(test,va=7:10,Fst=0.3,Ri=c(22,27,195),Rp=c("A","C","E"),ci=0)
# TracePst(test,va="QM1",Ri=c(3,7:17),Pw=c("C","D"),pts=20)

## The function is currently defined as
function (data, va = 0, ci = 1, boot = 1000, pe = 0.95, Fst = -1, 
    Pw = 0, Rp = 0, Ri = 0, xm = 2, pts = 30) 
{
    nonNa.clm <- function(data, clm) {
        nb.ind = dim(data)[1]
        nb.na = 0
        for (i in 1:nb.ind) if (is.na(data[i, clm])) 
            nb.na = nb.na + 1
        return(nb.ind - nb.na)
    }
    dat.fra.prep <- function(data) {
        nb.var = dim(data)[2] - 1
        data = as.data.frame(data)
        data[, 1] = as.character(data[, 1])
        for (i in 1:nb.var) {
            if (is.numeric(data[, i + 1]) == FALSE) 
                data[, i + 1] = as.numeric(as.character(data[, 
                  i + 1]))
        }
        dat.sta <- function(dat) {
            nb.vari = dim(dat)[2] - 1
            st.dev = rep(0, nb.vari)
            mea = rep(0, nb.vari)
            for (i in 1:nb.vari) {
                nna.clm = nonNa.clm(dat, i + 1)
                st.dev[i] = sqrt((nna.clm - 1)/nna.clm) * sd(dat[, 
                  i + 1], na.rm = TRUE)
                mea[i] = mean(dat[, i + 1], na.rm = TRUE)
            }
            for (j in 1:nb.vari) dat[, j + 1] = (dat[, j + 1] - 
                mea[j])/st.dev[j]
            return(dat)
        }
        data = dat.sta(data)
        return(data)
    }
    dat.rem.ind.pop <- function(data, ind = 0, pop = 0) {
        data = as.data.frame(data)
        dat.rem.ind <- function(dat, ind) {
            nb.rem.ind = length(ind)
            nb.ind = dim(dat)[1]
            for (i in 1:nb.rem.ind) dat = dat[row.names(dat)[1:(nb.ind - 
                i + 1)] != ind[i], ]
            return(dat)
        }
        dat.rem.pop <- function(dat, pop) {
            nb.rem.pop = length(pop)
            for (i in 1:nb.rem.pop) dat = dat[dat[, 1] != pop[i], 
                ]
            return(dat)
        }
        if (ind[1] != 0) 
            data = dat.rem.ind(data, ind)
        if (pop[1] != 0) 
            data = dat.rem.pop(data, pop)
        return(data)
    }
    dat.pw <- function(data, pw = 0) {
        if (pw[1] == 0) 
            return(data)
        else {
            data = data[data[, 1] == pw[1] | data[, 1] == pw[2], 
                ]
            return(data)
        }
    }
    nb.pop <- function(data) {
        data = data[order(data[, 1]), ]
        nb.ind = dim(data)[1]
        nb.pop = 1
        for (i in 1:(nb.ind - 1)) if (data[i, 1] != data[i + 
            1, 1]) 
            nb.pop = nb.pop + 1
        return(nb.pop)
    }
    pop.freq <- function(data) {
        data = data[order(data[, 1]), ]
        nb.ind = dim(data)[1]
        dat.fra = as.data.frame(data)
        nb.pop = 1
        for (i in 1:(nb.ind - 1)) if (data[i, 1] != data[i + 
            1, 1]) 
            nb.pop = nb.pop + 1
        pop.freq.vec = rep(1, nb.pop)
        name = rep(0, nb.pop)
        k = 1
        name[1] = as.character(dat.fra[1, 1])
        for (i in 2:nb.ind) if (dat.fra[i - 1, 1] == dat.fra[i, 
            1]) 
            pop.freq.vec[k] = pop.freq.vec[k] + 1
        else {
            k = k + 1
            name[k] = as.character(dat.fra[i, 1])
        }
        names(pop.freq.vec) = name
        return(pop.freq.vec)
    }
    Pst.val <- function(data, csh = 1) {
        nbpop = nb.pop(data)
        nb.var = dim(data)[2] - 1
        data = data[order(data[, 1]), ]
        if (nbpop == 1) 
            return(rep(0, nb.var))
        else {
            pop.freq = pop.freq(data)
            Pst.clm <- function(dat, clm) {
                mea = mean(dat[, clm], na.rm = TRUE)
                nna.clm = nonNa.clm(dat, clm)
                SSTotal = (nna.clm - 1) * var(dat[, clm], na.rm = TRUE)
                mea.pop = rep(0, nbpop)
                nna.pop.freq = rep(0, nbpop)
                nna.pop.freq[1] = nonNa.clm(dat[1:(pop.freq[1]), 
                  ], clm)
                nb.allna.pop = 0
                if (nna.pop.freq[1] == 0) 
                  nb.allna.pop = 1
                else mea.pop[1] = mean(dat[1:(pop.freq[1]), clm], 
                  na.rm = TRUE)
                for (i in 2:nbpop) {
                  nna.pop.freq[i] = nonNa.clm(dat[(sum(pop.freq[1:(i - 
                    1)]) + 1):(sum(pop.freq[1:i])), ], clm)
                  if (nna.pop.freq[i] != 0) 
                    mea.pop[i] = mean(dat[(sum(pop.freq[1:(i - 
                      1)]) + 1):(sum(pop.freq[1:i])), clm], na.rm = TRUE)
                  else nb.allna.pop = nb.allna.pop + 1
                }
                SSBetween = sum(nna.pop.freq * (mea.pop - mea)^2)
                SSWithin = SSTotal - SSBetween
                if ((nna.clm - nbpop + nb.allna.pop) * (nbpop - 
                  nb.allna.pop - 1) != 0) {
                  MSWithin = SSWithin/(nna.clm - nbpop + nb.allna.pop)
                  MSBetween = SSBetween/(nbpop - nb.allna.pop - 
                    1)
                  return(csh * MSBetween/(csh * MSBetween + 2 * 
                    MSWithin))
                }
                else {
                  if ((nna.clm - nbpop + nb.allna.pop) == 0) 
                    return(1)
                  else return(0)
                }
            }
            pst.val = rep(0, nb.var)
            for (j in 1:nb.var) pst.val[j] = Pst.clm(data, j + 
                1)
            return(pst.val)
        }
    }
    boot.pst.va <- function(data, csh, boot, clm) {
        nb.ind = dim(data)[1]
        dat = data[, c(1, clm)]
        boot.val = rep(0, boot)
        for (i in 1:boot) {
            da = dat[sample(1:nb.ind, nb.ind, T), ]
            boot.val[i] = Pst.val(da, csh)
        }
        return(boot.val)
    }
    ConInt.pst.va <- function(data, csh, boot, clm, per) {
        boot.pst.val = boot.pst.va(data = data, csh = csh, boot = boot, 
            clm = clm)
        boot.pst.val = sort(boot.pst.val)
        return(c(boot.pst.val[floor(boot * (1 - per)/2 + 1)], 
            boot.pst.val[ceiling(boot * (per + 1)/2)]))
    }
    Trace <- function(data, pts, boot, Fst, xm, ci) {
        tra.pst.fst.va <- function(data, pts, clm, Fst, xmax, 
            lab.pos) {
            data = data[, c(1, clm)]
            points <- function(nb.pts) {
                pst.va = Pst.val(data, 0)
                for (i in 1:nb.pts) pst.va = c(pst.va, Pst.val(data, 
                  xmax * i/nb.pts))
                return(pst.va)
            }
            pst.val = points(nb.pts = pts)
            csh.val = xm * c(0:pts)/pts
            plot(pst.val ~ csh.val, type = "l", xlab = "c/h^2", 
                ylab = "Pst", main = c("Pst variations:", names(data)[2]), 
                ylim = c(0, 1), col = "firebrick1")
            if (Fst != -1) {
                abline(h = Fst, col = "green", lty = 4)
                text(0.05 * lab.pos - 0.06, Fst + 0.04 * lab.pos - 
                  0.01, "Fst", col = "green")
            }
        }
        tra.confint.va <- function(clm) {
            point <- function(nb.pt) {
                ci.pst.va = ConInt.pst.va(data, csh = 0, boot = boot, 
                  clm = clm, per = pe)
                upbnd.val = ci.pst.va[2]
                lowbnd.val = ci.pst.va[1]
                for (i in 1:nb.pt) ci.pst.va = c(ci.pst.va, ConInt.pst.va(data, 
                  csh = xm * i/nb.pt, boot = boot, clm = clm, 
                  per = pe))
                for (i in 1:nb.pt) upbnd.val = c(ci.pst.va[2 + 
                  2 * i], upbnd.val)
                for (i in 1:nb.pt) lowbnd.val = c(lowbnd.val, 
                  ci.pst.va[1 + 2 * i])
                return(c(upbnd.val, lowbnd.val))
            }
            ci.pst.val = point(nb.pt = pts)
            csh.val = xm * c(0:pts)/pts
            rev.csh.val = rev(csh.val)
            plot(ci.pst.val ~ c(rev.csh.val, csh.val), type = "l", 
                xlab = "c/h^2", ylab = "Pst", main = c("Pst variations:", 
                  names(data)[clm]), ylim = c(0, 1), col = "chocolate4", 
                lty = 2)
        }
        nb.var = dim(data)[2] - 1
        nb.gra.lon = ceiling(sqrt(nb.var))
        par(mfrow = c(nb.gra.lon, nb.gra.lon))
        for (i in 1:nb.var) {
            tra.pst.fst.va(data, pts = pts, Fst = Fst, clm = i + 
                1, xmax = xm, lab.pos = nb.gra.lon)
            if (ci == 1) {
                par(new = TRUE)
                tra.confint.va(clm = i + 1)
            }
        }
    }
    if (va[1] == 0) {
        nb.var = dim(data)[2] - 1
        va = 1:nb.var
    }
    else {
        nb.var = length(va)
        for (i in 1:nb.var) {
            for (j in 2:dim(data)[2]) {
                if (names(data)[j] == va[i]) 
                  va[i] = j - 1
            }
        }
        va = as.numeric(va)
        if (is.na(sum(va)) == TRUE) 
            return("va is not valid!")
    }
    data = dat.fra.prep(data)
    data = dat.rem.ind.pop(data, ind = Ri, pop = Rp)
    data = dat.pw(data, Pw)
    print("Populations sizes are:")
    print(pop.freq(data))
    data = data[, c(1, va + 1)]
    dev.new()
    Trace(data, pts = pts, boot = boot, Fst = Fst, xm = xm, ci = ci)
  }

[Package Pstat version 1.2 Index]