\documentclass[11pt,notitlepage]{article}\usepackage[]{graphicx}\usepackage[]{color}
%% maxwidth is the original width if it is less than linewidth
%% otherwise use linewidth (to make sure the graphics do not exceed the margin)
\makeatletter
\def\maxwidth{ %
  \ifdim\Gin@nat@width>\linewidth
    \linewidth
  \else
    \Gin@nat@width
  \fi
}
\makeatother

\definecolor{fgcolor}{rgb}{0.345, 0.345, 0.345}
\newcommand{\hlnum}[1]{\textcolor[rgb]{0.686,0.059,0.569}{#1}}%
\newcommand{\hlstr}[1]{\textcolor[rgb]{0.192,0.494,0.8}{#1}}%
\newcommand{\hlcom}[1]{\textcolor[rgb]{0.678,0.584,0.686}{\textit{#1}}}%
\newcommand{\hlopt}[1]{\textcolor[rgb]{0,0,0}{#1}}%
\newcommand{\hlstd}[1]{\textcolor[rgb]{0.345,0.345,0.345}{#1}}%
\newcommand{\hlkwa}[1]{\textcolor[rgb]{0.161,0.373,0.58}{\textbf{#1}}}%
\newcommand{\hlkwb}[1]{\textcolor[rgb]{0.69,0.353,0.396}{#1}}%
\newcommand{\hlkwc}[1]{\textcolor[rgb]{0.333,0.667,0.333}{#1}}%
\newcommand{\hlkwd}[1]{\textcolor[rgb]{0.737,0.353,0.396}{\textbf{#1}}}%

\usepackage{framed}
\makeatletter
\newenvironment{kframe}{%
 \def\at@end@of@kframe{}%
 \ifinner\ifhmode%
  \def\at@end@of@kframe{\end{minipage}}%
  \begin{minipage}{\columnwidth}%
 \fi\fi%
 \def\FrameCommand##1{\hskip\@totalleftmargin \hskip-\fboxsep
 \colorbox{shadecolor}{##1}\hskip-\fboxsep
     % There is no \\@totalrightmargin, so:
     \hskip-\linewidth \hskip-\@totalleftmargin \hskip\columnwidth}%
 \MakeFramed {\advance\hsize-\width
   \@totalleftmargin\z@ \linewidth\hsize
   \@setminipage}}%
 {\par\unskip\endMakeFramed%
 \at@end@of@kframe}
\makeatother

\definecolor{shadecolor}{rgb}{.97, .97, .97}
\definecolor{messagecolor}{rgb}{0, 0, 0}
\definecolor{warningcolor}{rgb}{1, 0, 1}
\definecolor{errorcolor}{rgb}{1, 0, 0}
\newenvironment{knitrout}{}{} % an empty environment to be redefined in TeX

\usepackage{alltt}

\usepackage{indentfirst}
\usepackage{amsmath}
\usepackage{amssymb}
\usepackage{natbib}
\usepackage{url}
\usepackage{graphicx}
\usepackage{tikz-cd}
\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}}
\IfFileExists{upquote.sty}{\usepackage{upquote}}{}
\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 cut-and-pasted 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]{aster-package}.  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{aster2-package} 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{kingsolver-et-al}.
That study estimated selection on each fitness component separately,
and therefore could not quantify how selection operates over the lifecycle,
nor identify potential trade-offs 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}


This document was processed with R, version 3.1.3, using the CRAN
packages \texttt{aster2}, version 0.2.1,
and \texttt{knitr}, version 1.10.5.
\begin{knitrout}
\definecolor{shadecolor}{rgb}{0.969, 0.969, 0.969}\color{fgcolor}\begin{kframe}
\begin{alltt}
\hlkwd{library}\hlstd{(aster2)}
\end{alltt}


{\ttfamily\noindent\itshape\color{messagecolor}{\#\# This is beta software.\\\#\# Unless you need to do aster models with dependence groups,\\\#\#\ \ \ \  use package "{}aster"{} instead.\\\#\# See help(aster2-package) for differences from package "{}aster"{}\\\#\#\ \ \ \  and examples.}}\end{kframe}
\end{knitrout}

\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).
\begin{knitrout}
\definecolor{shadecolor}{rgb}{0.969, 0.969, 0.969}\color{fgcolor}\begin{kframe}
\begin{alltt}
\hlkwd{data}\hlstd{(hornworm)}
\hlstd{data} \hlkwb{<-} \hlstd{hornworm}
\end{alltt}
\end{kframe}
\end{knitrout}

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 pre-eclosion, 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:lambda-weights}.  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:aster-transform}
   \varphi_j = \theta_j -
   \sum_{\substack{G \in \mathcal{G} \\ q(G) = j}} c_G(\theta_G)
\end{equation}
and
\begin{equation} \label{eq:cumulant-function}
   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:aster-transform} 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:cumulant-function} 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 right-hand 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]{geyer-gdor}.
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 zero-truncated 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{aster-aphid}, which had to use
a conditional submodel because of time-dependent 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}.

\begin{knitrout}
\definecolor{shadecolor}{rgb}{0.969, 0.969, 0.969}\color{fgcolor}\begin{kframe}
\begin{alltt}
\hlstd{offspring} \hlkwb{<-} \hlkwd{as.numeric}\hlstd{(}\hlkwd{grepl}\hlstd{(}\hlstr{"B"}\hlstd{,data}\hlopt{$}\hlstd{redata}\hlopt{$}\hlstd{varb))}
\hlstd{modmat.alt} \hlkwb{<-} \hlkwd{model.matrix}\hlstd{(resp} \hlopt{~} \hlstd{varb} \hlopt{+}
  \hlstd{offspring}\hlopt{:}\hlstd{(Mass_Repro} \hlopt{+} \hlstd{Mass_2nd} \hlopt{+}
  \hlkwd{I}\hlstd{(Mass_Repro}\hlopt{^}\hlnum{2}\hlstd{)} \hlopt{+} \hlkwd{I}\hlstd{(Mass_2nd}\hlopt{^}\hlnum{2}\hlstd{)} \hlopt{+} \hlkwd{I}\hlstd{(Mass_2nd}\hlopt{*}\hlstd{Mass_Repro)} \hlopt{+}
  \hlstd{Time_2nd} \hlopt{+} \hlkwd{I}\hlstd{(Time_2nd}\hlopt{^}\hlnum{2}\hlstd{)} \hlopt{+} \hlkwd{I}\hlstd{(Mass_2nd}\hlopt{*}\hlstd{Time_2nd)} \hlopt{+}
  \hlkwd{I}\hlstd{(Time_2nd}\hlopt{*}\hlstd{Mass_Repro)),} \hlkwc{data} \hlstd{= data}\hlopt{$}\hlstd{redata)}
\hlstd{tau.alt} \hlkwb{<-} \hlkwd{crossprod}\hlstd{(modmat.alt, data}\hlopt{$}\hlstd{redata}\hlopt{$}\hlstd{resp)}
\end{alltt}
\end{kframe}
\end{knitrout}

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.
\begin{knitrout}
\definecolor{shadecolor}{rgb}{0.969, 0.969, 0.969}\color{fgcolor}\begin{kframe}
\begin{alltt}
\hlstd{modmat.null} \hlkwb{<-} \hlkwd{model.matrix}\hlstd{(resp} \hlopt{~} \hlstd{varb} \hlopt{+}
  \hlstd{offspring}\hlopt{:}\hlstd{(Mass_Repro} \hlopt{+} \hlstd{Mass_2nd} \hlopt{+}
  \hlkwd{I}\hlstd{(Mass_Repro}\hlopt{^}\hlnum{2}\hlstd{)} \hlopt{+} \hlkwd{I}\hlstd{(Mass_2nd}\hlopt{^}\hlnum{2}\hlstd{)} \hlopt{+} \hlkwd{I}\hlstd{(Mass_2nd}\hlopt{*}\hlstd{Mass_Repro)} \hlopt{+}
  \hlstd{Time_2nd),} \hlkwc{data} \hlstd{= data}\hlopt{$}\hlstd{redata)}
\hlstd{tau} \hlkwb{<-} \hlkwd{crossprod}\hlstd{(modmat.null, data}\hlopt{$}\hlstd{redata}\hlopt{$}\hlstd{resp)}
\hlstd{beta} \hlkwb{<-} \hlstd{beta.null} \hlkwb{<-} \hlkwd{transformUnconditional}\hlstd{(tau, modmat.null,}
  \hlstd{data,} \hlkwc{from} \hlstd{=} \hlstr{"tau"}\hlstd{,} \hlkwc{to} \hlstd{=} \hlstr{"beta"}\hlstd{,} \hlkwc{tolerance} \hlstd{=} \hlnum{1e-100}\hlstd{)}
\end{alltt}
\end{kframe}
\end{knitrout}


The term on the right in (\ref{score}) is calculated using the
\texttt{transformUnconditional} function.


\begin{knitrout}
\definecolor{shadecolor}{rgb}{0.969, 0.969, 0.969}\color{fgcolor}\begin{kframe}
\begin{alltt}
\hlstd{tau.null} \hlkwb{<-} \hlkwd{transformUnconditional}\hlstd{(}\hlkwd{c}\hlstd{(beta.null,}\hlkwd{rep}\hlstd{(}\hlnum{0}\hlstd{,}\hlnum{3}\hlstd{)),}
  \hlstd{modmat.alt, data,} \hlkwc{from} \hlstd{=} \hlstr{"beta"}\hlstd{,} \hlkwc{to} \hlstd{=} \hlstr{"tau"}\hlstd{)}
\hlstd{score} \hlkwb{<-} \hlstd{tau.alt} \hlopt{-} \hlstd{tau.null}
\end{alltt}
\end{kframe}
\end{knitrout}


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).

\begin{knitrout}
\definecolor{shadecolor}{rgb}{0.969, 0.969, 0.969}\color{fgcolor}\begin{kframe}
\begin{alltt}
\hlstd{Fisher.null} \hlkwb{<-} \hlkwd{jacobian}\hlstd{(}\hlkwd{c}\hlstd{(beta.null,}\hlkwd{rep}\hlstd{(}\hlnum{0}\hlstd{,}\hlnum{3}\hlstd{)), data,}
  \hlkwc{transform} \hlstd{=} \hlstr{"unconditional"}\hlstd{,} \hlkwc{from} \hlstd{=} \hlstr{"beta"}\hlstd{,}
  \hlkwc{to} \hlstd{=} \hlstr{"tau"}\hlstd{,} \hlkwc{modmat} \hlstd{= modmat.alt)}
\end{alltt}
\end{kframe}
\end{knitrout}


We now construct the Rao test statistic. The reference distribution for this test is $\chi^2_3$ and the p-value suggests that we fail to reject the null hypothesis when testing at the 0.05 significance level.


\begin{knitrout}
\definecolor{shadecolor}{rgb}{0.969, 0.969, 0.969}\color{fgcolor}\begin{kframe}
\begin{alltt}
\hlstd{Rao} \hlkwb{<-} \hlkwd{t}\hlstd{(score)} \hlopt{%*%} \hlkwd{solve}\hlstd{(Fisher.null,} \hlkwc{tol} \hlstd{=} \hlnum{1e-100}\hlstd{)} \hlopt{%*%} \hlstd{score}
\hlkwd{pchisq}\hlstd{(Rao,} \hlkwc{df} \hlstd{=} \hlnum{3}\hlstd{,} \hlkwc{lower} \hlstd{=} \hlnum{FALSE}\hlstd{)}
\end{alltt}
\begin{verbatim}
##           [,1]
## [1,] 0.1048548
\end{verbatim}
\end{kframe}
\end{knitrout}


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 p-value 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.


\begin{knitrout}
\definecolor{shadecolor}{rgb}{0.969, 0.969, 0.969}\color{fgcolor}\begin{kframe}
\begin{alltt}
\hlstd{modmat.null} \hlkwb{<-} \hlkwd{model.matrix}\hlstd{(resp} \hlopt{~} \hlstd{varb} \hlopt{+}
  \hlstd{offspring}\hlopt{:}\hlstd{(Mass_Repro} \hlopt{+} \hlkwd{I}\hlstd{(Mass_Repro}\hlopt{^}\hlnum{2}\hlstd{)} \hlopt{+} \hlstd{Time_2nd} \hlopt{+}
    \hlstd{Mass_2nd),}
  \hlkwc{data} \hlstd{= data}\hlopt{$}\hlstd{redata)}
\hlstd{tau} \hlkwb{<-} \hlkwd{crossprod}\hlstd{(modmat.null, data}\hlopt{$}\hlstd{redata}\hlopt{$}\hlstd{resp)}
\hlstd{beta.null} \hlkwb{<-} \hlkwd{transformUnconditional}\hlstd{(tau, modmat.null,}
  \hlstd{data,} \hlkwc{from} \hlstd{=} \hlstr{"tau"}\hlstd{,} \hlkwc{to} \hlstd{=} \hlstr{"beta"}\hlstd{,} \hlkwc{tolerance} \hlstd{=} \hlnum{1e-100}\hlstd{)}

\hlstd{modmat.alt} \hlkwb{<-} \hlkwd{model.matrix}\hlstd{(resp} \hlopt{~} \hlstd{varb} \hlopt{+}
  \hlstd{offspring}\hlopt{:}\hlstd{(Mass_Repro} \hlopt{+} \hlkwd{I}\hlstd{(Mass_Repro}\hlopt{^}\hlnum{2}\hlstd{)} \hlopt{+} \hlstd{Time_2nd} \hlopt{+}
  \hlstd{Mass_2nd} \hlopt{+} \hlkwd{I}\hlstd{(Mass_Repro}\hlopt{*}\hlstd{Mass_2nd)} \hlopt{+} \hlkwd{I}\hlstd{(Mass_2nd}\hlopt{^}\hlnum{2}\hlstd{)),}
  \hlkwc{data} \hlstd{= data}\hlopt{$}\hlstd{redata)}
\hlstd{tau.alt} \hlkwb{<-} \hlkwd{crossprod}\hlstd{(modmat.alt, data}\hlopt{$}\hlstd{redata}\hlopt{$}\hlstd{resp)}

\hlstd{tau.null} \hlkwb{<-} \hlkwd{transformUnconditional}\hlstd{(}\hlkwd{c}\hlstd{(beta.null,}\hlkwd{rep}\hlstd{(}\hlnum{0}\hlstd{,}\hlnum{2}\hlstd{)),}
  \hlstd{modmat.alt, data,} \hlkwc{from} \hlstd{=} \hlstr{"beta"}\hlstd{,} \hlkwc{to} \hlstd{=} \hlstr{"tau"}\hlstd{)}
\hlstd{score} \hlkwb{<-} \hlstd{tau.alt} \hlopt{-} \hlstd{tau.null}

\hlstd{Fisher.null} \hlkwb{<-} \hlkwd{jacobian}\hlstd{(}\hlkwd{c}\hlstd{(beta.null,}\hlkwd{rep}\hlstd{(}\hlnum{0}\hlstd{,}\hlnum{2}\hlstd{)), data,}
  \hlkwc{transform} \hlstd{=} \hlstr{"unconditional"}\hlstd{,} \hlkwc{from} \hlstd{=} \hlstr{"beta"}\hlstd{,}
  \hlkwc{to} \hlstd{=} \hlstr{"tau"}\hlstd{,} \hlkwc{modmat} \hlstd{= modmat.alt)}
\hlstd{Rao} \hlkwb{<-} \hlkwd{t}\hlstd{(score)} \hlopt{%*%} \hlkwd{solve}\hlstd{(Fisher.null,} \hlkwc{tol} \hlstd{=} \hlnum{1e-100}\hlstd{)} \hlopt{%*%} \hlstd{score}
\hlkwd{pchisq}\hlstd{(Rao,} \hlkwc{df} \hlstd{=} \hlnum{2}\hlstd{,} \hlkwc{lower} \hlstd{=} \hlnum{FALSE}\hlstd{)}
\end{alltt}
\begin{verbatim}
##              [,1]
## [1,] 3.371453e-08
\end{verbatim}
\end{kframe}
\end{knitrout}


We now consider removing the quadratic terms for mass at eclosion. The p-value is small enough to reject this smaller model at any reasonable significance level.


\begin{knitrout}
\definecolor{shadecolor}{rgb}{0.969, 0.969, 0.969}\color{fgcolor}\begin{kframe}
\begin{alltt}
\hlstd{modmat.null} \hlkwb{<-} \hlkwd{model.matrix}\hlstd{(resp} \hlopt{~} \hlstd{varb} \hlopt{+}
  \hlstd{offspring}\hlopt{:}\hlstd{(Mass_2nd} \hlopt{+} \hlkwd{I}\hlstd{(Mass_2nd}\hlopt{^}\hlnum{2}\hlstd{)} \hlopt{+} \hlstd{Time_2nd} \hlopt{+}
    \hlstd{Mass_Repro),}
  \hlkwc{data} \hlstd{= data}\hlopt{$}\hlstd{redata)}
\hlstd{tau} \hlkwb{<-} \hlkwd{crossprod}\hlstd{(modmat.null, data}\hlopt{$}\hlstd{redata}\hlopt{$}\hlstd{resp)}
\hlstd{beta.null} \hlkwb{<-} \hlkwd{transformUnconditional}\hlstd{(tau, modmat.null,}
  \hlstd{data,} \hlkwc{from} \hlstd{=} \hlstr{"tau"}\hlstd{,} \hlkwc{to} \hlstd{=} \hlstr{"beta"}\hlstd{,} \hlkwc{tolerance} \hlstd{=} \hlnum{1e-100}\hlstd{)}

\hlstd{modmat.alt} \hlkwb{<-} \hlkwd{model.matrix}\hlstd{(resp} \hlopt{~} \hlstd{varb} \hlopt{+}
  \hlstd{offspring}\hlopt{:}\hlstd{(Mass_2nd} \hlopt{+} \hlkwd{I}\hlstd{(Mass_2nd}\hlopt{^}\hlnum{2}\hlstd{)} \hlopt{+} \hlstd{Time_2nd} \hlopt{+}
    \hlstd{Mass_Repro} \hlopt{+} \hlkwd{I}\hlstd{(Mass_Repro}\hlopt{^}\hlnum{2}\hlstd{)} \hlopt{+} \hlkwd{I}\hlstd{(Mass_Repro}\hlopt{*}\hlstd{Mass_2nd)),}
  \hlkwc{data} \hlstd{= data}\hlopt{$}\hlstd{redata)}
\hlstd{tau.alt} \hlkwb{<-} \hlkwd{crossprod}\hlstd{(modmat.alt, data}\hlopt{$}\hlstd{redata}\hlopt{$}\hlstd{resp)}

\hlstd{tau.null} \hlkwb{<-} \hlkwd{transformUnconditional}\hlstd{(}\hlkwd{c}\hlstd{(beta.null,}\hlkwd{rep}\hlstd{(}\hlnum{0}\hlstd{,}\hlnum{2}\hlstd{)),}
  \hlstd{modmat.alt, data,} \hlkwc{from} \hlstd{=} \hlstr{"beta"}\hlstd{,} \hlkwc{to} \hlstd{=} \hlstr{"tau"}\hlstd{)}
\hlstd{score} \hlkwb{<-} \hlstd{tau.alt} \hlopt{-} \hlstd{tau.null}

\hlstd{Fisher.null} \hlkwb{<-} \hlkwd{jacobian}\hlstd{(}\hlkwd{c}\hlstd{(beta.null,}\hlkwd{rep}\hlstd{(}\hlnum{0}\hlstd{,}\hlnum{2}\hlstd{)), data,}
  \hlkwc{transform} \hlstd{=} \hlstr{"unconditional"}\hlstd{,} \hlkwc{from} \hlstd{=} \hlstr{"beta"}\hlstd{,}
  \hlkwc{to} \hlstd{=} \hlstr{"tau"}\hlstd{,} \hlkwc{modmat} \hlstd{= modmat.alt)}
\hlstd{Rao} \hlkwb{<-} \hlkwd{t}\hlstd{(score)} \hlopt{%*%} \hlkwd{solve}\hlstd{(Fisher.null,} \hlkwc{tol} \hlstd{=} \hlnum{1e-100}\hlstd{)} \hlopt{%*%} \hlstd{score}
\hlkwd{pchisq}\hlstd{(Rao,} \hlkwc{df} \hlstd{=} \hlnum{2}\hlstd{,} \hlkwc{lower} \hlstd{=} \hlnum{FALSE}\hlstd{)}
\end{alltt}
\begin{verbatim}
##              [,1]
## [1,] 6.739063e-89
\end{verbatim}
\end{kframe}
\end{knitrout}


Finally, we consider removing the linear term for age that individuals reach their second instar stage. The p-value is small enough to reject this smaller model at any reasonable significance level.


\begin{knitrout}
\definecolor{shadecolor}{rgb}{0.969, 0.969, 0.969}\color{fgcolor}\begin{kframe}
\begin{alltt}
\hlstd{modmat.null} \hlkwb{<-} \hlkwd{model.matrix}\hlstd{(resp} \hlopt{~} \hlstd{varb} \hlopt{+}
  \hlstd{offspring}\hlopt{:}\hlstd{(Mass_2nd} \hlopt{+} \hlkwd{I}\hlstd{(Mass_2nd}\hlopt{^}\hlnum{2}\hlstd{)} \hlopt{+} \hlstd{Mass_Repro} \hlopt{+}
  \hlkwd{I}\hlstd{(Mass_Repro}\hlopt{^}\hlnum{2}\hlstd{)} \hlopt{+} \hlkwd{I}\hlstd{(Mass_Repro}\hlopt{*}\hlstd{Mass_2nd)),}
  \hlkwc{data} \hlstd{= data}\hlopt{$}\hlstd{redata)}
\hlstd{tau} \hlkwb{<-} \hlkwd{crossprod}\hlstd{(modmat.null, data}\hlopt{$}\hlstd{redata}\hlopt{$}\hlstd{resp)}
\hlstd{beta.null} \hlkwb{<-} \hlkwd{transformUnconditional}\hlstd{(tau, modmat.null,}
  \hlstd{data,} \hlkwc{from} \hlstd{=} \hlstr{"tau"}\hlstd{,} \hlkwc{to} \hlstd{=} \hlstr{"beta"}\hlstd{,} \hlkwc{tolerance} \hlstd{=} \hlnum{1e-100}\hlstd{)}

\hlstd{modmat.alt} \hlkwb{<-} \hlkwd{model.matrix}\hlstd{(resp} \hlopt{~} \hlstd{varb} \hlopt{+}
  \hlstd{offspring}\hlopt{:}\hlstd{(Mass_2nd} \hlopt{+} \hlkwd{I}\hlstd{(Mass_2nd}\hlopt{^}\hlnum{2}\hlstd{)} \hlopt{+} \hlstd{Mass_Repro} \hlopt{+}
  \hlkwd{I}\hlstd{(Mass_Repro}\hlopt{^}\hlnum{2}\hlstd{)} \hlopt{+} \hlkwd{I}\hlstd{(Mass_Repro}\hlopt{*}\hlstd{Mass_2nd)} \hlopt{+} \hlstd{Time_2nd),}
  \hlkwc{data} \hlstd{= data}\hlopt{$}\hlstd{redata)}
\hlstd{tau.alt} \hlkwb{<-} \hlkwd{crossprod}\hlstd{(modmat.alt, data}\hlopt{$}\hlstd{redata}\hlopt{$}\hlstd{resp)}

\hlstd{tau.null} \hlkwb{<-} \hlkwd{transformUnconditional}\hlstd{(}\hlkwd{c}\hlstd{(beta.null,}\hlnum{0}\hlstd{),}
  \hlstd{modmat.alt, data,} \hlkwc{from} \hlstd{=} \hlstr{"beta"}\hlstd{,} \hlkwc{to} \hlstd{=} \hlstr{"tau"}\hlstd{)}
\hlstd{score} \hlkwb{<-} \hlstd{tau.alt} \hlopt{-} \hlstd{tau.null}

\hlstd{Fisher.null} \hlkwb{<-} \hlkwd{jacobian}\hlstd{(}\hlkwd{c}\hlstd{(beta.null,}\hlnum{0}\hlstd{), data,}
  \hlkwc{transform} \hlstd{=} \hlstr{"unconditional"}\hlstd{,} \hlkwc{from} \hlstd{=} \hlstr{"beta"}\hlstd{,}
  \hlkwc{to} \hlstd{=} \hlstr{"tau"}\hlstd{,} \hlkwc{modmat} \hlstd{= modmat.alt)}
\hlstd{Rao} \hlkwb{<-} \hlkwd{t}\hlstd{(score)} \hlopt{%*%} \hlkwd{solve}\hlstd{(Fisher.null,} \hlkwc{tol} \hlstd{=} \hlnum{1e-100}\hlstd{)} \hlopt{%*%} \hlstd{score}
\hlkwd{pchisq}\hlstd{(Rao,} \hlkwc{df} \hlstd{=} \hlnum{1}\hlstd{,} \hlkwc{lower} \hlstd{=} \hlnum{FALSE}\hlstd{)}
\end{alltt}
\begin{verbatim}
##              [,1]
## [1,] 7.878988e-05
\end{verbatim}
\end{kframe}
\end{knitrout}








%<<echo = 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 mean-value 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 mean-value 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 mean-value parameter from the linear model in both mass traits is now calculated. We transform our MLE for the submodel canonical parameter to the mean-value 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)
%@

%<<eval = FALSE, echo = FALSE>>=
%foot <- crossprod(modmat.quad, data$redata$resp)
%foot[c(29:31)] <- 0
%blah <- transformUnconditional(foot, tolerance = 1e-30,
%  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 p-value where the reference distribution for this test is $\chi^2$ with three degrees of freedom. The p-value 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 = 1e-50) %*% 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 p-value 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 = 1e-50) %*% score
%Rao

%df <- length(beta.full) - length(beta.quad)
%pval <- pchisq(Rao, df = df, lower.tail = FALSE)
%pval
%@



%<<echo = 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)
%@



%<<>>=
%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 = 1e-50) %*% 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 four-dimensional graph (fitness versus three covariates), it cannot be
visualized, and we make one three-dimensional 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.
\begin{knitrout}
\definecolor{shadecolor}{rgb}{0.969, 0.969, 0.969}\color{fgcolor}\begin{kframe}
\begin{alltt}
\hlstd{vars} \hlkwb{<-} \hlkwd{levels}\hlstd{(hornworm}\hlopt{$}\hlstd{redata}\hlopt{$}\hlstd{varb)}
\hlstd{pred} \hlkwb{<-} \hlstd{hornworm}\hlopt{$}\hlstd{pred}
\hlstd{group} \hlkwb{<-} \hlstd{hornworm}\hlopt{$}\hlstd{group}
\hlstd{code} \hlkwb{<-} \hlstd{hornworm}\hlopt{$}\hlstd{code}
\hlstd{families} \hlkwb{<-} \hlstd{hornworm}\hlopt{$}\hlstd{families}
\end{alltt}
\end{kframe}
\end{knitrout}

Now the $101^2$ hypothetical individuals with unique mass traits reaching
their second instar larval stage at age 2 are generated

\begin{knitrout}
\definecolor{shadecolor}{rgb}{0.969, 0.969, 0.969}\color{fgcolor}\begin{kframe}
\begin{alltt}
\hlstd{x} \hlkwb{<-} \hlkwd{seq}\hlstd{(}\hlkwc{from} \hlstd{=} \hlnum{0}\hlstd{,} \hlkwc{to} \hlstd{=} \hlnum{0.016}\hlstd{,} \hlkwc{by} \hlstd{=} \hlnum{0.016}\hlopt{/}\hlnum{100}\hlstd{)}
\hlstd{y} \hlkwb{<-} \hlkwd{seq}\hlstd{(}\hlkwc{from} \hlstd{=} \hlnum{0}\hlstd{,} \hlkwc{to} \hlstd{=} \hlnum{2.3}\hlstd{,} \hlkwc{by} \hlstd{=} \hlnum{2.3}\hlopt{/}\hlnum{100}\hlstd{)}
\hlstd{days} \hlkwb{<-} \hlnum{2}\hlopt{:}\hlnum{6} \hlopt{/} \hlnum{7}
\hlstd{long.sex} \hlkwb{<-} \hlkwd{factor}\hlstd{(}\hlkwd{rep}\hlstd{(}\hlstr{"F"}\hlstd{,} \hlnum{101}\hlopt{^}\hlnum{2}\hlstd{),} \hlkwc{levels} \hlstd{=} \hlkwd{c}\hlstd{(}\hlstr{"F"}\hlstd{,} \hlstr{"M"}\hlstd{,} \hlstr{"U"}\hlstd{))}
\hlstd{mnew} \hlkwb{<-} \hlkwd{data.frame}\hlstd{(long.sex,} \hlkwd{rep}\hlstd{(days[}\hlnum{1}\hlstd{],} \hlnum{101}\hlopt{^}\hlnum{2}\hlstd{),}
    \hlkwd{rep}\hlstd{(x,} \hlkwc{times} \hlstd{=} \hlnum{101}\hlstd{),} \hlkwd{rep}\hlstd{(y,} \hlkwc{each} \hlstd{=} \hlnum{101}\hlstd{))}
\hlstd{mnew} \hlkwb{<-} \hlkwd{cbind}\hlstd{(mnew,} \hlkwd{matrix}\hlstd{(}\hlkwd{rep}\hlstd{(}\hlnum{0}\hlstd{,} \hlnum{26}\hlopt{*}\hlnum{101}\hlopt{^}\hlnum{2}\hlstd{),} \hlkwc{ncol} \hlstd{=} \hlnum{26}\hlstd{))}
\hlstd{mnames} \hlkwb{<-} \hlkwd{c}\hlstd{(}\hlstr{"Sex"}\hlstd{,} \hlstr{"Time_2nd"}\hlstd{,} \hlstr{"Mass_2nd"}\hlstd{,} \hlstr{"Mass_Repro"}\hlstd{, vars)}
\hlkwd{names}\hlstd{(mnew)} \hlkwb{<-} \hlstd{mnames}
\end{alltt}
\end{kframe}
\end{knitrout}


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.

\begin{knitrout}
\definecolor{shadecolor}{rgb}{0.969, 0.969, 0.969}\color{fgcolor}\begin{kframe}
\begin{alltt}
\hlstd{data.mnew} \hlkwb{<-} \hlkwd{asterdata}\hlstd{(mnew,} \hlkwc{vars} \hlstd{= vars,} \hlkwc{pred} \hlstd{= pred,}
  \hlkwc{group} \hlstd{= group,} \hlkwc{code} \hlstd{= code,} \hlkwc{families} \hlstd{= families)}
\hlstd{foo} \hlkwb{<-} \hlkwd{as.character}\hlstd{(data.mnew}\hlopt{$}\hlstd{redata}\hlopt{$}\hlstd{varb)}
\hlstd{offspring.mnew} \hlkwb{<-} \hlkwd{as.numeric}\hlstd{(}\hlkwd{grepl}\hlstd{(}\hlstr{"B"}\hlstd{, foo))}
\hlstd{data.mnew}\hlopt{$}\hlstd{redata} \hlkwb{<-} \hlkwd{transform}\hlstd{(data.mnew}\hlopt{$}\hlstd{redata,}
  \hlkwc{offspring.mnew} \hlstd{= offspring.mnew)}

\hlstd{modmat.mnew} \hlkwb{<-} \hlkwd{model.matrix}\hlstd{(resp} \hlopt{~} \hlstd{varb} \hlopt{+}
  \hlstd{offspring.mnew}\hlopt{:}\hlstd{(Mass_Repro} \hlopt{+} \hlstd{Mass_2nd} \hlopt{+}
  \hlkwd{I}\hlstd{(Mass_Repro}\hlopt{^}\hlnum{2}\hlstd{)} \hlopt{+} \hlkwd{I}\hlstd{(Mass_2nd}\hlopt{^}\hlnum{2}\hlstd{)} \hlopt{+} \hlkwd{I}\hlstd{(Mass_2nd}\hlopt{*}\hlstd{Mass_Repro)} \hlopt{+}
  \hlstd{Time_2nd),} \hlkwc{data} \hlstd{= data.mnew}\hlopt{$}\hlstd{redata)}
\end{alltt}
\end{kframe}
\end{knitrout}


This model matrix is then used to find the MLE of the saturated model
unconditional canonical parameter $\varphi$.

\begin{knitrout}
\definecolor{shadecolor}{rgb}{0.969, 0.969, 0.969}\color{fgcolor}\begin{kframe}
\begin{alltt}
\hlstd{phi} \hlkwb{<-} \hlstd{modmat.mnew} \hlopt{%*%} \hlstd{beta}
\end{alltt}
\end{kframe}
\end{knitrout}

Estimates on the canonical scale are not directly interpretable. We need estimates for the mean-value 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.

\begin{knitrout}
\definecolor{shadecolor}{rgb}{0.969, 0.969, 0.969}\color{fgcolor}\begin{kframe}
\begin{alltt}
\hlstd{mus2} \hlkwb{<-} \hlkwd{transformSaturated}\hlstd{(phi, data.mnew,} \hlkwc{from} \hlstd{=} \hlstr{"phi"}\hlstd{,}
  \hlkwc{to} \hlstd{=} \hlstr{"mu"}\hlstd{)}
\hlstd{ind} \hlkwb{<-} \hlkwd{which}\hlstd{(}\hlkwd{grepl}\hlstd{(}\hlstr{"B"}\hlstd{,}\hlkwd{rownames}\hlstd{(}\hlkwd{data.frame}\hlstd{(modmat.mnew))))}
\hlstd{qux2} \hlkwb{<-} \hlkwd{matrix}\hlstd{(mus2[ind],} \hlkwc{nrow} \hlstd{=} \hlnum{101}\hlopt{^}\hlnum{2}\hlstd{,} \hlkwc{ncol} \hlstd{=} \hlnum{8}\hlstd{)}
\hlkwd{colnames}\hlstd{(qux2)} \hlkwb{<-} \hlkwd{paste}\hlstd{(}\hlstr{"c"}\hlstd{,} \hlnum{1}\hlopt{:}\hlnum{8}\hlstd{,} \hlkwc{sep} \hlstd{=} \hlstr{""}\hlstd{)}
\hlkwd{rownames}\hlstd{(qux2)} \hlkwb{<-} \hlkwd{c}\hlstd{(}\hlnum{1}\hlopt{:}\hlnum{101}\hlopt{^}\hlnum{2}\hlstd{)}
\hlstd{sums2} \hlkwb{<-} \hlkwd{as.numeric}\hlstd{(}\hlkwd{apply}\hlstd{(qux2,} \hlnum{1}\hlstd{, sum))}
\hlstd{zday2} \hlkwb{<-} \hlkwd{matrix}\hlstd{(sums2,} \hlkwc{nrow} \hlstd{=} \hlnum{101}\hlstd{)}
\end{alltt}
\end{kframe}
\end{knitrout}

The \texttt{maximum} function, defined below,
is used to find the two mass values yielding the highest expected fitness.
\begin{knitrout}
\definecolor{shadecolor}{rgb}{0.969, 0.969, 0.969}\color{fgcolor}\begin{kframe}
\begin{alltt}
\hlstd{maximum} \hlkwb{<-} \hlkwa{function}\hlstd{(}\hlkwc{sums}\hlstd{)\{}
  \hlstd{max.ind} \hlkwb{<-} \hlkwd{which}\hlstd{(sums} \hlopt{==} \hlkwd{max}\hlstd{(sums))}
  \hlstd{column} \hlkwb{<-} \hlkwd{floor}\hlstd{(max.ind} \hlopt{/} \hlnum{101}\hlstd{)}
  \hlstd{row} \hlkwb{<-} \hlstd{max.ind} \hlopt{-} \hlnum{101} \hlopt{*} \hlstd{column}
  \hlstd{point} \hlkwb{<-} \hlkwd{c}\hlstd{(x[row], y[column])}
  \hlkwd{return}\hlstd{(point)}
\hlstd{\}}
\end{alltt}
\end{kframe}
\end{knitrout}

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 so-called 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.
\begin{knitrout}
\definecolor{shadecolor}{rgb}{0.969, 0.969, 0.969}\color{fgcolor}\begin{kframe}
\begin{alltt}
\hlstd{mass.second} \hlkwb{<-} \hlstd{hornworm}\hlopt{$}\hlstd{redata}\hlopt{$}\hlstd{Mass_2nd}
\hlstd{mass.reprod} \hlkwb{<-} \hlstd{hornworm}\hlopt{$}\hlstd{redata}\hlopt{$}\hlstd{Mass_Repro}
\hlstd{sex} \hlkwb{<-} \hlkwd{as.character}\hlstd{(hornworm}\hlopt{$}\hlstd{redata}\hlopt{$}\hlstd{Sex)}
\hlkwd{unique}\hlstd{(sex)}
\end{alltt}
\begin{verbatim}
## [1] "F" "M" "U"
\end{verbatim}
\begin{alltt}
\hlstd{time.second} \hlkwb{<-} \hlstd{hornworm}\hlopt{$}\hlstd{redata}\hlopt{$}\hlstd{Time_2nd}
\hlstd{id} \hlkwb{<-} \hlstd{data}\hlopt{$}\hlstd{redata}\hlopt{$}\hlstd{id}
\hlstd{u} \hlkwb{<-} \hlkwd{unique}\hlstd{(id)}
\hlstd{idx} \hlkwb{<-} \hlkwd{match}\hlstd{(u, id)}
\hlstd{mass.second} \hlkwb{<-} \hlstd{mass.second[idx]}
\hlstd{mass.reprod} \hlkwb{<-} \hlstd{mass.reprod[idx]}
\hlstd{sex} \hlkwb{<-} \hlstd{sex[idx]}
\hlstd{time.second} \hlkwb{<-} \hlstd{time.second[idx]}
\hlstd{mass.second} \hlkwb{<-} \hlstd{mass.second[sex} \hlopt{==} \hlstr{"F"}\hlstd{]}
\hlstd{mass.reprod} \hlkwb{<-} \hlstd{mass.reprod[sex} \hlopt{==} \hlstr{"F"}\hlstd{]}
\end{alltt}
\end{kframe}
\end{knitrout}

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}
\begin{knitrout}
\definecolor{shadecolor}{rgb}{0.969, 0.969, 0.969}\color{fgcolor}\begin{kframe}
\begin{alltt}
\hlkwd{plot}\hlstd{(mass.second,} \hlkwc{pch} \hlstd{=} \hlnum{20}\hlstd{,}
  \hlstd{mass.reprod,} \hlkwc{xlab} \hlstd{=} \hlstr{"Mass at 
  2nd instar"}\hlstd{,} \hlkwc{ylab} \hlstd{=} \hlstr{"Mass at Reproduction"}\hlstd{,}
  \hlkwc{main} \hlstd{=} \hlstr{"Fitness Landscape for Ovariole Counts vs. Mass"}\hlstd{)}
\hlkwd{points}\hlstd{(}\hlkwd{maximum}\hlstd{(sums2)[}\hlnum{1}\hlstd{],} \hlkwd{maximum}\hlstd{(sums2)[}\hlnum{2}\hlstd{],} \hlkwc{col} \hlstd{=} \hlstr{"red"}\hlstd{,}
  \hlkwc{pch} \hlstd{=} \hlnum{19}\hlstd{)}
\hlkwd{contour}\hlstd{(x,y,zday2,} \hlkwc{add} \hlstd{=} \hlnum{TRUE}\hlstd{,} \hlkwc{nlevels} \hlstd{=} \hlnum{8}\hlstd{,} \hlkwc{levels} \hlstd{=}
  \hlkwd{c}\hlstd{(}\hlnum{100}\hlstd{,}\hlnum{150}\hlstd{,}\hlnum{200}\hlstd{,}\hlnum{250}\hlstd{,}\hlnum{300}\hlstd{,}\hlnum{325}\hlstd{,}\hlnum{350}\hlstd{,}\hlnum{360}\hlstd{) )}
\end{alltt}
\end{kframe}
\includegraphics[width=\maxwidth]{figure/unnamed-chunk-7-1} 

\end{knitrout}
\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.


\begin{knitrout}
\definecolor{shadecolor}{rgb}{0.969, 0.969, 0.969}\color{fgcolor}\begin{kframe}
\begin{alltt}
\hlstd{mnew} \hlkwb{<-} \hlkwd{data.frame}\hlstd{(long.sex,} \hlkwd{rep}\hlstd{(days[}\hlnum{2}\hlstd{],} \hlnum{101}\hlopt{^}\hlnum{2}\hlstd{),}
    \hlkwd{rep}\hlstd{(x,} \hlkwc{times} \hlstd{=} \hlnum{101}\hlstd{),} \hlkwd{rep}\hlstd{(y,} \hlkwc{each} \hlstd{=} \hlnum{101}\hlstd{))}
\hlstd{mnew} \hlkwb{<-} \hlkwd{cbind}\hlstd{(mnew,} \hlkwd{matrix}\hlstd{(}\hlkwd{rep}\hlstd{(}\hlnum{0}\hlstd{,} \hlnum{26}\hlopt{*}\hlnum{101}\hlopt{^}\hlnum{2}\hlstd{),} \hlkwc{ncol} \hlstd{=} \hlnum{26}\hlstd{))}
\hlstd{mnames} \hlkwb{<-} \hlkwd{c}\hlstd{(}\hlstr{"Sex"}\hlstd{,} \hlstr{"Time_2nd"}\hlstd{,} \hlstr{"Mass_2nd"}\hlstd{,} \hlstr{"Mass_Repro"}\hlstd{, vars)}
\hlkwd{names}\hlstd{(mnew)} \hlkwb{<-} \hlstd{mnames}
\hlstd{data.mnew} \hlkwb{<-} \hlkwd{asterdata}\hlstd{(mnew,} \hlkwc{vars} \hlstd{= vars,} \hlkwc{pred} \hlstd{= pred,}
  \hlkwc{group} \hlstd{= group,} \hlkwc{code} \hlstd{= code,} \hlkwc{families} \hlstd{= families)}
\hlstd{foo} \hlkwb{<-} \hlkwd{as.character}\hlstd{(data.mnew}\hlopt{$}\hlstd{redata}\hlopt{$}\hlstd{varb)}
\hlstd{offspring.mnew} \hlkwb{<-} \hlkwd{as.numeric}\hlstd{(}\hlkwd{grepl}\hlstd{(}\hlstr{"B"}\hlstd{, foo))}
\hlstd{data.mnew}\hlopt{$}\hlstd{redata} \hlkwb{<-} \hlkwd{transform}\hlstd{(data.mnew}\hlopt{$}\hlstd{redata,}
  \hlkwc{offspring.mnew} \hlstd{= offspring.mnew)}

\hlstd{modmat.mnew} \hlkwb{<-} \hlkwd{model.matrix}\hlstd{(resp} \hlopt{~} \hlstd{varb} \hlopt{+}
  \hlstd{offspring.mnew}\hlopt{:}\hlstd{(Mass_Repro} \hlopt{+} \hlstd{Mass_2nd} \hlopt{+}
  \hlkwd{I}\hlstd{(Mass_Repro}\hlopt{^}\hlnum{2}\hlstd{)} \hlopt{+} \hlkwd{I}\hlstd{(Mass_2nd}\hlopt{^}\hlnum{2}\hlstd{)} \hlopt{+} \hlkwd{I}\hlstd{(Mass_2nd}\hlopt{*}\hlstd{Mass_Repro)} \hlopt{+}
  \hlstd{Time_2nd),} \hlkwc{data} \hlstd{= data.mnew}\hlopt{$}\hlstd{redata)}
\hlstd{phi} \hlkwb{<-} \hlkwd{crossprod}\hlstd{(}\hlkwd{t}\hlstd{(modmat.mnew), beta)}
\hlstd{mus3} \hlkwb{<-} \hlkwd{transformSaturated}\hlstd{(phi, data.mnew,} \hlkwc{from} \hlstd{=} \hlstr{"phi"}\hlstd{,}
  \hlkwc{to} \hlstd{=} \hlstr{"mu"}\hlstd{)}

\hlstd{ind} \hlkwb{<-} \hlkwd{which}\hlstd{(}\hlkwd{grepl}\hlstd{(}\hlstr{"B"}\hlstd{,}\hlkwd{rownames}\hlstd{(}\hlkwd{data.frame}\hlstd{(modmat.mnew))))}
\hlstd{qux3} \hlkwb{<-} \hlkwd{matrix}\hlstd{(mus3[ind],} \hlkwc{nrow} \hlstd{=} \hlnum{10201}\hlstd{,} \hlkwc{ncol} \hlstd{=} \hlnum{8}\hlstd{)}
\hlkwd{colnames}\hlstd{(qux3)} \hlkwb{<-} \hlkwd{paste}\hlstd{(}\hlstr{"c"}\hlstd{,} \hlnum{1}\hlopt{:}\hlnum{8}\hlstd{,} \hlkwc{sep} \hlstd{=} \hlstr{""}\hlstd{)}
\hlkwd{rownames}\hlstd{(qux3)} \hlkwb{<-} \hlkwd{c}\hlstd{(}\hlnum{1}\hlopt{:}\hlnum{10201}\hlstd{)}
\hlstd{sums3} \hlkwb{<-} \hlkwd{as.numeric}\hlstd{(}\hlkwd{apply}\hlstd{(qux3,} \hlnum{1}\hlstd{, sum))}
\hlstd{zday3} \hlkwb{<-} \hlkwd{matrix}\hlstd{(sums3,} \hlkwc{nrow} \hlstd{=} \hlnum{101}\hlstd{)}

\hlstd{mnew} \hlkwb{<-} \hlkwd{data.frame}\hlstd{(long.sex,} \hlkwd{rep}\hlstd{(days[}\hlnum{3}\hlstd{],} \hlnum{101}\hlopt{^}\hlnum{2}\hlstd{),}
    \hlkwd{rep}\hlstd{(x,} \hlkwc{times} \hlstd{=} \hlnum{101}\hlstd{),} \hlkwd{rep}\hlstd{(y,} \hlkwc{each} \hlstd{=} \hlnum{101}\hlstd{))}
\hlstd{mnew} \hlkwb{<-} \hlkwd{cbind}\hlstd{(mnew,} \hlkwd{matrix}\hlstd{(}\hlkwd{rep}\hlstd{(}\hlnum{0}\hlstd{,} \hlnum{26}\hlopt{*}\hlnum{101}\hlopt{^}\hlnum{2}\hlstd{),} \hlkwc{ncol} \hlstd{=} \hlnum{26}\hlstd{))}
\hlstd{mnames} \hlkwb{<-} \hlkwd{c}\hlstd{(}\hlstr{"Sex"}\hlstd{,} \hlstr{"Time_2nd"}\hlstd{,} \hlstr{"Mass_2nd"}\hlstd{,} \hlstr{"Mass_Repro"}\hlstd{, vars)}
\hlkwd{names}\hlstd{(mnew)} \hlkwb{<-} \hlstd{mnames}
\hlstd{data.mnew} \hlkwb{<-} \hlkwd{asterdata}\hlstd{(mnew,} \hlkwc{vars} \hlstd{= vars,} \hlkwc{pred} \hlstd{= pred,}
  \hlkwc{group} \hlstd{= group,} \hlkwc{code} \hlstd{= code,} \hlkwc{families} \hlstd{= families)}
\hlstd{foo} \hlkwb{<-} \hlkwd{as.character}\hlstd{(data.mnew}\hlopt{$}\hlstd{redata}\hlopt{$}\hlstd{varb)}
\hlstd{offspring.mnew} \hlkwb{<-} \hlkwd{as.numeric}\hlstd{(}\hlkwd{grepl}\hlstd{(}\hlstr{"B"}\hlstd{, foo))}
\hlstd{data.mnew}\hlopt{$}\hlstd{redata} \hlkwb{<-} \hlkwd{transform}\hlstd{(data.mnew}\hlopt{$}\hlstd{redata,}
  \hlkwc{offspring.mnew} \hlstd{= offspring.mnew)}

\hlstd{modmat.mnew} \hlkwb{<-} \hlkwd{model.matrix}\hlstd{(resp} \hlopt{~} \hlstd{varb} \hlopt{+}
  \hlstd{offspring.mnew}\hlopt{:}\hlstd{(Mass_Repro} \hlopt{+} \hlstd{Mass_2nd} \hlopt{+}
  \hlkwd{I}\hlstd{(Mass_Repro}\hlopt{^}\hlnum{2}\hlstd{)} \hlopt{+} \hlkwd{I}\hlstd{(Mass_2nd}\hlopt{^}\hlnum{2}\hlstd{)} \hlopt{+} \hlkwd{I}\hlstd{(Mass_2nd}\hlopt{*}\hlstd{Mass_Repro)} \hlopt{+}
  \hlstd{Time_2nd),} \hlkwc{data} \hlstd{= data.mnew}\hlopt{$}\hlstd{redata)}
\hlstd{phi} \hlkwb{<-} \hlkwd{crossprod}\hlstd{(}\hlkwd{t}\hlstd{(modmat.mnew), beta)}
\hlstd{mus4} \hlkwb{<-} \hlkwd{transformSaturated}\hlstd{(phi, data.mnew,} \hlkwc{from} \hlstd{=} \hlstr{"phi"}\hlstd{,}
  \hlkwc{to} \hlstd{=} \hlstr{"mu"}\hlstd{)}

\hlstd{ind} \hlkwb{<-} \hlkwd{which}\hlstd{(}\hlkwd{grepl}\hlstd{(}\hlstr{"B"}\hlstd{,}\hlkwd{rownames}\hlstd{(}\hlkwd{data.frame}\hlstd{(modmat.mnew))))}
\hlstd{qux4} \hlkwb{<-} \hlkwd{matrix}\hlstd{(mus4[ind],} \hlkwc{nrow} \hlstd{=} \hlnum{10201}\hlstd{,} \hlkwc{ncol} \hlstd{=} \hlnum{8}\hlstd{)}
\hlkwd{colnames}\hlstd{(qux4)} \hlkwb{<-} \hlkwd{paste}\hlstd{(}\hlstr{"c"}\hlstd{,} \hlnum{1}\hlopt{:}\hlnum{8}\hlstd{,} \hlkwc{sep} \hlstd{=} \hlstr{""}\hlstd{)}
\hlkwd{rownames}\hlstd{(qux4)} \hlkwb{<-} \hlkwd{c}\hlstd{(}\hlnum{1}\hlopt{:}\hlnum{10201}\hlstd{)}
\hlstd{sums4} \hlkwb{<-} \hlkwd{as.numeric}\hlstd{(}\hlkwd{apply}\hlstd{(qux4,} \hlnum{1}\hlstd{, sum))}
\hlstd{zday4} \hlkwb{<-} \hlkwd{matrix}\hlstd{(sums4,} \hlkwc{nrow} \hlstd{=} \hlnum{101}\hlstd{)}

\hlstd{mnew} \hlkwb{<-} \hlkwd{data.frame}\hlstd{(long.sex,} \hlkwd{rep}\hlstd{(days[}\hlnum{4}\hlstd{],} \hlnum{101}\hlopt{^}\hlnum{2}\hlstd{),}
    \hlkwd{rep}\hlstd{(x,} \hlkwc{times} \hlstd{=} \hlnum{101}\hlstd{),} \hlkwd{rep}\hlstd{(y,} \hlkwc{each} \hlstd{=} \hlnum{101}\hlstd{))}
\hlstd{mnew} \hlkwb{<-} \hlkwd{cbind}\hlstd{(mnew,} \hlkwd{matrix}\hlstd{(}\hlkwd{rep}\hlstd{(}\hlnum{0}\hlstd{,} \hlnum{26}\hlopt{*}\hlnum{101}\hlopt{^}\hlnum{2}\hlstd{),} \hlkwc{ncol} \hlstd{=} \hlnum{26}\hlstd{))}
\hlstd{mnames} \hlkwb{<-} \hlkwd{c}\hlstd{(}\hlstr{"Sex"}\hlstd{,} \hlstr{"Time_2nd"}\hlstd{,} \hlstr{"Mass_2nd"}\hlstd{,} \hlstr{"Mass_Repro"}\hlstd{, vars)}
\hlkwd{names}\hlstd{(mnew)} \hlkwb{<-} \hlstd{mnames}
\hlstd{data.mnew} \hlkwb{<-} \hlkwd{asterdata}\hlstd{(mnew,} \hlkwc{vars} \hlstd{= vars,} \hlkwc{pred} \hlstd{= pred,}
  \hlkwc{group} \hlstd{= group,} \hlkwc{code} \hlstd{= code,} \hlkwc{families} \hlstd{= families)}
\hlstd{foo} \hlkwb{<-} \hlkwd{as.character}\hlstd{(data.mnew}\hlopt{$}\hlstd{redata}\hlopt{$}\hlstd{varb)}
\hlstd{offspring.mnew} \hlkwb{<-} \hlkwd{as.numeric}\hlstd{(}\hlkwd{grepl}\hlstd{(}\hlstr{"B"}\hlstd{, foo))}
\hlstd{data.mnew}\hlopt{$}\hlstd{redata} \hlkwb{<-} \hlkwd{transform}\hlstd{(data.mnew}\hlopt{$}\hlstd{redata,}
  \hlkwc{offspring.mnew} \hlstd{= offspring.mnew)}

\hlstd{modmat.mnew} \hlkwb{<-} \hlkwd{model.matrix}\hlstd{(resp} \hlopt{~} \hlstd{varb} \hlopt{+}
  \hlstd{offspring.mnew}\hlopt{:}\hlstd{(Mass_Repro} \hlopt{+} \hlstd{Mass_2nd} \hlopt{+}
  \hlkwd{I}\hlstd{(Mass_Repro}\hlopt{^}\hlnum{2}\hlstd{)} \hlopt{+} \hlkwd{I}\hlstd{(Mass_2nd}\hlopt{^}\hlnum{2}\hlstd{)} \hlopt{+} \hlkwd{I}\hlstd{(Mass_2nd}\hlopt{*}\hlstd{Mass_Repro)} \hlopt{+}
  \hlstd{Time_2nd),} \hlkwc{data} \hlstd{= data.mnew}\hlopt{$}\hlstd{redata)}
\hlstd{phi} \hlkwb{<-} \hlkwd{crossprod}\hlstd{(}\hlkwd{t}\hlstd{(modmat.mnew), beta)}
\hlstd{mus5} \hlkwb{<-} \hlkwd{transformSaturated}\hlstd{(phi, data.mnew,} \hlkwc{from} \hlstd{=} \hlstr{"phi"}\hlstd{,}
  \hlkwc{to} \hlstd{=} \hlstr{"mu"}\hlstd{)}

\hlstd{ind} \hlkwb{<-} \hlkwd{which}\hlstd{(}\hlkwd{grepl}\hlstd{(}\hlstr{"B"}\hlstd{,}\hlkwd{rownames}\hlstd{(}\hlkwd{data.frame}\hlstd{(modmat.mnew))))}
\hlstd{qux5} \hlkwb{<-} \hlkwd{matrix}\hlstd{(mus5[ind],} \hlkwc{nrow} \hlstd{=} \hlnum{10201}\hlstd{,} \hlkwc{ncol} \hlstd{=} \hlnum{8}\hlstd{)}
\hlkwd{colnames}\hlstd{(qux5)} \hlkwb{<-} \hlkwd{paste}\hlstd{(}\hlstr{"c"}\hlstd{,} \hlnum{1}\hlopt{:}\hlnum{8}\hlstd{,} \hlkwc{sep} \hlstd{=} \hlstr{""}\hlstd{)}
\hlkwd{rownames}\hlstd{(qux5)} \hlkwb{<-} \hlkwd{c}\hlstd{(}\hlnum{1}\hlopt{:}\hlnum{10201}\hlstd{)}
\hlstd{sums5} \hlkwb{<-} \hlkwd{as.numeric}\hlstd{(}\hlkwd{apply}\hlstd{(qux5,} \hlnum{1}\hlstd{, sum))}
\hlstd{zday5} \hlkwb{<-} \hlkwd{matrix}\hlstd{(sums5,} \hlkwc{nrow} \hlstd{=} \hlnum{101}\hlstd{)}
\end{alltt}
\end{kframe}
\end{knitrout}



%\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) )
%@







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}.



\begin{knitrout}
\definecolor{shadecolor}{rgb}{0.969, 0.969, 0.969}\color{fgcolor}\begin{kframe}
\begin{alltt}
\hlkwd{par}\hlstd{(}\hlkwc{mfrow} \hlstd{=} \hlkwd{c}\hlstd{(}\hlnum{2}\hlstd{,} \hlnum{2}\hlstd{),} \hlkwc{mar} \hlstd{=} \hlkwd{rep}\hlstd{(}\hlnum{0.5}\hlstd{,} \hlnum{4}\hlstd{),}
  \hlkwc{oma} \hlstd{=} \hlkwd{c}\hlstd{(}\hlnum{4}\hlstd{,} \hlnum{4}\hlstd{,} \hlnum{0}\hlstd{,} \hlnum{0}\hlstd{)} \hlopt{+} \hlnum{0.1}\hlstd{)}
\hlcom{# day 2}
\hlkwd{plot}\hlstd{(mass.second,}
  \hlstd{mass.reprod,}
  \hlkwc{axes} \hlstd{=} \hlnum{FALSE}\hlstd{,} \hlkwc{xlab} \hlstd{=} \hlstr{""}\hlstd{,} \hlkwc{ylab} \hlstd{=} \hlstr{""}\hlstd{,} \hlkwc{pch} \hlstd{=} \hlnum{20}\hlstd{)}
\hlkwd{box}\hlstd{()}
\hlkwd{axis}\hlstd{(}\hlkwc{side} \hlstd{=} \hlnum{2}\hlstd{,} \hlkwc{outer} \hlstd{=} \hlnum{TRUE}\hlstd{)}
\hlkwd{mtext}\hlstd{(}\hlstr{"mass at eclosion"}\hlstd{,} \hlkwc{outer} \hlstd{=} \hlnum{TRUE}\hlstd{,} \hlkwc{line} \hlstd{=} \hlnum{3}\hlstd{,}
  \hlkwc{side} \hlstd{=} \hlnum{2}\hlstd{,} \hlkwc{at} \hlstd{=} \hlnum{0.50}\hlstd{)}
\hlkwd{mtext}\hlstd{(}\hlstr{"mass at 2nd instar stage"}\hlstd{,} \hlkwc{outer} \hlstd{=} \hlnum{TRUE}\hlstd{,} \hlkwc{line} \hlstd{=} \hlnum{3}\hlstd{,}
  \hlkwc{side} \hlstd{=} \hlnum{1}\hlstd{,} \hlkwc{at} \hlstd{=} \hlnum{0.50}\hlstd{)}
\hlkwd{points}\hlstd{(}\hlkwd{maximum}\hlstd{(sums2)[}\hlnum{1}\hlstd{],} \hlkwd{maximum}\hlstd{(sums2)[}\hlnum{2}\hlstd{],} \hlkwc{pch} \hlstd{=} \hlnum{0}\hlstd{)}
\hlstd{levels} \hlkwb{<-} \hlkwd{seq}\hlstd{(}\hlnum{50}\hlstd{,} \hlnum{350}\hlstd{,} \hlkwc{by} \hlstd{=} \hlnum{50}\hlstd{)}
\hlkwd{contour}\hlstd{(x, y, zday2,} \hlkwc{add} \hlstd{=} \hlnum{TRUE}\hlstd{,} \hlkwc{levels} \hlstd{= levels)}
\hlcom{# day 3}
\hlkwd{plot}\hlstd{(mass.second,}
  \hlstd{mass.reprod,}
  \hlkwc{axes} \hlstd{=} \hlnum{FALSE}\hlstd{,} \hlkwc{xlab} \hlstd{=} \hlstr{""}\hlstd{,} \hlkwc{ylab} \hlstd{=} \hlstr{""}\hlstd{,} \hlkwc{pch} \hlstd{=} \hlnum{20}\hlstd{)}
\hlkwd{box}\hlstd{()}
\hlkwd{points}\hlstd{(}\hlkwd{maximum}\hlstd{(sums3)[}\hlnum{1}\hlstd{],} \hlkwd{maximum}\hlstd{(sums3)[}\hlnum{2}\hlstd{],} \hlkwc{pch} \hlstd{=} \hlnum{0}\hlstd{)}
\hlkwd{contour}\hlstd{(x, y, zday3,} \hlkwc{add} \hlstd{=} \hlnum{TRUE}\hlstd{,} \hlkwc{levels} \hlstd{= levels)}
\hlcom{# day 4}
\hlkwd{plot}\hlstd{(mass.second,}
  \hlstd{mass.reprod,}
  \hlkwc{axes} \hlstd{=} \hlnum{FALSE}\hlstd{,} \hlkwc{xlab} \hlstd{=} \hlstr{""}\hlstd{,} \hlkwc{ylab} \hlstd{=} \hlstr{""}\hlstd{,} \hlkwc{pch} \hlstd{=} \hlnum{20}\hlstd{)}
\hlkwd{box}\hlstd{()}
\hlkwd{axis}\hlstd{(}\hlkwc{side} \hlstd{=} \hlnum{1}\hlstd{,} \hlkwc{outer} \hlstd{=} \hlnum{TRUE}\hlstd{)}
\hlkwd{axis}\hlstd{(}\hlkwc{side} \hlstd{=} \hlnum{2}\hlstd{,} \hlkwc{outer} \hlstd{=} \hlnum{TRUE}\hlstd{)}
\hlkwd{points}\hlstd{(}\hlkwd{maximum}\hlstd{(sums4)[}\hlnum{1}\hlstd{],} \hlkwd{maximum}\hlstd{(sums4)[}\hlnum{2}\hlstd{],} \hlkwc{pch} \hlstd{=} \hlnum{0}\hlstd{)}
\hlkwd{contour}\hlstd{(x, y, zday4,} \hlkwc{add} \hlstd{=} \hlnum{TRUE}\hlstd{,} \hlkwc{levels} \hlstd{= levels)}
\hlcom{# day 5}
\hlkwd{plot}\hlstd{(mass.second,}
  \hlstd{mass.reprod,}
  \hlkwc{axes} \hlstd{=} \hlnum{FALSE}\hlstd{,} \hlkwc{xlab} \hlstd{=} \hlstr{""}\hlstd{,} \hlkwc{ylab} \hlstd{=} \hlstr{""}\hlstd{,} \hlkwc{pch} \hlstd{=} \hlnum{20}\hlstd{)}
\hlkwd{box}\hlstd{()}
\hlkwd{axis}\hlstd{(}\hlkwc{side} \hlstd{=} \hlnum{1}\hlstd{,} \hlkwc{outer} \hlstd{=} \hlnum{TRUE}\hlstd{)}
\hlkwd{points}\hlstd{(}\hlkwd{maximum}\hlstd{(sums5)[}\hlnum{1}\hlstd{],} \hlkwd{maximum}\hlstd{(sums5)[}\hlnum{2}\hlstd{],} \hlkwc{pch} \hlstd{=} \hlnum{0}\hlstd{)}
\hlkwd{contour}\hlstd{(x, y, zday5,} \hlkwc{add} \hlstd{=} \hlnum{TRUE}\hlstd{,} \hlkwc{levels} \hlstd{= levels)}
\end{alltt}
\end{kframe}
\end{knitrout}






\begin{figure}
\begin{center}
%%%%% DON'T REPEAT YOURSELF (the DRY rule)
%%%%% Quote chunks, don't repeat chunks.
\begin{knitrout}
\definecolor{shadecolor}{rgb}{0.969, 0.969, 0.969}\color{fgcolor}
\includegraphics[width=\maxwidth]{figure/unnamed-chunk-9-1} 

\end{knitrout}
\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
363.6 (top left),
342.4 (top right),
322.4 (bottom left),
and 303.6 (bottom-right). }
\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. 


\begin{knitrout}
\definecolor{shadecolor}{rgb}{0.969, 0.969, 0.969}\color{fgcolor}\begin{kframe}
\begin{alltt}
\hlstd{mu.star} \hlkwb{<-} \hlstd{tau.alt[}\hlkwd{grepl}\hlstd{(}\hlstr{"varbB"}\hlstd{,} \hlkwd{rownames}\hlstd{(tau.alt))]}
\end{alltt}
\end{kframe}
\end{knitrout}


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.


\begin{knitrout}
\definecolor{shadecolor}{rgb}{0.969, 0.969, 0.969}\color{fgcolor}\begin{kframe}
\begin{alltt}
\hlstd{nind} \hlkwb{<-} \hlkwd{length}\hlstd{(}\hlkwd{unique}\hlstd{(hornworm}\hlopt{$}\hlstd{redata}\hlopt{$}\hlstd{id))}
\hlstd{nind}
\end{alltt}
\begin{verbatim}
## [1] 162
\end{verbatim}
\begin{alltt}
\hlstd{n} \hlkwb{<-} \hlstd{nind}
\hlstd{gr} \hlkwb{<-} \hlkwa{function}\hlstd{(}\hlkwc{lam}\hlstd{)\{}
  \hlstd{n} \hlopt{-} \hlkwd{sum}\hlstd{(mu.star} \hlopt{*} \hlkwd{exp}\hlstd{(}\hlopt{-}\hlstd{(lam} \hlopt{*} \hlstd{(}\hlnum{32} \hlopt{+} \hlkwd{seq}\hlstd{(}\hlkwc{along} \hlstd{= mu.star)))))}
\hlstd{\}}
\hlstd{lout} \hlkwb{<-} \hlkwd{uniroot}\hlstd{(gr,} \hlkwc{lower} \hlstd{=} \hlopt{-}\hlnum{105}\hlstd{,} \hlkwc{upper} \hlstd{=} \hlnum{105}\hlstd{)}
\hlstd{lambda.hat} \hlkwb{<-} \hlstd{lout}\hlopt{$}\hlstd{root}
\hlstd{lambda.hat}
\end{alltt}
\begin{verbatim}
## [1] 0.1215653
\end{verbatim}
\end{kframe}
\end{knitrout}

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} = 0.122$. 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}) = 1.129$ 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:stable-reweight}
   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


\begin{knitrout}
\definecolor{shadecolor}{rgb}{0.969, 0.969, 0.969}\color{fgcolor}\begin{kframe}
\begin{alltt}
\hlstd{weight} \hlkwb{<-} \hlkwa{function}\hlstd{(}\hlkwc{t}\hlstd{)} \hlkwd{exp}\hlstd{(}\hlopt{-}\hlstd{lambda.hat} \hlopt{*} \hlstd{(}\hlnum{32} \hlopt{+} \hlstd{t))}
\hlkwd{weight}\hlstd{(}\hlkwd{c}\hlstd{(}\hlnum{1}\hlopt{:}\hlnum{8}\hlstd{))}
\end{alltt}
\begin{verbatim}
## [1] 0.018103433 0.016031192 0.014196154 0.012571166 0.011132186 0.009857921
## [7] 0.008729517 0.007730277
\end{verbatim}
\end{kframe}
\end{knitrout}


%\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:lambda-weights}
   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)
\begin{knitrout}
\definecolor{shadecolor}{rgb}{0.969, 0.969, 0.969}\color{fgcolor}\begin{kframe}
\begin{alltt}
\hlstd{foo} \hlkwb{<-} \hlstd{data}\hlopt{$}\hlstd{redata}\hlopt{$}\hlstd{varb}
\hlstd{offspring} \hlkwb{<-} \hlkwd{as.numeric}\hlstd{(}\hlkwd{grepl}\hlstd{(}\hlstr{"B"}\hlstd{, foo))}
\hlstd{bar} \hlkwb{<-} \hlkwd{grepl}\hlstd{(}\hlstr{"B"}\hlstd{, foo)}
\hlstd{baz} \hlkwb{<-} \hlstd{offspring}
\hlstd{baz[bar]} \hlkwb{<-} \hlkwd{sub}\hlstd{(}\hlstr{"B"}\hlstd{,} \hlstr{""}\hlstd{, foo)[bar]}
\hlstd{baz} \hlkwb{<-} \hlkwd{as.numeric}\hlstd{(baz)}
\hlstd{baz[bar]} \hlkwb{<-} \hlkwd{weight}\hlstd{(baz[bar]} \hlopt{-} \hlnum{32}\hlstd{)}
\hlstd{offspring} \hlkwb{<-} \hlstd{baz}
\hlstd{offspring} \hlkwb{<-} \hlstd{offspring} \hlopt{*} \hlnum{100}
\end{alltt}
\end{kframe}
\end{knitrout}


A new model matrix is created that incorporates reweighted fitness.


\begin{knitrout}
\definecolor{shadecolor}{rgb}{0.969, 0.969, 0.969}\color{fgcolor}\begin{kframe}
\begin{alltt}
\hlstd{modmat.full} \hlkwb{<-} \hlkwd{model.matrix}\hlstd{(resp} \hlopt{~} \hlstd{varb} \hlopt{+}
  \hlstd{offspring}\hlopt{:}\hlstd{(Mass_Repro} \hlopt{+} \hlstd{Mass_2nd} \hlopt{+}
  \hlkwd{I}\hlstd{(Mass_Repro}\hlopt{^}\hlnum{2}\hlstd{)} \hlopt{+} \hlkwd{I}\hlstd{(Mass_2nd}\hlopt{^}\hlnum{2}\hlstd{)} \hlopt{+} \hlkwd{I}\hlstd{(Mass_2nd}\hlopt{*}\hlstd{Mass_Repro)} \hlopt{+}
  \hlstd{Time_2nd),} \hlkwc{data} \hlstd{= data}\hlopt{$}\hlstd{redata)}
\end{alltt}
\end{kframe}
\end{knitrout}


We now obtain the new aster submodel mean-value parameter $\tau$ associated with reweighted fitness. The \texttt{transformUnconditional} function is used to find the corresponding submodel canonical parameter vector $\beta$.


\begin{knitrout}
\definecolor{shadecolor}{rgb}{0.969, 0.969, 0.969}\color{fgcolor}\begin{kframe}
\begin{alltt}
\hlstd{tau.full} \hlkwb{<-} \hlkwd{crossprod}\hlstd{(modmat.full, data}\hlopt{$}\hlstd{redata}\hlopt{$}\hlstd{resp)}
\hlstd{beta.full} \hlkwb{<-} \hlkwd{transformUnconditional}\hlstd{(tau.full,}
  \hlstd{modmat.full, data,} \hlkwc{from} \hlstd{=} \hlstr{"tau"}\hlstd{,} \hlkwc{to} \hlstd{=} \hlstr{"beta"}\hlstd{)}
\end{alltt}
\end{kframe}
\end{knitrout}


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
\begin{knitrout}
\definecolor{shadecolor}{rgb}{0.969, 0.969, 0.969}\color{fgcolor}\begin{kframe}
\begin{alltt}
\hlstd{x} \hlkwb{<-} \hlkwd{seq}\hlstd{(}\hlkwc{from} \hlstd{=} \hlnum{0}\hlstd{,} \hlkwc{to} \hlstd{=} \hlnum{0.016}\hlstd{,} \hlkwc{by} \hlstd{=} \hlnum{0.016}\hlopt{/}\hlnum{100}\hlstd{)}
\hlstd{y} \hlkwb{<-} \hlkwd{seq}\hlstd{(}\hlkwc{from} \hlstd{=} \hlnum{0}\hlstd{,} \hlkwc{to} \hlstd{=} \hlnum{2.3}\hlstd{,} \hlkwc{by} \hlstd{=} \hlnum{2.3}\hlopt{/}\hlnum{100}\hlstd{)}
\hlstd{days} \hlkwb{<-} \hlnum{2}\hlopt{:}\hlnum{6} \hlopt{/} \hlnum{7}
\hlstd{mnew} \hlkwb{<-} \hlkwd{data.frame}\hlstd{(long.sex,} \hlkwd{rep}\hlstd{(days[}\hlnum{1}\hlstd{],} \hlnum{101}\hlopt{^}\hlnum{2}\hlstd{),}
    \hlkwd{rep}\hlstd{(x,} \hlkwc{times} \hlstd{=} \hlnum{101}\hlstd{),} \hlkwd{rep}\hlstd{(y,} \hlkwc{each} \hlstd{=} \hlnum{101}\hlstd{))}
\hlstd{mnew} \hlkwb{<-} \hlkwd{cbind}\hlstd{(mnew,} \hlkwd{matrix}\hlstd{(}\hlkwd{rep}\hlstd{(}\hlnum{0}\hlstd{,} \hlnum{26}\hlopt{*}\hlnum{101}\hlopt{^}\hlnum{2}\hlstd{),} \hlkwc{ncol} \hlstd{=} \hlnum{26}\hlstd{))}
\hlstd{mnames} \hlkwb{<-} \hlkwd{c}\hlstd{(}\hlstr{"Sex"}\hlstd{,} \hlstr{"Time_2nd"}\hlstd{,} \hlstr{"Mass_2nd"}\hlstd{,} \hlstr{"Mass_Repro"}\hlstd{, vars)}
\hlkwd{names}\hlstd{(mnew)} \hlkwb{<-} \hlstd{mnames}
\end{alltt}
\end{kframe}
\end{knitrout}


This model matrix is converted into an \texttt{asterdata} object in order to find estimates of fitness using the \texttt{aster2} package.


\begin{knitrout}
\definecolor{shadecolor}{rgb}{0.969, 0.969, 0.969}\color{fgcolor}\begin{kframe}
\begin{alltt}
\hlstd{data.mnew} \hlkwb{<-} \hlkwd{asterdata}\hlstd{(mnew,} \hlkwc{vars} \hlstd{= vars,} \hlkwc{pred} \hlstd{= pred,}
  \hlkwc{group} \hlstd{= group,} \hlkwc{code} \hlstd{= code,} \hlkwc{families} \hlstd{= families)}
\end{alltt}
\end{kframe}
\end{knitrout}


The population growth rate is now incorporated into the model matrix for the hypothetical individuals.


\begin{knitrout}
\definecolor{shadecolor}{rgb}{0.969, 0.969, 0.969}\color{fgcolor}\begin{kframe}
\begin{alltt}
\hlstd{foo} \hlkwb{<-} \hlkwd{as.character}\hlstd{(data.mnew}\hlopt{$}\hlstd{redata}\hlopt{$}\hlstd{varb)}
\hlstd{bar} \hlkwb{<-} \hlkwd{grepl}\hlstd{(}\hlstr{"B"}\hlstd{, foo)}
\hlstd{offspring.mnew} \hlkwb{<-} \hlkwd{as.numeric}\hlstd{(bar)}

\hlstd{baz} \hlkwb{<-} \hlstd{offspring.mnew}
\hlstd{baz[bar]} \hlkwb{<-} \hlkwd{sub}\hlstd{(}\hlstr{"B"}\hlstd{,} \hlstr{""}\hlstd{, foo)[bar]}
\hlstd{baz} \hlkwb{<-} \hlkwd{as.numeric}\hlstd{(baz)}
\hlstd{baz[bar]} \hlkwb{<-} \hlkwd{weight}\hlstd{(baz[bar]} \hlopt{-} \hlnum{32}\hlstd{)}
\hlstd{offspring.mnew} \hlkwb{<-} \hlstd{baz}
\hlstd{offspring.mnew} \hlkwb{<-} \hlnum{100} \hlopt{*} \hlstd{offspring.mnew}
\hlstd{data.mnew}\hlopt{$}\hlstd{redata} \hlkwb{<-} \hlkwd{transform}\hlstd{(data.mnew}\hlopt{$}\hlstd{redata,}
  \hlkwc{offspring.mnew} \hlstd{= offspring.mnew)}
\end{alltt}
\end{kframe}
\end{knitrout}


$\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.


\begin{knitrout}
\definecolor{shadecolor}{rgb}{0.969, 0.969, 0.969}\color{fgcolor}\begin{kframe}
\begin{alltt}
\hlstd{modmat.mnew} \hlkwb{<-} \hlkwd{model.matrix}\hlstd{(resp} \hlopt{~} \hlstd{varb} \hlopt{+}
  \hlstd{offspring.mnew}\hlopt{:}\hlstd{(Mass_Repro} \hlopt{+} \hlstd{Mass_2nd} \hlopt{+}
  \hlkwd{I}\hlstd{(Mass_Repro}\hlopt{^}\hlnum{2}\hlstd{)} \hlopt{+} \hlkwd{I}\hlstd{(Mass_2nd}\hlopt{^}\hlnum{2}\hlstd{)} \hlopt{+} \hlkwd{I}\hlstd{(Mass_2nd}\hlopt{*}\hlstd{Mass_Repro)} \hlopt{+}
  \hlstd{Time_2nd),} \hlkwc{data} \hlstd{= data.mnew}\hlopt{$}\hlstd{redata)}
\hlstd{phi} \hlkwb{<-} \hlkwd{crossprod}\hlstd{(}\hlkwd{t}\hlstd{(modmat.mnew), beta.full)}
\hlstd{mus.fit2} \hlkwb{<-} \hlkwd{transformSaturated}\hlstd{(phi, data.mnew,} \hlkwc{from} \hlstd{=} \hlstr{"phi"}\hlstd{,}
  \hlkwc{to} \hlstd{=} \hlstr{"mu"}\hlstd{)}
\end{alltt}
\end{kframe}
\end{knitrout}


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
\begin{knitrout}
\definecolor{shadecolor}{rgb}{0.969, 0.969, 0.969}\color{fgcolor}\begin{kframe}
\begin{alltt}
\hlstd{ind} \hlkwb{<-} \hlkwd{which}\hlstd{(}\hlkwd{grepl}\hlstd{(}\hlstr{"B"}\hlstd{,}\hlkwd{rownames}\hlstd{(}\hlkwd{data.frame}\hlstd{(modmat.mnew))))}
\hlstd{qux.fit2} \hlkwb{<-} \hlkwd{matrix}\hlstd{(mus.fit2[ind],} \hlkwc{nrow} \hlstd{=} \hlnum{10201}\hlstd{,} \hlkwc{ncol} \hlstd{=} \hlnum{8}\hlstd{)}
\hlstd{sumslam.fit2} \hlkwb{<-} \hlkwd{as.numeric}\hlstd{(}\hlkwd{apply}\hlstd{(qux.fit2,} \hlnum{1}\hlstd{, sum))}
\hlstd{zday.fit2} \hlkwb{<-} \hlkwd{matrix}\hlstd{(sumslam.fit2,} \hlkwc{nrow} \hlstd{=} \hlnum{101}\hlstd{)}
\end{alltt}
\end{kframe}
\end{knitrout}




\begin{figure}
\begin{center}
\begin{knitrout}
\definecolor{shadecolor}{rgb}{0.969, 0.969, 0.969}\color{fgcolor}\begin{kframe}
\begin{alltt}
\hlkwd{plot}\hlstd{(mass.second,} \hlkwc{pch} \hlstd{=} \hlnum{20}\hlstd{,}
  \hlstd{mass.reprod,} \hlkwc{xlab} \hlstd{=} \hlstr{"Mass at 
  2nd instar"}\hlstd{,} \hlkwc{ylab} \hlstd{=} \hlstr{"Mass at Reproduction"}\hlstd{,}
  \hlkwc{main} \hlstd{=} \hlstr{"Fitness Landscape for Ovariole Counts vs. Mass"}\hlstd{)}
\hlkwd{points}\hlstd{(}\hlkwd{maximum}\hlstd{(sums2)[}\hlnum{1}\hlstd{],} \hlkwd{maximum}\hlstd{(sums2)[}\hlnum{2}\hlstd{],} \hlkwc{col} \hlstd{=} \hlstr{"red"}\hlstd{,}
  \hlkwc{pch} \hlstd{=} \hlnum{19}\hlstd{)}
\hlkwd{contour}\hlstd{(x,y,zday2,} \hlkwc{add} \hlstd{=} \hlnum{TRUE}\hlstd{,} \hlkwc{nlevels} \hlstd{=} \hlnum{8}\hlstd{,} \hlkwc{levels} \hlstd{=}
  \hlkwd{c}\hlstd{(}\hlnum{100}\hlstd{,}\hlnum{150}\hlstd{,}\hlnum{200}\hlstd{,}\hlnum{250}\hlstd{,}\hlnum{300}\hlstd{,}\hlnum{325}\hlstd{,}\hlnum{350}\hlstd{,}\hlnum{360}\hlstd{) )}


\hlkwd{plot}\hlstd{(mass.second,} \hlkwc{pch} \hlstd{=} \hlnum{20}\hlstd{,}
  \hlstd{mass.reprod,} \hlkwc{xlab} \hlstd{=} \hlstr{"Mass at 
  2nd instar"}\hlstd{,} \hlkwc{ylab} \hlstd{=} \hlstr{"Mass at Reproduction"}\hlstd{,}
  \hlkwc{main} \hlstd{=} \hlstr{"Fitness Landscape for Ovariole Counts vs. Mass"}\hlstd{)}
\hlkwd{points}\hlstd{(}\hlkwd{maximum}\hlstd{(sumslam.fit2)[}\hlnum{1}\hlstd{],} \hlkwd{maximum}\hlstd{(sumslam.fit2)[}\hlnum{2}\hlstd{],}
  \hlkwc{col} \hlstd{=} \hlstr{"red"}\hlstd{,} \hlkwc{pch} \hlstd{=} \hlnum{19}\hlstd{)}
\hlkwd{contour}\hlstd{(x,y,zday.fit2,} \hlkwc{add} \hlstd{=} \hlnum{TRUE}\hlstd{,} \hlkwc{nlevels} \hlstd{=} \hlnum{8}\hlstd{)}
\end{alltt}
\end{kframe}
\includegraphics[width=6cm,height=6cm]{figure/comparison-landscapes-1} 
\includegraphics[width=6cm,height=6cm]{figure/comparison-landscapes-2} 

\end{knitrout}
\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
\begin{knitrout}
\definecolor{shadecolor}{rgb}{0.969, 0.969, 0.969}\color{fgcolor}\begin{kframe}
\begin{alltt}
\hlstd{mymu.full} \hlkwb{<-} \hlkwd{transformUnconditional}\hlstd{(beta.full,}
  \hlstd{modmat.full, data,} \hlkwc{from} \hlstd{=} \hlstr{"beta"}\hlstd{,} \hlkwc{to} \hlstd{=} \hlstr{"mu"}\hlstd{)}
\hlkwd{names}\hlstd{(mymu.full)} \hlkwb{<-} \hlkwd{rownames}\hlstd{(modmat.full)}
\hlstd{mymu.full.births} \hlkwb{<-} \hlstd{mymu.full[}\hlkwd{grepl}\hlstd{(}\hlstr{"B"}\hlstd{,} \hlkwd{names}\hlstd{(mymu.full))]}
\hlstd{mean.fitness} \hlkwb{<-} \hlkwd{mean}\hlstd{(mymu.full.births)}
\hlstd{mean.fitness} \hlkwb{<-} \hlstd{mean.fitness} \hlopt{*} \hlnum{8}
\hlstd{zday.fit2} \hlkwb{<-} \hlstd{zday.fit2} \hlopt{/} \hlstd{mean.fitness}
\end{alltt}
\end{kframe}
\end{knitrout}


The same routine is performed for the 30603 individuals that reach their second instar life stage at ages 3, 4, and 5.


\begin{knitrout}
\definecolor{shadecolor}{rgb}{0.969, 0.969, 0.969}\color{fgcolor}\begin{kframe}
\begin{alltt}
\hlstd{mnew} \hlkwb{<-} \hlkwd{data.frame}\hlstd{(long.sex,} \hlkwd{rep}\hlstd{(days[}\hlnum{2}\hlstd{],} \hlnum{101}\hlopt{^}\hlnum{2}\hlstd{),}
    \hlkwd{rep}\hlstd{(x,} \hlkwc{times} \hlstd{=} \hlnum{101}\hlstd{),} \hlkwd{rep}\hlstd{(y,} \hlkwc{each} \hlstd{=} \hlnum{101}\hlstd{))}
\hlstd{mnew} \hlkwb{<-} \hlkwd{cbind}\hlstd{(mnew,} \hlkwd{matrix}\hlstd{(}\hlkwd{rep}\hlstd{(}\hlnum{0}\hlstd{,} \hlnum{26}\hlopt{*}\hlnum{101}\hlopt{^}\hlnum{2}\hlstd{),} \hlkwc{ncol} \hlstd{=} \hlnum{26}\hlstd{))}
\hlstd{mnames} \hlkwb{<-} \hlkwd{c}\hlstd{(}\hlstr{"Sex"}\hlstd{,} \hlstr{"Time_2nd"}\hlstd{,} \hlstr{"Mass_2nd"}\hlstd{,} \hlstr{"Mass_Repro"}\hlstd{, vars)}
\hlkwd{names}\hlstd{(mnew)} \hlkwb{<-} \hlstd{mnames}
\hlstd{data.mnew} \hlkwb{<-} \hlkwd{asterdata}\hlstd{(mnew,} \hlkwc{vars} \hlstd{= vars,} \hlkwc{pred} \hlstd{= pred,}
  \hlkwc{group} \hlstd{= group,} \hlkwc{code} \hlstd{= code,} \hlkwc{families} \hlstd{= families)}
\hlstd{foo} \hlkwb{<-} \hlkwd{as.character}\hlstd{(data.mnew}\hlopt{$}\hlstd{redata}\hlopt{$}\hlstd{varb)}
\hlstd{offspring.mnew} \hlkwb{<-} \hlkwd{as.numeric}\hlstd{(bar)}

\hlstd{baz} \hlkwb{<-} \hlstd{offspring.mnew}
\hlstd{baz[bar]} \hlkwb{<-} \hlkwd{sub}\hlstd{(}\hlstr{"B"}\hlstd{,} \hlstr{""}\hlstd{, foo)[bar]}
\hlstd{baz} \hlkwb{<-} \hlkwd{as.numeric}\hlstd{(baz)}
\hlstd{baz[bar]} \hlkwb{<-} \hlkwd{weight}\hlstd{(baz[bar]} \hlopt{-} \hlnum{32}\hlstd{)}
\hlstd{offspring.mnew} \hlkwb{<-} \hlstd{baz}
\hlstd{offspring.mnew} \hlkwb{<-} \hlnum{100} \hlopt{*} \hlstd{offspring.mnew}
\hlstd{data.mnew}\hlopt{$}\hlstd{redata} \hlkwb{<-} \hlkwd{transform}\hlstd{(data.mnew}\hlopt{$}\hlstd{redata,}
  \hlkwc{offspring.mnew} \hlstd{= offspring.mnew)}

\hlstd{modmat.mnew} \hlkwb{<-} \hlkwd{model.matrix}\hlstd{(resp} \hlopt{~} \hlstd{varb} \hlopt{+}
  \hlstd{offspring.mnew}\hlopt{:}\hlstd{(Mass_Repro} \hlopt{+} \hlstd{Mass_2nd} \hlopt{+}
  \hlkwd{I}\hlstd{(Mass_Repro}\hlopt{^}\hlnum{2}\hlstd{)} \hlopt{+} \hlkwd{I}\hlstd{(Mass_2nd}\hlopt{^}\hlnum{2}\hlstd{)} \hlopt{+} \hlkwd{I}\hlstd{(Mass_2nd}\hlopt{*}\hlstd{Mass_Repro)} \hlopt{+}
  \hlstd{Time_2nd),} \hlkwc{data} \hlstd{= data.mnew}\hlopt{$}\hlstd{redata)}
\hlstd{phi} \hlkwb{<-} \hlkwd{crossprod}\hlstd{(}\hlkwd{t}\hlstd{(modmat.mnew), beta.full)}
\hlstd{mus.fit3} \hlkwb{<-} \hlkwd{transformSaturated}\hlstd{(phi, data.mnew,} \hlkwc{from} \hlstd{=} \hlstr{"phi"}\hlstd{,}
  \hlkwc{to} \hlstd{=} \hlstr{"mu"}\hlstd{)}

\hlstd{births.fit3} \hlkwb{<-} \hlstd{mus.fit3[ind]}
\hlstd{quxlam.fit3} \hlkwb{<-} \hlkwd{matrix}\hlstd{(births.fit3,} \hlkwc{nrow} \hlstd{=} \hlnum{10201}\hlstd{)}
\hlstd{sumslam.fit3} \hlkwb{<-} \hlkwd{apply}\hlstd{(quxlam.fit3,} \hlnum{1}\hlstd{, sum)}
\hlstd{zday.fit3} \hlkwb{<-} \hlkwd{matrix}\hlstd{(sumslam.fit3,} \hlkwc{nrow} \hlstd{=} \hlnum{101}\hlstd{)}
\hlstd{zday.fit3} \hlkwb{<-} \hlstd{zday.fit3} \hlopt{/} \hlstd{mean.fitness}


\hlstd{mnew} \hlkwb{<-} \hlkwd{data.frame}\hlstd{(long.sex,} \hlkwd{rep}\hlstd{(days[}\hlnum{3}\hlstd{],} \hlnum{101}\hlopt{^}\hlnum{2}\hlstd{),}
    \hlkwd{rep}\hlstd{(x,} \hlkwc{times} \hlstd{=} \hlnum{101}\hlstd{),} \hlkwd{rep}\hlstd{(y,} \hlkwc{each} \hlstd{=} \hlnum{101}\hlstd{))}
\hlstd{mnew} \hlkwb{<-} \hlkwd{cbind}\hlstd{(mnew,} \hlkwd{matrix}\hlstd{(}\hlkwd{rep}\hlstd{(}\hlnum{0}\hlstd{,} \hlnum{26}\hlopt{*}\hlnum{101}\hlopt{^}\hlnum{2}\hlstd{),} \hlkwc{ncol} \hlstd{=} \hlnum{26}\hlstd{))}
\hlstd{mnames} \hlkwb{<-} \hlkwd{c}\hlstd{(}\hlstr{"Sex"}\hlstd{,} \hlstr{"Time_2nd"}\hlstd{,} \hlstr{"Mass_2nd"}\hlstd{,} \hlstr{"Mass_Repro"}\hlstd{, vars)}
\hlkwd{names}\hlstd{(mnew)} \hlkwb{<-} \hlstd{mnames}
\hlstd{data.mnew} \hlkwb{<-} \hlkwd{asterdata}\hlstd{(mnew,} \hlkwc{vars} \hlstd{= vars,} \hlkwc{pred} \hlstd{= pred,}
  \hlkwc{group} \hlstd{= group,} \hlkwc{code} \hlstd{= code,} \hlkwc{families} \hlstd{= families)}
\hlstd{foo} \hlkwb{<-} \hlkwd{as.character}\hlstd{(data.mnew}\hlopt{$}\hlstd{redata}\hlopt{$}\hlstd{varb)}
\hlstd{offspring.mnew} \hlkwb{<-} \hlkwd{as.numeric}\hlstd{(bar)}

\hlstd{baz} \hlkwb{<-} \hlstd{offspring.mnew}
\hlstd{baz[bar]} \hlkwb{<-} \hlkwd{sub}\hlstd{(}\hlstr{"B"}\hlstd{,} \hlstr{""}\hlstd{, foo)[bar]}
\hlstd{baz} \hlkwb{<-} \hlkwd{as.numeric}\hlstd{(baz)}
\hlstd{baz[bar]} \hlkwb{<-} \hlkwd{weight}\hlstd{(baz[bar]} \hlopt{-} \hlnum{32}\hlstd{)}
\hlstd{offspring.mnew} \hlkwb{<-} \hlstd{baz}
\hlstd{offspring.mnew} \hlkwb{<-} \hlnum{100} \hlopt{*} \hlstd{offspring.mnew}
\hlstd{data.mnew}\hlopt{$}\hlstd{redata} \hlkwb{<-} \hlkwd{transform}\hlstd{(data.mnew}\hlopt{$}\hlstd{redata,}
  \hlkwc{offspring.mnew} \hlstd{= offspring.mnew)}

\hlstd{modmat.mnew} \hlkwb{<-} \hlkwd{model.matrix}\hlstd{(resp} \hlopt{~} \hlstd{varb} \hlopt{+}
  \hlstd{offspring.mnew}\hlopt{:}\hlstd{(Mass_Repro} \hlopt{+} \hlstd{Mass_2nd} \hlopt{+}
  \hlkwd{I}\hlstd{(Mass_Repro}\hlopt{^}\hlnum{2}\hlstd{)} \hlopt{+} \hlkwd{I}\hlstd{(Mass_2nd}\hlopt{^}\hlnum{2}\hlstd{)} \hlopt{+} \hlkwd{I}\hlstd{(Mass_2nd}\hlopt{*}\hlstd{Mass_Repro)} \hlopt{+}
  \hlstd{Time_2nd),} \hlkwc{data} \hlstd{= data.mnew}\hlopt{$}\hlstd{redata)}
\hlstd{phi} \hlkwb{<-} \hlkwd{crossprod}\hlstd{(}\hlkwd{t}\hlstd{(modmat.mnew), beta.full)}
\hlstd{mus.fit4} \hlkwb{<-} \hlkwd{transformSaturated}\hlstd{(phi, data.mnew,} \hlkwc{from} \hlstd{=} \hlstr{"phi"}\hlstd{,}
  \hlkwc{to} \hlstd{=} \hlstr{"mu"}\hlstd{)}

\hlstd{births.fit4} \hlkwb{<-} \hlstd{mus.fit4[ind]}
\hlstd{quxlam.fit4} \hlkwb{<-} \hlkwd{matrix}\hlstd{(births.fit4,} \hlkwc{nrow} \hlstd{=} \hlnum{10201}\hlstd{)}
\hlstd{sumslam.fit4} \hlkwb{<-} \hlkwd{apply}\hlstd{(quxlam.fit4,} \hlnum{1}\hlstd{, sum)}
\hlstd{zday.fit4} \hlkwb{<-} \hlkwd{matrix}\hlstd{(sumslam.fit4,} \hlkwc{nrow} \hlstd{=} \hlnum{101}\hlstd{)}
\hlstd{zday.fit4} \hlkwb{<-} \hlstd{zday.fit4} \hlopt{/} \hlstd{mean.fitness}


\hlstd{mnew} \hlkwb{<-} \hlkwd{data.frame}\hlstd{(long.sex,} \hlkwd{rep}\hlstd{(days[}\hlnum{4}\hlstd{],} \hlnum{101}\hlopt{^}\hlnum{2}\hlstd{),}
    \hlkwd{rep}\hlstd{(x,} \hlkwc{times} \hlstd{=} \hlnum{101}\hlstd{),} \hlkwd{rep}\hlstd{(y,} \hlkwc{each} \hlstd{=} \hlnum{101}\hlstd{))}
\hlstd{mnew} \hlkwb{<-} \hlkwd{cbind}\hlstd{(mnew,} \hlkwd{matrix}\hlstd{(}\hlkwd{rep}\hlstd{(}\hlnum{0}\hlstd{,} \hlnum{26}\hlopt{*}\hlnum{101}\hlopt{^}\hlnum{2}\hlstd{),} \hlkwc{ncol} \hlstd{=} \hlnum{26}\hlstd{))}
\hlstd{mnames} \hlkwb{<-} \hlkwd{c}\hlstd{(}\hlstr{"Sex"}\hlstd{,} \hlstr{"Time_2nd"}\hlstd{,} \hlstr{"Mass_2nd"}\hlstd{,} \hlstr{"Mass_Repro"}\hlstd{, vars)}
\hlkwd{names}\hlstd{(mnew)} \hlkwb{<-} \hlstd{mnames}
\hlstd{data.mnew} \hlkwb{<-} \hlkwd{asterdata}\hlstd{(mnew,} \hlkwc{vars} \hlstd{= vars,} \hlkwc{pred} \hlstd{= pred,}
  \hlkwc{group} \hlstd{= group,} \hlkwc{code} \hlstd{= code,} \hlkwc{families} \hlstd{= families)}
\hlstd{foo} \hlkwb{<-} \hlkwd{as.character}\hlstd{(data.mnew}\hlopt{$}\hlstd{redata}\hlopt{$}\hlstd{varb)}
\hlstd{offspring.mnew} \hlkwb{<-} \hlkwd{as.numeric}\hlstd{(bar)}

\hlstd{baz} \hlkwb{<-} \hlstd{offspring.mnew}
\hlstd{baz[bar]} \hlkwb{<-} \hlkwd{sub}\hlstd{(}\hlstr{"B"}\hlstd{,} \hlstr{""}\hlstd{, foo)[bar]}
\hlstd{baz} \hlkwb{<-} \hlkwd{as.numeric}\hlstd{(baz)}
\hlstd{baz[bar]} \hlkwb{<-} \hlkwd{weight}\hlstd{(baz[bar]} \hlopt{-} \hlnum{32}\hlstd{)}
\hlstd{offspring.mnew} \hlkwb{<-} \hlstd{baz}
\hlstd{offspring.mnew} \hlkwb{<-} \hlnum{100} \hlopt{*} \hlstd{offspring.mnew}
\hlstd{data.mnew}\hlopt{$}\hlstd{redata} \hlkwb{<-} \hlkwd{transform}\hlstd{(data.mnew}\hlopt{$}\hlstd{redata,}
  \hlkwc{offspring.mnew} \hlstd{= offspring.mnew)}

\hlstd{modmat.mnew} \hlkwb{<-} \hlkwd{model.matrix}\hlstd{(resp} \hlopt{~} \hlstd{varb} \hlopt{+}
  \hlstd{offspring.mnew}\hlopt{:}\hlstd{(Mass_Repro} \hlopt{+} \hlstd{Mass_2nd} \hlopt{+}
  \hlkwd{I}\hlstd{(Mass_Repro}\hlopt{^}\hlnum{2}\hlstd{)} \hlopt{+} \hlkwd{I}\hlstd{(Mass_2nd}\hlopt{^}\hlnum{2}\hlstd{)} \hlopt{+} \hlkwd{I}\hlstd{(Mass_2nd}\hlopt{*}\hlstd{Mass_Repro)} \hlopt{+}
  \hlstd{Time_2nd),} \hlkwc{data} \hlstd{= data.mnew}\hlopt{$}\hlstd{redata)}
\hlstd{phi} \hlkwb{<-} \hlkwd{crossprod}\hlstd{(}\hlkwd{t}\hlstd{(modmat.mnew), beta.full)}
\hlstd{mus.fit5} \hlkwb{<-} \hlkwd{transformSaturated}\hlstd{(phi, data.mnew,} \hlkwc{from} \hlstd{=} \hlstr{"phi"}\hlstd{,}
  \hlkwc{to} \hlstd{=} \hlstr{"mu"}\hlstd{)}

\hlstd{births.fit5} \hlkwb{<-} \hlstd{mus.fit5[ind]}
\hlstd{quxlam.fit5} \hlkwb{<-} \hlkwd{matrix}\hlstd{(births.fit5,} \hlkwc{nrow} \hlstd{=} \hlnum{10201}\hlstd{)}
\hlstd{sumslam.fit5} \hlkwb{<-} \hlkwd{apply}\hlstd{(quxlam.fit5,} \hlnum{1}\hlstd{, sum)}
\hlstd{zday.fit5} \hlkwb{<-} \hlkwd{matrix}\hlstd{(sumslam.fit5,} \hlkwc{nrow} \hlstd{=} \hlnum{101}\hlstd{)}
\hlstd{zday.fit5} \hlkwb{<-} \hlstd{zday.fit5} \hlopt{/} \hlstd{mean.fitness}
\end{alltt}
\end{kframe}
\end{knitrout}



%We now make the four-panelled fitness landscape where fitness is reweighted to include the population growth rate. 


%<<eval = 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(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}
%<<fitnesslandscape-reweighted, 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(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 four-panelled 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.





%<<eval=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)
%@






%\begin{figure}
%\begin{center}
%<<4by2-part1,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
%round(max(zday2),1) (top left),
%round(max(zday3),1) (bottom left),
%round(max(zday.fit2),2) (top right),
%and round(max(zday.fit3),2) (bottom-right). 
%}
%\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) %%%%%%%%%%%%
%<<eval = FALSE>>=
%# 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}
%<<4by2-part2,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
%round(max(zday4),1) (top left),
%round(max(zday5),1) (bottom left),
%round(max(zday.fit4),2) (top right),
%and round(max(zday.fit5),2) (bottom-right). 
%}
%\label{fig:bar}
%\end{figure}


We now make the eight-paneled 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.


\begin{knitrout}
\definecolor{shadecolor}{rgb}{0.969, 0.969, 0.969}\color{fgcolor}\begin{kframe}
\begin{alltt}
\hlstd{xlim} \hlkwb{<-} \hlkwd{range}\hlstd{(x)}
\hlstd{ylim} \hlkwb{<-} \hlkwd{range}\hlstd{(y)}

\hlkwd{par}\hlstd{(}\hlkwc{mfrow} \hlstd{=} \hlkwd{c}\hlstd{(}\hlnum{4}\hlstd{,} \hlnum{2}\hlstd{),} \hlkwc{mar} \hlstd{=} \hlkwd{rep}\hlstd{(}\hlnum{0.5}\hlstd{,} \hlnum{4}\hlstd{),}
  \hlkwc{oma} \hlstd{=} \hlkwd{c}\hlstd{(}\hlnum{4}\hlstd{,} \hlnum{4}\hlstd{,} \hlnum{0}\hlstd{,} \hlnum{4}\hlstd{)} \hlopt{+} \hlnum{0.1}\hlstd{)}
\hlcom{# day 2}
\hlkwd{plot}\hlstd{(mass.second,}
    \hlstd{mass.reprod,}
    \hlkwc{axes} \hlstd{=} \hlnum{FALSE}\hlstd{,} \hlkwc{xlab} \hlstd{=} \hlstr{""}\hlstd{,} \hlkwc{ylab} \hlstd{=} \hlstr{""}\hlstd{,} \hlkwc{pch} \hlstd{=} \hlnum{20}\hlstd{,}
    \hlkwc{xlim} \hlstd{= xlim,} \hlkwc{ylim} \hlstd{= ylim)}
\hlkwd{box}\hlstd{()}
\hlkwd{axis}\hlstd{(}\hlkwc{side} \hlstd{=} \hlnum{2}\hlstd{,} \hlkwc{outer} \hlstd{=} \hlnum{TRUE}\hlstd{)}
\hlkwd{mtext}\hlstd{(}\hlstr{"mass at eclosion"}\hlstd{,} \hlkwc{outer} \hlstd{=} \hlnum{TRUE}\hlstd{,} \hlkwc{line} \hlstd{=} \hlnum{3}\hlstd{,} \hlkwc{side} \hlstd{=} \hlnum{2}\hlstd{,}
  \hlkwc{at} \hlstd{=} \hlnum{0.50}\hlstd{)}
\hlkwd{mtext}\hlstd{(}\hlstr{"mass at 2nd instar stage"}\hlstd{,} \hlkwc{outer} \hlstd{=} \hlnum{TRUE}\hlstd{,} \hlkwc{line} \hlstd{=} \hlnum{3}\hlstd{,}
  \hlkwc{side} \hlstd{=} \hlnum{1}\hlstd{,} \hlkwc{at} \hlstd{=} \hlnum{0.50}\hlstd{)}
\hlkwd{mtext}\hlstd{(}\hlstr{"age second instar stage reached"}\hlstd{,} \hlkwc{outer} \hlstd{=} \hlnum{TRUE}\hlstd{,}
  \hlkwc{line} \hlstd{=} \hlnum{1}\hlstd{,} \hlkwc{side} \hlstd{=} \hlnum{4}\hlstd{,} \hlkwc{at} \hlstd{=} \hlnum{0.50}\hlstd{)}
\hlkwd{mtext}\hlstd{(}\hlkwd{as.character}\hlstd{(}\hlnum{2}\hlopt{:}\hlnum{5}\hlstd{),} \hlkwc{outer} \hlstd{=} \hlnum{TRUE}\hlstd{,} \hlkwc{line} \hlstd{=} \hlnum{3}\hlstd{,}
  \hlkwc{side} \hlstd{=} \hlnum{4}\hlstd{,} \hlkwc{at} \hlstd{=} \hlkwd{c}\hlstd{(}\hlnum{7}\hlstd{,} \hlnum{5}\hlstd{,} \hlnum{3}\hlstd{,} \hlnum{1}\hlstd{)} \hlopt{/} \hlnum{8}\hlstd{)}
\hlkwd{points}\hlstd{(}\hlkwd{maximum}\hlstd{(sums2)[}\hlnum{1}\hlstd{],} \hlkwd{maximum}\hlstd{(sums2)[}\hlnum{2}\hlstd{],} \hlkwc{pch} \hlstd{=} \hlnum{0}\hlstd{)}
\hlstd{levels} \hlkwb{<-} \hlkwd{seq}\hlstd{(}\hlnum{50}\hlstd{,} \hlnum{350}\hlstd{,} \hlkwc{by} \hlstd{=} \hlnum{50}\hlstd{)}
\hlkwd{contour}\hlstd{(x, y, zday2,} \hlkwc{add} \hlstd{=} \hlnum{TRUE}\hlstd{,} \hlkwc{levels} \hlstd{= levels)}


\hlkwd{plot}\hlstd{(mass.second,}
    \hlstd{mass.reprod,}
    \hlkwc{axes} \hlstd{=} \hlnum{FALSE}\hlstd{,} \hlkwc{xlab} \hlstd{=} \hlstr{""}\hlstd{,} \hlkwc{ylab} \hlstd{=} \hlstr{""}\hlstd{,} \hlkwc{pch} \hlstd{=} \hlnum{20}\hlstd{,}
    \hlkwc{xlim} \hlstd{= xlim,} \hlkwc{ylim} \hlstd{= ylim)}
\hlkwd{box}\hlstd{()}
\hlkwd{points}\hlstd{(}\hlkwd{maximum}\hlstd{(sumslam.fit2)[}\hlnum{1}\hlstd{],} \hlkwd{maximum}\hlstd{(sumslam.fit2)[}\hlnum{2}\hlstd{],}
  \hlkwc{pch} \hlstd{=} \hlnum{0}\hlstd{)}
\hlkwd{contour}\hlstd{(x, y, zday.fit2,} \hlkwc{nlevels} \hlstd{=} \hlnum{7}\hlstd{,} \hlkwc{add} \hlstd{=} \hlnum{TRUE}\hlstd{)}


\hlcom{# day 3}
\hlkwd{plot}\hlstd{(mass.second,}
    \hlstd{mass.reprod,}
    \hlkwc{axes} \hlstd{=} \hlnum{FALSE}\hlstd{,} \hlkwc{xlab} \hlstd{=} \hlstr{""}\hlstd{,} \hlkwc{ylab} \hlstd{=} \hlstr{""}\hlstd{,} \hlkwc{pch} \hlstd{=} \hlnum{20}\hlstd{,}
    \hlkwc{xlim} \hlstd{= xlim,} \hlkwc{ylim} \hlstd{= ylim)}
\hlkwd{box}\hlstd{()}
\hlkwd{axis}\hlstd{(}\hlkwc{side} \hlstd{=} \hlnum{2}\hlstd{,} \hlkwc{outer} \hlstd{=} \hlnum{TRUE}\hlstd{)}
\hlkwd{mtext}\hlstd{(}\hlstr{"mass at eclosion"}\hlstd{,} \hlkwc{outer} \hlstd{=} \hlnum{TRUE}\hlstd{,} \hlkwc{line} \hlstd{=} \hlnum{3}\hlstd{,}
  \hlkwc{side} \hlstd{=} \hlnum{2}\hlstd{,} \hlkwc{at} \hlstd{=} \hlnum{0.50}\hlstd{)}
\hlkwd{points}\hlstd{(}\hlkwd{maximum}\hlstd{(sums3)[}\hlnum{1}\hlstd{],} \hlkwd{maximum}\hlstd{(sums3)[}\hlnum{2}\hlstd{],} \hlkwc{pch} \hlstd{=} \hlnum{0}\hlstd{)}
\hlkwd{contour}\hlstd{(x, y, zday3,} \hlkwc{add} \hlstd{=} \hlnum{TRUE}\hlstd{,} \hlkwc{levels} \hlstd{= levels)}


\hlkwd{plot}\hlstd{(mass.second,}
    \hlstd{mass.reprod,}
    \hlkwc{axes} \hlstd{=} \hlnum{FALSE}\hlstd{,} \hlkwc{xlab} \hlstd{=} \hlstr{""}\hlstd{,} \hlkwc{ylab} \hlstd{=} \hlstr{""}\hlstd{,} \hlkwc{pch} \hlstd{=} \hlnum{20}\hlstd{,}
    \hlkwc{xlim} \hlstd{= xlim,} \hlkwc{ylim} \hlstd{= ylim)}
\hlkwd{box}\hlstd{()}
\hlkwd{points}\hlstd{(}\hlkwd{maximum}\hlstd{(sumslam.fit3)[}\hlnum{1}\hlstd{],} \hlkwd{maximum}\hlstd{(sumslam.fit3)[}\hlnum{2}\hlstd{],}
  \hlkwc{pch} \hlstd{=} \hlnum{0}\hlstd{)}
\hlkwd{contour}\hlstd{(x, y, zday.fit3,} \hlkwc{nlevels} \hlstd{=} \hlnum{7}\hlstd{,} \hlkwc{add} \hlstd{=} \hlnum{TRUE}\hlstd{)}


\hlcom{# day 4}
\hlkwd{plot}\hlstd{(mass.second,}
    \hlstd{mass.reprod,}
    \hlkwc{axes} \hlstd{=} \hlnum{FALSE}\hlstd{,} \hlkwc{xlab} \hlstd{=} \hlstr{""}\hlstd{,} \hlkwc{ylab} \hlstd{=} \hlstr{""}\hlstd{,} \hlkwc{pch} \hlstd{=} \hlnum{20}\hlstd{,}
    \hlkwc{xlim} \hlstd{= xlim,} \hlkwc{ylim} \hlstd{= ylim)}
\hlkwd{box}\hlstd{()}
\hlkwd{axis}\hlstd{(}\hlkwc{side} \hlstd{=} \hlnum{2}\hlstd{,} \hlkwc{outer} \hlstd{=} \hlnum{TRUE}\hlstd{)}
\hlkwd{mtext}\hlstd{(}\hlstr{"mass at eclosion"}\hlstd{,} \hlkwc{outer} \hlstd{=} \hlnum{TRUE}\hlstd{,} \hlkwc{line} \hlstd{=} \hlnum{3}\hlstd{,}
  \hlkwc{side} \hlstd{=} \hlnum{2}\hlstd{,} \hlkwc{at} \hlstd{=} \hlnum{0.50}\hlstd{)}
\hlkwd{points}\hlstd{(}\hlkwd{maximum}\hlstd{(sums4)[}\hlnum{1}\hlstd{],} \hlkwd{maximum}\hlstd{(sums4)[}\hlnum{2}\hlstd{],} \hlkwc{pch} \hlstd{=} \hlnum{0}\hlstd{)}
\hlkwd{contour}\hlstd{(x, y, zday4,} \hlkwc{add} \hlstd{=} \hlnum{TRUE}\hlstd{,} \hlkwc{levels} \hlstd{= levels)}


\hlkwd{plot}\hlstd{(mass.second,}
    \hlstd{mass.reprod,}
    \hlkwc{axes} \hlstd{=} \hlnum{FALSE}\hlstd{,} \hlkwc{xlab} \hlstd{=} \hlstr{""}\hlstd{,} \hlkwc{ylab} \hlstd{=} \hlstr{""}\hlstd{,} \hlkwc{pch} \hlstd{=} \hlnum{20}\hlstd{,}
    \hlkwc{xlim} \hlstd{= xlim,} \hlkwc{ylim} \hlstd{= ylim)}
\hlkwd{box}\hlstd{()}
\hlkwd{points}\hlstd{(}\hlkwd{maximum}\hlstd{(sumslam.fit4)[}\hlnum{1}\hlstd{],} \hlkwd{maximum}\hlstd{(sumslam.fit4)[}\hlnum{2}\hlstd{],}
  \hlkwc{pch} \hlstd{=} \hlnum{0}\hlstd{)}
\hlkwd{contour}\hlstd{(x, y, zday.fit4,} \hlkwc{nlevels} \hlstd{=} \hlnum{7}\hlstd{,} \hlkwc{add} \hlstd{=} \hlnum{TRUE}\hlstd{)}


\hlcom{# day 5}
\hlkwd{plot}\hlstd{(mass.second,}
    \hlstd{mass.reprod,}
    \hlkwc{axes} \hlstd{=} \hlnum{FALSE}\hlstd{,} \hlkwc{xlab} \hlstd{=} \hlstr{""}\hlstd{,} \hlkwc{ylab} \hlstd{=} \hlstr{""}\hlstd{,} \hlkwc{pch} \hlstd{=} \hlnum{20}\hlstd{,}
    \hlkwc{xlim} \hlstd{= xlim,} \hlkwc{ylim} \hlstd{= ylim)}
\hlkwd{box}\hlstd{()}
\hlkwd{axis}\hlstd{(}\hlkwc{side} \hlstd{=} \hlnum{1}\hlstd{,} \hlkwc{outer} \hlstd{=} \hlnum{TRUE}\hlstd{)}
\hlkwd{axis}\hlstd{(}\hlkwc{side} \hlstd{=} \hlnum{2}\hlstd{,} \hlkwc{outer} \hlstd{=} \hlnum{TRUE}\hlstd{)}
\hlkwd{points}\hlstd{(}\hlkwd{maximum}\hlstd{(sums5)[}\hlnum{1}\hlstd{],} \hlkwd{maximum}\hlstd{(sums5)[}\hlnum{2}\hlstd{],} \hlkwc{pch} \hlstd{=} \hlnum{0}\hlstd{)}
\hlkwd{contour}\hlstd{(x, y, zday5,} \hlkwc{add} \hlstd{=} \hlnum{TRUE}\hlstd{,} \hlkwc{levels} \hlstd{= levels)}


\hlkwd{plot}\hlstd{(mass.second,}
    \hlstd{mass.reprod,}
    \hlkwc{axes} \hlstd{=} \hlnum{FALSE}\hlstd{,} \hlkwc{xlab} \hlstd{=} \hlstr{""}\hlstd{,} \hlkwc{ylab} \hlstd{=} \hlstr{""}\hlstd{,} \hlkwc{pch} \hlstd{=} \hlnum{20}\hlstd{,}
    \hlkwc{xlim} \hlstd{= xlim,} \hlkwc{ylim} \hlstd{= ylim)}
\hlkwd{box}\hlstd{()}
\hlkwd{axis}\hlstd{(}\hlkwc{side} \hlstd{=} \hlnum{1}\hlstd{,} \hlkwc{outer} \hlstd{=} \hlnum{TRUE}\hlstd{)}
\hlkwd{points}\hlstd{(}\hlkwd{maximum}\hlstd{(sumslam.fit5)[}\hlnum{1}\hlstd{],} \hlkwd{maximum}\hlstd{(sumslam.fit5)[}\hlnum{2}\hlstd{],}
  \hlkwc{pch} \hlstd{=} \hlnum{0}\hlstd{)}
\hlkwd{contour}\hlstd{(x, y, zday.fit5,} \hlkwc{nlevels} \hlstd{=} \hlnum{7}\hlstd{,} \hlkwc{add} \hlstd{=} \hlnum{TRUE}\hlstd{)}
\end{alltt}
\end{kframe}
\end{knitrout}










\begin{figure}[p]
\begin{center}
%%%%% DON'T REPEAT YOURSELF (the DRY rule)
%%%%% Quote chunks, don't repeat chunks.
\begin{knitrout}
\definecolor{shadecolor}{rgb}{0.969, 0.969, 0.969}\color{fgcolor}
\includegraphics[width=\maxwidth]{figure/new4by2-1} 

\end{knitrout}
\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)
363.6,
342.4,
322.4,
303.6,
(right column, from top to bottom)
1.39,
1.28,
1.20,
1.14.}
% 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 %%%%%
%<<eval = 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,
%  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}
%<<4by2-tech, 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[((x-1)*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[((x-1)*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[((x-1)*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[((x-1)*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:stable-reweight} 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 = 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}
%<<label=fig4, 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(sumslam2)[1], maximum(sumslam2)[2], pch = 0)
%contour(x, y, zlam2, nlevels = 7, add = TRUE)
%# day 3
%plot(mass.second, 
%  mass.reprod,
%  axes = FALSE, xlab = "", ylab = "", pch = 20)
%box()
%points(maximum(sumslam3)[1], maximum(sumslam3)[2], pch = 0)
%contour(x, y, zlam3, 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(sumslam4)[1], maximum(sumslam4)[2], pch = 0)
%contour(x, y, zlam4, 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(sumslam5)[1], maximum(sumslam5)[2], pch = 0)
%contour(x, y, zlam5, nlevels = 7, add = TRUE)
%@
%\end{center}
%\caption{Relative fitness landscapes accounting for $\lambda$.
%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. 
%The maximum values are round(max(zlam2),2) (top left),
%round(max(zlam3),2) (top right),
%round(max(zlam4),2) (bottom left),
%and round(max(zlam5),2) (bottom right).}
%\label{fig:four}
%\end{figure}



%\subsection{Fitness landscape weighted with respect to a hypothetically declining population}

%Changing the mean population growth rate of the population, $\lambda$, 
%%%and reweighting the fitness landscapes according to the new value 
%can, in principle, substantially change the landscape. Here we consider an alternate value of $\lambda$ to be $\ensuremath{-0.122}$ which is the negative of the estimate of the actual population growth rate. This new value corresponds to a population that is decreasing in size. The new fitness landscapes for this hypothetical, declining population (Fig.~\ref{fig:five}) differ in scale but not in the shape of landscapes (cf.\ Fig.~\ref{fig:two} and Fig.~\ref{fig:four}). 
%For these data, mean population growth rate has little effect on the patterns of phenotypic selection on these traits. \\

%To construct Fig.~\ref{fig:five} we first obtain the new $\lambda$ value corresponding to the hypothetically declining population. Estimated mean expected fitness is obtained the same way as before. 


%<<>>=
%lambda.hat2 <- - lambda.hat 
%myweights <- exp(- lambda.hat2 * mydays)
%mean.fitness.neglam <- mean(mymu.full.births * myweights) * 8
%@


%The weights with respect to this new population growth rate are


%<<>>=
%weight.neg <- function(t) exp(-lambda.hat2 * (32 + t))
%weight.neg(c(1:8))
%@


%We perform the same calculations as before to produce the relative fitness landscapes reweighted with respect to the population growth rate for this hypothetically declining population. The 40804 hypothetical individuals used are the same individuals that produced the previous fitness landscapes.


%<<>>=
%neglam2 <- function(x) weight.neg(x) * births2[((x-1)*10201+1):(x*10201)]
%quxneglam2 <- unlist(apply( matrix(c(1:8)),1, FUN = neglam2))
%sumsneglam2 <- apply(quxneglam2, 1, sum)
%zneglam2 <- matrix(sumsneglam2, nrow = 101)

%# compute relative fitness
%zneglam2 <- zneglam2 / mean.fitness.neglam

%neglam3 <- function(x) weight.neg(x) * births3[((x-1)*10201+1):(x*10201)]
%quxneglam3 <- unlist(apply( matrix(c(1:8)),1, FUN = neglam3))
%sumsneglam3 <- apply(quxneglam3, 1, sum)
%zneglam3 <- matrix(sumsneglam3, nrow = 101)

%# compute relative fitness
%zneglam3 <- zneglam3 / mean.fitness.neglam

%neglam4 <- function(x) weight.neg(x) * births4[((x-1)*10201+1):(x*10201)]
%quxneglam4 <- unlist(apply( matrix(c(1:8)),1, FUN = neglam4))
%sumsneglam4 <- apply(quxneglam4, 1, sum)
%zneglam4 <- matrix(sumsneglam4, nrow = 101)

%# compute relative fitness
%zneglam4 <- zneglam4 / mean.fitness.neglam

%neglam5 <- function(x) weight.neg(x) * births5[((x-1)*10201+1):(x*10201)]
%quxneglam5 <- unlist(apply( matrix(c(1:8)),1, FUN = neglam5))
%sumsneglam5 <- apply(quxneglam5, 1, sum)
%zneglam5 <- matrix(sumsneglam5, nrow = 101)

%# compute relative fitness
%zneglam5 <- zneglam5 / mean.fitness.neglam
%@


%\begin{figure}[h!]
%\begin{center}
%<<label=neglam, echo = FALSE>>=
%par(mfrow = c(2, 2), mar = rep(0.5, 4), oma = c(4, 4, 0, 0) + 0.1)

%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(sumsneglam2)[1], maximum(sumsneglam2)[2], pch = 0)
%contour(x,y,zneglam2, add = TRUE)

%plot(mass.second, 
%  mass.reprod, 
%  axes = FALSE, xlab = "", ylab = "", pch = 20)
%box()
%points(maximum(sumsneglam3)[1], maximum(sumsneglam3)[2], pch = 0)
%contour(x,y,zneglam3, add = TRUE)

%plot(mass.second, 
%  mass.reprod, 
%  axes = FALSE, xlab = "", ylab = "", pch = 20)
%box()
%axis(side = 2, outer = TRUE)
%axis(side = 1, outer = TRUE)
%points(maximum(sumsneglam4)[1], maximum(sumsneglam4)[2], pch = 0)
%contour(x,y,zneglam4, add = TRUE)

%plot(mass.second, 
%  mass.reprod,
%  axes = FALSE, xlab = "", ylab = "", pch = 20)
%box()
%axis(side = 1, outer = TRUE)
%points(maximum(sumsneglam5)[1], maximum(sumsneglam5)[2], pch = 0)
%contour(x,y,zneglam5, add = TRUE)
%@
%\end{center}
%\caption{fitness landscapes accounting for $\lambda$ of
%a hypothetically declining population.
%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. 
%The maximum values are
%1.5 (top-left),
%1.42 (top-right),
%1.33 (bottom-left),
%and 1.25 (bottom-right).}
%\label{fig:five}
%\end{figure}




%\begin{figure}[h!]
%\begin{center}
%<<out.width='6cm', out.height='6cm',fig.show='hold'>>=
%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(sumslam2)[1], maximum(sumslam2)[2], pch = 0)
%contour(x, y, zlam2, nlevels = 7, add = TRUE)

%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)
%@
%\end{center}
%\caption{The left panel is the fitness landscape where fitness is reweighted according to the population growth rate. The right panel is the fitness landscape with no reweighting. The contours appear to be no different.}
%\label{fig:six}
%\end{figure}

































%\newpage
%Redo the analysis with fitness reweighted to incorporate the population growth rate.


%<<>>=
%foo <- data$redata$varb
%bar <- grepl("B",data$redata$varb)
%offspring <- as.numeric(bar)
%baz <- offspring
%baz[bar] <- sub("B", "", foo)[bar]
%baz <- as.numeric(baz)
%baz[bar] <- weight(baz[bar] - 32)
%offspring <- baz
%data$redata <- transform(data$redata, 
%  offspring = offspring)

%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")
%@


%<<>>=
%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 mean-value parameter from the linear model in both mass traits is now calculated. We transform our MLE for the submodel canonical parameter to the mean-value 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)
%@

%<<eval = FALSE, echo = FALSE>>=
%foot <- crossprod(modmat.quad, data$redata$resp)
%foot[c(29:31)] <- 0
%blah <- transformUnconditional(foot, tolerance = 1e-30,
%  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 p-value where the reference distribution for this test is $\chi^2$ with three degrees of freedom. The p-value 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 = 1e-50) %*% 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 p-value 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 = 1e-50) %*% score
%Rao

%df <- length(beta.full) - length(beta.quad)
%pval <- pchisq(Rao, df = df, lower.tail = FALSE)
%pval
%@



%<<echo = 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)
%@


%All more 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 ovariole count (expected Darwinian fitness). For example, the model including the interaction between mass at reproduction and age to the second instar stage is seen to be not significant.



%<<>>=
%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 = 1e-50) %*% score
%Rao

%df <- length(beta.bigger) - length(beta.full)
%pval <- pchisq(Rao, df = df, lower.tail = FALSE)
%pval
%@






\appendix

\section{Make the Data}

In this appendix we show how the \texttt{hornworm} dataset in the \astertwo
package was initially created starting from the data for the original analysis
\citep{kingsolver-et-al}, which was a comma separated values (CSV) file.
This file can be found at the same location at the University of Minnesota
Digital Conservancy (\url{http://conservancy.umn.edu/}) where this technical
report is found.
\begin{knitrout}
\definecolor{shadecolor}{rgb}{0.969, 0.969, 0.969}\color{fgcolor}\begin{kframe}
\begin{alltt}
\hlstd{MSexta} \hlkwb{<-} \hlkwd{read.csv}\hlstd{(}\hlstr{"MSexta.Aster.csv"}\hlstd{)}
\end{alltt}


{\ttfamily\noindent\color{warningcolor}{\#\# Warning in file(file, "{}rt"{}): cannot open file 'MSexta.Aster.csv': No such file or directory}}

{\ttfamily\noindent\bfseries\color{errorcolor}{\#\# Error in file(file, "{}rt"{}): cannot open the connection}}\begin{alltt}
\hlkwd{names}\hlstd{(MSexta)}
\end{alltt}


{\ttfamily\noindent\bfseries\color{errorcolor}{\#\# Error in eval(expr, envir, enclos): object 'MSexta' not found}}\end{kframe}
\end{knitrout}

Some of the column names are changed for convenience.
Any column name that ends in \texttt{Repro} means that the data in that column
pertains to the eclosion stage of \emph{M. sexta}.

\begin{knitrout}
\definecolor{shadecolor}{rgb}{0.969, 0.969, 0.969}\color{fgcolor}\begin{kframe}
\begin{alltt}
\hlkwd{names}\hlstd{(MSexta)[}\hlkwd{c}\hlstd{(}\hlnum{3}\hlstd{,}\hlnum{6}\hlstd{,}\hlnum{7}\hlstd{,}\hlnum{8}\hlstd{)]} \hlkwb{<-} \hlkwd{c}\hlstd{(}\hlstr{"Surv_Repro"}\hlstd{,}\hlstr{"Time_Repro"}\hlstd{,}
    \hlstr{"Mass_Repro"}\hlstd{,}\hlstr{"Eggs"}\hlstd{)}
\end{alltt}


{\ttfamily\noindent\bfseries\color{errorcolor}{\#\# Error in names(MSexta)[c(3, 6, 7, 8)] <- c("{}Surv\_Repro"{}, "{}Time\_Repro"{}, : object 'MSexta' not found}}\begin{alltt}
\hlkwd{names}\hlstd{(MSexta)}
\end{alltt}


{\ttfamily\noindent\bfseries\color{errorcolor}{\#\# Error in eval(expr, envir, enclos): object 'MSexta' not found}}\end{kframe}
\end{knitrout}

%#names(MSexta)
%#sapply(MSexta, class)
%#levels(MSexta$id)
%#levels(MSexta$Sex)
%#sapply(MSexta, function(x) any(is.na(x)))

%#MSexta[MSexta$Surv_Repro == 0, ]
%# note that time to reproduction and mass at reproduction are undefined
%# if individual did not survive to reproduction
%#MSexta[MSexta$Surv_Repro == 1 & MSexta$Sex == "M", ]
%# note that for all males ovariole count is undefined, have to trim this
%# node off of graph for males.  But also have some other undefined values.
%#MSexta[MSexta$Surv_Repro == 1 & MSexta$Sex == "F", ]
%# note that for some females ovariole count is also undefined, have to trim this
%# node off of graph for these females.

This dataset has two individuals that have an ovariole count
but have sex declared as \texttt{NA}.  We change sex from \texttt{NA}
to female for these individuals.

\begin{knitrout}
\definecolor{shadecolor}{rgb}{0.969, 0.969, 0.969}\color{fgcolor}\begin{kframe}
\begin{alltt}
\hlstd{sex.change} \hlkwb{<-} \hlstd{(}\hlopt{!} \hlkwd{is.na}\hlstd{(MSexta}\hlopt{$}\hlstd{Eggs))} \hlopt{&} \hlkwd{is.na}\hlstd{(MSexta}\hlopt{$}\hlstd{Sex)}
\end{alltt}


{\ttfamily\noindent\bfseries\color{errorcolor}{\#\# Error in eval(expr, envir, enclos): object 'MSexta' not found}}\begin{alltt}
\hlkwd{sum}\hlstd{(sex.change)}
\end{alltt}


{\ttfamily\noindent\bfseries\color{errorcolor}{\#\# Error in eval(expr, envir, enclos): object 'sex.change' not found}}\begin{alltt}
\hlstd{MSexta[sex.change,} \hlstr{"Sex"}\hlstd{]} \hlkwb{<-} \hlstr{"F"}
\end{alltt}


{\ttfamily\noindent\bfseries\color{errorcolor}{\#\# Error in MSexta[sex.change, "{}Sex"{}] <- "{}F"{}: object 'MSexta' not found}}\end{kframe}
\end{knitrout}


Sex cannot be determined for \emph{M. sexta} that die before pupation
so the dataset at this point has \texttt{NA} values recorded for the
sex of such individuals.  The \texttt{aster} and \texttt{aster2} packages
do not allow \texttt{NA} values in covariates,
so change sex from \texttt{NA} to \texttt{U} for unobservable.

\begin{knitrout}
\definecolor{shadecolor}{rgb}{0.969, 0.969, 0.969}\color{fgcolor}\begin{kframe}
\begin{alltt}
\hlstd{fred} \hlkwb{<-} \hlkwd{as.character}\hlstd{(MSexta}\hlopt{$}\hlstd{Sex)}
\end{alltt}


{\ttfamily\noindent\bfseries\color{errorcolor}{\#\# Error in eval(expr, envir, enclos): object 'MSexta' not found}}\begin{alltt}
\hlstd{fred[}\hlkwd{is.na}\hlstd{(fred)]} \hlkwb{<-} \hlstr{"U"}
\end{alltt}


{\ttfamily\noindent\bfseries\color{errorcolor}{\#\# Error in fred[is.na(fred)] <- "{}U"{}: object 'fred' not found}}\begin{alltt}
\hlstd{MSexta}\hlopt{$}\hlstd{Sex} \hlkwb{<-} \hlkwd{as.factor}\hlstd{(fred)}
\end{alltt}


{\ttfamily\noindent\bfseries\color{errorcolor}{\#\# Error in is.factor(x): object 'fred' not found}}\end{kframe}
\end{knitrout}

There is some data on female \emph{M. sexta} that reached eclosion but
had \texttt{NA} recorded for ovariole count.  Since the fitness of such
individuals is unknown (missing data), they provide no information about
fitness, so we remove them from the data.

\begin{knitrout}
\definecolor{shadecolor}{rgb}{0.969, 0.969, 0.969}\color{fgcolor}\begin{kframe}
\begin{alltt}
\hlstd{condition} \hlkwb{<-} \hlstd{(MSexta}\hlopt{$}\hlstd{Sex} \hlopt{==} \hlstr{"F"}\hlstd{)} \hlopt{&} \hlkwd{is.na}\hlstd{(MSexta}\hlopt{$}\hlstd{Eggs)} \hlopt{&}
    \hlstd{(MSexta}\hlopt{$}\hlstd{Surv_Repro} \hlopt{==} \hlnum{1}\hlstd{)}
\end{alltt}


{\ttfamily\noindent\bfseries\color{errorcolor}{\#\# Error in eval(expr, envir, enclos): object 'MSexta' not found}}\begin{alltt}
\hlkwd{sum}\hlstd{(condition)}
\end{alltt}


{\ttfamily\noindent\bfseries\color{errorcolor}{\#\# Error in eval(expr, envir, enclos): object 'condition' not found}}\begin{alltt}
\hlstd{MSexta} \hlkwb{<-} \hlstd{MSexta[}\hlopt{!}\hlstd{condition,]}
\end{alltt}


{\ttfamily\noindent\bfseries\color{errorcolor}{\#\# Error in eval(expr, envir, enclos): object 'MSexta' not found}}\end{kframe}
\end{knitrout}


We now add the survival to pupation variable, called $P$
in the aster graph (p.~\pageref{graph}).
This variable is constructed using the \texttt{Sex} variable.
Every individual that has sex recorded as male or female survived to pupation,
and every individual that has no sex recorded died before pupation. 


\begin{knitrout}
\definecolor{shadecolor}{rgb}{0.969, 0.969, 0.969}\color{fgcolor}\begin{kframe}
\begin{alltt}
\hlstd{P} \hlkwb{<-} \hlkwd{as.numeric}\hlstd{(MSexta}\hlopt{$}\hlstd{Sex} \hlopt{!=} \hlstr{"U"}\hlstd{)}
\end{alltt}


{\ttfamily\noindent\bfseries\color{errorcolor}{\#\# Error in eval(expr, envir, enclos): object 'MSexta' not found}}\begin{alltt}
\hlstd{MSexta} \hlkwb{<-} \hlkwd{data.frame}\hlstd{(MSexta,} \hlkwc{P} \hlstd{= P)}
\end{alltt}


{\ttfamily\noindent\bfseries\color{errorcolor}{\#\# Error in data.frame(MSexta, P = P): object 'MSexta' not found}}\end{kframe}
\end{knitrout}

Other nodes in the graph  are constructed using
the survival to reproduction variables \verb@Surv_Repro@ and \verb@Time_Repro@.
\begin{knitrout}
\definecolor{shadecolor}{rgb}{0.969, 0.969, 0.969}\color{fgcolor}\begin{kframe}
\begin{alltt}
\hlstd{surv.repro} \hlkwb{<-} \hlstd{MSexta}\hlopt{$}\hlstd{Surv_Repro}
\end{alltt}


{\ttfamily\noindent\bfseries\color{errorcolor}{\#\# Error in eval(expr, envir, enclos): object 'MSexta' not found}}\begin{alltt}
\hlstd{time.repro} \hlkwb{<-} \hlstd{MSexta}\hlopt{$}\hlstd{Time_Repro}
\end{alltt}


{\ttfamily\noindent\bfseries\color{errorcolor}{\#\# Error in eval(expr, envir, enclos): object 'MSexta' not found}}\end{kframe}
\end{knitrout}
Time of death is not recorded for individuals that do not reach eclosion.
\begin{knitrout}
\definecolor{shadecolor}{rgb}{0.969, 0.969, 0.969}\color{fgcolor}\begin{kframe}
\begin{alltt}
\hlkwd{all}\hlstd{(}\hlkwd{is.na}\hlstd{(time.repro)} \hlopt{==} \hlstd{(surv.repro} \hlopt{==} \hlnum{0}\hlstd{))}
\end{alltt}


{\ttfamily\noindent\bfseries\color{errorcolor}{\#\# Error in eval(expr, envir, enclos): object 'time.repro' not found}}\end{kframe}
\end{knitrout}

So that \texttt{NA} values do not propagate into the variables we are creating,
we eliminate these \texttt{NA} values.
\begin{knitrout}
\definecolor{shadecolor}{rgb}{0.969, 0.969, 0.969}\color{fgcolor}\begin{kframe}
\begin{alltt}
\hlstd{time.repro[}\hlkwd{is.na}\hlstd{(time.repro)]} \hlkwb{<-} \hlnum{0}
\end{alltt}


{\ttfamily\noindent\bfseries\color{errorcolor}{\#\# Error in time.repro[is.na(time.repro)] <- 0: object 'time.repro' not found}}\end{kframe}
\end{knitrout}

Since the day of death for individuals that did not survive to eclosion
was not recorded,
we cannot say on what day death occurred for these individuals and hence
record the death at age 33 days.  For all other individuals, we record
when they reached eclosion by creating all the indicator variables in
the graph.  For each of the ``$T$ with subscripts'' variables
in the graph, we create an R variable of that name except we
flatten the subscripts ($T_{35_1}$ becomes \texttt{T351}).  Because
all deaths after pupation are recorded on day 33, that day is special.
There is a \texttt{T330} variable but no $\texttt{T}x\texttt{0}$ variable for
$x > 33$.  There are $\texttt{T}x\texttt{1}$ and $\texttt{T}x\texttt{2}$
variables for all $x$ between 33 and 40.
\begin{knitrout}
\definecolor{shadecolor}{rgb}{0.969, 0.969, 0.969}\color{fgcolor}\begin{kframe}
\begin{alltt}
\hlstd{T330} \hlkwb{<-} \hlkwd{as.numeric}\hlstd{((MSexta}\hlopt{$}\hlstd{Surv_Repro} \hlopt{==} \hlnum{0}\hlstd{)} \hlopt{&} \hlstd{(MSexta}\hlopt{$}\hlstd{P} \hlopt{==} \hlnum{1}\hlstd{))}
\end{alltt}


{\ttfamily\noindent\bfseries\color{errorcolor}{\#\# Error in eval(expr, envir, enclos): object 'MSexta' not found}}\begin{alltt}
\hlkwa{for} \hlstd{(j} \hlkwa{in} \hlnum{33}\hlopt{:}\hlnum{40}\hlstd{) \{}
    \hlstd{varname1} \hlkwb{<-} \hlkwd{paste}\hlstd{(}\hlstr{"T"}\hlstd{, j,} \hlstr{"1"}\hlstd{,} \hlkwc{sep} \hlstd{=} \hlstr{""}\hlstd{)}
    \hlstd{varname2} \hlkwb{<-} \hlkwd{paste}\hlstd{(}\hlstr{"T"}\hlstd{, j,} \hlstr{"2"}\hlstd{,} \hlkwc{sep} \hlstd{=} \hlstr{""}\hlstd{)}
    \hlkwd{assign}\hlstd{(varname1,} \hlkwd{as.numeric}\hlstd{(time.repro} \hlopt{>} \hlstd{j))}
    \hlkwd{assign}\hlstd{(varname2,} \hlkwd{as.numeric}\hlstd{(time.repro} \hlopt{==} \hlstd{j))}
\hlstd{\}}
\end{alltt}


{\ttfamily\noindent\bfseries\color{errorcolor}{\#\# Error in assign(varname1, as.numeric(time.repro > j)): object 'time.repro' not found}}\end{kframe}
\end{knitrout}

Now we also need the ovariole count variables, and we change \texttt{NA}
to zero there too, which is harmless because all females that survive
to eclosion now have egg counts.

We have another problem that will not arise until we use
the \texttt{asterdata} function, but it is easiest to fix here.
The graph for females is different from the graph for males and unobserved.
Females have egg counts, which are the $\texttt{B}x$ nodes of the graph.
Others don't.  But the R function \texttt{asterdata} only deals with
data such that every individual has the same graph.  The idea is do that first,
and fix it later using the \texttt{subset} function to delete any unwanted
nodes.  But if we do that,
the males and unknown must have valid egg counts (even though they are
bogus and will be deleted later) or the \texttt{asterdata} function will
complain.  Since we are using
the zero-truncated Poisson distribution for egg counts, every male
who reaches eclosion must have at least one egg (even though this is
bogus and we are going to delete his egg count node later).  The unknowns
do not reach eclosion, so no fix-up is needed for them.
\begin{knitrout}
\definecolor{shadecolor}{rgb}{0.969, 0.969, 0.969}\color{fgcolor}\begin{kframe}
\begin{alltt}
\hlstd{eggs} \hlkwb{<-} \hlstd{MSexta}\hlopt{$}\hlstd{Eggs}
\end{alltt}


{\ttfamily\noindent\bfseries\color{errorcolor}{\#\# Error in eval(expr, envir, enclos): object 'MSexta' not found}}\begin{alltt}
\hlstd{sex} \hlkwb{<-} \hlkwd{as.character}\hlstd{(MSexta}\hlopt{$}\hlstd{Sex)}
\end{alltt}


{\ttfamily\noindent\bfseries\color{errorcolor}{\#\# Error in eval(expr, envir, enclos): object 'MSexta' not found}}\begin{alltt}
\hlstd{eggs[sex} \hlopt{==} \hlstr{"M"}\hlstd{]} \hlkwb{<-} \hlnum{NA}
\end{alltt}


{\ttfamily\noindent\bfseries\color{errorcolor}{\#\# Error in eggs[sex == "{}M"{}] <- NA: object 'eggs' not found}}\begin{alltt}
\hlkwd{all}\hlstd{(}\hlkwd{is.na}\hlstd{(eggs)} \hlopt{==} \hlstd{((surv.repro} \hlopt{==} \hlnum{0}\hlstd{)} \hlopt{|} \hlstd{(sex} \hlopt{!=} \hlstr{"F"}\hlstd{)))}
\end{alltt}


{\ttfamily\noindent\bfseries\color{errorcolor}{\#\# Error in eval(expr, envir, enclos): object 'eggs' not found}}\begin{alltt}
\hlstd{eggs[}\hlkwd{is.na}\hlstd{(eggs)]} \hlkwb{<-} \hlnum{0}
\end{alltt}


{\ttfamily\noindent\bfseries\color{errorcolor}{\#\# Error in eggs[is.na(eggs)] <- 0: object 'eggs' not found}}\begin{alltt}
\hlstd{eggs[(sex} \hlopt{==} \hlstr{"M"}\hlstd{)} \hlopt{&} \hlstd{(surv.repro} \hlopt{==} \hlnum{1}\hlstd{)]} \hlkwb{<-} \hlnum{1}
\end{alltt}


{\ttfamily\noindent\bfseries\color{errorcolor}{\#\# Error in eggs[(sex == "{}M"{}) \& (surv.repro == 1)] <- 1: object 'eggs' not found}}\end{kframe}
\end{knitrout}
The variable $\texttt{B}x$ is equal to the egg count for that individual
if the individual reached eclosion on day $x$ and is zero otherwise.
\begin{knitrout}
\definecolor{shadecolor}{rgb}{0.969, 0.969, 0.969}\color{fgcolor}\begin{kframe}
\begin{alltt}
\hlkwa{for} \hlstd{(j} \hlkwa{in} \hlnum{33}\hlopt{:}\hlnum{40}\hlstd{) \{}
    \hlstd{varname} \hlkwb{<-} \hlkwd{paste}\hlstd{(}\hlstr{"B"}\hlstd{, j,} \hlkwc{sep} \hlstd{=} \hlstr{""}\hlstd{)}
    \hlkwd{assign}\hlstd{(varname,} \hlkwd{ifelse}\hlstd{(time.repro} \hlopt{==} \hlstd{j, eggs,} \hlnum{0}\hlstd{))}
\hlstd{\}}
\end{alltt}


{\ttfamily\noindent\bfseries\color{errorcolor}{\#\# Error in ifelse(time.repro == j, eggs, 0): object 'time.repro' not found}}\end{kframe}
\end{knitrout}

Finally we need to fix the one covariate that still has \texttt{NA} values.
\begin{knitrout}
\definecolor{shadecolor}{rgb}{0.969, 0.969, 0.969}\color{fgcolor}\begin{kframe}
\begin{alltt}
\hlkwd{all}\hlstd{(}\hlkwd{is.na}\hlstd{(MSexta}\hlopt{$}\hlstd{Mass_Repro)} \hlopt{==} \hlstd{(MSexta}\hlopt{$}\hlstd{Surv_Repro} \hlopt{==} \hlnum{0}\hlstd{))}
\end{alltt}


{\ttfamily\noindent\bfseries\color{errorcolor}{\#\# Error in eval(expr, envir, enclos): object 'MSexta' not found}}\begin{alltt}
\hlstd{MSexta}\hlopt{$}\hlstd{Mass_Repro[}\hlkwd{is.na}\hlstd{(MSexta}\hlopt{$}\hlstd{Mass_Repro)]} \hlkwb{<-} \hlnum{0}
\end{alltt}


{\ttfamily\noindent\bfseries\color{errorcolor}{\#\# Error in MSexta\$Mass\_Repro[is.na(MSexta\$Mass\_Repro)] <- 0: object 'MSexta' not found}}\end{kframe}
\end{knitrout}

We now stuff all these variables into a data frame, in the process changing
the units on some of the variables: the units for \verb@Time_2nd@ change
from days to weeks, and the units for \verb@Mass_2nd@ and \verb@Mass_Repro@
change from milligrams to grams.  These changes of units were perhaps
unnecessary, but were done with the idea that they might make computations
more stable.

\begin{knitrout}
\definecolor{shadecolor}{rgb}{0.969, 0.969, 0.969}\color{fgcolor}\begin{kframe}
\begin{alltt}
\hlstd{Msextadata} \hlkwb{<-} \hlkwd{data.frame}\hlstd{(}\hlkwc{LarvaID} \hlstd{= MSexta}\hlopt{$}\hlstd{LarvaID,} \hlkwc{Sex} \hlstd{= MSexta}\hlopt{$}\hlstd{Sex,}
    \hlkwc{Time_2nd} \hlstd{= MSexta}\hlopt{$}\hlstd{Time_2nd} \hlopt{/} \hlnum{7}\hlstd{,} \hlkwc{Mass_2nd} \hlstd{= MSexta}\hlopt{$}\hlstd{Mass_2nd} \hlopt{/} \hlnum{1000}\hlstd{,}
    \hlkwc{Mass_Repro} \hlstd{= MSexta}\hlopt{$}\hlstd{Mass_Repro} \hlopt{/} \hlnum{1000}\hlstd{, P, T330, T331, T332, B33,}
    \hlstd{T341, T342, B34, T351, T352, B35, T361, T362, B36, T371, T372,}
    \hlstd{B37, T381, T382, B38, T391, T392, B39, T401, T402, B40)}
\end{alltt}


{\ttfamily\noindent\bfseries\color{errorcolor}{\#\# Error in data.frame(LarvaID = MSexta\$LarvaID, Sex = MSexta\$Sex, Time\_2nd = MSexta\$Time\_2nd/7, : object 'MSexta' not found}}\end{kframe}
\end{knitrout}

And check for no \texttt{NA} values.

\begin{knitrout}
\definecolor{shadecolor}{rgb}{0.969, 0.969, 0.969}\color{fgcolor}\begin{kframe}
\begin{alltt}
\hlopt{!} \hlkwd{any}\hlstd{(}\hlkwd{is.na}\hlstd{(Msextadata))}
\end{alltt}


{\ttfamily\noindent\bfseries\color{errorcolor}{\#\# Error in eval(expr, envir, enclos): object 'Msextadata' not found}}\end{kframe}
\end{knitrout}

Now we need to specify the rest of the graph shown in Figure~1.
\begin{knitrout}
\definecolor{shadecolor}{rgb}{0.969, 0.969, 0.969}\color{fgcolor}\begin{kframe}
\begin{alltt}
\hlstd{vars} \hlkwb{<-} \hlkwd{c}\hlstd{(}\hlstr{"P"}\hlstd{,} \hlstr{"T330"}\hlstd{,} \hlstr{"T331"}\hlstd{,} \hlstr{"T332"}\hlstd{,} \hlstr{"B33"}\hlstd{,} \hlstr{"T341"}\hlstd{,} \hlstr{"T342"}\hlstd{,}
    \hlstr{"B34"}\hlstd{,} \hlstr{"T351"}\hlstd{,} \hlstr{"T352"}\hlstd{,} \hlstr{"B35"}\hlstd{,} \hlstr{"T361"}\hlstd{,} \hlstr{"T362"}\hlstd{,} \hlstr{"B36"}\hlstd{,}
    \hlstr{"T371"}\hlstd{,} \hlstr{"T372"}\hlstd{,} \hlstr{"B37"}\hlstd{,} \hlstr{"T381"}\hlstd{,} \hlstr{"T382"}\hlstd{,} \hlstr{"B38"}\hlstd{,} \hlstr{"T391"}\hlstd{,}
    \hlstr{"T392"}\hlstd{,} \hlstr{"B39"}\hlstd{,} \hlstr{"T401"}\hlstd{,} \hlstr{"T402"}\hlstd{,} \hlstr{"B40"}\hlstd{)}
\hlstd{pred} \hlkwb{<-} \hlkwd{c}\hlstd{(}\hlnum{0}\hlstd{,} \hlnum{1}\hlstd{,} \hlnum{1}\hlstd{,} \hlnum{1}\hlstd{,} \hlnum{4}\hlstd{,} \hlnum{3}\hlstd{,} \hlnum{3}\hlstd{,} \hlnum{7}\hlstd{,} \hlnum{6}\hlstd{,} \hlnum{6}\hlstd{,} \hlnum{10}\hlstd{,} \hlnum{9}\hlstd{,} \hlnum{9}\hlstd{,} \hlnum{13}\hlstd{,} \hlnum{12}\hlstd{,} \hlnum{12}\hlstd{,}
    \hlnum{16}\hlstd{,} \hlnum{15}\hlstd{,} \hlnum{15}\hlstd{,} \hlnum{19}\hlstd{,} \hlnum{18}\hlstd{,} \hlnum{18}\hlstd{,} \hlnum{22}\hlstd{,} \hlnum{21}\hlstd{,} \hlnum{21}\hlstd{,} \hlnum{25}\hlstd{)}
\hlstd{group} \hlkwb{<-} \hlkwd{c}\hlstd{(}\hlnum{0}\hlstd{,} \hlnum{0}\hlstd{,} \hlnum{2}\hlstd{,} \hlnum{3}\hlstd{,} \hlnum{0}\hlstd{,} \hlnum{0}\hlstd{,} \hlnum{6}\hlstd{,} \hlnum{0}\hlstd{,} \hlnum{0}\hlstd{,} \hlnum{9}\hlstd{,} \hlnum{0}\hlstd{,} \hlnum{0}\hlstd{,} \hlnum{12}\hlstd{,} \hlnum{0}\hlstd{,} \hlnum{0}\hlstd{,} \hlnum{15}\hlstd{,}
    \hlnum{0}\hlstd{,} \hlnum{0}\hlstd{,} \hlnum{18}\hlstd{,} \hlnum{0}\hlstd{,} \hlnum{0}\hlstd{,} \hlnum{21}\hlstd{,} \hlnum{0}\hlstd{,} \hlnum{0}\hlstd{,} \hlnum{24}\hlstd{,} \hlnum{0}\hlstd{)}
\hlstd{families} \hlkwb{<-} \hlkwd{list}\hlstd{(}\hlstr{"bernoulli"}\hlstd{,} \hlkwd{fam.multinomial}\hlstd{(}\hlnum{3}\hlstd{),}
    \hlkwd{fam.zero.truncated.poisson}\hlstd{(),} \hlkwd{fam.multinomial}\hlstd{(}\hlnum{2}\hlstd{))}
\hlstd{code} \hlkwb{<-} \hlkwd{c}\hlstd{(}\hlnum{1}\hlstd{,} \hlnum{2}\hlstd{,} \hlnum{2}\hlstd{,} \hlnum{2}\hlstd{,} \hlnum{3}\hlstd{,} \hlnum{4}\hlstd{,} \hlnum{4}\hlstd{,} \hlnum{3}\hlstd{,} \hlnum{4}\hlstd{,} \hlnum{4}\hlstd{,} \hlnum{3}\hlstd{,} \hlnum{4}\hlstd{,} \hlnum{4}\hlstd{,} \hlnum{3}\hlstd{,} \hlnum{4}\hlstd{,} \hlnum{4}\hlstd{,} \hlnum{3}\hlstd{,}
    \hlnum{4}\hlstd{,} \hlnum{4}\hlstd{,} \hlnum{3}\hlstd{,} \hlnum{4}\hlstd{,} \hlnum{4}\hlstd{,} \hlnum{3}\hlstd{,} \hlnum{4}\hlstd{,} \hlnum{4}\hlstd{,} \hlnum{3}\hlstd{)}
\end{alltt}
\end{kframe}
\end{knitrout}

One of these variables that specify the graph is obvious.
\texttt{vars} names the variables that are pasted together to form
the response vector (all the variables at all the nodes of the graph).
The others are less obvious.  \texttt{pred} specifies the arrows of the
graph, and \texttt{group} specifies the lines.  If \texttt{pred[j]} is
not zero, then there is an arrow from node \texttt{pred[j]} to node \texttt{j}.
Otherwise, there is a line from the initial node
(marked 1 in the aster graph) to node \texttt{j}.
Similarly, if \texttt{group[j]} is nonzero, then there is a line from
node \texttt{group[j]} to node \texttt{j}.
\texttt{code} is an index vector into the families vector;
\verb@families[code[j]]@ is the family for the conditional distribution
of the dependence group containing node \texttt{j} given its predecessor
node.  Since these are so hard to follow, we check that these variables
are set right.
\begin{knitrout}
\definecolor{shadecolor}{rgb}{0.969, 0.969, 0.969}\color{fgcolor}\begin{kframe}
\begin{alltt}
\hlstd{charpred} \hlkwb{<-} \hlkwd{c}\hlstd{(}\hlstr{"Initial"}\hlstd{, vars)[pred} \hlopt{+} \hlnum{1}\hlstd{]}
\hlstd{groupidx} \hlkwb{<-} \hlkwd{seq}\hlstd{(}\hlkwc{along} \hlstd{= group)}
\hlkwa{for} \hlstd{(j} \hlkwa{in} \hlkwd{seq}\hlstd{(}\hlkwc{along} \hlstd{= groupidx))}
    \hlkwa{if} \hlstd{(group[j]} \hlopt{!=} \hlnum{0}\hlstd{)}
        \hlstd{groupidx[j]} \hlkwb{<-} \hlstd{groupidx[group[j]]}
\hlkwd{data.frame}\hlstd{(}\hlkwc{predecessor} \hlstd{= charpred,} \hlkwc{successor} \hlstd{= vars,}
     \hlkwc{group} \hlstd{= groupidx, code)}
\end{alltt}
\begin{verbatim}
##    predecessor successor group code
## 1      Initial         P     1    1
## 2            P      T330     2    2
## 3            P      T331     2    2
## 4            P      T332     2    2
## 5         T332       B33     5    3
## 6         T331      T341     6    4
## 7         T331      T342     6    4
## 8         T342       B34     8    3
## 9         T341      T351     9    4
## 10        T341      T352     9    4
## 11        T352       B35    11    3
## 12        T351      T361    12    4
## 13        T351      T362    12    4
## 14        T362       B36    14    3
## 15        T361      T371    15    4
## 16        T361      T372    15    4
## 17        T372       B37    17    3
## 18        T371      T381    18    4
## 19        T371      T382    18    4
## 20        T382       B38    20    3
## 21        T381      T391    21    4
## 22        T381      T392    21    4
## 23        T392       B39    23    3
## 24        T391      T401    24    4
## 25        T391      T402    24    4
## 26        T402       B40    26    3
\end{verbatim}
\end{kframe}
\end{knitrout}

In each line here, \texttt{predecessor} is the predecessor
(variable at the tail of an arrow) and
\texttt{successor} is the successor
(variable at the head of an arrow)
for some arrow in the graph,
\texttt{group} is now an arbitrary number that is the same for successor
nodes in the same dependence group and different for successor nodes
in different dependence groups, and \texttt{code} is the same as before.
We see that only $\{ \texttt{T330}, \texttt{T331}, \texttt{T332} \}$
form a three-node dependence group, and their distribution is (same
for all three, as it must be, because this is their joint distribution)
three-dimensional multinomial.
We see that $\{ \texttt{T}x\texttt{1}, \texttt{T}x\texttt{2} \}$
form a two-node dependence group for $x > 33$, and their distribution is (same
for both) two-dimensional multinomial.
We see that $P$ is a dependence group all by itself, and its distribution
is Bernoulli.
We see that $\texttt{B}x$ is a dependence group all by itself for $x \ge 33$,
and its distribution is zero-truncated Poisson.
We see that the predecessor of $P$ is the initial node,
the predecessor of the $\texttt{T33}$ dependence group is $P$,
the predecessor of the $\texttt{T}x$ dependence group is
the $\texttt{T}w\texttt{1}$ node where $w = x - 1$
(the predecessor of $T_{x_y}$ is $T_{w_1}$),
and the predecessor of $\texttt{B}x$ is $\texttt{T}x\texttt{2}$ for $x \ge 33$.
It all checks.

Now we make the asterdata object.
\begin{knitrout}
\definecolor{shadecolor}{rgb}{0.969, 0.969, 0.969}\color{fgcolor}\begin{kframe}
\begin{alltt}
\hlstd{fred} \hlkwb{<-} \hlkwd{asterdata}\hlstd{(Msextadata,} \hlkwc{vars} \hlstd{= vars,} \hlkwc{pred} \hlstd{= pred,}
    \hlkwc{group} \hlstd{= group,} \hlkwc{code} \hlstd{= code,} \hlkwc{families} \hlstd{= families)}
\end{alltt}


{\ttfamily\noindent\bfseries\color{errorcolor}{\#\# Error in is.data.frame(data): object 'Msextadata' not found}}\begin{alltt}
\hlkwd{dim}\hlstd{(fred}\hlopt{$}\hlstd{redata)}
\end{alltt}


{\ttfamily\noindent\bfseries\color{errorcolor}{\#\# Error in eval(expr, envir, enclos): object 'fred' not found}}\begin{alltt}
\hlkwd{nrow}\hlstd{(Msextadata)} \hlopt{*} \hlkwd{length}\hlstd{(vars)}
\end{alltt}


{\ttfamily\noindent\bfseries\color{errorcolor}{\#\# Error in nrow(Msextadata): object 'Msextadata' not found}}\begin{alltt}
\hlkwd{names}\hlstd{(fred}\hlopt{$}\hlstd{redata)}
\end{alltt}


{\ttfamily\noindent\bfseries\color{errorcolor}{\#\# Error in eval(expr, envir, enclos): object 'fred' not found}}\end{kframe}
\end{knitrout}
And then we remove nodes of the graph that are not supposed to be there:
ovariole count nodes for males and unknown.
\begin{knitrout}
\definecolor{shadecolor}{rgb}{0.969, 0.969, 0.969}\color{fgcolor}\begin{kframe}
\begin{alltt}
\hlstd{condition2} \hlkwb{<-} \hlstd{(fred}\hlopt{$}\hlstd{redata}\hlopt{$}\hlstd{Sex} \hlopt{==} \hlstr{"F"}\hlstd{)} \hlopt{|} \hlstd{(}\hlopt{!} \hlkwd{grepl}\hlstd{(}\hlstr{"B"}\hlstd{, fred}\hlopt{$}\hlstd{redata}\hlopt{$}\hlstd{varb))}
\end{alltt}


{\ttfamily\noindent\bfseries\color{errorcolor}{\#\# Error in eval(expr, envir, enclos): object 'fred' not found}}\begin{alltt}
\hlstd{fred} \hlkwb{<-} \hlkwd{subset}\hlstd{(fred,} \hlkwc{subset} \hlstd{= condition2,} \hlkwc{successors} \hlstd{=} \hlnum{FALSE}\hlstd{)}
\end{alltt}


{\ttfamily\noindent\bfseries\color{errorcolor}{\#\# Error in subset(fred, subset = condition2, successors = FALSE): object 'fred' not found}}\begin{alltt}
\hlkwd{dim}\hlstd{(fred}\hlopt{$}\hlstd{redata)}
\end{alltt}


{\ttfamily\noindent\bfseries\color{errorcolor}{\#\# Error in eval(expr, envir, enclos): object 'fred' not found}}\begin{alltt}
\hlkwd{nrow}\hlstd{(Msextadata)} \hlopt{*} \hlkwd{length}\hlstd{(vars)} \hlopt{-}
    \hlkwd{sum}\hlstd{(Msextadata}\hlopt{$}\hlstd{Sex} \hlopt{!=} \hlstr{"F"}\hlstd{)} \hlopt{*} \hlkwd{sum}\hlstd{(}\hlkwd{grepl}\hlstd{(}\hlstr{"B"}\hlstd{, vars))}
\end{alltt}


{\ttfamily\noindent\bfseries\color{errorcolor}{\#\# Error in nrow(Msextadata): object 'Msextadata' not found}}\begin{alltt}
\hlkwd{names}\hlstd{(fred}\hlopt{$}\hlstd{redata)}
\end{alltt}


{\ttfamily\noindent\bfseries\color{errorcolor}{\#\# Error in eval(expr, envir, enclos): object 'fred' not found}}\end{kframe}
\end{knitrout}

Now we check this against the \texttt{hornworm} dataset from
the \texttt{aster2} package.
\begin{knitrout}
\definecolor{shadecolor}{rgb}{0.969, 0.969, 0.969}\color{fgcolor}\begin{kframe}
\begin{alltt}
\hlkwd{identical}\hlstd{(fred, hornworm)}
\end{alltt}


{\ttfamily\noindent\bfseries\color{errorcolor}{\#\# Error in identical(fred, hornworm): object 'fred' not found}}\end{kframe}
\end{knitrout}
If \texttt{FALSE}, this is i18n insanity, which is not really R's fault.
The details section of \texttt{help(Comparison)} has an explicit disclaimer
that sort order is not guaranteed to work the same way on different machines.
\begin{knitrout}
\definecolor{shadecolor}{rgb}{0.969, 0.969, 0.969}\color{fgcolor}\begin{kframe}
\begin{alltt}
\hlstd{foo} \hlkwb{<-} \hlstd{fred}\hlopt{$}\hlstd{redata}\hlopt{$}\hlstd{LarvaID}
\end{alltt}


{\ttfamily\noindent\bfseries\color{errorcolor}{\#\# Error in eval(expr, envir, enclos): object 'fred' not found}}\begin{alltt}
\hlstd{bar} \hlkwb{<-} \hlstd{hornworm}\hlopt{$}\hlstd{redata}\hlopt{$}\hlstd{LarvaID}
\hlkwd{identical}\hlstd{(foo, bar)}
\end{alltt}
\begin{verbatim}
## [1] FALSE
\end{verbatim}
\begin{alltt}
\hlstd{foo} \hlkwb{<-} \hlkwd{as.character}\hlstd{(foo)}
\hlstd{bar} \hlkwb{<-} \hlkwd{as.character}\hlstd{(bar)}
\hlkwd{identical}\hlstd{(foo, bar)}
\end{alltt}
\begin{verbatim}
## [1] FALSE
\end{verbatim}
\begin{alltt}
\hlstd{baz} \hlkwb{<-} \hlkwd{levels}\hlstd{(hornworm}\hlopt{$}\hlstd{redata}\hlopt{$}\hlstd{LarvaID)}
\hlkwd{all}\hlstd{(foo} \hlopt{%in%} \hlstd{baz)}
\end{alltt}
\begin{verbatim}
## [1] FALSE
\end{verbatim}
\begin{alltt}
\hlstd{qux} \hlkwb{<-} \hlkwd{factor}\hlstd{(foo,} \hlkwc{levels} \hlstd{= baz)}
\hlstd{fred}\hlopt{$}\hlstd{redata}\hlopt{$}\hlstd{LarvaID} \hlkwb{<-} \hlstd{qux}
\end{alltt}


{\ttfamily\noindent\bfseries\color{errorcolor}{\#\# Error in fred\$redata\$LarvaID <- qux: object 'fred' not found}}\begin{alltt}
\hlkwd{identical}\hlstd{(fred, hornworm)}
\end{alltt}


{\ttfamily\noindent\bfseries\color{errorcolor}{\#\# Error in identical(fred, hornworm): object 'fred' not found}}\end{kframe}
\end{knitrout}

\begin{knitrout}
\definecolor{shadecolor}{rgb}{0.969, 0.969, 0.969}\color{fgcolor}\begin{kframe}


{\ttfamily\noindent\bfseries\color{errorcolor}{\#\# Error in eval(expr, envir, enclos): object 'fred' not found}}\end{kframe}
\end{knitrout}












\newpage
\begin{thebibliography}{}

% \bibitem[Arnold and Wade(1984)]{arnold-wade}
% Arnold, S.~J., and M.~J. Wade. 1984.
% \newblock On the measurement of natural and sexual selection: theory.
% \newblock \emph{Evolution} 38:709--719.

%\bibitem[Blows and Brooks(2003)]{blows-brooks}
%Blows, M.W., and R. Brooks. 2003.
%\newblock Measuring nonlinear selection.
%\newblock \emph{Am. Nat.} 162:815--820.

\bibitem[Charlesworth(1980)]{charles}
Charlesworth, B. 1980.
\newblock \emph{Evolution in age-structured populations}.
\newblock Cambridge Univ. Press, Cambridge, U.~K.

%\bibitem[Charmantier and Gienapp(2014)]{charmantier-gienapp}
%Charmantier, A., and P. Gienapp. 2014.
%\newblock Climate change and timing of avian breeding and migration:
%    evolutionary versus plastic changes.
%\newblock \emph{Evol. Appl.} 7:15-28.

%\bibitem[Clutton-Brock and Sheldon(2010)]{clutton-brock-sheldon}
%Clutton-Brock, T., and B.~C. Sheldon. 2010.
%\newblock Individuals and populations: the role of long-term,
%    individual-based studies of animals in ecology and evolutionary biology.
%\newblock \emph{Trends Ecol. Evol.} 25:562-573.

% \bibitem[Diamond, et al.(2010)Diamond, Hawkins, Nijhout, and Kingsolver]
%     {diamond-et-al}
% Diamond, S.~E., S.~D. Hawkins, H.~F. Nijhout, and J.~G. Kingsolver. 2010.
% \newblock Evolutionary divergence of field and laboratory populations of
%     \emph{Manduca sexta} in response to host-plant quality.
% \newblock Ecological Entomology 35:166--174.

%\bibitem[Diamond and Kingsolver(2010a)]{diamond-kingsolver-too}
%Diamond, S.~E., and J.~G. Kingsolver. 2010a.
%\newblock Environmental dependence of thermal reaction norms: Host plant
%    quality can reverse the temperature-size rule.
%\newblock \emph{Am. Nat.} 175:1--10.

%\bibitem[Diamond and Kingsolver(2010b)]{diamond-kingsolver}
%Diamond, S.~E., and J.~G. Kingsolver. 2010b.
%\newblock Fitness consequences of host plant choice: a field experiment.
%\newblock \emph{Oikos} 119:542--550.

\bibitem[Fisher(1930)]{fisher}
Fisher, R.~A. 1930.
\newblock \emph{The genetical theory of natural selection}.
\newblock Clarendon Press, Oxford, U.~K.

\bibitem[Geyer(2009)]{geyer-gdor}
Geyer, C.~J. (2009).
\newblock Likelihood inference in exponential families and directions
    of recession.
\newblock \emph{Electronic Journal of Statistics}, \textbf{3}, 259--289.

\bibitem[Geyer(2010)Geyer]{phil}
Geyer, C.~J. (2010).
\newblock A philosophical look at aster models.
\newblock Technical Report No. 676. School of Statistics,
    University of Minnesota.
\newblock \url{http://purl.umn.edu/57163}.

\bibitem[Geyer(2014)]{aster-package}
Geyer, C.~J. 2014.
\newblock R package \texttt{aster} (aster models), version 0.8-30.
\newblock \url{http://cran.r-project.org/package=aster}.

\bibitem[Geyer(2015)]{aster2-package}
Geyer, C.~J. 2015.
\newblock R package aster2 (aster models), version 0.2-1.
\newblock \url{http://cran.r-project.org/package=aster2}.

\bibitem[Geyer, et al.(2007)Geyer, Wagenius, and Shaw]{aster1}
Geyer, C.~J., S. Wagenius, and R.~G. Shaw. 2007.
\newblock Aster models for life history analysis.
\newblock \emph{Biometrika} 94:415--426.

%\bibitem[Kingsolver and Diamond(2011)]{kingsolver-diamond}
%Kingsolver, J.~G., and S.~E. Diamond. 2011.
%\newblock Phenotypic selection in natural populations: What limits
%    directional selection?
%\newblock \emph{Am. Nat.} 177:346--357.

\bibitem[Kingsolver, et al.(2012)Kingsolver, Diamond, Seiter,
    and Higgins]{kingsolver-et-al}
Kingsolver, J.~G., S.~E. Diamond, S.~A. Seiter, and J.~K. Higgins. 2012.
\newblock Direct and indirect phenotypic selection on developmental
    trajectories in \emph{Manduca sexta}.
\newblock \emph{Funct. Ecol.} 26:598--607.

%\bibitem[Kingsolver, et al.(2001)Kingsolver, Hoekstra, Hoekstra, Berrigan,
%    Vignieri, Hill, Hoang, Gilbert, and Beerli]{kingsolver-et-al-2001}
%Kingsolver, J.~G., H.~E. Hoekstra, J.~M. Hoekstra, D. Berrigan, S.~N. Vignieri,
%    C.~E. Hill, A. Hoang, P. Gilbert, and P. Beerli. 2001.
%\newblock The strength of phenotypic selection in natural populations.
%\newblock \emph{Am. Nat.} 157:245-261.

%\bibitem[Klein, et al.(2011)Klein, Carpentier, Oddou-Muratorio]{klein}
%Klein, E.~K., Carpentier, F.~H., Oddou-Muratorio, S. 2011.
%\newblock Estimating the variance of male fecundity from genotypes of 
%    progeny arrays: evaluation of the Bayesian forward approach.
%\newblock \emph{Methods Ecol. Evol.}, \textbf{2}, 349-361.

%\bibitem[Lande and Arnold(1983)]{lande}
%Lande, R., and S. Arnold. 1983.
%\newblock The measurement of selection on correlated characters.
%\newblock \emph{Evolution} 37:1210--1226.

%\bibitem[Lenski and Service(1982)]{lenski-service}
%Lenski, R.~E., and P.~M. Service. 1982.
%\newblock The statistical analysis of population growth rates calculated
%    from schedules of survivorship and fecundity.
%\newblock \emph{Ecology} 63:655–662.

%\bibitem[Madden and Chamberlin(1945)]{madden-chamberlin}
%Madden, A.~H., and F.~S. Chamberlin. 1945.
%\newblock Biology of the tobacco hornworm in the southern cigar-tobacco
%    district.
%\newblock \emph{U. S. Dep. Agric. Tech. Bull.} 896.
%\newblock \url{http://purl.umn.edu/169886}.

% \bibitem[Mechaber and Hildebrand(2000)]{mechaber-hildebrand}
% Mechaber, W.~L., and J.~G. Hildebrand. 2000.
% \newblock Novel, non-solanaceous hostplant record for \emph{Manduca sexta}
%     (Lepidoptera: Sphingidae) in the southwestern United States.
% \newblock \emph{Ann. Entomol. Soc. Am.} 93:447--451.

% \bibitem[Mira and Bernays(2002)]{mira-bernays}
% Mira, A., and E.~A. Bernays. 2002.
% \newblock Trade-offs in host use by \emph{Manduca sexta:} plant characters
%     vs natural enemies.
% \newblock \emph{Oikos} 97:387--397.

%\bibitem[Nijhout and Davidowitz(2009)]{nijhout-davidowitz}
%Nijhout, H.~F., and G. Davidowitz. (2009).
%\newblock The developmental-physiological basis of phenotypic plasticity.
%\newblock Pp. 589--608 \emph{in} D.~W. Whitman and T.~N. Ananthakrishnan, eds.
%\newblock \emph{Phenotypic plasticity of insects: mechanisms and consequences}.
%\newblock Science Publishers, Enfield, NH.

%\bibitem[Ozgul et al.(2010)Ozgul, Childs, Oli, Armitage, Blumstein, Olson,
%    Tuljapurkar, and Coulson]{ozgul-et-al-2010}
%Ozgul, A., D.~Z. Childs, M.~K. Oli, K.~B. Armitage, D.~T. Blumstein,
%   L.~E. Olson, S. Tuljapurkar, and T. Coulson. 2010.
%\newblock Coupled dynamics of body mass and population growth in response
%    to environmental change.
%\newblock \emph{Nature} 466:482--485.

%\bibitem[Ozgul et al.(2009)Ozgul, Tuljapurkar, Benton, Pemberton,
%    Clutton-Brock, and Coulson]{ozgul-et-al-2009}
%Ozgul, A., S. Tuljapurkar, T.~G. Benton, J.~M. Pemberton, T.~H. Clutton-Brock,
%    and T. Coulson. 2009.
%\newblock The dynamics of phenotypic change and the shrinking sheep
%    of St. Kilda.
%\newblock \emph{Science} 325:464-467.

\bibitem[R Development Core Team, 2014]{rcore}
R Development Core Team. 2014.
\newblock R: A language and environment for statistical computing.
\newblock R Foundation for Statistical Computing, Vienna, Austria.
\newblock \url{http://www.R-project.org}.

%\bibitem[Roff(2002)]{roff}
%Roff, D.~A. 2002.
%\newblock \emph{Life history evolution}.
%\newblock Sinauer Associates, Sunderland MA.

\bibitem[Shaw and Geyer(2010)]{aster3}
Shaw, R.~G., and C.~J. Geyer. 2010.
\newblock Inferring fitness landscapes.
\newblock \emph{Evolution} 64:2510--2520.

\bibitem[Shaw, et al.(2008)Shaw, Geyer, Wagenius, Hangelbroek, and
    Etterson]{aster2}
Shaw, R.~G., C.~J. Geyer, S. Wagenius, H. Hangelbroek, and J.~R. Etterson. 2008.
\newblock Unifying life-history analyses for 
    inference of fitness and population growth.
\newblock \emph{Am. Nat.} 172:E35--E47.

\bibitem[Shaw, et~al.(2015)Shaw, Wagenius, and Geyer]{aster-aphid}
Shaw, R.~G., Wagenius, S., and Geyer, C.~J. (2015).
\newblock The susceptibility of \emph{Echinacea angustifolia} to a
    specialist aphid: eco-evolutionary perspective on genotypic variation
    and demographic consequences.
\newblock \emph{Journal of Ecology}, \textbf{103}, 809--818.

%\bibitem[Siepielski, et al.(2009)Siepielski, DiBattista,
%    and Carlson]{siepielski-et-al}
%Siepielski, A.~M., J.~D. DiBattista, and S.~M. Carlson. 2009.
%\newblock It's about time: the temporal dynamics of phenotypic selection
%    in the wild.
%\newblock \emph{Ecol. Lett.} 12:1261--1276.

%\bibitem[Yamamoto(1969)]{yamamoto}
%Yamamoto, R.~T. 1969.
%\newblock Mass rearing of the tobacco hornworm. II. Larval rearing
%    and pupation.
%\newblock \emph{J. Econ. Entomol.} 62:1427--1437.

%\bibitem[Yamauchi and Yoshitake(1984)]{yamauchi-yoshitake}
%Yamauchi, H., and N. Yoshitake. 1984.
%\newblock Developmental stages of ovarian follicles of the silkworm,
%    \emph{Bombyx mori} L.
%\newblock \emph{J. Morphol.} 179:21--31.

\end{thebibliography}





\end{document}
