| 3. Mechanistic models {sra} | R Documentation |
Descriptive models of artificial-selection responses: quantitative genetics models
Description
The sra functions for mechanistic model are wrappers for the maximum-likelihood optimization routine mle. They implement classical quantitative genetics models, fit them to artificial-selection time series, and provide estimates of e.g. the effective population size, the mutational variance, the strength of genetic / environmental canalization, the directionality and strength of epistasis, etc., given some assumptions about the properties of the genetic architecture.
Usage
sraCstvar(sradata, start=NULL, fixed=NULL, macroE=FALSE, Bulmer=TRUE, ...)
sraDrift(sradata, start=NULL, fixed=NULL, macroE=FALSE, Bulmer=TRUE, ...)
sraMutation(sradata, start=NULL, fixed=NULL, macroE=FALSE, Bulmer=TRUE, ...)
sraCanalization(sradata, start=NULL, fixed=NULL, macroE=FALSE, Bulmer=TRUE, ...)
sraCanalizationOpt(sradata, start=NULL, fixed=NULL, macroE=FALSE, Bulmer=TRUE, ...)
sraSelection(sradata, start=NULL, fixed=NULL, macroE=FALSE, Bulmer=TRUE, ...)
sraDirepistasis(sradata, start=NULL, fixed=NULL, macroE=FALSE, ...)
Arguments
sradata |
A data object generated by |
start |
A named list of starting values for the convergence algorithm. |
fixed |
A named list of the parameters that have to be kept constant. |
macroE |
Whether or not macro-environmental effects (random deviation of the phenotypic mean each generation) should be included. This might have some side effects. |
Bulmer |
Whether or not the impact of linkage disequilibrium (Bulmer effect) due to selection on variance should be accounted for. |
... |
Additional parameters to be passed to |
Details
All functions (except sraDirepistasis) rely on the same underlying model, and thus simply provide convenient shortcuts for different sets of parameters to fit.
mu0 is the initial phenotype of the population.
logvarA0 is the logarithm of the initial additive variance in the population.
logvarE0 is the logarithm of the initial environmental variance in the population.
logNe is the logarithm of the effective population size.
logn is the logarithm of the effective number of loci.
logvarM is the logarithm of the mutational variance.
kc and kg are the strength of environmental and genetic canalization, respectively.
o corresponds to the 'optimum' phenotype. When fixed and set to NA, o is identical to mu0. For convenience, the same optimum is used for environmental canalization, genetic canalization and natural stabilizing selection.
s corresponds to the strength of natural selection.
logvarepsilon is the logarithm of the variance of the epistatic coefficient (\varepsilon). This parameter is fixed by default, since it is unlikely that realistic data sets contain enough information to estimate it properly.
The dynamic model that is fitted (except of the directional epistasis, detailed in the next paragraph) is:
\mu(t+1) = \mu(t) + VarA(t) * (\beta(t) + s \delta(t))
VarA(t+1) = VarM + VarA(t) * (1 - 1/(2*N_e)) * e^{kg*(|\delta(t+1)| - |\delta(t)|)}
VarE(t+1) = VarE(t) * e^{kc*|\delta(t)|}
\mu(1), VarA(1) and varE(1) are parameters of the model, \beta(t) is the selection gradient, calculated for each generation from the data set, and \delta(t) = \mu(t) - o.
The directional epistasis model has its own setting:
\mu(t+1) = \mu(t) + varA(t) * \beta(t)
VarA(t+1) = VarA(t) + 2*\beta(t)*\varepsilon*VarA(t)^2
VarE(t+1) = VarE(1) + (\varepsilon^2 + Var\varepsilon)*VarA(t)^2
Where epsilon is a key parameter, representing the directionality of epistasis (Hansen and Wagner 2001, Carter et al. 2005). The properties of the likelihood function when epsilon varies makes numerical convergence tricky, and the function SRAdirepistasis actually performs two model fits: one for positive epsilons (the estimated parameter being logepsilon) and one for negative epsilons (estimating logminusepsilon). The likelihood of both models is then compared, and the best model is returned, providing either logepsilon or logminusepsilon. Although part of the model, the parameter logvarepsilon appeared to affect environmental variance only weakly, and its estimation is problematic in most cases.
These models are extensively described in le Rouzic et al 2010.
Value
The functions return objects of class srafit, a list containing information about the model, the data, and the parameter estimates. Some standard R functions can be applied to the object, including AIC (AIC.srafit), logLik (logLik.srafit), vcov (vcov.srafit), coef (coef.srafit) confint (confint.srafit), and plot (plot.srafit).
Author(s)
Arnaud Le Rouzic
References
Carter, A.J.R., Hermisson, J. and Hansen, T.F. (2005) The role of epistatic genetic interactions in the response to selection and the evolution of evolvability. Theor. Pop. Biol. 68, 179-196.
Hansen, T.F. and Wagner, G.P. (2001) Modelling genetic architecture: a multilinear theory of gene interaction. Theor. Pop. biol. 59, 61-86.
Le Rouzic, A., Houle, D., and Hansen, T.F. (2010) A modelling framework for the analysis of artificial selection-response time series. in prep.
See Also
sraAutoreg, sraAutoregHerit and sraAutoregEvolv for phenomenological models, sraAutoregTsMinuslogL and sraAutoregTimeseries for some details about the internal functions, AIC.srafit, logLik.srafit,vcov.srafit, coef.srafit, confint.srafit, plot.srafit for the analysis of the results.
Examples
########### Generating a dummy dataset ################
m <- c(12,11,12,14,18,17,19,22,20,19)
v <- c(53,47,97,155,150,102,65,144,179,126)
s <- c(15,14,14,17,21,20,22,25,24,NA)
n <- c(100,80,120,60,100,90,110,80,60,100)
########## Making a sra data set #######################
data <- sraData(phen.mean=m, phen.var=v, phen.sel=s, N=n)
#################### Data Analysis ####################
cstvar <- sraCstvar(data)
drift <- sraDrift(data)
direpi <- sraDirepistasis(data)
# In case of convergence problems, better starting values can be provided:
direpi <- sraDirepistasis(data, start=list(mu0=10, logvarA0=log(20), logvarE0=NA),
fixed=list(logNe=log(50)))
plot(cstvar)
AIC(direpi)