The test data that comes with mlgenomeu is included as part of this package, we
can load it via data()
library(dfe) data('ma_control') data('ma_experimental')
If these data really come from a mutation accumulation experiment it's an odd one, as the among-line variance is lower in the experimental data than the control:
c(Control=var(ma_control), Experiment=var(ma_experimental))
d_exp <- density(ma_experimental) plot(density(ma_control), xlab="Some fitness measure", ylim=c(0,0.8), main="Emperical density kernels") lines(d_exp$x, d_exp$y, col='red')
The functions in dfe
take the difference in fitness between ancestral and
mutation-accumulation lines and input. We can do that, and see the mean fitness
of the MA lines is sligthly lower than that of the controls (as it makes
calculation much mroe straight forward, we treat positive values as deleterious
mutations):
w <- ma_experimental - mean(ma_control) mean(w)
The classic Bateman-Mukai estimators are privded by the function BM()
which
takes the dchange in mean fitness ($\Delta M$) and the variance of the ancestral
lines subjected to the same fitness assay ($V_e$) as inputs.
Unfortunately, the low variane in the MA lines is going to make for a
non-sensical estimate of the mutation rate ($U$) and effect size ($E[a]$):
BM(w, var(ma_control))
Negative mutations of very large effect... We'll see a sensible use of these estimators later.
If you are familar with mlgenomeu there are couple of things to note about our implimentation of this model. First, we use a paramaterisation of the Gamma distribution with shape ($\alpha$) and rate ($\beta$) parmaters, meaning the mean fitness-effect of mutations ($E[a]$ in mlgenomu) is given by $\frac{\alpha}{\beta}$, and variance in fitness is $\frac{\alpha}{\beta^2}$. Second, we take the mean difference in fitness for the MA lines as input.
In the documentation for mlgenomeu, you are told to set some specific starting
values and fixed-values to start performing a profile likelihood search. We
have dveloped a method of moments estimator to find good starting values,
depending only on a suggestd mutation rate. For now we need to first explicitly
calculate these starting values using mom_ma_gamma()
mom_ma_gamma(obs=w, Ut=0.5, Ve=var(ma_control))
Remember, the mean effect of mutations is $\frac{\alpha}{\beta}$, which, given our estimate isi $\frac{0.003}{0.046} \approx 0.0652$
We can use these values to calculate the log liklihood of a specific
distribution of fitness effects and mutation rate. The likelihood functions in
dfe
all start with dma_
which is short for density of a mutation
accumulation model. To calculate the model we have to provide values for the
shape and rate that specify the Gamma distribution, $Ve$ and $Ut$ the
per-experiment mutation rate:
dma_gamma(obs=w, shape=0.003, rate=0.046, Ve=var(ma_control), Ut=0.5)
The function fit_ma_gamma
can perform a maximum liklihood search. This
function wraps the R function stats4::mle
, which in turn using the base R
function optim
to perform the search. We just need to provide a lists of fixed
and starting values. Let's follow the mlgenomeu example and fix two paramaters
(Ve
and the shape
) and estimate the rest:
fit <- fit_ma_gamma(obs=w, fixed=list(shape=0.003, Ve=var(ma_control)), start=list(Ut=0.5, rate=0.046)) fit
The function throws some warnings, which are the result of the very low number we end up estimating. Ignoring the warnings for now, we can see our ML estimates of the shape paramater and mutation rate are getitng close to zero. In fact, we can show that a model with zero mutations gives approximately as good a fit:
fit dma_gamma(obs=w, shape=0.003, rate=0.046, Ve=var(ma_control), Ut=0)
This also gives us a change to confirm the density function is working as expected, since in the case of zero mutations the likelihood should be the same as that caclulated from a normal distribution with $mean = 0$ and $\sigma = Ve$:
log(prod(dnorm(w, mean=0, sd=sqrt(var(ma_control)))))
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.