%\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. <>= install.packages("OnAge") @ <>= 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(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. <>= ## 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: <>= search.range <- c(6, 12) search.range.TF <- search.range.CH <- search.range @ <>= res.tf <- onset.test(ll.real, dataFTF, dataMTF, search.range.TF, do.plot=TRUE) @ <>= res.ch <- onset.test(ll.real, dataFCH, dataMCH, search.range.CH, do.plot=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}. <>= 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. <>= # data('pre-computed-simulation') load(file = "pre-computed-simulation.RData") @ <>= ## 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}: <>= 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) <>= 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. <>= failed.simulation <- which(!cvg.tf)[1] failed.f.tf <- data.f.tf[[failed.simulation]] failed.m.tf <- data.m.tf[[failed.simulation]] @ <>= 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. <>= 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$. <>= 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. <>= ## 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() @ \bibliographystyle{plainnat} \bibliography{bibli} \end{document}