| 2. Phenomenological models {sra} | R Documentation |
Descriptive models of artificial-selection responses: auto-regressive models
Description
The sraAutoreg functions are wrappers for the maximum-likelihood optimization function mle. They propose descriptive models for the dynamics of genetic architectures of different complexities based on an auto-regressive framework, additional parameters corresponding to different generation lags. The model can also be fit considering logatithmic, "heritability" and "evolvability" scales.
Usage
sraAutoreg(sradata, active = c(FALSE, TRUE, FALSE, FALSE),
start = NULL, fixed = NULL, negative.k = FALSE,
rand = 0, rep = 1, ...)
sraAutoregLog(sradata, active = c(FALSE, TRUE, FALSE, FALSE),
start = NULL, fixed = NULL, negative.k = FALSE,
rand = 0, rep = 1, ...)
sraAutoregHerit(sradata, active = c(FALSE, TRUE, FALSE, FALSE),
start = NULL, fixed = NULL, negative.k = FALSE,
rand = 0, rep = 1, ...)
sraAutoregEvolv(sradata, active = c(FALSE, TRUE, FALSE, FALSE),
start = NULL, fixed = NULL, negative.k = FALSE,
rand = 0, rep = 1, ...)
Arguments
sradata |
A data object generated by |
active |
A vector of four booleans, corresponding to the active lags for both genetic and environmental variances (see Details). By default, only lag 1 is active, corresponding to an exponential change of the variances. |
start |
A named list of starting values for the convergence algorithm. |
fixed |
A named list of the parameters that have to be kept constant. |
negative.k |
Whether or not the k parameters can take negative values. Negtive values for k lead to more complex likelihood functions, and the resulting dynamics may display cyclic patterns. |
rand |
Amount of randomness for the starting values. Useful in case of convergence issues. Although this variable can take any positive value, reasonable figures should not exceed 0.2. |
rep |
Number of convergence attempts. When the likelihood function is complex, which is often the case when the number of parameters exceeds 6 to 8, convergence may often fail or end up on a local maximum. When |
... |
Additional parameters to be passed to |
Details
Model
The following summarizes the models developed in Le Rouzic et al. 2010.
The mean of the population \mu changes according to the Lande equation (Lande and Arnold 1983):
\mu(t+1) = \mu(t) + VarA(t)*\beta(t),
where \beta(t) is the selection gradient at generation t.
The genetic architecture models predict the dynamics of a parameter P as:
P(t+1) = k0 + k1*P(t) + k2*P(t-1) + k3*P(t-2)
Models with time lags > 3 could be easily implemented, but convergence issues increase with the number of parameters. The first time points are calculated as if P(t<1)=P(1), e.g. P(3) = k0 + k1*P(2) + k2*P(1) + k3*P(1).
Each model considers the dynamics of two independent parameters, one related to the additive genetic variance (VarA), one related to the residual (environmental) variance VarE (which actually also accounts for all non-additive genetic variance).
The default model sraAutoreg considers directly the dynamics of VarA (parameters: kA0, kA1, kA2, and kA3) and the dynamics of VarE (parameters kE0, kE1, kE2, and kE3).
The log scale turns a multiplicative trait into an additive one, and is particularly relevant for ratio-scale traits (i.e. most quantitative traits such as size, fertility, etc). The original data is transformed assuming log-normality, and the likelihood is computed based on the log-normal density function.
The "heritability" model sraAutoregHerit focuses on the dynamics of h^2=VarA/(VarA+VarE), described by the parameters kA0, kA1, kA2, and kA3, and considers that kE0, kE1, kE2, and kE3 describe the dynamics of the phenotypic variance VarP=VarA+VarE. Therefore, VarE is constrained both by the dynamics of VarP and the independent dynamics of h^2.
The "evolvability" model considers that kA0, kA1, kA2, and kA3 describe the dynamics of IA(t) = VarA(t)/(\mu(t)^2), and kE0, kE1, kE2, and kE3 the dynamics of IE(t) = VarE(t)/(\mu(t)^2).
Shortcut for active and inactive parameters
The user will often have to fit models of different complexity. This can be achieve by manipulating the active vector. c(FALSE,FALSE,FALSE,FALSE) corresponds to a constant-variance model (no dynamic parameter), c(TRUE,FALSE,FALSE,FALSE) to a case in which only kA0 and kE0 are active, c(TRUE,TRUE,FALSE,FALSE) to active parameters for lags 0 and 1 only, etc. The total number of parameters in the model will be 3+2*x, where x is the number of TRUE in the vector active.
To bypass the constrains of this shortcut, it is possible to specify the active and inactive parameters through the list of starting values. A combination such as active=c(TRUE,FALSE,TRUE,FALSE), start=list(kA1=0,kE3=NA), fixed=list(kE2=1) will lead to a model with 8 active parameters (mu0, varA0, varE0, kA0, kE0, kA1 (which starting value will be 0), kA2, and kE3 (which starting value, specified as NA, will be determined via the function sraStartingvalues. All other parameters are fixed.
Parameterization
The models thus involve up to 11 parameters: three initial values (\mu(1), VarA(1) and VarE(1)), four parameters to describe the dynamics of the additive variance (or relative variable such as IA or h^2) (kA0, kA1, kA2, and kA3), and four parameters for the environmental variance (or IE, or h^2): kE0, kE1, kE2, and kE3. To make numerical convergence more efficient, the following parameterization was implemented: parameters mu0, logvarA0 and logvarE0 correspond to the estimates of the initial values of, respectively, the population mean, the logarithm of the additive variance, and the logarithm of the environmental variance. The parameters kA0 and kE0 are calculated as relative to the initial values of the dynamic variable, e.g. relativekA0 = k0A/VarA(1) (so that relativekA0 has the same unit and the same order of magnitude as kA1, kA2 and kA3).
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
Le Rouzic, A., Houle, D., and Hansen, T.F. (2011) A modelling framework for the analysis of artificial selection-response time series. Genetics Research.
Lande, R., and Arnold, S. (1983) The measurement of selection on correlated characters. Evolution 37:1210-1226.
See Also
sraCstvar, sraDrift and all other mechanistic 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
# Making the example reproducible
########### 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 ####################
# Autoregressive models
autor <- sraAutoreg(data)
# Details of the model:
AIC(autor)
coef(autor)
plot(autor)
plot(autor, var=TRUE)
# Alternative scales
autor.log <- sraAutoregLog(data)
autor.herit <- sraAutoregHerit(data)
autor.evolv <- sraAutoregEvolv(data)
# Changes in the complexity of the model:
autor0 <- sraAutoreg(data, active=c(TRUE,TRUE,FALSE,FALSE))
# In case of convergence issues
autor1 <- sraAutoreg(data, active=c(TRUE,TRUE,TRUE,TRUE), rep=2, rand=0.1)