%\VignetteIndexEntry{Using OnAge}
%\VignettePackage{OnAge}
\documentclass[11pt]{article}
\usepackage{times}
\usepackage{hyperref}
\usepackage{geometry}
\usepackage{natbib}
\usepackage[pdftex]{graphicx}
\usepackage{url}
\usepackage[utf8]{inputenc}
\SweaveOpts{keep.source=TRUE,prefix=TRUE}
\newcommand{\HH}[1]{{\bf H_#1}}
% R part
\newcommand{\R}[1]{{\textsf{#1}}}
\newcommand{\Rfunction}[1]{{\texttt{#1}}}
\newcommand{\Robject}[1]{{\texttt{#1}}}
\newcommand{\Rpackage}[1]{{\textit{#1}}}
\newcommand{\Rclass}[1]{{\textit{#1}}}
\newcommand{\Metas}[1]{{\texttt{#1}}}
\begin{document}
\title{\Rpackage{OnAge}: Test of between-group differences in the onset of senescence}
\author{Laurent Jacob \and Fréderic Douhard \and Jean-François
Lemaître \and Jean-Michel Gaillard \and Aurélie Siberchicot}
\maketitle
\begin{abstract}
\Rpackage{OnAge} implements the procedure used
in~\cite{Douhard2017cost} to test whether the onset of senescence in
an individual trait is the same between two groups (\emph{e.g.},
females vs males or two populations). The procedure is a likelihood
ratio test comparing a model with single onset for both groups to a
model allowing different onsets.
\end{abstract}
\section{Introduction}
In the living world, the process of senescence is the rule rather than
the exception~\citep{Nussey2013Senescence}. A decrease in the
probability to survive (i.e. actuarial senescence), reproduce
(i.e. reproductive senescence), in body condition or other phenotypic
measurements (e.g. level of immune performance) with increasing age
has now been documented in a tremendous number of species
(e.g.~\citep{Nussey2013Senescence, Gaillard2017Senescence,
Cheynel2017Immunosenescence}). However, senescence patterns can be
highly variable within and among
species~\citep{Jones2008Senescence}. For many years, most studies have
tried to understand the evolutionary roots of this variability by
comparing the rate of a senescence of a given trait between two (or
more) species or populations. These rates of senescence were measured
by fitting, on age-specific data, a mathematical model (e.g. a
Gompertz or a Weibull in the case of survival, a quadratic or a
threshold mode in the case of body mass) starting from the age at
first reproduction. The rational of fitting such models from the age
at first reproduction takes its origin in two pioneer contributions of
the evolutionary biology of
senescence~\citep{Williams1957Pleiotropy,Hamilton1966moulding}
following Williams' initial assumption that the force of selection
starts to decline from the age of reproductive
maturity~\citep{Williams1957Pleiotropy}. However, recent and detailed
investigations of age-specific changes in several life history traits
have revealed that the decline in performance rarely starts from the
age at first reproduction~
\citep{Jones2008Senescence,Hayward2015Asynchrony, Peron2010Senescence}. In
addition, there is increasing evidence that ecological
and biological differences among populations can impact the onset
rather than the rate of senescence~\citep{Tidiere2015Does}. The study
of the onset of senescence is still at its infancy and there is
currently a great need to develop methods that allow detecting and
statistically assessing the difference in age at the onset of senescence
in a given trait across different populations. The package \Rpackage{OnAge}
fulfills these objectives.
\section{Software features}
\Rpackage{OnAge} exports a single function
\Rfunction{onset.test}. This functions takes as input a log likelihood
function \Rfunction{ll}, two data frames \Robject{data1} and
\Robject{data2} containing the data of each group, and a vector
describing over which range of onset of senescence the log likelihood
should be maximized. The first argument of \Rfunction{ll} should be
the age at the onset of senescence, its second argument should be a data
frame. The function should return the log-likelihood of a model which
depends on this onset of senescence. In the example of this vignette, we
use a mixed effect linear model. The (fixed) effect of age is linear
after the onset and absent before. The two data frames should contain
all relevant fields for the model in the two groups between which
the differential onset is tested. They will be passed as second
arguments to \Rfunction{ll}.
\Rfunction{onset.test} returns an asymptotic p-value for the null
hypothesis that the age at the onset of senescence is the same in both
groups. This p-value relies on a log likelihood ratio. More precisely,
if $L(\theta_1, \theta_2)$ is the likelihood parameterized by the ages
at the onset of senescence $(\theta_1, \theta_2)$ in two groups ---
and possibly other parameters ---, we test the null hypothesis
$\HH{0}: \theta_1 = \theta_2$. We form the log ratio statistic
comparing the maximum likelihood under $\HH{0}$ and a complementary
$\HH{1}$ where $\theta_1$ and $\theta_2$ are allowed to take on
different optimal values. Wilk’s theorem~\citep[Thm 6.5
p.432]{Shao2003Mathematical} provides that under $\HH{0}$, the log
likelihood ratio statistic converges towards a $\chi_2^r$ random
variable, where $r$ is the difference between the dimensionality of
the parameter under $\HH{1}$ and $\HH{0}$. In our case, $r=1$ since we
only relax the assumption that the two ages at the onset of senescence
are equal.
\Rfunction{onset.test} additionally returns maximum likelihood
estimates for the age at the onset of senescence in both groups and in
the population obtained by merging them, and optional confidence
intervals for these three ages at the onset of senescence. It also
provides the maximized log-likelihood under $\HH{0}$ and under
$\HH{1}$, the log-likelihood ratio statistic and a binary flag
indicating whether the optimization failed to converge at any point.
\section{Case study}
We use the dataset of~\cite{Douhard2017cost}, containing the age, body
mass and 14 other attributes of 454 roe deer from Chizé and Trois
Fontaines. We illustrate how the package can be used to test whether
the onset of body mass senescence varies between males and females.
We first illustrate how to use our test on the body mass measured
in~\cite{Douhard2017cost}. In a second step, we apply our test to body
mass simulated from a mixed model. This allows us to evaluate:
\begin{enumerate}
\item whether our asymptotic p-values are correctly calibrated,
\emph{i.e.}, whether for all $\alpha\in [0, 1]$, a proportion
$\alpha$ of the experiments simulated under $\HH{0}$ yield a p-value
smaller than $\alpha$ --- meaning that the p-value actually
translates into a false positive rate.
\item as a sanity check, whether our test has some power to detect
a differential onset.
\item whether the optimization process runs into numerical problems,
how these problems can be diagnosed and how they affect the result.
\end{enumerate}
\subsection{Loading the library and the data}
We install and load the \Rpackage{OnAge} package by typing or pasting the
following codes in R command line.
<<lib, echo=TRUE, eval = FALSE>>=
install.packages("OnAge")
@
<<lib, echo=TRUE>>=
library(OnAge)
@
We then load the dataset and split it into four differents groups:
females Chizé (FCH), males Chizé (MCH), females Trois Fontaines (FTF),
males Trois Fontaines (MTF).
<<data, echo=TRUE>>=
data(RoeDeerMassData)
RoeDeerMassData$ID <- factor(RoeDeerMassData$ID)
RoeDeerMassData$cohort <- factor(RoeDeerMassData$cohort)
dataFCH <- RoeDeerMassData[RoeDeerMassData$sex%in%"F"&
RoeDeerMassData$population%in%"CH", ]
dataMCH <- RoeDeerMassData[RoeDeerMassData$sex%in%"M"&
RoeDeerMassData$population%in%"CH", ]
dataFTF <- RoeDeerMassData[RoeDeerMassData$sex%in%"F"&
RoeDeerMassData$population%in%"TF", ]
dataMTF <- RoeDeerMassData[RoeDeerMassData$sex%in%"M"&
RoeDeerMassData$population%in%"TF", ]
@
\subsection{Defining the model}
We define the model in which we test for differential onset. Here, we
use a mixed effect linear model representing the body mass as a linear
combination of fixed effects for factors \Robject{b1(age, thr)}
representing the age of the individual, age at last capture, last year
of last capture and random effects for factors \Robject{ID}
(individual) and cohort. Function \Rfunction{b1} transforms its input
\Robject{age} such that a linear function with slope $\alpha$ of the
transformed \Robject{age} is a piecewise linear function of the
original \Robject{age}, with effect zero before the threshold
\Robject{bp} and slope $\alpha$ after. The likelihood ratio test will
compare the log-likelihood of a joint model where all individuals
follow the same distribution to one allowing different onset
\Robject{thr} for males and females.
<<model, echo=TRUE>>=
## b1: function for piecewise regression (transforms x into 0 before bp)
b1 <- function(x, bp) ifelse(x < bp, 0, x - bp)
## Use this function to define the model in which the differential
## onset hypothesis is tested.
ll.real <- function(thr, dataIn){
logLik(lme4::lmer(body.mass ~ b1(age, thr) + age.at.last.capture +
last.year.of.capture + (1|ID) + (1|cohort),
data=dataIn, REML="FALSE"))
}
## Same model using simulated body mass
ll.sim <- function(thr, dataIn){
logLik(lme4::lmer(body.mass.sim ~ b1(age, thr) + age.at.last.capture +
last.year.of.capture + (1|ID) + (1|cohort),
data=dataIn, REML="FALSE"))
}
@
\subsection{Testing both populations for sex-related differential age
at the onset of senescence}
Once the model has been defined, testing $\HH{0}$ within each
population using \Rfunction{onset.test} is rather straightforward:
<<test, echo=TRUE>>=
search.range <- c(6, 12)
search.range.TF <- search.range.CH <- search.range
@
<<testTF, echo=TRUE, include=FALSE, fig=TRUE, width=12, height=4>>=
res.tf <- onset.test(ll.real, dataFTF, dataMTF, search.range.TF,
do.plot=TRUE)
@
<<testCH, echo=TRUE, include=FALSE, fig=TRUE, width=12, height=4>>=
res.ch <- onset.test(ll.real, dataFCH, dataMCH, search.range.CH,
do.plot=TRUE)
@
<<test-result, echo=TRUE>>=
cat(sprintf("p-value for differential age at onset is %g in
Trois Fontaines, %g in Chizé", res.tf$pv, res.ch$pv))
@
\Rfunction{onset.test} also returns a maximum likelihood estimate of
the onset within each population and optionally computes a confidence
interval for each onset. In addition, it can provide a plot of the
log-likelihood against the onset, showing the maximum likelihood
estimate $\hat{\theta}_{ML}$ as a red star and the confidence interval
as vertical dotted lines, as illustrated on
Figure~\ref{fig:llprof-real}.
\begin{figure}[ht]
\begin{center}
\includegraphics[width=0.9\textwidth]{RoeDeerMass-simulation-testTF}
\includegraphics[width=0.9\textwidth]{RoeDeerMass-simulation-testCH}
\caption{\label{fig:llprof-real}Log-likelihood profiles on the data
of~\cite{Douhard2017cost}: Trois Fontaines (top) and Chizé
(bottom) cohorts. Left: female individuals, middle: male
individuals, right: overall population. The red star denotes the
maximum likelihood estimate of the onset. The two vertical dotted
lines denote the lower and upper bounds of the computed $95\%$
confidence interval. The horizontal dotted line represents the
log-likelihood above which a value of the onset belongs to the
confidence interval.}
\end{center}
\end{figure}
The confidence interval is computed by inverting the acceptance region
of a likelihood ratio test of the null hypothesis $\theta = \theta_0$:
it is formed by all onsets $\theta_0$ for which this null hypothesis
is not rejected at level $1-\alpha$, where $\alpha$ is the desired
confidence level. The method looks for roots of the function
$2(L(\hat{\theta}_{ML}) - L(\theta_0)) - q_{0.95}$ where $q_{0.95}$ is
the 95-th percentile of the $\chi_2^1$ distribution. A caveat is that
this method assumes that the log-likelihood is monotonically
decreasing, \emph{i.e.}, that the log-likelihood ratio statistic is
monotonically increasing around the maximum likelihood estimator of
onset. If this is not the case, the true confidence interval may be
formed by a union of intervals and the method will return a set formed
by some of these intervals and all intervals in between. This happens
on the Trois Fontaines female subset (top left), where the
log-likelihood curve goes below the horizontal line between the
vertical ones. The actual confidence interval is the set of onsets
whose log-likelihood is above the horizontal line, \emph{i.e},
$[7, 9.8] \cup [10.4, 10.7]$, whereas the method returns $[7, 10.7]$.
Fortunately, it is easy to check this visually by setting
\texttt{do.plot=TRUE} in the \Rfunction{onset.test} function
(Figure~\ref{fig:llprof-real}).
\subsection{Running the simulation}
We now simulate body mass based on the same model used in
\Rfunction{ll.sim}.
<<simulation, echo=TRUE>>=
sfunc <- function(dataIn, b1, thr, bw.offset, age.effect,
age.last.effect, last.effect, id.s2, coh.s2, s2){
id.dum <- sapply(levels(dataIn$ID),
FUN=function(ll) 1*(levels(dataIn$ID)[dataIn$ID] == ll))
id.effect <- matrix(rnorm(ncol(id.dum), sd=sqrt(id.s2)), ncol=1)
coh.dum <- sapply(levels(dataIn$cohort),
FUN=function(ll) 1*(levels(dataIn$cohort)[dataIn$cohort] == ll))
coh.effect <- matrix(rnorm(ncol(coh.dum), sd=sqrt(coh.s2)), ncol=1)
dataIn$body.mass.sim <- bw.offset +
(b1(dataIn$age, thr) * age.effect +
dataIn$age.at.last.capture * age.last.effect +
dataIn$last.year.of.capture * last.effect +
id.dum %*% id.effect +
coh.dum %*% coh.effect +
rnorm(nrow(dataIn), sd=sqrt(s2)))
return(dataIn)
}
## Noise level in the linear model
s2 <- 10
## True onset for males and females under H0
thr.m <- thr.f <- 7
## Difference of onsets under H1
thr.delta <- 2
## Baseline body mass
bw.offset <- 10
## Fixed effects on bw
age.effect <- -1
age.last.effect <- 1
last.effect <- 1
## Variance of random effects on bw
id.s2 <- 1
coh.s2 <- 1
@
We run\footnote{The simulations are not actually run to produce the
vignette, we just load pre-computed results to save time.} $5000$
simulations under $\HH{0}$ and $5000$ others under $\HH{1}$. For each
simulation, we generate one set of body mass from the Chizé data and
another from the Trois Fontaines data and apply the
\Rfunction{onset.test} function within each generated dataset to
compare the age at the onset of senescence between males and females.
We also re-run likelihood maximization to store two indicators of
numerical optimization quality: the number of evaluations and the norm
of the gradient at the optimum.
<<load-simu, echo=FALSE, results=hide>>=
# data('pre-computed-simulation')
load(file = "pre-computed-simulation.RData")
@
<<run, echo=TRUE, eval=FALSE, results=hide>>=
## Number of simulations we want to run under H0 and H1.
n.h0 <- 5000
n.h1 <- 5000
n.rep <- n.h0 + n.h1
pv.tf <- pv.ch <- llr.tf <- llr.ch <-
lh1.tf <- lh1.ch <- lh0.tf <- lh0.ch <- rep(NA, n.rep)
cvg.tf <- cvg.ch <- rep(NA, n.rep)
warn.tf <- warn.ch <- rep(FALSE, n.rep)
data.f.ch <- data.m.ch <- data.f.tf <- data.m.tf <- list()
ftf.grad <- ftf.feval <- ftf.joint.grad <- ftf.joint.feval <- rep(NA, n.rep)
mtf.grad <- mtf.feval <- mtf.joint.grad <- mtf.joint.feval <- rep(NA, n.rep)
fch.grad <- fch.feval <- fch.joint.grad <- fch.joint.feval <- rep(NA, n.rep)
mch.grad <- mch.feval <- mch.joint.grad <- mch.joint.feval <- rep(NA, n.rep)
## Range over which we optimize the onset. In this example, going up
## to 17 yo makes the simulation unstable: the loglikelihood under H1
## has a (suboptimal) local maximum in large values (for which we have
## few samples), leading to inaccurate (and sometimes negative)
## loglikelihood ratio statistics.
search.range <- c(6, 12) # data not available before 6 years old
search.range.TF <- search.range.CH <- search.range
## Main loop for simulations
for(rr in 1:n.rep){
print(rr)
## Simulate data from the Chizé population
data.f.ch[[rr]] <- sfunc(dataFCH, b1, thr.f, bw.offset, age.effect,
age.last.effect, last.effect, id.s2, coh.s2, s2)
data.m.ch[[rr]] <- sfunc(dataMCH, b1, thr.m, bw.offset, age.effect,
age.last.effect, last.effect, id.s2, coh.s2, s2)
## Simulate data from the Trois Fontaines population
data.f.tf[[rr]] <- sfunc(dataFTF, b1, thr.f, bw.offset, age.effect,
age.last.effect, last.effect, id.s2, coh.s2, s2)
data.m.tf[[rr]] <- sfunc(dataMTF, b1, thr.m, bw.offset, age.effect,
age.last.effect, last.effect, id.s2, coh.s2, s2)
## Compute the likelihood ratio test for this Trois Fontaines simulation
test.TF <- tryCatch({res=onset.test(ll.sim, data.f.tf[[rr]],
data.m.tf[[rr]], search.range.TF, CI.lvl=NA)
res$warn=FALSE
res},
warning=function(w) {
res <- onset.test(ll.sim, data.f.tf[[rr]],
data.m.tf[[rr]], search.range.TF, CI.lvl=NA)
res$warn <- TRUE
return(res)
})
llr.tf[rr] <- test.TF$llr
lh1.tf[rr] <- test.TF$lh1
lh0.tf[rr] <- test.TF$lh0
pv.tf[rr] <- test.TF$pv
cvg.tf[rr] <- test.TF$cvg.ok
warn.tf[rr] <- test.TF$warn
## Compute the likelihood ratio test for this Chizé simulation
test.CH <- tryCatch({res=onset.test(ll.sim, data.f.ch[[rr]],
data.m.ch[[rr]], search.range.CH, CI.lvl=NA)
res$warn=FALSE
res},
warning=function(w) {
res <- onset.test(ll.sim, data.f.ch[[rr]],
data.m.ch[[rr]], search.range.CH, CI.lvl=NA)
res$warn <- TRUE
return(res)
})
llr.ch[rr] <- test.CH$llr
lh1.ch[rr] <- test.CH$lh1
lh0.ch[rr] <- test.CH$lh0
pv.ch[rr] <- test.CH$pv
cvg.ch[rr] <- test.CH$cvg.ok
warn.ch[rr] <- test.CH$warn
## Optimality check: amplitude of the gradient (should be close to
## 0) and number of function evaluation (if equal to the largest
## allowed value, it is likely that the optimization did not
## converge).
ftf.lm <- lmer(body.mass.sim ~ b1(age, test.TF$est.1) +
age.at.last.capture + last.year.of.capture +
(1|ID) + (1|CohortF), data=data.f.tf[[rr]], REML='FALSE')
mtf.lm <- lmer(body.mass.sim ~ b1(age, test.TF$est.2) +
age.at.last.capture + last.year.of.capture +
(1|ID) + (1|CohortF), data=data.m.tf[[rr]], REML='FALSE')
ftf.joint.lm <- lmer(body.mass.sim ~ b1(age, test.TF$est.joint) +
age.at.last.capture + last.year.of.capture +
(1|ID) + (1|CohortF), data=data.f.tf[[rr]], REML='FALSE')
mtf.joint.lm <- lmer(body.mass.sim ~ b1(age, test.TF$est.joint) +
age.at.last.capture + last.year.of.capture +
(1|ID) + (1|CohortF), data=data.m.tf[[rr]], REML='FALSE')
fch.lm <- lmer(body.mass.sim ~ b1(age, test.CH$est.1) +
age.at.last.capture + last.year.of.capture +
(1|ID) + (1|CohortF), data=data.f.ch[[rr]], REML='FALSE')
mch.lm <- lmer(body.mass.sim ~ b1(age, test.CH$est.2) +
age.at.last.capture + last.year.of.capture +
(1|ID) + (1|CohortF), data=data.m.ch[[rr]], REML='FALSE')
fch.joint.lm <- lmer(body.mass.sim ~ b1(age, test.CH$est.joint) +
age.at.last.capture + last.year.of.capture +
(1|ID) + (1|CohortF), data=data.f.ch[[rr]], REML='FALSE')
mch.joint.lm <- lmer(body.mass.sim ~ b1(age, test.CH$est.joint) +
age.at.last.capture + last.year.of.capture +
(1|ID) + (1|CohortF), data=data.m.ch[[rr]], REML='FALSE')
ftf.grad[rr] <- max(abs(ftf.lm@optinfo$derivs$gradient))
ftf.feval[rr] <- ftf.lm@optinfo$feval
mtf.grad[rr] <- max(abs(mtf.lm@optinfo$derivs$gradient))
mtf.feval[rr] <- mtf.lm@optinfo$feval
ftf.joint.grad[rr] <- max(abs(ftf.joint.lm@optinfo$derivs$gradient))
ftf.joint.feval[rr] <- ftf.joint.lm@optinfo$feval
mtf.joint.grad[rr] <- max(abs(mtf.joint.lm@optinfo$derivs$gradient))
mtf.joint.feval[rr] <- mtf.joint.lm@optinfo$feval
fch.grad[rr] <- max(abs(fch.lm@optinfo$derivs$gradient))
fch.feval[rr] <- fch.lm@optinfo$feval
mch.grad[rr] <- max(abs(mch.lm@optinfo$derivs$gradient))
mch.feval[rr] <- mch.lm@optinfo$feval
fch.joint.grad[rr] <- max(abs(fch.joint.lm@optinfo$derivs$gradient))
fch.joint.feval[rr] <- fch.joint.lm@optinfo$feval
mch.joint.grad[rr] <- max(abs(mch.joint.lm@optinfo$derivs$gradient))
mch.joint.feval[rr] <- mch.joint.lm@optinfo$feval
if(rr == n.h0){
thr.f = thr.m + thr.delta
}
}
@
\subsection{Numerical optimization: troubleshooting}
Our test statistic relies on the numerical optimization of the
likelihood defined by \Rfunction{ll.sim}. Numerical optimization may fail
to converge, sometimes leading to negative log-likelihood
ratios. \Rfunction{onset.test} detects such cases and sets a
\texttt{cvg.ok} flag to \texttt{FALSE}:
<<cvg, echo=TRUE>>=
cat(sprintf('Negative log-likelihood obtained in proportion %g of the Trois
Fontaine and %g of the Chizé simulations', mean(!cvg.tf), mean(!cvg.ch)))
@
In our case, most parameters are optimized using descent methods as
implemented in the \Rfunction{lmer} function of \Rpackage{lme4}. The
age at the onset of senescence is optimized by the
\Rfunction{optimize} function of the \Rpackage{stat} package. Both can
fail and cause problems to the test computation. We use our simulation
to show how to detect and fix these problems.
\subsubsection{Descent methods in \Rfunction{lmer}}
Looking at the \Rfunction{summary} of the gradient $\ell_\infty$ norms
and number of function evaluation can help detect problems in
numerical optimization:
% save(file='../data/pre-computed-simulation.RData', failed.f.tf,
% failed.m.tf, ftf.grad, mtf.grad, ftf.joint.grad, mtf.joint.grad,
% fch.grad, mch.grad, fch.joint.grad, mch.joint.grad, ftf.feval,
% mtf.feval, ftf.joint.feval, mtf.joint.feval, fch.feval, mch.feval,
% fch.joint.feval, mch.joint.feval, search.range.TF, cvg.tf, cvg.ch,
% llr.ch, llr.tf, pv.ch, pv.tf, n.h0, res.tf, res.ch)
<<optim, echo=TRUE>>=
summary(ftf.grad)
summary(mtf.grad)
summary(ftf.joint.grad)
summary(mtf.joint.grad)
summary(fch.grad)
summary(mch.grad)
summary(fch.joint.grad)
summary(mch.joint.grad)
summary(ftf.feval)
summary(mtf.feval)
summary(ftf.joint.feval)
summary(mtf.joint.feval)
summary(fch.feval)
summary(mch.feval)
summary(fch.joint.feval)
summary(mch.joint.feval)
@
In our case, all gradients have small norm and optimization always
stopped long before reaching the maximum number of evaluations
(default $10000$). If you detect a large gradient norm or number of
evaluation, the main workarounds would be to
\begin{itemize}
\item Change the stopping conditions of the descent method
(\emph{e.g.} increase the number of maximal iterations).
\item Change the model, which may be overdetermined and lead to
numerically unstable results.
\end{itemize}
\subsubsection{\Rfunction{optimize} function}
In our experience, the negative log-likelihoods do not occur because
of the descent methods: both the gradient norm and the number of
function evaluations took reasonable values. However, plotting the
log-likelihood profile with respect to the age at the onset of
senescence reveals that these problems correspond to cases where the
log-likelihood is non-concave with respect to this variable,
suggesting that \Rfunction{optimize} returns a suboptimal local
maximum.
<<load-failed-run, echo=TRUE, eval=FALSE>>=
failed.simulation <- which(!cvg.tf)[1]
failed.f.tf <- data.f.tf[[failed.simulation]]
failed.m.tf <- data.m.tf[[failed.simulation]]
@
<<profile-simu-plot, echo=TRUE, include=FALSE, fig=TRUE, width=12, height=4>>=
res.failed.tf <- onset.test(ll.sim, failed.f.tf, failed.m.tf,
search.range.TF, do.plot=TRUE)
@
\begin{figure}[ht]
\begin{center}
\includegraphics[width=\textwidth]{RoeDeerMass-simulation-profile-simu-plot}
\caption{\label{fig:llprof-simu}Log-likelihood profiles on a simulation
producing a negative log-likelihood ratio statistic. For male individuals
(middle plot), the \Rfunction{optimize} function returned a local maximum.}
\end{center}
\end{figure}
In this particular case, the log-likelihood profile of the male
population has two local maxima (Figure~\ref{fig:llprof-simu}, middle
panel). Log-likelihood profile visualization can easily be used on any
dataset to ensure \Rfunction{optimize} does not return a suboptimal
solution --- by setting \texttt{do.plot=TRUE} in the
\Rfunction{onset.test} function. If it does, it it possible to re-run
it after adjusting the search range to remove non-global
maxima. Figure~\ref{fig:llprof-real} shows that the data
of~\cite{Douhard2017cost} lead to regular log-likelihood profiles. The
profiles sometimes show smalls local maxima (\emph{e.g.}, for females
in Trois Fontaines) but \Rfunction{optimize} which finds the global
maximum --- represented by the red star.
\subsection{Results}
Figure~\ref{fig:qqplot} shows a Q-Q plot of the empirical quantiles of
the log-likelihood ratio statistics obtained in our simulations under
$\HH{0}$ against the quantiles of a $\chi_2^1$ random variable,
showing that the statistic is approximately
$\chi_2^1$-distributed.
<<QQPlot, echo=TRUE, results=hide, fig=TRUE, include=FALSE, width=12, heigh=6>>=
oldpar <- par()
par(mfrow=c(1, 2), pty="s", mar=c(2, 6, 1, 1)+0.1)
qqplot(rchisq(1e4, 1), llr.ch[1:n.h0], main='Chizé', pch=20,
xlab=expression(chi[2](1) ~ 'quantiles'), ylab=expression(
'Empirical quantiles of the log-likelihood ratio test statistic under'~H[0]))
abline(a=0, b=1, col='red')
qqplot(rchisq(1e4, 1), llr.tf[1:n.h0], main='Trois Fontaines', pch=20,
xlab=expression(chi[2](1) ~ 'quantiles'), ylab=expression(
'Empirical quantiles of the log-likelihood ratio test statistic under'~H[0]))
abline(a=0, b=1, col='red')
par(oldpar)
@
\begin{figure}[ht]
\begin{center}
\includegraphics[width=\textwidth]{RoeDeerMass-simulation-QQPlot}
\caption{\label{fig:qqplot}Q-Q plot of the empirical quantiles of the log-likelihood
ratio test statistic under $\HH{0}$ against the quantiles of a
$\chi_2^1$ variable.}
\end{center}
\end{figure}
Accordingly, Figure~\ref{fig:calibration}
shows that for all $\alpha\in [0, 0.1]$, the proportion of the
experiments simulated under $\HH{0}$ yielding a p-value smaller than
$\alpha$ is close to $\alpha$.
<<calibration, echo=TRUE, results=hide, fig=TRUE, include=FALSE, width=6, heigh=6>>=
ch.thr <- unique(sort(pv.ch))
tf.thr <- unique(sort(pv.tf))
ch.lvl <- ch.pwr <- rep(-1, length(ch.thr))
tf.lvl <- tf.pwr <- rep(-1, length(tf.thr))
for(tt in 1:length(ch.thr)){
ch.lvl[tt] <- mean((pv.ch[1:n.h0] <= ch.thr[tt]))
ch.pwr[tt] <- mean((pv.ch[-(1:n.h0)] <= ch.thr[tt]))
}
for(tt in 1:length(tf.thr)){
tf.lvl[tt] <- mean((pv.tf[1:n.h0] <= tf.thr[tt]))
tf.pwr[tt] <- mean((pv.tf[-(1:n.h0)] <= tf.thr[tt]))
}
par(pty="s", mar=c(2, 4, 0, 2) + 0.1)
plot(ch.thr, ch.lvl, cex.lab=1.5,
xlab='P-value threshold',
ylab='False positive rate',
col='blue', type="s", lwd=2,
xlim=c(0, 0.1), ylim=c(0, 0.1))
lines(tf.thr, tf.lvl, col='red', type='l',lwd=2)
abline(a=0, b=1)
legend("bottomright", c('Chizé', 'Trois Fontaines'), lwd=2,
col=c('blue', 'red'))
@
\begin{figure}[ht]
\begin{center}
\includegraphics[width=0.4\textwidth]{RoeDeerMass-simulation-calibration}
\caption{\label{fig:calibration}Calibration plots of the proportion of
simulation under $\HH{0}$ yielding a p-value below a threshold against
this threshold.}
\end{center}
\end{figure}
Finally, Figure~\ref{fig:roc} shows that our test has some power to
detect changes in the age at the onset of senescence --- obviously
depending on the signal to noise ratio used in the simulation setting.
<<ROC, echo=TRUE, results=hide, fig=TRUE, include=FALSE, width=6, heigh=6>>=
## Plot calibration and ROC curve
par(pty="s", mar=c(2, 4, 0, 2) + 0.1)
plot(ch.lvl, ch.pwr, cex.lab=1.5,
xlab='False positive rate',
ylab='True positive rate',
col='blue', type="s", lwd=2)
lines(tf.lvl, tf.pwr, col='red', type='l', lwd=2)
abline(a=0, b=1)
legend("bottomright", c('Chizé', 'Trois Fontaines'), lwd=2,
col=c('blue', 'red'))
@
\begin{figure}[ht]
\begin{center}
\includegraphics[width=0.4\textwidth]{RoeDeerMass-simulation-ROC}
\caption{\label{fig:roc}ROC curves: true positive rate versus false
positive rate.}
\end{center}
\end{figure}
\section{Session Information}
<<sessionInfo, echo=FALSE>>=
sessionInfo()
@
\bibliographystyle{plainnat}
\bibliography{bibli}
\end{document}