\documentclass[11pt,notitlepage]{article}
\usepackage{indentfirst}
\usepackage{amsmath}
\usepackage{amssymb}
\usepackage{natbib}
\usepackage{url}
\usepackage{graphicx}
\usepackage{tikzcd}
\usepackage[utf8]{inputenc}
\usepackage{xspace}
% remove room for headers, add to textheight
%\addtolength{\textheight}{\headheight}
%\addtolength{\textheight}{\headsep}
%\setlength{\headheight}{0 pt}
%\setlength{\headsep}{0 pt}
% adjust so right and left hand pages have different margins
% not sure why these particular numbers were picked or if they
% still make sense
% \showthe\evensidemargin
% \showthe\oddsidemargin
% \showthe\textwidth
% \evensidemargin 15 pt
% \oddsidemargin 61.5 pt
% \textwidth 392.5 pt
%\evensidemargin 28.75 pt
%\oddsidemargin 75.25 pt
%\textwidth 365 pt
%%%%% NEW go with 1.25 inch margin on all sides
%\setlength{\textheight}{\paperheight}
%\addtolength{\textheight}{ 2 in}
%\setlength{\topmargin}{0.25 pt}
%\setlength{\headheight}{0 pt}
%\setlength{\headsep}{0 pt}
%\addtolength{\textheight}{ \topmargin}
%\addtolength{\textheight}{ \topmargin}
%\addtolength{\textheight}{ \footskip}
%\setlength{\oddsidemargin}{0.25 in}
%\setlength{\evensidemargin}{0.25 in}
%\addtolength{\textwidth}{ \oddsidemargin}
%\addtolength{\textwidth}{ \evensidemargin}
\newcommand{\R}{\mathbb{R}}
\newcommand{\E}{\text{E}}
\newcommand{\Var}{\text{Var}}
\newcommand{\Fi}{\Sigma^{1}}
\newcommand{\eFi}{\hat{\Sigma}^{1}}
\newcommand{\Rsoft}{\texttt{R}\;}
\newcommand{\astertwo}{\texttt{aster2}\xspace}
\newcommand{\asterone}{\texttt{aster}\xspace}
\newcommand{\betahat}{\hat{\beta}}
\newcommand{\tauhat}{\hat{\tau}}
\begin{document}
\vspace*{0.9375in}
\begin{center}
{\bfseries Supporting Data Analysis for ``An Integrated Analysis \\
of Phenotypic Selection on Insect Body Size \\
and Development Time''} \\
By \\
Daniel J. Eck, Ruth G. Shaw, \\
Charles J. Geyer, and Joel G. Kingsolver \\
Technical Report No.~698 \\
School of Statistics \\
University of Minnesota \\
May 13, 2015 \\
Revised \today
\end{center}
\thispagestyle{empty}
\cleardoublepage
\setcounter{page}{1}
\thispagestyle{empty}
\begin{abstract}
This technical report (TR) gives details of the data analysis backing
up a paper having the same authors as this TR and having the title that
is quoted in the title of this TR. It uses the R package \texttt{knitr}
to process its source file (Rnw file) containing \LaTeX\ and R,
so that all R that
appears in this document is not cutandpasted but actually run and
actually produces the results that it appears to produce. This document
reproduces all tables and figures from the paper and hence the entire data
analysis reported in the paper. The data for the analysis are in the
R object \texttt{hornworm} in the R package \texttt{aster2}, which is a
CRAN package. An appendix of this technical report shows how this R
object of class \texttt{"asterdata"} was created from the raw data,
a comma separated values (CSV) file.
\end{abstract}
\section{Introduction}
In the paper that this technical report supplements, we describe
a new generation of aster models \citep*{aster1} that incorporate age of
first reproduction in the model for fitness.
Aster models are implemented in the R statistical computing language
\citep{rcore} via the CRAN package \asterone \citep[originally
available in 2005]{asterpackage}. This package does not allow
accounting for variation in generation time in the expression of fitness.
To allow for this component of fitness, the CRAN package \astertwo
\citep{aster2package} has
the added capability of specifying
``dependence groups,'' such that
life history stages, e.~g., insect larval instars, pupa, and adult,
and, in particular, variation in the age at which individuals reach these stages, are included in the model for fitness.
Here we use \texttt{aster2} to take these stages into account
in conducting selection analyses for
a dataset from a previous study of phenotypic selection on body size
and age at different developmental stages in \emph{Manduca sexta}
\citep{kingsolveretal}.
That study estimated selection on each fitness component separately,
and therefore could not quantify how selection operates over the lifecycle,
nor identify potential tradeoffs in selection among fitness components.
Our aster analyses of these data provide an integrated view of phenotypic
selection on size and age across development and fitness components in
this study system. We discuss how these methods can be applied to
selection analyses in other systems.
\section{R}
<>=
foo < getRversion()
bar < packageVersion("aster2")
qux < packageVersion("knitr")
@
This document was processed with R, version \Sexpr{foo}, using the CRAN
packages \texttt{aster2}, version \Sexpr{bar},
and \texttt{knitr}, version \Sexpr{qux}.
<>=
library(aster2)
@
\section{Data}
Data are in the R object \texttt{hornworm} in the \texttt{aster2} package.
We are also going to use \texttt{data} as a synonym of this data object
(because that was what we called it when we started this analysis,
the name \texttt{hornworm} came later).
<>=
data(hornworm)
data < hornworm
@
We show how this data set, an R object of class \texttt{"asterdata"},
was created in the appendix.
\section{Graph}
The aster graph for females in these data is shown below.
\label{graph}
\begin{center}
\setlength{\unitlength}{0.25 cm}
\begin{picture}(47,15)(10,0)
\scriptsize
\put(10,5){\makebox(0,0){$1$}}
\put(15,5){\makebox(0,0){$P$}}
\put(11.25,5){\vector(1,0){2.5}}
\put(20,0){\makebox(0,0){$T_{33_0}$}}
\put(20,5){\makebox(0,0){$T_{33_1}$}}
\put(20,10){\makebox(0,0){$T_{33_2}$}}
\put(20,6.25){\line(0,1){2.5}}
\put(20,3.75){\line(0,1){2.5}}
\put(16.25,3.75){\vector(1,1){2.5}}
\put(16.25,5){\vector(1,0){2.5}}
\put(16.25,6.25){\vector(1,1){2.5}}
\put(20,15){\makebox(0,0){$B_{33}$}}
\put(20,11.25){\vector(0,1){2.5}}
\put(25,5){\makebox(0,0){$T_{34_1}$}}
\put(25,10){\makebox(0,0){$T_{34_2}$}}
\put(25,6.25){\line(0,1){2.5}}
\put(21.25,5){\vector(1,0){2.5}}
\put(21.25,6.25){\vector(1,1){2.5}}
\put(25,15){\makebox(0,0){$B_{34}$}}
\put(25,11.25){\vector(0,1){2.5}}
\put(30,5){\makebox(0,0){$T_{35_1}$}}
\put(30,10){\makebox(0,0){$T_{35_2}$}}
\put(30,6.25){\line(0,1){2.5}}
\put(26.25,5){\vector(1,0){2.5}}
\put(26.25,6.25){\vector(1,1){2.5}}
\put(30,15){\makebox(0,0){$B_{35}$}}
\put(30,11.25){\vector(0,1){2.5}}
\put(35,5){\makebox(0,0){$T_{36_1}$}}
\put(35,10){\makebox(0,0){$T_{36_2}$}}
\put(35,6.25){\line(0,1){2.5}}
\put(31.25,5){\vector(1,0){2.5}}
\put(31.25,6.25){\vector(1,1){2.5}}
\put(35,15){\makebox(0,0){$B_{36}$}}
\put(35,11.25){\vector(0,1){2.5}}
\put(40,5){\makebox(0,0){$T_{37_1}$}}
\put(40,10){\makebox(0,0){$T_{37_2}$}}
\put(40,6.25){\line(0,1){2.5}}
\put(36.25,5){\vector(1,0){2.5}}
\put(36.25,6.25){\vector(1,1){2.5}}
\put(40,15){\makebox(0,0){$B_{37}$}}
\put(40,11.25){\vector(0,1){2.5}}
\put(45,5){\makebox(0,0){$T_{38_1}$}}
\put(45,10){\makebox(0,0){$T_{38_2}$}}
\put(45,6.25){\line(0,1){2.5}}
\put(41.25,5){\vector(1,0){2.5}}
\put(41.25,6.25){\vector(1,1){2.5}}
\put(45,15){\makebox(0,0){$B_{38}$}}
\put(45,11.25){\vector(0,1){2.5}}
\put(50,5){\makebox(0,0){$T_{39_1}$}}
\put(50,10){\makebox(0,0){$T_{39_2}$}}
\put(50,6.25){\line(0,1){2.5}}
\put(46.25,5){\vector(1,0){2.5}}
\put(46.25,6.25){\vector(1,1){2.5}}
\put(50,15){\makebox(0,0){$B_{39}$}}
\put(50,11.25){\vector(0,1){2.5}}
\put(55,5){\makebox(0,0){$T_{40_1}$}}
\put(55,10){\makebox(0,0){$T_{40_2}$}}
\put(55,6.25){\line(0,1){2.5}}
\put(51.25,5){\vector(1,0){2.5}}
\put(51.25,6.25){\vector(1,1){2.5}}
\put(55,15){\makebox(0,0){$B_{40}$}}
\put(55,11.5){\vector(0,1){2.5}}
\end{picture}
\end{center}
The aster graph for males and unknown (sex not determined because died
before pupation) is the same except there are no $B_x$ nodes.
Arrows go from predecessor nodes to successor nodes.
Lines (that are not arrows) link dependence groups.
Nodes are labeled by their associated variables.
$P$ node is pupation indicator, $T$ nodes are survival and eclosion indicators,
$B$ nodes are ovariole counts.
Subscripts indicate age (in days), subsubscripts indicate
variables in the same dependence group (0 = death,
1 = surviving but preeclosion, 2 = eclosion at this time).
All deaths before reproduction were modeled as occurring on day 33
because the day of death for individuals who died after pupation but
before eclosion was not recorded.
As in all aster graphs, predecessor equals zero implies successor equals zero.
If the variable at any node is zero, then the variables at all downstream
nodes are also zero (all nodes one gets to by following arrows from a node
where the variable is zero). Thus at most one variable of any $T_x$ dependence
group is nonzero, and the nonzero $T_{x_y}$ variables indicate the path
the life history for an individual takes. Similarly, at most one $B_x$
variable is nonzero (because at most one $T_{x_2}$ variable is nonzero),
the one where $x$ is the day the individual actually eclosed (if it did).
Our measure of Darwinian fitness for this graph is at first
(Section~\ref{sec:selection} below) total ovariole
count, the sum of all the $B_x$ nodes. Since at most one of these
nodes is nonzero (as explained above), this is just the number of
ovarioles the individual had if it survived to eclosion (and zero otherwise).
Later (Section~\ref{sec:grow} below) our measure of Darwinian fitness
is a weighted sum of ovariole counts, weighted according to
age and population growth rate, with weights given by
\eqref{eq:lambdaweights}. This is fitness adjusted for
population growth rate.
\section{A Digression on Aster Model Theory}
Because a referee of the paper this technical report accompanies asked for it,
we give a brief explanation of aster models and their theory.
This section can be skipped by readers who want to get on to the analysis.
The full aster graph has a node for every measured component of fitness
for every individual and also a constant node for every individual.
The figure on p.~\pageref{graph} is only the part of the full graph for
one female individual. The node marked 1 indicates that this part of
the graph is for one individual (in some data sets data on multiple
individuals is lumped together and then the constant node says how many
individuals that is).
The \emph{response vector} $y$ of an aster model
contains every measured component of fitness for every individual,
one component for every node in the full graph except for constant nodes.
Every component of $y$ belongs to exactly one dependence group $G$.
Some components are in dependence groups by themselves.
For example, in the figure on p.~\pageref{graph} the nodes
$T_{33_0}$, $T_{33_1}$, and $T_{33_2}$ form a dependence group of size 3,
the nodes $T_{34_0}$ and $T_{34_1}$ form a dependence group of size 2,
the node $B_{35}$ forms a dependence group of size 1 (a group by
itself), and the node $P$ is another group by itself.
Each dependence group $G$ has exactly one predecessor node $q(G)$.
In the figure on p.~\pageref{graph}
the predecessor of $\{ T_{33_0}, T_{33_1}, T_{33_2} \}$ is $P$,
the predecessor of $\{ T_{34_0}, T_{34_1} \}$ is $T_{33_1}$,
the predecessor of $B_{35}$ is $T_{35_2}$, and
the predecessor of $P$ is the constant node.
The \emph{saturated aster model} has one parameter for every component
of the response vector.
The part of the response vector $y_G$ containing the components in dependence
group $G$ is exponential family with sample size $q(G)$ and hence makes
contribution to the log likelihood
$$
y_G^T \theta_G  y_{q(G)} c_G(\theta_G)
$$
(\citealp[equation~(2),]{aster1} and \citealp[Section~2.4]{phil}).
The vector $\theta$ whose parts are $\theta_G$ is the
\emph{conditional canonical parameter vector}. The function $c_G$ is
the \emph{cumulant function} for dependence group $G$. The whole aster
model is also an exponential family (not just the conditional distributions
of its dependence groups). The whole log likelihood is the sum over
all dependence groups for all individuals
$$
l(\theta) = \sum_{G \in \mathcal{G}} \bigl[
y_G^T \theta_G  y_{q(G)} c_G(\theta_G) \bigr]
$$
where $\mathcal{G}$ is the family of all dependence groups. To see that
this has exponential family form, we rewrite it as
$$
l(\theta)
=
\sum_{j \in J} y_j \left[ \theta_j

\sum_{\substack{G \in \mathcal{G} \\ q(G) = j}} c_G(\theta_G) \right]

\sum_{\substack{G \in \mathcal{G} \\ j \notin J}} y_{q(G)} c_G(\theta_G)
$$
where $J = \bigcup \mathcal{G}$ the set of all nonconstant nodes of the
full graph, and match it to the exponential family form
$$
l(\varphi) = y^T \varphi  c(\varphi)
$$
seeing that this works with
\begin{equation} \label{eq:astertransform}
\varphi_j = \theta_j 
\sum_{\substack{G \in \mathcal{G} \\ q(G) = j}} c_G(\theta_G)
\end{equation}
and
\begin{equation} \label{eq:cumulantfunction}
c(\varphi) =
\sum_{\substack{G \in \mathcal{G} \\ j \notin J}} y_{q(G)} c_G(\theta_G)
\end{equation}
The vector $\varphi$ having components $\varphi_j$ is the \emph{unconditional
canonical parameter vector}, and the function $c$ is the cumulant function
of the saturated model. Equation~\eqref{eq:astertransform} is called the
\emph{aster transform}. It defines an invertible infinitely differentiable
change of parameter \citep[Section~2.3]{aster1}. It may seem odd that
\eqref{eq:cumulantfunction} purports to define a function of $\varphi$ by
giving it as a function of $\theta$, but this is valid because of the
invertibility of the aster transform.
Neither $\theta$ nor $\varphi$ has a simple connection with the components
of $y$. For that we want mean value parameters.
The \emph{unconditional mean value parameter} is the
vector $\mu$ having components
$$
\mu_j = E(y_j)
$$
and the \emph{conditional canonical parameter} is the vector $\xi$ having
components
$$
\xi_j = E(y_j  y_{p(j)} = 1)
$$
where $p(j) = q(G)$ for the unique $G$ that contains $j$, so $y_{p(j)}$ is
the predecessor of $y_j$ like $y_{q(G)}$ is the predecessor of $y_G$.
It is a fundamental property of aster models that means multiply
\begin{equation} \label{eq:fundamental}
\mu_j = \xi_j \mu_{p(j)}
\end{equation}
\citep[eq.~(12)]{aster1}, so applying \eqref{eq:fundamental} recursively
we have
\begin{align*}
\mu_j & = \xi_j \xi_{p(j)} \mu_{p(p(j))}
\\
\mu_j & = \xi_j \xi_{p(j)} \xi_{p(p(j))} \mu_{p(p(p(j)))}
\\
\mu_j & = \xi_j \xi_{p(j)} \xi_{p(p(j))} \xi_{p(p(p(j)))} \mu_{p(p(p(p(j))))}
\end{align*}
and so forth going all the way back to when the $\mu$ on the righthand side
is at a constant node and hence is just the constant (the expectation of a
constant is that constant). In short, unconditional means are products
of conditional means.
We now have four parameterizations of the saturated model $\theta$, $\varphi$,
$\xi$, and $\mu$. All have their purposes. Any may be needed in some
application. In most aster models all of the transformations between any
of these parameters are invertible and infinitely differentiable and can
be calculated by the computer. For models having multinomial dependence
groups like our hornworm data, the $\theta$ and $\varphi$ parameterizations
are not identifiable because each part of the response vector $y_G$ for
a multinomial dependence group satisfies the linear constraint
$$
\sum_{j \in G} y_j = y_{q(G)}
$$
and any linear constraint satisfied by the canonical statistic of
an exponential family causes nonidentifiability
\citep[Theorem~1 and the following discussion]{geyergdor}.
This nonidentifiability is broken by fixing one component of $\theta_G$
for each multinomial dependence group $G$ if we are using the $\theta$
parameterization and similarly for $\varphi$. With that proviso
all of the parameter transformations are invertible and infinitely
differentiable.
Of course, saturated models are uninteresting for data analysis because
they have too many parameters. They are just the largest model within
which we seek parsimonious submodels. Modeling mean values as linear
functions of means, like linear models (LM) do, is nonsense because of the
constraints on means (conditional means for both Bernoulli and multinomial
are between zero and one, those for zerotruncated Poisson are nonnegative).
Even generalized linear models (GLM) don't do that. Instead GLM that are
exponential families (logistic regression and Poisson regression with log
link) model canonical parameters linearly. Aster models do the same.
The question is which canonical parameter vector $\theta$ or $\varphi$?
The best answer is $\varphi$ because that is the parameter that makes
the whole aster model an exponential family. An \emph{unconditional
canonical affine submodel} has
$$
\varphi = a + M \beta
$$
where $a$ is a known vector and $M$ a known matrix, called
the \emph{offset vector} and \emph{model matrix} in the terminology
of the R functions \texttt{lm} which fits LM and \texttt{glm} which fits GLM
and of the R package \texttt{aster2} (the R package \texttt{aster}
says ``origin'' instead of ``offset''), and where $\beta$ is the submodel
parameter vector. This makes the submodel
a full exponential family. The submodel log likelihood is
$$
l(\beta) = y^T (a + M \beta)  c(a + M \beta)
= y^T a + y^T M \beta  c(a + M \beta)
$$
and the term not containing the parameter $\beta$ can be dropped giving
$$
l(\beta) = y^T M \beta  c(a + M \beta)
= (M^T y)^T \beta  c_\text{sub}(\beta)
$$
and this shows we have another exponential family with canonical statistic
$M^T y$, canonical parameter $\beta$ and cumulant function $c_\text{sub}$.
The submodel mean value parameter is the expectation of the canonical statistic
$$
\tau = E(M^T y) = M^T \mu
$$
With the proviso about fixing one component of $\varphi$ or $\theta$ in
each multinomial dependence group and with the additional proviso that
$M$ is full rank, we have six identifiable parameterizations
$\theta$, $\varphi$, $\xi$, $\mu$, $\beta$, and $\tau$ and all of the
transformations between them are invertible and infinitely differentiable.
The following diagram shows these parameters and some of their transformations.
\begin{center}
\begin{tikzpicture}[descr/.style={fill=white},text height=1.5ex, text depth=0.25ex]
\node (a) at (0,0) {$\theta$};
\node (b) at (4,0) {$\varphi$};
\node (c) at (8,0) {$\beta$};
\node (d) at (0,3) {$\xi$};
\node (e) at (4,3) {$\mu$};
\node (f) at (8,3) {$\tau$};
\path[>,font=\scriptsize]
([yshift= 5pt]a.east) edge node[above] {\text{aster transform}} ([yshift= 5pt]b.west)
([yshift= 5pt]b.west) edge node[below] {\text{inverse aster transform}} ([yshift= 5pt]a.east)
([yshift= 5pt]b.east) edge node[above] {} ([yshift= 5pt]c.west)
([yshift= 5pt]c.west) edge node[below] {$\varphi = M\beta$} ([yshift= 5pt]b.east)
([xshift= 5pt]a.south) edge node[right] {$\xi_G = \nabla c_G(\theta_G)$} ([xshift= 5pt]d.north)
([xshift= 5pt]d.north) edge node[right] {} ([xshift= 5pt]a.south)
([xshift= 5pt]b.south) edge node[right] {$\mu = \nabla c(\varphi)$} ([xshift= 5pt]e.north)
([xshift= 5pt]e.north) edge node[right] {} ([xshift= 5pt]b.south)
([xshift= 5pt]c.south) edge node[right] {$\tau = \nabla c_\text{sub}(\beta)$} ([xshift= 5pt]f.north)
([xshift= 5pt]f.north) edge node[right] {} ([xshift= 5pt]c.south)
([yshift= 5pt]e.east) edge node[above] {$\tau = M^T\mu$} ([yshift= 5pt]f.west)
([yshift= 5pt]f.west) edge node[below] {} ([yshift= 5pt]e.east)
([yshift= 5pt]d.east) edge node[above] {\text{multiplication}} ([yshift= 5pt]e.west)
([yshift= 5pt]e.west) edge node[below] {\text{division}} ([yshift= 5pt]d.east);
\end{tikzpicture}
\label{diagram}
\end{center}
As the diagram shows, transformations from canonical to mean value parameters
are given by differentiating cumulant functions. The inverses of these
transformations have no closed form expression but can be done by solving
optimization problems, they are essentially equivalent to doing maximum
likelihood.
Maximum likelihood estimation for the submodel differentiates
the submodel log likelihood
$$
\nabla l(\beta) = M^T y  \nabla c_\text{sub}(\beta)
= M^T y  M^T E_\beta(y)
$$
sets the derivative equal to zero and solves for $\beta$. Thus the
maximum likelihood estimate (MLE) for $\beta$ satisfies
$$
M^T y = \nabla c_\text{sub}(\hat{\beta}) = M^T E_{\hat{\beta}}(y)
$$
By invariance of maximum likelihood the MLE for $\tau$ is what the MLE
for $\beta$ maps to under the parameter transformation
$$
\hat{\tau} = \nabla c_\text{sub}(\hat{\beta}) = M^T y
$$
Everything in this paragraph holds for every full exponential family
(nothing special about aster models here). MLE for them satisfy
the ``observed equals expected'' property (the MLE of the mean value parameter
is the observed value of the canonical statistic, in this case
$\hat{\tau} = M^T y$).
Exponential families give aster unconditional canonical affine submodels
many more good properties: the multivariate monotonicity property of maps
from canonical to mean value parameters (\citealp[sec.~2.9]{phil} and
\citealp[appendix]{aster3}), the sufficient dimension reduction property
(all MLE are sufficient statistics, \citealp[sec.~2.10]{phil}), and the
maximum entropy property \citep[sec.~2.12]{phil}. In particular, the
multivariate monotonicity property of the maps
$\varphi \longleftrightarrow \mu$ enables much important theory
\citep[appendix]{aster3}. All of these properties and more
\citep[discussion]{phil} make unconditional canonical affine submodels
by far the most useful.
Conditional canonical affine submodels model $\theta$ affinely
$$
\theta = a + M \beta
$$
These are not full exponential families. They are, of course curved
exponential families because they are smooth submodels of the saturated model.
But that doesn't give them any of the desirable properties of the unconditional
canonical affine submodels. Generally, conditional submodels are a lot less
parsimonious than unconditional submodels \citep{aster1,phil}. As far as
we know there are only two published examples of conditional submodels:
Example~1 in \citet{aster2}, which arguably could have been done with
an unconditional submodel, and \citet{asteraphid}, which had to use
a conditional submodel because of timedependent covariates (the
first essential use of such models). We ignore conditional submodels for the
rest of this technical report because the R package \texttt{aster2}
currently has no way to fit them.
\section{Model Fitting}
Fitting an aster model using the \astertwo package is not as easy as using
the \asterone package. The \astertwo package has no model fitting
function that takes a formula and fits a model, like the \texttt{aster}
function in the \asterone package. Nor does it have generic
\texttt{summary}, \texttt{anova}, and \texttt{predict} functions that support
doing hypothesis tests and confidence intervals.
% The reason is that when the \astertwo package was considered about half done,
% a new line of research (random effects aster models) opened up, and support
% for it was added to the \asterone package (because it was more advanced).
% So now the \astertwo package is less than a third done (it does not
% have random effects and needs them too). But research on random effects aster
% models continues, still delaying the completion of the \astertwo package.
The \astertwo package is the only way to do aster models with dependence
groups. Currently, for everything else, use the \asterone package.
The \astertwo package has parameter
transformation functions that do all transformations between
all parameterizations of aster models shown in the diagram
on p.~\pageref{diagram}.
In being able to do all of these transformations, the \texttt{aster2} package
is superior to the \texttt{aster} package, which only does some of them
(using the functions \texttt{predict} and \texttt{astertransform}).
But it does require more work to do maximum likelihood using the \astertwo
package than the \asterone package. The observed equals expected property
applied to aster models says that the MLE of the submodel mean value
parameter $\hat{\tau}$ is equal to the submodel sufficient statistic $M^T y$.
Then by invariance of maximum likelihood, the MLE of the submodel canonical
parameter $\hat{\beta}$ is the transformation of $\hat{\tau}$.
Thus maximum likelihood in the \astertwo package is done as follows.
\begin{enumerate}
\item Make the model matrix $M$.
The R function \texttt{model.matrix} can be used to do this.
\item Make the submodel sufficient statistic $M^T y = \hat{\tau}$.
\item Transform $\hat{\tau}$ to $\hat{\beta}$ using the R function
\texttt{transformUnconditional}.
\end{enumerate}
Then, the MLE having been obtained, there is no support for hypothesis
tests or confidence intervals. But the \astertwo package can calculate Fisher
information for any of these parameters. The Fisher information matrix
for $\beta$ is
$$
I(\beta) = \nabla^2 c_\text{sub}(\beta) = \nabla h(\beta)
$$
where $h$ is the transformation $\beta \longrightarrow \tau$.
The Fisher information matrix for $\tau$ is
$$
I(\tau) = \bigl[ \nabla^2 c_\text{sub}(\beta) \bigr]^{1}
= \nabla h^{1}(\tau)
$$
And the Jacobian matrices $\nabla h(\beta)$ and $\nabla h^{ 1}(\tau)$
are computed by the R function \texttt{jacobian}.
Then from Fisher information one can compute various hypothesis tests and
confidence intervals. In this document we only do Rao tests.
The \astertwo package does not have a function that calculates log likelihoods
(in the \asterone package both R functions \texttt{aster} and \texttt{mlogl}
do this), so we have no easy way to calculate likelihood ratios and likelihood
ratio test statistics. That is why we do not do likelihood ratio tests.
\section{Model selection} \label{sec:selection}
Our analysis selects a best model via backwards selection.
We start with a model that is full quadratic in all three covariates
(age at second instar, mass at second instar, and mass at eclosion).
We then show by doing a Rao test of model comparison that we can
drop quadratic terms involving age.
We then show by doing more Rao tests of model comparison that we cannot
drop any other terms.
To conduct such a test, the score function and inverse observed Fisher
information need to be calculated under the null model. The score function is
\begin{equation} \label{score}
M_\text{alter}^TY  \nabla c(M_\text{alter}\betahat_\text{null})
\end{equation}
where $M_\text{alter}$ and $\betahat_\text{null}$ are the model matrix for
the full quadratic model and the MLE for $\beta$ using the null model,
respectively (we do not have an offset vector in any of our models).
We first calculate the term on the left in \eqref{score}.
<<>>=
offspring < as.numeric(grepl("B",data$redata$varb))
modmat.alt < model.matrix(resp ~ varb +
offspring:(Mass_Repro + Mass_2nd +
I(Mass_Repro^2) + I(Mass_2nd^2) + I(Mass_2nd*Mass_Repro) +
Time_2nd + I(Time_2nd^2) + I(Mass_2nd*Time_2nd) +
I(Time_2nd*Mass_Repro)), data = data$redata)
tau.alt < crossprod(modmat.alt, data$redata$resp)
@
We now calculate $\betahat_\text{null}$ by mapping $\hat{\tau}_\text{null}$
to the $\beta$ parameterization using the \texttt{transformUnconditional}
function, this null hypothesis being the one in which \verb@Time_Repro@
is linear and the other two covariates are quadratic.
%The model matrix constructed here corresponds to an aster model with two linear terms for mass at the second instar larval stage and reproductive stage. The maximum likelihood estimate $\tauhat$ is $M^TY$. We obtain $\betahat$ by using the \texttt{transformUnconditional} function in the \astertwo package. $\betahat$ is needed indirectly for model comparison.
<>=
modmat.null < model.matrix(resp ~ varb +
offspring:(Mass_Repro + Mass_2nd +
I(Mass_Repro^2) + I(Mass_2nd^2) + I(Mass_2nd*Mass_Repro) +
Time_2nd), data = data$redata)
tau < crossprod(modmat.null, data$redata$resp)
beta < beta.null < transformUnconditional(tau, modmat.null,
data, from = "tau", to = "beta", tolerance = 1e100)
@
The term on the right in (\ref{score}) is calculated using the
\texttt{transformUnconditional} function.
<>=
tau.null < transformUnconditional(c(beta.null,rep(0,3)),
modmat.alt, data, from = "beta", to = "tau")
score < tau.alt  tau.null
@
Observed Fisher information is calculated using the \texttt{jacobian} function. Three zeros are added to $\hat{\beta}_\text{null}$ in order to calculate
observed Fisher information for the alternative model at the MLE for the
null model. The three zeros are for the terms
(the quadratic term for age at second instar and the two crossproduct terms
of age at second instar with one of the mass covariates) that are present
in the alternative but absent in the null (which is the same as being
set to zero).
<<>>=
Fisher.null < jacobian(c(beta.null,rep(0,3)), data,
transform = "unconditional", from = "beta",
to = "tau", modmat = modmat.alt)
@
We now construct the Rao test statistic. The reference distribution for this test is $\chi^2_3$ and the pvalue suggests that we fail to reject the null hypothesis when testing at the 0.05 significance level.
<<>>=
Rao < t(score) %*% solve(Fisher.null, tol = 1e100) %*% score
pchisq(Rao, df = 3, lower = FALSE)
@
All less complicated models are found to be insignificant when testing at the 0.05 significance level. Our final model includes the full quadratic structure between the two mass terms and a linear term for the age at which an individual \emph{M. sexta} reaches its second instar larval stage as useful predictors of unconditional expected ovariole counts. First, we consider removing the quadratic terms for mass at second instar. The pvalue is small enough to reject this smaller model at any reasonable significance level. The steps to conduct this test are similar to those that conducted our first hypothesis test.
<>=
modmat.null < model.matrix(resp ~ varb +
offspring:(Mass_Repro + I(Mass_Repro^2) + Time_2nd +
Mass_2nd),
data = data$redata)
tau < crossprod(modmat.null, data$redata$resp)
beta.null < transformUnconditional(tau, modmat.null,
data, from = "tau", to = "beta", tolerance = 1e100)
modmat.alt < model.matrix(resp ~ varb +
offspring:(Mass_Repro + I(Mass_Repro^2) + Time_2nd +
Mass_2nd + I(Mass_Repro*Mass_2nd) + I(Mass_2nd^2)),
data = data$redata)
tau.alt < crossprod(modmat.alt, data$redata$resp)
tau.null < transformUnconditional(c(beta.null,rep(0,2)),
modmat.alt, data, from = "beta", to = "tau")
score < tau.alt  tau.null
Fisher.null < jacobian(c(beta.null,rep(0,2)), data,
transform = "unconditional", from = "beta",
to = "tau", modmat = modmat.alt)
Rao < t(score) %*% solve(Fisher.null, tol = 1e100) %*% score
pchisq(Rao, df = 2, lower = FALSE)
@
We now consider removing the quadratic terms for mass at eclosion. The pvalue is small enough to reject this smaller model at any reasonable significance level.
<>=
modmat.null < model.matrix(resp ~ varb +
offspring:(Mass_2nd + I(Mass_2nd^2) + Time_2nd +
Mass_Repro),
data = data$redata)
tau < crossprod(modmat.null, data$redata$resp)
beta.null < transformUnconditional(tau, modmat.null,
data, from = "tau", to = "beta", tolerance = 1e100)
modmat.alt < model.matrix(resp ~ varb +
offspring:(Mass_2nd + I(Mass_2nd^2) + Time_2nd +
Mass_Repro + I(Mass_Repro^2) + I(Mass_Repro*Mass_2nd)),
data = data$redata)
tau.alt < crossprod(modmat.alt, data$redata$resp)
tau.null < transformUnconditional(c(beta.null,rep(0,2)),
modmat.alt, data, from = "beta", to = "tau")
score < tau.alt  tau.null
Fisher.null < jacobian(c(beta.null,rep(0,2)), data,
transform = "unconditional", from = "beta",
to = "tau", modmat = modmat.alt)
Rao < t(score) %*% solve(Fisher.null, tol = 1e100) %*% score
pchisq(Rao, df = 2, lower = FALSE)
@
Finally, we consider removing the linear term for age that individuals reach their second instar stage. The pvalue is small enough to reject this smaller model at any reasonable significance level.
<>=
modmat.null < model.matrix(resp ~ varb +
offspring:(Mass_2nd + I(Mass_2nd^2) + Mass_Repro +
I(Mass_Repro^2) + I(Mass_Repro*Mass_2nd)),
data = data$redata)
tau < crossprod(modmat.null, data$redata$resp)
beta.null < transformUnconditional(tau, modmat.null,
data, from = "tau", to = "beta", tolerance = 1e100)
modmat.alt < model.matrix(resp ~ varb +
offspring:(Mass_2nd + I(Mass_2nd^2) + Mass_Repro +
I(Mass_Repro^2) + I(Mass_Repro*Mass_2nd) + Time_2nd),
data = data$redata)
tau.alt < crossprod(modmat.alt, data$redata$resp)
tau.null < transformUnconditional(c(beta.null,0),
modmat.alt, data, from = "beta", to = "tau")
score < tau.alt  tau.null
Fisher.null < jacobian(c(beta.null,0), data,
transform = "unconditional", from = "beta",
to = "tau", modmat = modmat.alt)
Rao < t(score) %*% solve(Fisher.null, tol = 1e100) %*% score
pchisq(Rao, df = 1, lower = FALSE)
@
%<>=
%mymu.full < transformUnconditional(beta.full,
% modmat.full, data, from = "beta", to = "mu")
%mytau.full < crossprod(modmat.full, mymu.full)
%#all.equal(mytau.full, tau.full)
%names(mymu.full) < rownames(modmat.full)
%@
%<<>>=
%offspring < as.numeric(grepl("B",data$redata$varb))
%modmat < model.matrix(resp ~ varb + offspring:(Mass_2nd +
% Mass_Repro), data = data$redata)
%tau < crossprod(modmat, data$redata$resp)
%beta < transformUnconditional(tau, modmat,
% data, from = "tau", to = "beta")
%@
%We now consider a full quadratic model in both mass traits. There is no conventional software in the \astertwo software to compare models. The likelihood ratio test, which is used in \texttt{anova.aster}, can not be implemented since no \texttt{mlogl} function exists in the \astertwo package. Thus, a Rao test is conducted ``by hand" to see if either the linear model or the quadratic model is preferred. To conduct this test, the score function and inverse observed Fisher information need to be calculated under the null model, in this case the linear model. The score function is $M_{quad}^TY  \nabla c_{sub}(M_{quad}\betahat_\text{null})$ where $M_{quad}$ and $\betahat_\text{null}$ are the model matrix for the quadratic model and the MLE for $\beta$ using the linear model respectively. The score is calculated using the \texttt{transformUnconditional} function. This function maps $\beta$ to $\tau$ where $\tau$ is the aster submodel meanvalue parameter. Observed Fisher information is calculated using the \texttt{jacobian} function. Three zeros are added to the null model estimate of $\beta$ in order to calculate observed Fisher information. This is result of constraining the two quadratic terms and the single interaction term to be zero when maximizing the likelihood of the null model. These steps are explained in more detail below.
%We first calculate the MLE of the submodel meanvalue parameter from the full quadratic model.
%<<>>=
%modmat.quad < model.matrix(resp ~ varb +
% offspring:(Mass_2nd + Mass_Repro + I(Mass_2nd^2) +
% I(Mass_Repro^2) + I(Mass_2nd*Mass_Repro)),
% data = data$redata)
%tau.quad < crossprod(modmat.quad, data$redata$resp)
%@
%The MLE of the submodel meanvalue parameter from the linear model in both mass traits is now calculated. We transform our MLE for the submodel canonical parameter to the meanvalue parameterization under the null hypothesis that the quadratic terms are not present in our model. The quadratic terms are the last three elements of our submodel canonical parameter vector for the full model. Therefore we add three zeros to the end of \texttt{beta} before calling \texttt{transformUnconditional}. The score function needed for the Rao test is then calculated.
%<<>>=
%tau.null < transformUnconditional(c(beta,rep(0,3)),
% modmat.quad, data, from = "beta", to = "tau")
%score < tau.quad  tau.null
%@
%The \texttt{jacobian} function is then used to find the Jacobian of the map from $\beta$ to $\tau$ under the null hypothesis.
%<<>>=
%Fisher.null < jacobian(c(beta,rep(0,3)), data,
% transform = "unconditional", from = "beta",
% to = "tau", modmat = modmat.quad)
%@
%<>=
%foot < crossprod(modmat.quad, data$redata$resp)
%foot[c(29:31)] < 0
%blah < transformUnconditional(foot, tolerance = 1e30,
% modmat.quad, data, from = "tau", to = "beta")
%Fisher.null2 < jacobian(blah, data,
% transform = "unconditional", from = "beta",
% to = "tau", modmat = modmat.quad)
%score2 < foot  tau.quad
%Rao2 < t(score2) %*% Fisher.null2 %*% score2
%@
%We build the Rao test statistic and then calculated the pvalue where the reference distribution for this test is $\chi^2$ with three degrees of freedom. The pvalue from the Rao test is lower than any reasonable significance level. The full quadratic model is preferred over the simpler linear model.
%<<>>=
%Rao < t(score) %*% solve(Fisher.null, tol = 1e50) %*% score
%Rao
%df < 3
%pval < pchisq(Rao, df = df, lower.tail = FALSE)
%pval
%@
%We now wish to compare the full quadratic model to the model that is full quadratic in the two mass terms and has a linear term for the time it takes individuals to reach their second instar life stage. The MLE of the submodel canonical parameter vector for the full quadratic model is obtained using the \texttt{transformUnconditional} function. This is now an estimator of the submodel canonical parameter vector corresponding to the null hypothesis of our model selection procedure.
%<<>>=
%beta.quad < transformUnconditional(tau.quad,
% modmat.quad, data, from = "tau", to = "beta")
%@
%The exact same steps are executed to perform the Rao test that compares the two models of interest. The pvalue from the Rao test indicates that the larger model is preferred at any reasonable choice of significance level. This model includes a linear term for the time it takes individuals to reach their second instar life stage to the full quadratic model selected in the previous test.
%<<>>=
%modmat.full < model.matrix(resp ~ varb +
% offspring:(Mass_2nd + Mass_Repro + I(Mass_2nd^2) +
% I(Mass_Repro^2) + I(Mass_2nd*Mass_Repro) + Time_2nd),
% data = data$redata)
%tau.full < crossprod(modmat.full, data$redata$resp)
%beta.full < transformUnconditional(tau.full,
% modmat.full, data, from = "tau", to = "beta")
%tau.null < transformUnconditional(c(beta.quad,0),
% modmat.full, data, from = "beta", to = "tau")
%Fisher.null < jacobian(c(beta.quad,0), data, transform =
% "unconditional", from = "beta", to = "tau",
% modmat = modmat.full)
%score < tau.full  tau.null
%Rao < t(score) %*% solve(Fisher.null, tol = 1e50) %*% score
%Rao
%df < length(beta.full)  length(beta.quad)
%pval < pchisq(Rao, df = df, lower.tail = FALSE)
%pval
%@
%<>=
%mymu.full < transformUnconditional(beta.full,
% modmat.full, data, from = "beta", to = "mu")
%mytau.full < crossprod(modmat.full, mymu.full)
%#all.equal(mytau.full, tau.full)
%names(mymu.full) < rownames(modmat.full)
%@
%<<>>=
%modmat.bigger < model.matrix(resp ~ varb +
% offspring:(Mass_2nd + Mass_Repro + I(Mass_2nd^2) +
% I(Mass_Repro^2) + I(Mass_2nd*Mass_Repro) + Time_2nd +
% I(Time_2nd*Mass_Repro)),
% data = data$redata)
%canonical.bigger < crossprod(modmat.bigger,
% data$redata$resp)
%beta.bigger < transformUnconditional(canonical.bigger,
% modmat.bigger, data, from = "tau", to = "beta")
%tau.null < transformUnconditional(c(beta.full,0),
% modmat.bigger, data, from = "beta", to = "tau")
%Fisher.null < jacobian(c(beta.full,0), data,
% transform = "unconditional", from = "beta",
% to = "tau", modmat = modmat.bigger)
%score < canonical.bigger  tau.null
%Rao < t(score) %*% solve(Fisher.null, tol = 1e50) %*% score
%Rao
%df < length(beta.bigger)  length(beta.full)
%pval < pchisq(Rao, df = df, lower.tail = FALSE)
%pval
%@
The results of the tests that compare our final model to smaller models are concisely summarized in Table~\ref{tab:one}.
\begin{table}
\begin{tabular}{lcc}
\hline
\hline
%\hline
%\multicolumn{3}{c}{Rao tests for smaller models} \\
null model & df & $P$value \\
\hline
removes quadratic terms for mass at second instar & 2 & $3.37 \times 10^{8}$ \\
removes quadratic terms for mass at eclosion & 2 & $< 10^{10}$\\
removes linear term for age at second instar & 1 & $7.88 \times 10^{5}$\\
\hline
\end{tabular}
\caption{Rao tests for smaller models. $P$values and degrees of freedom
for Rao tests of three smaller models against the larger model that
includes linear, quadratic, and interaction term for the two mass traits and
a linear term for age at second larval instar stage.}
\label{tab:one}
\end{table}
\section{Fitness landscapes}
Now that we have selected the best model, we want to plot fitness landscape,
unconditional expected fitness as a function of covariates. We do this by
predicting, on the basis of the model, values for fitness for all observed
values of age at second instar and samples of 101 values
of mass at second instar and 101 values of mass at eclosion. Since this is
a fourdimensional graph (fitness versus three covariates), it cannot be
visualized, and we make one threedimensional graph for each age at
2nd instar using
the R function \texttt{contour}.
Covariate values that maximize expected fitness are also of interest.
Before we
do this we need to extract some useful information from the \texttt{hornworm}
dataset.
<>=
vars < levels(hornworm$redata$varb)
pred < hornworm$pred
group < hornworm$group
code < hornworm$code
families < hornworm$families
@
Now the $101^2$ hypothetical individuals with unique mass traits reaching
their second instar larval stage at age 2 are generated
<<>>=
x < seq(from = 0, to = 0.016, by = 0.016/100)
y < seq(from = 0, to = 2.3, by = 2.3/100)
days < 2:6 / 7
long.sex < factor(rep("F", 101^2), levels = c("F", "M", "U"))
mnew < data.frame(long.sex, rep(days[1], 101^2),
rep(x, times = 101), rep(y, each = 101))
mnew < cbind(mnew, matrix(rep(0, 26*101^2), ncol = 26))
mnames < c("Sex", "Time_2nd", "Mass_2nd", "Mass_Repro", vars)
names(mnew) < mnames
@
We then make the object of class \texttt{asterdata} for this new data
and the model matrix for our best model and these new data.
<>=
data.mnew < asterdata(mnew, vars = vars, pred = pred,
group = group, code = code, families = families)
foo < as.character(data.mnew$redata$varb)
offspring.mnew < as.numeric(grepl("B", foo))
data.mnew$redata < transform(data.mnew$redata,
offspring.mnew = offspring.mnew)
modmat.mnew < model.matrix(resp ~ varb +
offspring.mnew:(Mass_Repro + Mass_2nd +
I(Mass_Repro^2) + I(Mass_2nd^2) + I(Mass_2nd*Mass_Repro) +
Time_2nd), data = data.mnew$redata)
@
This model matrix is then used to find the MLE of the saturated model
unconditional canonical parameter $\varphi$.
<<>>=
phi < modmat.mnew %*% beta
@
Estimates on the canonical scale are not directly interpretable. We need estimates for the meanvalue parameters of the unconditional model ($\mu$). These are found using the \texttt{transformSaturated} function. Only the estimates corresponding to ovariole count nodes are of interest. For every hypothetical individual there are eight predictions of ovariole count, one for every age that the individual can reach its reproduction stage. These predictions are summed to arrive at an estimate of expected ovariole count for each of the $101^2$ hypothetical individuals.
Of course, in this biological system observed fitness for a single individual is nonzero at only one
age, the age at which the individual reproduces (at most one $Bx$ variable
in the graph on p.~\pageref{graph} is nonzero). But expected fitness is
nonzero at all ages (all components of the unconditional mean value parameter
vector $\mu$ are nonzero). Adding the components of $\mu$ for all the $Bx$
nodes of the graph adds the contribution to expected fitness of reproduction
at all ages.
<>=
mus2 < transformSaturated(phi, data.mnew, from = "phi",
to = "mu")
ind < which(grepl("B",rownames(data.frame(modmat.mnew))))
qux2 < matrix(mus2[ind], nrow = 101^2, ncol = 8)
colnames(qux2) < paste("c", 1:8, sep = "")
rownames(qux2) < c(1:101^2)
sums2 < as.numeric(apply(qux2, 1, sum))
zday2 < matrix(sums2, nrow = 101)
@
The \texttt{maximum} function, defined below,
is used to find the two mass values yielding the highest expected fitness.
<<>>=
maximum < function(sums){
max.ind < which(sums == max(sums))
column < floor(max.ind / 101)
row < max.ind  101 * column
point < c(x[row], y[column])
return(point)
}
@
To make the plot we need to pull more things out of the \texttt{hornworm}
dataset.
The data are in the \texttt{redata} component of the \texttt{asterdata} object.
We pull it out. This is an object in socalled long format (in the terminology
of the R function \texttt{reshape}) so all of the covariates are repeated
as many times as there are nodes in the graph. We get rid of this repetition
because we are only interested in covariates here. And then we keep the
data for females only.
<>=
mass.second < hornworm$redata$Mass_2nd
mass.reprod < hornworm$redata$Mass_Repro
sex < as.character(hornworm$redata$Sex)
unique(sex)
time.second < hornworm$redata$Time_2nd
id < data$redata$id
u < unique(id)
idx < match(u, id)
mass.second < mass.second[idx]
mass.reprod < mass.reprod[idx]
sex < sex[idx]
time.second < time.second[idx]
mass.second < mass.second[sex == "F"]
mass.reprod < mass.reprod[sex == "F"]
@
The fitness landscape for these $101^2$ hypothetical individuals
reaching the second instar stage at age 2 is plotted in Figure~\ref{fig:day2}.
\begin{figure}
\begin{center}
<<>>=
plot(mass.second, pch = 20,
mass.reprod, xlab = "Mass at
2nd instar", ylab = "Mass at Reproduction",
main = "Fitness Landscape for Ovariole Counts vs. Mass")
points(maximum(sums2)[1], maximum(sums2)[2], col = "red",
pch = 19)
contour(x,y,zday2, add = TRUE, nlevels = 8, levels =
c(100,150,200,250,300,325,350,360) )
@
\end{center}
\caption{Fitness Landscape for Those Reaching Second Instar at Day 2.}
\label{fig:day2}
\end{figure}
The code below constructs the fitness landscape plots for hypothetical
individuals reaching the second instar larval stage at ages 3, 4, and 5.
<>=
mnew < data.frame(long.sex, rep(days[2], 101^2),
rep(x, times = 101), rep(y, each = 101))
mnew < cbind(mnew, matrix(rep(0, 26*101^2), ncol = 26))
mnames < c("Sex", "Time_2nd", "Mass_2nd", "Mass_Repro", vars)
names(mnew) < mnames
data.mnew < asterdata(mnew, vars = vars, pred = pred,
group = group, code = code, families = families)
foo < as.character(data.mnew$redata$varb)
offspring.mnew < as.numeric(grepl("B", foo))
data.mnew$redata < transform(data.mnew$redata,
offspring.mnew = offspring.mnew)
modmat.mnew < model.matrix(resp ~ varb +
offspring.mnew:(Mass_Repro + Mass_2nd +
I(Mass_Repro^2) + I(Mass_2nd^2) + I(Mass_2nd*Mass_Repro) +
Time_2nd), data = data.mnew$redata)
phi < crossprod(t(modmat.mnew), beta)
mus3 < transformSaturated(phi, data.mnew, from = "phi",
to = "mu")
ind < which(grepl("B",rownames(data.frame(modmat.mnew))))
qux3 < matrix(mus3[ind], nrow = 10201, ncol = 8)
colnames(qux3) < paste("c", 1:8, sep = "")
rownames(qux3) < c(1:10201)
sums3 < as.numeric(apply(qux3, 1, sum))
zday3 < matrix(sums3, nrow = 101)
mnew < data.frame(long.sex, rep(days[3], 101^2),
rep(x, times = 101), rep(y, each = 101))
mnew < cbind(mnew, matrix(rep(0, 26*101^2), ncol = 26))
mnames < c("Sex", "Time_2nd", "Mass_2nd", "Mass_Repro", vars)
names(mnew) < mnames
data.mnew < asterdata(mnew, vars = vars, pred = pred,
group = group, code = code, families = families)
foo < as.character(data.mnew$redata$varb)
offspring.mnew < as.numeric(grepl("B", foo))
data.mnew$redata < transform(data.mnew$redata,
offspring.mnew = offspring.mnew)
modmat.mnew < model.matrix(resp ~ varb +
offspring.mnew:(Mass_Repro + Mass_2nd +
I(Mass_Repro^2) + I(Mass_2nd^2) + I(Mass_2nd*Mass_Repro) +
Time_2nd), data = data.mnew$redata)
phi < crossprod(t(modmat.mnew), beta)
mus4 < transformSaturated(phi, data.mnew, from = "phi",
to = "mu")
ind < which(grepl("B",rownames(data.frame(modmat.mnew))))
qux4 < matrix(mus4[ind], nrow = 10201, ncol = 8)
colnames(qux4) < paste("c", 1:8, sep = "")
rownames(qux4) < c(1:10201)
sums4 < as.numeric(apply(qux4, 1, sum))
zday4 < matrix(sums4, nrow = 101)
mnew < data.frame(long.sex, rep(days[4], 101^2),
rep(x, times = 101), rep(y, each = 101))
mnew < cbind(mnew, matrix(rep(0, 26*101^2), ncol = 26))
mnames < c("Sex", "Time_2nd", "Mass_2nd", "Mass_Repro", vars)
names(mnew) < mnames
data.mnew < asterdata(mnew, vars = vars, pred = pred,
group = group, code = code, families = families)
foo < as.character(data.mnew$redata$varb)
offspring.mnew < as.numeric(grepl("B", foo))
data.mnew$redata < transform(data.mnew$redata,
offspring.mnew = offspring.mnew)
modmat.mnew < model.matrix(resp ~ varb +
offspring.mnew:(Mass_Repro + Mass_2nd +
I(Mass_Repro^2) + I(Mass_2nd^2) + I(Mass_2nd*Mass_Repro) +
Time_2nd), data = data.mnew$redata)
phi < crossprod(t(modmat.mnew), beta)
mus5 < transformSaturated(phi, data.mnew, from = "phi",
to = "mu")
ind < which(grepl("B",rownames(data.frame(modmat.mnew))))
qux5 < matrix(mus5[ind], nrow = 10201, ncol = 8)
colnames(qux5) < paste("c", 1:8, sep = "")
rownames(qux5) < c(1:10201)
sums5 < as.numeric(apply(qux5, 1, sum))
zday5 < matrix(sums5, nrow = 101)
@
%\begin{figure}[h!]
%\begin{center}
%<<>>=
%plot(mass.second, mass.reprod, pch = 20,
% xlab = "Mass at 2nd instar", ylab = "Mass at Reproduction",
% main = "Fitness Landscape for ovariole counts vs. Mass")
%points(maximum(sums3)[1], maximum(sums3)[2], col = "red", pch = 19)
%contour(x,y,zday3, add = TRUE, nlevels = 8, levels =
% c(50,100,150,200,250,300,350,400) )
%@
%\end{center}
%\end{figure}
%\begin{figure}[h!]
%\begin{center}
%<<>>=
%plot(mass.second, mass.reprod, pch = 20,
% xlab = "Mass at 2nd instar", ylab = "Mass at Reproduction",
% main = "Fitness Landscape for ovariole counts vs. Mass")
%points(maximum(sums4)[1], maximum(sums4)[2], col = "red", pch = 19)
%contour(x,y,zday4, add = TRUE, nlevels = 8, levels =
% c(50,100,150,200,250,300,350,400) )
%@
%\end{center}
%\end{figure}
%<<>>=
%plot(mass.second, mass.reprod, pch = 20,
% xlab = "Mass at 2nd instar", ylab = "Mass at Reproduction",
% main = "Fitness Landscape for ovariole counts vs. Mass")
%points(maximum(sums5)[1], maximum(sums5)[2], col = "red", pch = 19)
%contour(x,y,zday5, add = TRUE, nlevels = 8, levels =
% c(50,100,150,200,250,300,350,400) )
%@
<>=
ind.max < which(zday2 == max(zday2))
col.max < ceiling(ind.max / 101)
mass2nd.max3 < zday3[, col.max]
mass2nd.max2 < zday2[, col.max]
mass2nd.max4 < zday4[, col.max]
mass2nd.max5 < zday5[, col.max]
@
The code below builds the fitness landscapes. These fitness landscapes appear in the left hand column of Figure 2 in the paper, which is called Figure~\ref{fig:paper} in this technical report. In the paper, the fitness landscapes appear smaller than those on display in Figure~\ref{fig:two}.
<>=
par(mfrow = c(2, 2), mar = rep(0.5, 4),
oma = c(4, 4, 0, 0) + 0.1)
# day 2
plot(mass.second,
mass.reprod,
axes = FALSE, xlab = "", ylab = "", pch = 20)
box()
axis(side = 2, outer = TRUE)
mtext("mass at eclosion", outer = TRUE, line = 3,
side = 2, at = 0.50)
mtext("mass at 2nd instar stage", outer = TRUE, line = 3,
side = 1, at = 0.50)
points(maximum(sums2)[1], maximum(sums2)[2], pch = 0)
levels < seq(50, 350, by = 50)
contour(x, y, zday2, add = TRUE, levels = levels)
# day 3
plot(mass.second,
mass.reprod,
axes = FALSE, xlab = "", ylab = "", pch = 20)
box()
points(maximum(sums3)[1], maximum(sums3)[2], pch = 0)
contour(x, y, zday3, add = TRUE, levels = levels)
# day 4
plot(mass.second,
mass.reprod,
axes = FALSE, xlab = "", ylab = "", pch = 20)
box()
axis(side = 1, outer = TRUE)
axis(side = 2, outer = TRUE)
points(maximum(sums4)[1], maximum(sums4)[2], pch = 0)
contour(x, y, zday4, add = TRUE, levels = levels)
# day 5
plot(mass.second,
mass.reprod,
axes = FALSE, xlab = "", ylab = "", pch = 20)
box()
axis(side = 1, outer = TRUE)
points(maximum(sums5)[1], maximum(sums5)[2], pch = 0)
contour(x, y, zday5, add = TRUE, levels = levels)
@
\begin{figure}
\begin{center}
%%%%% DON'T REPEAT YOURSELF (the DRY rule)
%%%%% Quote chunks, don't repeat chunks.
<>=
<>
@
\end{center}
\caption{Fitness landscapes for expected unconditional ovariole counts vs.\ mass
at eclosion and mass at 2nd instar stage. Different panels are for different
ages (in days since hatching) at which individuals reached the second instar
larval stage (age 2 is top left, age 3 is top right, age 4 is bottom left,
and age 5 is bottom right). Mass at second instar stage is on the $x$axis,
and mass at eclosion is on the $y$axis. The boxes denote the maxima.
%%%Plots are for individuals reaching the second instar stage at
%%day 2 (top left), day 3 (top right), day 4 (bottom left), and day 5 (bottom
%%right). Squares denote locations of maxima.
The maximum values are
\Sexpr{round(max(zday2),1)} (top left),
\Sexpr{round(max(zday3),1)} (top right),
\Sexpr{round(max(zday4),1)} (bottom left),
and \Sexpr{round(max(zday5),1)} (bottomright). }
\label{fig:two}
\end{figure}
Figure~\ref{fig:two} displays the entire fitness landscape for all hypothetical individuals. The fitness landscapes generated show that expected unconditional ovariole count is predominantly explained by an individual's mass at eclosion.
%%%fitness holding age at eclosion constant is predominantly explained by an individual's mass at eclosion.
As mass at eclosion increases, expected fitness increases with a maximum fitness found for female \emph{M. sexta} weighing roughly 2 grams at eclosion. Expected fitness is also influenced by age at which individuals reach the second instar larval stage. The contours show that individuals reaching the second instar larval stage earlier have higher expected numbers of offspring. The contours suggest that fitness depends only weakly on mass at the second instar stage. However, the formal hypothesis test shows that this relationship is highly significant.
In addition, expected unconditional ovariole count declines with increasing age to 2nd instar (Fig.~\ref{fig:two}). For example, maximum expected ovariole count declines by $16 \%$ as age to 2nd instar increases from 2 (upper left panel) to 5 (lower right panel) days. This effect is largely due to the effects of age at 2nd instar on survival to eclosion: slower development (later age at 2nd instar) is associated with lower survival.
\section{Population growth rate} \label{sec:grow}
The preceding analysis accounts for survival and ovariole count as components of fitness, but does not take into account the role of variation in timing of reproduction in fitness variation. To incorporate this effect we must consider the expected population growth rate. The population growth rate parameter ($\lambda$) for the observed population of \emph{M. sexta} is estimated from the stable age equation discussed in \citep*[p.~26]{fisher}, as the basis of accounting for individuals' age at reproduction in its lifetime fitness. In our context, this is
\begin{equation} \label{eq:stable}
1 = \frac{1}{n} \sum_{i = 1}^n \sum_{x = 33}^{40} \mu_{i x} e^{\lambda x}
\end{equation}
where $\mu_{i x}$ is the unconditional expected ovariole count for individual $i$ at day $x$ which is given below.
<<>>=
mu.star < tau.alt[grepl("varbB", rownames(tau.alt))]
@
Having $e^{\lambda x}$ instead of $\rho^x$ in \eqref{eq:stable} follows \citet{charles}. The population growth rate parameter $\lambda$ is calculated here using the \texttt{uniroot} function.
<<>>=
nind < length(unique(hornworm$redata$id))
nind
n < nind
gr < function(lam){
n  sum(mu.star * exp((lam * (32 + seq(along = mu.star)))))
}
lout < uniroot(gr, lower = 105, upper = 105)
lambda.hat < lout$root
lambda.hat
@
In most treatments of the stable age equation,
starting with Fisher, the term $\mu_{i x}$ in \eqref{eq:stable}
is written as the product of probability
of survival to age $x$ and the conditional expectation of number of offspring
at age $x$ given survival to age $x$, but we do not do that because
$\mu_{i x}$ is calculated directly by the aster software.
Most treatments of the stable age equation,
starting with Fisher, do not average over all individuals in the data,
the $(1 / n) \sum_{i = 1}^n$ in \eqref{eq:stable}. That is because those
treatments do not allow for variation among individuals. Consequently, they
use the same model for all individuals and apply the stable age
equation to one individual (and hence to all because all are the same according
to the model). Here, where $\mu_{i x}$ is different for different individuals
because of the covariates in the model (mass at second instar,
mass at eclosion, and age at second instar), we replace the means
for a typical individual with the average over all individuals in the data.
% this is not used anywhere
%<<>>=
%mean.mu.star < mean(mu.star * exp((lambda.hat *
% (32 + seq(along = mu.star)))))
%blah < mu.star * exp((lambda.hat * (32 +
% seq(along = mu.star)))) / mean.mu.star
%@
From the $\hat{\mu}_{i x}$ produced by the aster model and \eqref{eq:stable} we obtain the estimate $\hat{\lambda} = \Sexpr{round(lambda.hat, 3)}$. The positive value of $\lambda$ indicates a growing population. This very large value of $\hat{\lambda}$, indicating that the population grows by a factor of $\exp(\hat{\lambda}) = \Sexpr{round(exp(lambda.hat), 3)}$ per day, results from the fact that many sources of mortality in natural populations were excluded from the experiment (e.~g., low densities of larvae, the netting to exclude predation by birds, and lack of predation before larvae were moved from the lab to the field garden). Such overestimates of population growth rate are typical of experiments that do not have all sources of natural mortality and all sources of failure to reproduce.
We examine the effects that $\lambda$ has on expected fitness by reweighting the fitness landscape according to
\begin{equation} \label{eq:stablereweight}
w(z) = \sum_{x = 33}^{40} \mu_x(z) e^{\lambda x},
\end{equation}
where $\mu_x(z)$ is now expected reproduction at age $x$ for a hypothetical individual having trait values given by $z$ \citep[p. 134]{charles}. The weights according to the population growth rate used are
<<>>=
weight < function(t) exp(lambda.hat * (32 + t))
weight(c(1:8))
@
%\section{Under construction from here on out}
Fitness is now reweighted to incorporate the population growth rate.
% After the population growth rate $\lambda$ has been determined
We refit the model using weights
\begin{equation} \label{eq:lambdaweights}
w_j = f_j e^{ \lambda t_j}
\end{equation}
where $f_j$ are the zero or one weights indicating nodes that contribute directly to fitness and $t_j$ is the age of the individual at node $j$. This weighting accounts for population growth rate \citep[p.~134]{charles}. We incorporate the weights using the code below.
% changes (re weight fitness)
<<>>=
foo < data$redata$varb
offspring < as.numeric(grepl("B", foo))
bar < grepl("B", foo)
baz < offspring
baz[bar] < sub("B", "", foo)[bar]
baz < as.numeric(baz)
baz[bar] < weight(baz[bar]  32)
offspring < baz
offspring < offspring * 100
@
A new model matrix is created that incorporates reweighted fitness.
<<>>=
modmat.full < model.matrix(resp ~ varb +
offspring:(Mass_Repro + Mass_2nd +
I(Mass_Repro^2) + I(Mass_2nd^2) + I(Mass_2nd*Mass_Repro) +
Time_2nd), data = data$redata)
@
We now obtain the new aster submodel meanvalue parameter $\tau$ associated with reweighted fitness. The \texttt{transformUnconditional} function is used to find the corresponding submodel canonical parameter vector $\beta$.
<>=
tau.full < crossprod(modmat.full, data$redata$resp)
beta.full < transformUnconditional(tau.full,
modmat.full, data, from = "tau", to = "beta")
@
The same hypothetical individuals used to build the fitness landscape before adjusting for the population growth rate are used to build the fitness landscape after we adjust for the population growth rate.
% create the same hypothetical individuals to
% generate a fitness landscape for individuals reaching
% their second instar life stage at age 2
<<>>=
x < seq(from = 0, to = 0.016, by = 0.016/100)
y < seq(from = 0, to = 2.3, by = 2.3/100)
days < 2:6 / 7
mnew < data.frame(long.sex, rep(days[1], 101^2),
rep(x, times = 101), rep(y, each = 101))
mnew < cbind(mnew, matrix(rep(0, 26*101^2), ncol = 26))
mnames < c("Sex", "Time_2nd", "Mass_2nd", "Mass_Repro", vars)
names(mnew) < mnames
@
This model matrix is converted into an \texttt{asterdata} object in order to find estimates of fitness using the \texttt{aster2} package.
<>=
data.mnew < asterdata(mnew, vars = vars, pred = pred,
group = group, code = code, families = families)
@
The population growth rate is now incorporated into the model matrix for the hypothetical individuals.
<<>>=
foo < as.character(data.mnew$redata$varb)
bar < grepl("B", foo)
offspring.mnew < as.numeric(bar)
baz < offspring.mnew
baz[bar] < sub("B", "", foo)[bar]
baz < as.numeric(baz)
baz[bar] < weight(baz[bar]  32)
offspring.mnew < baz
offspring.mnew < 100 * offspring.mnew
data.mnew$redata < transform(data.mnew$redata,
offspring.mnew = offspring.mnew)
@
$\mu_x(z)$ is calculated in two steps. We first transform from $\beta$ to $\varphi$ manually and then use the \texttt{transformSaturated} function to transform from $\varphi$ to $\mu$ as done previously.
<<>>=
modmat.mnew < model.matrix(resp ~ varb +
offspring.mnew:(Mass_Repro + Mass_2nd +
I(Mass_Repro^2) + I(Mass_2nd^2) + I(Mass_2nd*Mass_Repro) +
Time_2nd), data = data.mnew$redata)
phi < crossprod(t(modmat.mnew), beta.full)
mus.fit2 < transformSaturated(phi, data.mnew, from = "phi",
to = "mu")
@
We now build the fitness landscape after adjusting for the population growth rate. The fitness landscape calculated using the code below is only for individuals that reach their second instar stage at age 2. The figure generated compares the fitness landscape constructed here to the fitness landscape unadjusted for $\lambda$ for individuals that reach their second instar stage at age 2.
% make the fitness landscape for reweighted fitness for
% individuals reaching eclosion on age 2
<<>>=
ind < which(grepl("B",rownames(data.frame(modmat.mnew))))
qux.fit2 < matrix(mus.fit2[ind], nrow = 10201, ncol = 8)
sumslam.fit2 < as.numeric(apply(qux.fit2, 1, sum))
zday.fit2 < matrix(sumslam.fit2, nrow = 101)
@
\begin{figure}
\begin{center}
<>=
plot(mass.second, pch = 20,
mass.reprod, xlab = "Mass at
2nd instar", ylab = "Mass at Reproduction",
main = "Fitness Landscape for Ovariole Counts vs. Mass")
points(maximum(sums2)[1], maximum(sums2)[2], col = "red",
pch = 19)
contour(x,y,zday2, add = TRUE, nlevels = 8, levels =
c(100,150,200,250,300,325,350,360) )
plot(mass.second, pch = 20,
mass.reprod, xlab = "Mass at
2nd instar", ylab = "Mass at Reproduction",
main = "Fitness Landscape for Ovariole Counts vs. Mass")
points(maximum(sumslam.fit2)[1], maximum(sumslam.fit2)[2],
col = "red", pch = 19)
contour(x,y,zday.fit2, add = TRUE, nlevels = 8)
@
\end{center}
\caption{Fitness landscapes without (left panel) and with (right panel) adjustment for population
growth rate $\lambda$ for hypothetical individuals reaching their second instar stage at age 2.}
\label{fig:comp}
\end{figure}
We can see that the contours of the two fitness landscapes in Figure~\ref{fig:comp} differ. We now change the fitness landscape adjusted for the population growth rate to a relative fitness landscape. This is done by dividing estimated expected ovariole counts by the mean of estimated expected ovariole counts.
%mymu.full < transformUnconditional(beta.full,
% modmat.full, data, from = "beta", to = "mu")
%names(mymu.full) < rownames(modmat.full)
%mymu.full.births < mymu.full[grepl("B", names(mymu.full))]
%mydays < sub("^.*\\.B", "", names(mymu.full.births))
%mydays < as.numeric(mydays[grepl("B", names(mymu.full.births))])
%myweights.lam < exp( lambda.hat * mydays)
%mean.fitness < mean(mymu.full.births * myweights.lam)
%mean.fitness < mean.fitness * 8
<>=
mymu.full < transformUnconditional(beta.full,
modmat.full, data, from = "beta", to = "mu")
names(mymu.full) < rownames(modmat.full)
mymu.full.births < mymu.full[grepl("B", names(mymu.full))]
mean.fitness < mean(mymu.full.births)
mean.fitness < mean.fitness * 8
zday.fit2 < zday.fit2 / mean.fitness
@
The same routine is performed for the 30603 individuals that reach their second instar life stage at ages 3, 4, and 5.
<>=
mnew < data.frame(long.sex, rep(days[2], 101^2),
rep(x, times = 101), rep(y, each = 101))
mnew < cbind(mnew, matrix(rep(0, 26*101^2), ncol = 26))
mnames < c("Sex", "Time_2nd", "Mass_2nd", "Mass_Repro", vars)
names(mnew) < mnames
data.mnew < asterdata(mnew, vars = vars, pred = pred,
group = group, code = code, families = families)
foo < as.character(data.mnew$redata$varb)
offspring.mnew < as.numeric(bar)
baz < offspring.mnew
baz[bar] < sub("B", "", foo)[bar]
baz < as.numeric(baz)
baz[bar] < weight(baz[bar]  32)
offspring.mnew < baz
offspring.mnew < 100 * offspring.mnew
data.mnew$redata < transform(data.mnew$redata,
offspring.mnew = offspring.mnew)
modmat.mnew < model.matrix(resp ~ varb +
offspring.mnew:(Mass_Repro + Mass_2nd +
I(Mass_Repro^2) + I(Mass_2nd^2) + I(Mass_2nd*Mass_Repro) +
Time_2nd), data = data.mnew$redata)
phi < crossprod(t(modmat.mnew), beta.full)
mus.fit3 < transformSaturated(phi, data.mnew, from = "phi",
to = "mu")
births.fit3 < mus.fit3[ind]
quxlam.fit3 < matrix(births.fit3, nrow = 10201)
sumslam.fit3 < apply(quxlam.fit3, 1, sum)
zday.fit3 < matrix(sumslam.fit3, nrow = 101)
zday.fit3 < zday.fit3 / mean.fitness
mnew < data.frame(long.sex, rep(days[3], 101^2),
rep(x, times = 101), rep(y, each = 101))
mnew < cbind(mnew, matrix(rep(0, 26*101^2), ncol = 26))
mnames < c("Sex", "Time_2nd", "Mass_2nd", "Mass_Repro", vars)
names(mnew) < mnames
data.mnew < asterdata(mnew, vars = vars, pred = pred,
group = group, code = code, families = families)
foo < as.character(data.mnew$redata$varb)
offspring.mnew < as.numeric(bar)
baz < offspring.mnew
baz[bar] < sub("B", "", foo)[bar]
baz < as.numeric(baz)
baz[bar] < weight(baz[bar]  32)
offspring.mnew < baz
offspring.mnew < 100 * offspring.mnew
data.mnew$redata < transform(data.mnew$redata,
offspring.mnew = offspring.mnew)
modmat.mnew < model.matrix(resp ~ varb +
offspring.mnew:(Mass_Repro + Mass_2nd +
I(Mass_Repro^2) + I(Mass_2nd^2) + I(Mass_2nd*Mass_Repro) +
Time_2nd), data = data.mnew$redata)
phi < crossprod(t(modmat.mnew), beta.full)
mus.fit4 < transformSaturated(phi, data.mnew, from = "phi",
to = "mu")
births.fit4 < mus.fit4[ind]
quxlam.fit4 < matrix(births.fit4, nrow = 10201)
sumslam.fit4 < apply(quxlam.fit4, 1, sum)
zday.fit4 < matrix(sumslam.fit4, nrow = 101)
zday.fit4 < zday.fit4 / mean.fitness
mnew < data.frame(long.sex, rep(days[4], 101^2),
rep(x, times = 101), rep(y, each = 101))
mnew < cbind(mnew, matrix(rep(0, 26*101^2), ncol = 26))
mnames < c("Sex", "Time_2nd", "Mass_2nd", "Mass_Repro", vars)
names(mnew) < mnames
data.mnew < asterdata(mnew, vars = vars, pred = pred,
group = group, code = code, families = families)
foo < as.character(data.mnew$redata$varb)
offspring.mnew < as.numeric(bar)
baz < offspring.mnew
baz[bar] < sub("B", "", foo)[bar]
baz < as.numeric(baz)
baz[bar] < weight(baz[bar]  32)
offspring.mnew < baz
offspring.mnew < 100 * offspring.mnew
data.mnew$redata < transform(data.mnew$redata,
offspring.mnew = offspring.mnew)
modmat.mnew < model.matrix(resp ~ varb +
offspring.mnew:(Mass_Repro + Mass_2nd +
I(Mass_Repro^2) + I(Mass_2nd^2) + I(Mass_2nd*Mass_Repro) +
Time_2nd), data = data.mnew$redata)
phi < crossprod(t(modmat.mnew), beta.full)
mus.fit5 < transformSaturated(phi, data.mnew, from = "phi",
to = "mu")
births.fit5 < mus.fit5[ind]
quxlam.fit5 < matrix(births.fit5, nrow = 10201)
sumslam.fit5 < apply(quxlam.fit5, 1, sum)
zday.fit5 < matrix(sumslam.fit5, nrow = 101)
zday.fit5 < zday.fit5 / mean.fitness
@
%We now make the fourpanelled fitness landscape where fitness is reweighted to include the population growth rate.
%<>=
%par(mfrow = c(2, 2), mar = rep(0.5, 4), oma = c(4, 4, 0, 0) + 0.1)
%# day 2
%plot(mass.second,
% mass.reprod,
% axes = FALSE, xlab = "", ylab = "", pch = 20)
%box()
%axis(side = 2, outer = TRUE)
%mtext("mass at eclosion", outer = TRUE, line = 3, side = 2, at = 0.50)
%mtext("mass at 2nd instar stage", outer = TRUE, line = 3, side = 1, at = 0.50)
%points(maximum(sumslam.fit2)[1], maximum(sumslam.fit2)[2], pch = 0)
%contour(x, y, zday.fit2, nlevels = 7, add = TRUE)
%# day 3
%plot(mass.second,
% mass.reprod,
% axes = FALSE, xlab = "", ylab = "", pch = 20)
%box()
%points(maximum(sumslam.fit3)[1], maximum(sumslam.fit3)[2], pch = 0)
%contour(x, y, zday.fit3, nlevels = 7, add = TRUE)
%# day 4
%plot(mass.second,
% mass.reprod,
% axes = FALSE, xlab = "", ylab = "", pch = 20)
%box()
%axis(side = 2, outer = TRUE)
%axis(side = 1, outer = TRUE)
%points(maximum(sumslam.fit4)[1], maximum(sumslam.fit4)[2], pch = 0)
%contour(x, y, zday.fit4, nlevels = 7, add = TRUE)
%# day 5
%plot(mass.second,
% mass.reprod,
% axes = FALSE, xlab = "", ylab = "", pch = 20)
%box()
%axis(side = 1, outer = TRUE)
%points(maximum(sumslam.fit5)[1], maximum(sumslam.fit5)[2], pch = 0)
%contour(x, y, zday.fit5, nlevels = 7, add = TRUE)
%@
%\begin{figure}
%\begin{center}
%<>=
%par(mfrow = c(2, 2), mar = rep(0.5, 4), oma = c(4, 4, 0, 0) + 0.1)
%# day 2
%plot(mass.second,
% mass.reprod,
% axes = FALSE, xlab = "", ylab = "", pch = 20)
%box()
%axis(side = 2, outer = TRUE)
%mtext("mass at eclosion", outer = TRUE, line = 3, side = 2, at = 0.50)
%mtext("mass at 2nd instar stage", outer = TRUE, line = 3, side = 1, at = 0.50)
%points(maximum(sumslam.fit2)[1], maximum(sumslam.fit2)[2], pch = 0)
%contour(x, y, zday.fit2, nlevels = 7, add = TRUE)
%# day 3
%plot(mass.second,
% mass.reprod,
% axes = FALSE, xlab = "", ylab = "", pch = 20)
%box()
%points(maximum(sumslam.fit3)[1], maximum(sumslam.fit3)[2], pch = 0)
%contour(x, y, zday.fit3, nlevels = 7, add = TRUE)
%# day 4
%plot(mass.second,
% mass.reprod,
% axes = FALSE, xlab = "", ylab = "", pch = 20)
%box()
%axis(side = 2, outer = TRUE)
%axis(side = 1, outer = TRUE)
%points(maximum(sumslam.fit4)[1], maximum(sumslam.fit4)[2], pch = 0)
%contour(x, y, zday.fit4, nlevels = 7, add = TRUE)
%# day 5
%plot(mass.second,
% mass.reprod,
% axes = FALSE, xlab = "", ylab = "", pch = 20)
%box()
%axis(side = 1, outer = TRUE)
%points(maximum(sumslam.fit5)[1], maximum(sumslam.fit5)[2], pch = 0)
%contour(x, y, zday.fit5, nlevels = 7, add = TRUE)
%@
%\end{center}
%\caption{foo}
%\end{figure}
%%%%%%%%%%%% the 2 two by two plots, before the most recent change %%%%%%%%%%%%
%We now make the two fourpanelled fitness landscape that appears in the paper. These plots allow the reader to discern the differences between estimated expected Darwinian fitness and estimated expected Darwinian fitness adjusted for the population growth rate.
%<>=
%par(mfrow = c(2, 2), mar = rep(0.5, 4),
% oma = c(4, 4, 0, 0) + 0.1)
%# day 2
%plot(mass.second,
% mass.reprod,
% axes = FALSE, xlab = "", ylab = "", pch = 20)
%box()
%axis(side = 2, outer = TRUE)
%mtext("mass at eclosion", outer = TRUE, line = 3,
% side = 2, at = 0.50)
%mtext("mass at 2nd instar stage", outer = TRUE, line = 3,
% side = 1, at = 0.50)
%points(maximum(sums2)[1], maximum(sums2)[2], pch = 0)
%levels < seq(50, 350, by = 50)
%contour(x, y, zday2, add = TRUE, levels = levels)
%plot(mass.second,
% mass.reprod,
% axes = FALSE, xlab = "", ylab = "", pch = 20)
%box()
%points(maximum(sumslam.fit2)[1], maximum(sumslam.fit2)[2],
% pch = 0)
%contour(x, y, zday.fit2, nlevels = 7, add = TRUE)
%# day 3
%plot(mass.second,
% mass.reprod,
% axes = FALSE, xlab = "", ylab = "", pch = 20)
%box()
%axis(side = 1, outer = TRUE)
%axis(side = 2, outer = TRUE)
%points(maximum(sums3)[1], maximum(sums3)[2], pch = 0)
%contour(x, y, zday3, add = TRUE, levels = levels)
%plot(mass.second,
% mass.reprod,
% axes = FALSE, xlab = "", ylab = "", pch = 20)
%box()
%axis(side = 1, outer = TRUE)
%points(maximum(sumslam.fit3)[1], maximum(sumslam.fit3)[2],
% pch = 0)
%contour(x, y, zday.fit3, nlevels = 7, add = TRUE)
%@
%\begin{figure}
%\begin{center}
%<<4by2part1,echo=FALSE>>=
%par(mfrow = c(2, 2), mar = rep(0.5, 4),
% oma = c(4, 4, 0, 0) + 0.1)
%# day 2
%plot(mass.second,
% mass.reprod,
% axes = FALSE, xlab = "", ylab = "", pch = 20)
%box()
%axis(side = 2, outer = TRUE)
%mtext("mass at eclosion", outer = TRUE, line = 3, side = 2,
% at = 0.50)
%mtext("mass at 2nd instar stage", outer = TRUE, line = 3,
% side = 1, at = 0.50)
%points(maximum(sums2)[1], maximum(sums2)[2], pch = 0)
%levels < seq(50, 350, by = 50)
%contour(x, y, zday2, add = TRUE, levels = levels)
%plot(mass.second,
% mass.reprod,
% axes = FALSE, xlab = "", ylab = "", pch = 20)
%box()
%points(maximum(sumslam.fit2)[1], maximum(sumslam.fit2)[2],
% pch = 0)
%contour(x, y, zday.fit2, nlevels = 7, add = TRUE)
%# day 3
%plot(mass.second,
% mass.reprod,
% axes = FALSE, xlab = "", ylab = "", pch = 20)
%box()
%axis(side = 1, outer = TRUE)
%axis(side = 2, outer = TRUE)
%points(maximum(sums3)[1], maximum(sums3)[2], pch = 0)
%contour(x, y, zday3, add = TRUE, levels = levels)
%plot(mass.second,
% mass.reprod,
% axes = FALSE, xlab = "", ylab = "", pch = 20)
%box()
%axis(side = 1, outer = TRUE)
%points(maximum(sumslam.fit3)[1], maximum(sumslam.fit3)[2],
% pch = 0)
%contour(x, y, zday.fit3, nlevels = 7, add = TRUE)
%@
%\end{center}
%\caption{Fitness landscapes with and without adjustment for population
%growth rate $\lambda$. Left column not adjusted, right column adjusted.
%Top row 2nd instar stage reached at age two,
%bottom row 2nd instar stage reached at age three
%(Fig.~\ref{fig:bar} has other ages). Numbers on contours are absolute
%fitness (unconditional expected ovariole counts) in the left column
%and are relative fitness (absolute fitness divided its average
%over the population) in the right column.
%All plots plot fitness vs.\ mass at eclosion and mass at 2nd instar stage.
%Boxes denote location of maxima.
%The maximum values are
%\Sexpr{round(max(zday2),1)} (top left),
%\Sexpr{round(max(zday3),1)} (bottom left),
%\Sexpr{round(max(zday.fit2),2)} (top right),
%and \Sexpr{round(max(zday.fit3),2)} (bottomright).
%}
%\label{fig:foo}
%\end{figure}
% These 8 plots are fitness landscapes. Fitness landscapes appearing in the first column are for unconditional ovariole count while fitness landscapes appearing in the second column are for fitness adjusted for the population growth rate. The first row is for individuals that reach their second instar stage at age 2, the second row is for individuals that reach their second instar stage at 3, and so on.
%%%%%%%%%%%% the 2 two by two plots, before the most recent change (part 2) %%%%%%%%%%%%
%<>=
%# day 4
%plot(mass.second,
% mass.reprod,
% axes = FALSE, xlab = "", ylab = "", pch = 20)
%box()
%mtext("mass at eclosion", outer = TRUE, line = 3, side = 2,
% at = 0.50)
%mtext("mass at 2nd instar stage", outer = TRUE, line = 3,
% side = 1, at = 0.50)
%axis(side = 2, outer = TRUE)
%points(maximum(sums4)[1], maximum(sums4)[2], pch = 0)
%contour(x, y, zday4, add = TRUE, levels = levels)
%plot(mass.second,
% mass.reprod,
% axes = FALSE, xlab = "", ylab = "", pch = 20)
%box()
%points(maximum(sumslam.fit4)[1], maximum(sumslam.fit4)[2],
% pch = 0)
%contour(x, y, zday.fit4, nlevels = 7, add = TRUE)
%# day 5
%plot(mass.second,
% mass.reprod,
% axes = FALSE, xlab = "", ylab = "", pch = 20)
%box()
%axis(side = 1, outer = TRUE)
%axis(side = 2, outer = TRUE)
%points(maximum(sums5)[1], maximum(sums5)[2], pch = 0)
%contour(x, y, zday5, add = TRUE, levels = levels)
%plot(mass.second,
% mass.reprod,
% axes = FALSE, xlab = "", ylab = "", pch = 20)
%box()
%axis(side = 1, outer = TRUE)
%points(maximum(sumslam.fit5)[1], maximum(sumslam.fit5)[2],
% pch = 0)
%contour(x, y, zday.fit5, nlevels = 7, add = TRUE)
%@
%\begin{figure}
%\begin{center}
%<<4by2part2,echo=FALSE>>=
%par(mfrow = c(2, 2), mar = rep(0.5, 4), oma =
% c(4, 4, 0, 0) + 0.1)
%# day 4
%plot(mass.second,
% mass.reprod,
% axes = FALSE, xlab = "", ylab = "", pch = 20)
%box()
%mtext("mass at eclosion", outer = TRUE, line = 3, side = 2,
% at = 0.50)
%mtext("mass at 2nd instar stage", outer = TRUE, line = 3,
% side = 1, at = 0.50)
%axis(side = 2, outer = TRUE)
%points(maximum(sums4)[1], maximum(sums4)[2], pch = 0)
%contour(x, y, zday4, add = TRUE, levels = levels)
%plot(mass.second,
% mass.reprod,
% axes = FALSE, xlab = "", ylab = "", pch = 20)
%box()
%points(maximum(sumslam.fit4)[1], maximum(sumslam.fit4)[2],
% pch = 0)
%contour(x, y, zday.fit4, nlevels = 7, add = TRUE)
%# day 5
%plot(mass.second,
% mass.reprod,
% axes = FALSE, xlab = "", ylab = "", pch = 20)
%box()
%axis(side = 1, outer = TRUE)
%axis(side = 2, outer = TRUE)
%points(maximum(sums5)[1], maximum(sums5)[2], pch = 0)
%contour(x, y, zday5, add = TRUE, levels = levels)
%plot(mass.second,
% mass.reprod,
% axes = FALSE, xlab = "", ylab = "", pch = 20)
%box()
%axis(side = 1, outer = TRUE)
%points(maximum(sumslam.fit5)[1], maximum(sumslam.fit5)[2],
% pch = 0)
%contour(x, y, zday.fit5, nlevels = 7, add = TRUE)
%@
%\end{center}
%\caption{Fitness landscapes with and without adjustment for population
%growth rate $\lambda$. Left column not adjusted, right column adjusted.
%Top row 2nd instar stage reached at age four,
%bottom row 2nd instar stage reached at age five
%(Fig.~\ref{fig:foo} has other ages). Numbers on contours are absolute
%fitness (unconditional expected ovariole counts) in the left column
%and are relative fitness (absolute fitness divided its average
%over the population) in the right column.
%All plots plot fitness vs.\ mass at eclosion and mass at 2nd instar stage.
%Boxes denote location of maxima.
%The maximum values are
%\Sexpr{round(max(zday4),1)} (top left),
%\Sexpr{round(max(zday5),1)} (bottom left),
%\Sexpr{round(max(zday.fit4),2)} (top right),
%and \Sexpr{round(max(zday.fit5),2)} (bottomright).
%}
%\label{fig:bar}
%\end{figure}
We now make the eightpaneled plot that appears in the paper, see Figure~\ref{fig:paper}. The plots in these panels are slightly different than those that appeared in Figure~\ref{fig:comp}. We have zoomed out in order to capture more of the contour shapes. This allows one to easily see the differences in the contours before and after adjusting for the population growth rate.
<>=
xlim < range(x)
ylim < range(y)
par(mfrow = c(4, 2), mar = rep(0.5, 4),
oma = c(4, 4, 0, 4) + 0.1)
# day 2
plot(mass.second,
mass.reprod,
axes = FALSE, xlab = "", ylab = "", pch = 20,
xlim = xlim, ylim = ylim)
box()
axis(side = 2, outer = TRUE)
mtext("mass at eclosion", outer = TRUE, line = 3, side = 2,
at = 0.50)
mtext("mass at 2nd instar stage", outer = TRUE, line = 3,
side = 1, at = 0.50)
mtext("age second instar stage reached", outer = TRUE,
line = 1, side = 4, at = 0.50)
mtext(as.character(2:5), outer = TRUE, line = 3,
side = 4, at = c(7, 5, 3, 1) / 8)
points(maximum(sums2)[1], maximum(sums2)[2], pch = 0)
levels < seq(50, 350, by = 50)
contour(x, y, zday2, add = TRUE, levels = levels)
plot(mass.second,
mass.reprod,
axes = FALSE, xlab = "", ylab = "", pch = 20,
xlim = xlim, ylim = ylim)
box()
points(maximum(sumslam.fit2)[1], maximum(sumslam.fit2)[2],
pch = 0)
contour(x, y, zday.fit2, nlevels = 7, add = TRUE)
# day 3
plot(mass.second,
mass.reprod,
axes = FALSE, xlab = "", ylab = "", pch = 20,
xlim = xlim, ylim = ylim)
box()
axis(side = 2, outer = TRUE)
mtext("mass at eclosion", outer = TRUE, line = 3,
side = 2, at = 0.50)
points(maximum(sums3)[1], maximum(sums3)[2], pch = 0)
contour(x, y, zday3, add = TRUE, levels = levels)
plot(mass.second,
mass.reprod,
axes = FALSE, xlab = "", ylab = "", pch = 20,
xlim = xlim, ylim = ylim)
box()
points(maximum(sumslam.fit3)[1], maximum(sumslam.fit3)[2],
pch = 0)
contour(x, y, zday.fit3, nlevels = 7, add = TRUE)
# day 4
plot(mass.second,
mass.reprod,
axes = FALSE, xlab = "", ylab = "", pch = 20,
xlim = xlim, ylim = ylim)
box()
axis(side = 2, outer = TRUE)
mtext("mass at eclosion", outer = TRUE, line = 3,
side = 2, at = 0.50)
points(maximum(sums4)[1], maximum(sums4)[2], pch = 0)
contour(x, y, zday4, add = TRUE, levels = levels)
plot(mass.second,
mass.reprod,
axes = FALSE, xlab = "", ylab = "", pch = 20,
xlim = xlim, ylim = ylim)
box()
points(maximum(sumslam.fit4)[1], maximum(sumslam.fit4)[2],
pch = 0)
contour(x, y, zday.fit4, nlevels = 7, add = TRUE)
# day 5
plot(mass.second,
mass.reprod,
axes = FALSE, xlab = "", ylab = "", pch = 20,
xlim = xlim, ylim = ylim)
box()
axis(side = 1, outer = TRUE)
axis(side = 2, outer = TRUE)
points(maximum(sums5)[1], maximum(sums5)[2], pch = 0)
contour(x, y, zday5, add = TRUE, levels = levels)
plot(mass.second,
mass.reprod,
axes = FALSE, xlab = "", ylab = "", pch = 20,
xlim = xlim, ylim = ylim)
box()
axis(side = 1, outer = TRUE)
points(maximum(sumslam.fit5)[1], maximum(sumslam.fit5)[2],
pch = 0)
contour(x, y, zday.fit5, nlevels = 7, add = TRUE)
@
\begin{figure}[p]
\begin{center}
%%%%% DON'T REPEAT YOURSELF (the DRY rule)
%%%%% Quote chunks, don't repeat chunks.
<>=
<>
@
\end{center}
\caption{Fitness landscapes without (left column) and with (right column) adjustment for population
growth rate $\lambda$.
Rows top to bottom are 2nd instar stage
reached at age 2, 3, 4, and 5.
Numbers on contours are absolute
fitness (unconditional expected ovariole counts) in the left column
and are relative fitness (absolute fitness divided by its average
over the population) in the right column.
All plots display fitness as contours vs. mass at eclosion and mass at 2nd instar stage.
Boxes denote locations of maxima.
Maximum values are (left column, from top to bottom)
\Sexpr{formatC(max(sums2), digits = 1, format = "f")},
\Sexpr{formatC(max(sums3), digits = 1, format = "f")},
\Sexpr{formatC(max(sums4), digits = 1, format = "f")},
\Sexpr{formatC(max(sums5), digits = 1, format = "f")},
(right column, from top to bottom)
\Sexpr{formatC(max(zday.fit2), digits = 2, format = "f")},
\Sexpr{formatC(max(zday.fit3), digits = 2, format = "f")},
\Sexpr{formatC(max(zday.fit4), digits = 2, format = "f")},
\Sexpr{formatC(max(zday.fit5), digits = 2, format = "f")}.}
% The maximum values for unconditional expected ovariole counts for each age are
% 363.6 (age 2),
% 342.4 (age 3),
% 322.4 (age 4),
% and 303.6 (age 5).}
%. Left column not adjusted: contours are fitness (unconditional
%expected ovariole count). Right column adjusted: contours are relative fitness
%(fitness divided by mean fitness).
%Boxes denote locations of maxima, the values
%being 363.6, 342.4, 322.4, and 303.6 going down left column
%and 1.39, 1.28, 1.20, and 1.14 going down right column.}
\label{fig:paper}
\end{figure}
%%%%% All 8 contours with unlabeled axes %%%%%
%<>=
%par(mfrow = c(4, 2), mar = rep(0.5, 4), oma = c(4, 4, 0, 0) + 0.1)
%# day 2
%plot(mass.second,
% mass.reprod,
% axes = FALSE, xlab = "", ylab = "", pch = 20)
%box()
%axis(side = 2, outer = TRUE)
%mtext("mass at eclosion", outer = TRUE, line = 3, side = 2, at = 0.50)
%mtext("mass at 2nd instar stage", outer = TRUE, line = 3, side = 1, at = 0.50)
%points(maximum(sums2)[1], maximum(sums2)[2], pch = 0)
%levels < seq(50, 350, by = 50)
%contour(x, y, zday2, add = TRUE, levels = levels,
% drawlabels = FALSE)
%plot(mass.second,
% mass.reprod,
% axes = FALSE, xlab = "", ylab = "", pch = 20)
%box()
%points(maximum(sumslam.fit2)[1], maximum(sumslam.fit2)[2], pch = 0)
%contour(x, y, zday.fit2, nlevels = 7, add = TRUE,
% drawlabels = FALSE)
%# day 3
%plot(mass.second,
% mass.reprod,
% axes = FALSE, xlab = "", ylab = "", pch = 20)
%box()
%axis(side = 2, outer = TRUE)
%mtext("mass at eclosion", outer = TRUE, line = 3, side = 2, at = 0.50)
%points(maximum(sums3)[1], maximum(sums3)[2], pch = 0)
%contour(x, y, zday3, add = TRUE, levels = levels,
% drawlabels = FALSE)
%plot(mass.second,
% mass.reprod,
% axes = FALSE, xlab = "", ylab = "", pch = 20)
%box()
%points(maximum(sumslam.fit3)[1], maximum(sumslam.fit3)[2], pch = 0)
%#contour(x, y, zday.fit3, nlevels = 7, add = TRUE)
%contour(x, y, zday.fit3, nlevels = 7, add = TRUE,
% drawlabels = FALSE)
%# day 4
%plot(mass.second,
% mass.reprod,
% axes = FALSE, xlab = "", ylab = "", pch = 20)
%box()
%axis(side = 2, outer = TRUE)
%mtext("mass at eclosion", outer = TRUE, line = 3, side = 2, at = 0.50)
%points(maximum(sums4)[1], maximum(sums4)[2], pch = 0)
%contour(x, y, zday4, add = TRUE, levels = levels,
% drawlabels = FALSE)
%plot(mass.second,
% mass.reprod,
% axes = FALSE, xlab = "", ylab = "", pch = 20)
%box()
%points(maximum(sumslam.fit4)[1], maximum(sumslam.fit4)[2], pch = 0)
%contour(x, y, zday.fit4, nlevels = 7, add = TRUE,
% drawlabels = FALSE)
%# day 5
%plot(mass.second,
% mass.reprod,
% axes = FALSE, xlab = "", ylab = "", pch = 20)
%box()
%axis(side = 1, outer = TRUE)
%axis(side = 2, outer = TRUE)
%points(maximum(sums5)[1], maximum(sums5)[2], pch = 0)
%contour(x, y, zday5, add = TRUE, levels = levels,
% drawlabels = FALSE)
%plot(mass.second,
% mass.reprod,
% axes = FALSE, xlab = "", ylab = "", pch = 20)
%box()
%axis(side = 1, outer = TRUE)
%points(maximum(sumslam.fit5)[1], maximum(sumslam.fit5)[2], pch = 0)
%contour(x, y, zday.fit5, nlevels = 7, add = TRUE,
% drawlabels = FALSE)
%@
%\begin{figure}
%\begin{center}
%<<4by2tech, echo = FALSE>>=
%par(mfrow = c(4, 2), mar = rep(0.5, 4), oma = c(4, 4, 0, 0) + 0.1)
%# day 2
%plot(mass.second,
% mass.reprod,
% axes = FALSE, xlab = "", ylab = "", pch = 20)
%box()
%axis(side = 2, outer = TRUE)
%mtext("mass at eclosion", outer = TRUE, line = 3, side = 2, at = 0.50)
%mtext("mass at 2nd instar stage", outer = TRUE, line = 3, side = 1, at = 0.50)
%points(maximum(sums2)[1], maximum(sums2)[2], pch = 0)
%levels < seq(50, 350, by = 50)
%#contour(x, y, zday2, add = TRUE, levels = levels)
%contour(x, y, zday2, add = TRUE, levels = levels,
% drawlabels = FALSE)
%plot(mass.second,
% mass.reprod,
% axes = FALSE, xlab = "", ylab = "", pch = 20)
%box()
%points(maximum(sumslam.fit2)[1], maximum(sumslam.fit2)[2], pch = 0)
%#contour(x, y, zday.fit2, nlevels = 7, add = TRUE)
%contour(x, y, zday.fit2, nlevels = 7, add = TRUE,
% drawlabels = FALSE)
%# day 3
%plot(mass.second,
% mass.reprod,
% axes = FALSE, xlab = "", ylab = "", pch = 20)
%box()
%axis(side = 2, outer = TRUE)
%mtext("mass at eclosion", outer = TRUE, line = 3, side = 2, at = 0.50)
%points(maximum(sums3)[1], maximum(sums3)[2], pch = 0)
%#contour(x, y, zday3, add = TRUE, levels = levels)
%contour(x, y, zday3, add = TRUE, levels = levels,
% drawlabels = FALSE)
%plot(mass.second,
% mass.reprod,
% axes = FALSE, xlab = "", ylab = "", pch = 20)
%box()
%points(maximum(sumslam.fit3)[1], maximum(sumslam.fit3)[2], pch = 0)
%#contour(x, y, zday.fit3, nlevels = 7, add = TRUE)
%contour(x, y, zday.fit3, nlevels = 7, add = TRUE,
% drawlabels = FALSE)
%# day 4
%plot(mass.second,
% mass.reprod,
% axes = FALSE, xlab = "", ylab = "", pch = 20)
%box()
%axis(side = 2, outer = TRUE)
%mtext("mass at eclosion", outer = TRUE, line = 3, side = 2, at = 0.50)
%points(maximum(sums4)[1], maximum(sums4)[2], pch = 0)
%#contour(x, y, zday4, add = TRUE, levels = levels)
%contour(x, y, zday4, add = TRUE, levels = levels,
% drawlabels = FALSE)
%plot(mass.second,
% mass.reprod,
% axes = FALSE, xlab = "", ylab = "", pch = 20)
%box()
%points(maximum(sumslam.fit4)[1], maximum(sumslam.fit4)[2], pch = 0)
%#contour(x, y, zday.fit4, nlevels = 7, add = TRUE)
%contour(x, y, zday.fit4, nlevels = 7, add = TRUE,
% drawlabels = FALSE)
%# day 5
%plot(mass.second,
% mass.reprod,
% axes = FALSE, xlab = "", ylab = "", pch = 20)
%box()
%axis(side = 1, outer = TRUE)
%axis(side = 2, outer = TRUE)
%points(maximum(sums5)[1], maximum(sums5)[2], pch = 0)
%#contour(x, y, zday5, add = TRUE, levels = levels)
%contour(x, y, zday5, add = TRUE, levels = levels,
% drawlabels = FALSE)
%plot(mass.second,
% mass.reprod,
% axes = FALSE, xlab = "", ylab = "", pch = 20)
%box()
%axis(side = 1, outer = TRUE)
%points(maximum(sumslam.fit5)[1], maximum(sumslam.fit5)[2], pch = 0)
%#contour(x, y, zday.fit5, nlevels = 7, add = TRUE)
%contour(x, y, zday.fit5, nlevels = 7, add = TRUE,
% drawlabels = FALSE)
%@
%\end{center}
%\caption{These 8 plots are fitness landscapes. Fitness landscapes appearing in the first column are
%for unconditional ovariole count while fitness landscapes appearing in the second column are for fitness
%adjusted for the population growth rate. The first row is for individuals that reach their second instar
%stage at age 2, the second row is for individuals that reach their second instar stage at 3, and so on. }
%\end{figure}
%$\mu_x(z)$ is calculated using the \texttt{transformUnconditional} function. This calculation could be made by first transforming from $\beta$ to $\varphi$ manually and then using \texttt{transformSaturated} to transform from $\varphi$ to $\mu$ as done previously.
%<<>>=
%mymu.full < transformUnconditional(beta.full,
% modmat.full, data, from = "beta", to = "mu")
%names(mymu.full) < rownames(modmat.full)
%@
%Reweighted estimated expected ovariole count is then computed along with the mean of reweighted estimated expected ovariole count. This mean is taken over the estimated expected ovariole count of females within our observed dataset. There are 8 expected ovariole count values for each individual. These correspond to expected ovariole count for an individual at every possible day at which they may reach their reproductive stage. We therefore multiply our estimated mean expected ovariole count by 8 so that our average has a denominator corresponding to the number of individuals rather than number of nodes.
%<<>>= % first line used to be grepl("varbB", rownames(tau.full))
%mymu.full.births < mymu.full[grepl("B", names(mymu.full))]
%mydays < sub("^.*\\.B", "", names(mymu.full.births))
%mydays < as.numeric(mydays[grepl("B", names(mymu.full.births))])
%myweights.lam < exp( lambda.hat * mydays)
%mean.fitness < mean(mymu.full.births * myweights.lam)
%mean.fitness < mean.fitness * 8
%@
%\subsection{Fitness landscape weighted with respect to the population growth rate}
%The fitness landscape for expected ovariole count reweighted with respect to the population growth rate is constructed the same way as before with the exception that estimates of $\mu$ are weighted.
%<<>>=
%mnew < data.frame(cbind(rep(x,101),
% unlist(lapply(y, FUN = gary))))
%mnew < cbind(id, rep(1,101^2), rep(days[1],101^2),
% mnew, matrix( rep(0,26*101^2), nrow = 101^2, ncol = 26))
%names(mnew) < names(Msextadata)
%data.mnew < asterdata(mnew, vars = vars, pred = pred,
% group = group, code = code, families = families)
%foo < as.character(data.mnew$redata$varb)
%bar < grepl("B", foo)
%offspring.mnew < as.numeric(bar)
%baz < offspring.mnew
%baz[bar] < sub("B", "", foo)[bar]
%baz < as.numeric(baz)
%baz[bar] < weight(baz[bar]  32)
%offspring.mnew < baz
%data.mnew$redata < transform(data.mnew$redata,
% offspring.mnew = offspring.mnew)
%modmat.mnew < model.matrix(resp ~ varb +
% offspring.mnew:(Mass_2nd + Mass_Repro + I(Mass_2nd^2) +
% I(Mass_Repro^2) + I(Mass_2nd*Mass_Repro) + Time_2nd),
% data = data.mnew$redata)
%phi < crossprod(t(modmat.mnew), beta.full)
%mus2 < transformSaturated(phi, data.mnew, from = "phi",
% to = "mu")
%ind < which(grepl("B",rownames(data.frame(modmat.mnew))))
%births2 < mus2[ind]
%#%john2 < function(x) weight(x) * births2[((x1)*10201+1):(x*10201)]
%#%quxlam2 < unlist(apply( matrix(c(1:8)),1, FUN = john2))
%quxlam2 < matrix(births2, nrow = 10201)
%sumslam2 < apply(quxlam2, 1, sum)
%zlam2 < matrix(sumslam2, nrow = 101)
%@
%We now divide by the estimated mean expected ovariole count to produce relative fitness for our hypothetical individuals
%<<>>=
%#%zlam2 < zlam2 / (mean.fitness)
%@
%Our relative fitness landscape is seen below.
%<<>>=
%plot(mass.second,
% mass.reprod, xlab =
% "Mass at 2nd instar", ylab = "Mass at Reproduction", pch = 20,
% main = "Fitness Landscape for ovariole counts vs.Mass")
%contour(x,y,zlam2, add = TRUE)
%@
%We perform the same calculations to produce relative fitness landscapes for the other 30603 hypothetical individuals that reach their second instar life stage at ages 3, 4, and 5.
%<<>>=
%mnew < data.frame(cbind(rep(x,101),
% unlist(lapply(y, FUN = gary))))
%mnew < cbind(id, rep(1,101^2), rep(days[2],101^2),
% mnew, matrix( rep(0,26*101^2), nrow = 101^2, ncol = 26))
%names(mnew) < names(Msextadata)
%data.mnew < asterdata(mnew, vars = vars, pred = pred,
% group = group, code = code, families = families)
%foo < as.character(data.mnew$redata$varb)
%offspring.mnew < as.numeric(bar)
%baz < offspring.mnew
%baz[bar] < sub("B", "", foo)[bar]
%baz < as.numeric(baz)
%baz[bar] < weight(baz[bar]  32)
%offspring.mnew < baz
%data.mnew$redata < transform(data.mnew$redata,
% offspring.mnew = offspring.mnew)
%modmat.mnew < model.matrix(resp ~ varb +
% offspring.mnew:(Mass_2nd + Mass_Repro + I(Mass_2nd^2) +
% I(Mass_Repro^2) + I(Mass_2nd*Mass_Repro) + Time_2nd),
% data = data.mnew$redata)
%phi < crossprod(t(modmat.mnew), beta.full)
%mus3 < transformSaturated(phi, data.mnew, from = "phi",
% to = "mu")
%births3 < mus3[ind]
%#%john3 < function(x) weight(x) * births3[((x1)*10201+1):(x*10201)]
%#%quxlam3 < unlist(apply( matrix(c(1:8)),1, FUN = john3))
%quxlam3 < matrix(births3, nrow = 10201)
%sumslam3 < apply(quxlam3, 1, sum)
%zlam3 < matrix(sumslam3, nrow = 101)
%# compute relative fitness
%#zlam3 < zlam3 / (mean.fitness * 8)
%mnew < data.frame(cbind(rep(x,101),
% unlist(lapply(y, FUN = gary))))
%mnew < cbind(id, rep(1,101^2), rep(days[3],101^2),
% mnew, matrix( rep(0,26*101^2), nrow = 101^2, ncol = 26))
%names(mnew) < names(Msextadata)
%data.mnew < asterdata(mnew, vars = vars, pred = pred,
% group = group, code = code, families = families)
%foo < as.character(data.mnew$redata$varb)
%offspring.mnew < as.numeric(bar)
%baz < offspring.mnew
%baz[bar] < sub("B", "", foo)[bar]
%baz < as.numeric(baz)
%baz[bar] < weight(baz[bar]  32)
%offspring.mnew < baz
%data.mnew$redata < transform(data.mnew$redata,
% offspring.mnew = offspring.mnew)
%modmat.mnew < model.matrix(resp ~ varb +
% offspring.mnew:(Mass_2nd + Mass_Repro + I(Mass_2nd^2) +
% I(Mass_Repro^2) + I(Mass_2nd*Mass_Repro) + Time_2nd),
% data = data.mnew$redata)
%phi < crossprod(t(modmat.mnew), beta.full)
%mus4 < transformSaturated(phi, data.mnew, from = "phi",
% to = "mu")
%births4 < mus4[ind]
%#%john4 < function(x) weight(x) * births4[((x1)*10201+1):(x*10201)]
%#%quxlam4 < unlist(apply( matrix(c(1:8)),1, FUN = john4))
%quxlam4 < matrix(births4, nrow = 10201)
%sumslam4 < apply(quxlam4, 1, sum)
%zlam4 < matrix(sumslam4, nrow = 101)
%# compute relative fitness
%#zlam4 < zlam4 / (mean.fitness * 8)
%mnew < data.frame(cbind(rep(x,101),
% unlist(lapply(y, FUN = gary))))
%mnew < cbind(id, rep(1,101^2), rep(days[4],101^2),
% mnew, matrix( rep(0,26*101^2), nrow = 101^2, ncol = 26))
%names(mnew) < names(Msextadata)
%data.mnew < asterdata(mnew, vars = vars, pred = pred,
% group = group, code = code, families = families)
%foo < as.character(data.mnew$redata$varb)
%offspring.mnew < as.numeric(bar)
%baz < offspring.mnew
%baz[bar] < sub("B", "", foo)[bar]
%baz < as.numeric(baz)
%baz[bar] < weight(baz[bar]  32)
%offspring.mnew < baz
%data.mnew$redata < transform(data.mnew$redata,
% offspring.mnew = offspring.mnew)
%modmat.mnew < model.matrix(resp ~ varb +
% offspring.mnew:(Mass_2nd + Mass_Repro + I(Mass_2nd^2) +
% I(Mass_Repro^2) + I(Mass_2nd*Mass_Repro) + Time_2nd),
% data = data.mnew$redata)
%phi < crossprod(t(modmat.mnew), beta.full)
%mus5 < transformSaturated(phi, data.mnew, from = "phi",
% to = "mu")
%births5 < mus5[ind]
%#%john5 < function(x) weight(x) * births5[((x1)*10201+1):(x*10201)]
%#%quxlam5 < unlist(apply( matrix(c(1:8)),1, FUN = john5))
%quxlam5 < matrix(births5, nrow = 10201)
%sumslam5 < apply(quxlam5, 1, sum)
%zlam5 < matrix(sumslam5, nrow = 101)
%# compute relative fitness
%#zlam5 < zlam5 / (mean.fitness * 8)
%@
%<<>>=
%plot(mass.second, mass.reprod, xlab =
%"Mass at 2nd instar", ylab = "Mass at Reproduction",
%main = "Fitness Landscape for ovariole counts vs.Mass")
%contour(x,y,zlam3, add = TRUE)
%@
%<<>>=
%plot(mass.second, mass.reprod, xlab =
%"Mass at 2nd instar", ylab = "Mass at Reproduction",
%main = "Fitness Landscape for ovariole counts vs.Mass")
%contour(x,y,zlam4, add = TRUE)
%@
%<<>>=
%plot(mass.second, mass.reprod, xlab =
%"Mass at 2nd instar", ylab = "Mass at Reproduction",
%main = "Fitness Landscape for ovariole counts vs.Mass")
%contour(x,y,zlam5, add = TRUE)
%@
%Comparison of \eqref{eq:stable} and \eqref{eq:stablereweight} shows that
%if we replace $z$ for a hypothetical individual by the covariate vector
%$z_i$ for individual $i$ and average over all individuals in the data,
%we get 1. This indicates that $w(z)$ is \emph{relative fitness} (fitness divided by mean
%fitness) for a hypothetical individual with covariate vector $z$.
%Then computing $w(z)$ for $z$ values ranging over a grid
%and making a contour plot of these values gives a fitness landscape
%adjusted for population growth rate (Figs.~\ref{fig:four} and~\ref{fig:five}). \\
%%%find the population growth rate. The population growth rate parameter for the observed population of \emph{M. sexta} is calculated from the stable age equation
%%%\[ 1 = \sum_{x = 33}^{40} \mu(x) e^{\lambda x} \]
%%%as discussed in \citep*[p.~26]{fisher}. The unconditional expectation of fitness for day $x$, denoted by $\mu_x$, is equivalent to Fisher's $l_x b_x$ where in our model $\lambda$ is in days instead of years. From the stable age equation, $\lambda = \Sexpr{round(lambda.hat, 3)}$. The positive value of $\lambda$ indicates a growing population. We examine the effects that $\lambda$ has on expected fitness by reweighting the fitness landscape according to
%%%\[ w_{ijk} = \sum_{x = 33}^{40} \mu_{ijk}(x) e^{\lambda x}. \]
%%%The subscripts $i$, $j$, and $k$ represent the $i^{th}$ mass at the second instar stage, the $j^{th}$ mass at eclosion, and the $k^{th}$ day at which an individual reaches its second instar stage respectively \citep[p.~134]{charles}. $\mu_{ijk}(x)$ is thus the unconditional expected fitness for a female \emph{M. sexta} on day $x$ having the $i^{th}$ value for mass at second instar stage, the $j^{th}$ value for mass at eclosion, and the $k^{th}$ value for day when second instar stage is reached.
%The estimated fitness landscapes taking into account population growth and timing of reproduction (Fig.~\ref{fig:four})
%%%tells essentially the same story as the first fitness landscape.
%is qualitatively similar to that for unconditional ovariole count without accounting for these aspects (Fig.~\ref{fig:two}).
%An increase in mass at eclosion increases expected fitness for each fixed age. As before, expected fitness does not appear to vary much with changes in mass at the second instar stage. Similarly, expected fitness decreases with increasing age at 2nd instar; for example, maximum fitness decreases by $15 \%$ as age to 2nd instar increases from 2 to 5 days (Fig.~\ref{fig:four}). This suggests that, in this study, variation in generation time does not contribute importantly to patterns of phenotypic selection on these traits.
%\begin{figure}[h!]
%\begin{center}
%<