\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}}}%
\let\hlipl\hlkwb

\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{amsthm}
\usepackage{amsmath}
\usepackage{amssymb}
\usepackage{natbib}
\usepackage{url}
\usepackage{graphicx}
\usepackage[utf8]{inputenc}
\usepackage{tikz-cd}


\newcommand{\R}{\mathbb{R}}
\newcommand{\E}{\text{E}}
\newcommand{\Var}{\text{Var}}
\newcommand{\Fi}{\Sigma}
\newcommand{\iFi}{\Sigma^{-1}}
\newcommand{\eiFi}{\widehat{\Sigma}^{-1}}
\newcommand{\eFi}{\widehat{\Sigma}}
\newcommand{\Rsoft}{\texttt{R}}
\newcommand{\astertwo}{\texttt{aster2}\;}
\newcommand{\betahat}{\hat{\beta}}
\newcommand{\tauhat}{\hat{\tau}}
\newcommand{\betaenv}{\hat{\beta}_{\text{env}}}
\newcommand{\tauenv}{\hat{\tau}_{\text{env}}}
\newcommand{\upenv}{\hat{\upsilon}_{\text{env}}}
\newcommand{\B}{\mathrm{B}}
\newcommand{\T}{\mathrm{T}}
\newcommand{\G}{\mathcal{G}}
\newcommand{\D}{\mathcal{D}}
\newcommand{\Su}{\mathcal{S}}
\newcommand{\Menv}{M_{\text{env}}}
\newcommand{\genv}{\hat{g}_{\text{env}}}
\newcommand{\MD}{M_{\text{1D}}}
\newcommand{\upD}{\hat{\upsilon}_{\text{1D}}}

\newcommand{\inner}[1]{\langle #1 \rangle}

\usepackage{geometry}
\geometry{margin=1in}

%\setlength{\parindent}{0cm}
\newtheorem{lem}{Lemma} 
\newtheorem{thm}{Theorem} 
\newtheorem{prop}{Proposition} 
\newtheorem{cor}{Corollary} 
\newtheorem{defn}{Definition}
\allowdisplaybreaks
\IfFileExists{upquote.sty}{\usepackage{upquote}}{}
\begin{document}

  \vspace*{0.9375in}
%  \begin{center}
%    {\bfseries Web-based Supplementary Materials for } \\
%    {\bfseries ``An Application of Envelope Methodology and Aster Models''} \vspace{0.15cm} \\
%    by Daniel J. Eck, Charles J. Geyer, and R. Dennis Cook. \vspace{0.5cm} \\
%    % Technical Report No.~801 (I'm just making up this number)\\ % need new number
%    School of Statistics \\
%    University of Minnesota \\
%     \today
%  \end{center}
  \begin{center}
    {\bfseries Supporting Data Analysis for \\
    ``Combining Envelope Methodology and Aster Models for Variance Reduction 
      in Life History Analyses''} \\
    By \\
    Daniel J. Eck, Charles J. Geyer, and R. Dennis Cook \\
    Technical Report No.~699 \\ 
    School of Statistics \\
    University of Minnesota \\
     \today 
  \end{center}

  \thispagestyle{empty}
  \cleardoublepage
  \setcounter{page}{1}
  \thispagestyle{empty}



In this technical report, we reproduce the examples that appear in 
\citet*{eck} and show that the combination of envelope methodology with 
aster models estimates expected Darwinian fitness with lower variability 
than estimation conducted with aster models alone. The \texttt{envlpaster} 
package \citep{envlp-package} is developed to facilitate these calculations. 
This package implements the bootstrap algorithms discussed in the paper and 
contains the two datasets corresponding to the examples discussed in the 
paper. There are other functions in this package that we do not explain in 
this technical report. Those functions fall into two categories. They are 
either necessary for implementing the 1D algorithm or they implement other 
bootstrap techniques not discussed in \citet*{eck}. The other bootstrap 
techniques are explained in the help pages of the \texttt{envlpaster} 
package. 

The bootstraps and model selection calculations in this technical report take 
substantial time to run. As a result, these calculations are not performed 
using this \texttt{knitr} script; they are performed in five \texttt{.RData} 
files and the code that produced these calculations is given in this technical 
report. The five \texttt{.RData} files are: 
\begin{itemize}
\item[1.] \texttt{ex1-Efronboot.RData}
\item[2.] \texttt{ex1-secondboot.RData}
\item[3.] \texttt{ex1-AIC.RData}
%\item[4.] \texttt{ex2-Efronboot.RData}
%\item[5.] \texttt{ex2-secondboot.RData}
\item[4.] \texttt{techreport\_eigenboot\_fitboot.RData}
\item[5.] \texttt{main-all-terms-DavidLowrydata.RData}
%\item[7.] \texttt{ex2-Efronboot-1d.RData}
%\item[8.] \texttt{ex2-secondboot-1d.RData}
%\item[9.] \texttt{ex2-MLEboot.RData}
%\item[10.] \texttt{ex2-AIC.RData}
\end{itemize}

Daniel J. Eck (\url{daniel.eck@yale.edu}) is the current maintainer of these 
\texttt{.RData} files and they are available on request.


\section{Replicating Example 1}
A population of 3000 organisms was simulated to form the dataset used in this 
aster analysis. %We generated this data to include a true envelope model. 
The code generating this dataset is provided in Appendix B. The last line of 
Appendix B shows that this code reproduces the dataset included in the 
\texttt{envlpaster} package. The organisms in this dataset either die or 
remain living for ten years. We observe three random variables for each year 
$i$ of the study. These variables are $U_i, Vi,$ and $W_i$. The variables 
$U_i$ and $Vi$ are Bernoulli where $U_i$ models survival and $V_i$ models 
whether or not an individual reproduces conditional on survival. The variables 
$W_i$ count offspring produced in the $i$-th year conditional on reproduction. 
Thus, $W_i$ is modeled as zero-truncated Poisson given $V_i$. The graphical 
structure for one individual is given by the graph shown in Figure 
\ref{graph3}. In these data, the variable $\sum_i W_i$ is considered to be 
Darwinian fitness. There are two covariates $(z_1,z_2)$ thought to be 
associated with Darwinian fitness and the aster model selected by the 
likelihood ratio test (LRT) is a full quadratic model with respect to these 
covariates. A full aster analysis of this data and its construction can be 
seen in \citet*{geyer4}. Note that our data differs from that in 
\citet*{geyer4} since we have included a true envelope space. 


\begin{figure}
\begin{center}
\setlength{\unitlength}{0.45 in}
\thicklines
\begin{picture}(10.3,2.5)(-0.1,-2.3)
\put(0,0){\makebox(0,0){$1_{\hphantom{0}}$}}
\put(1,0){\makebox(0,0){$U_1$}}
\put(2,0){\makebox(0,0){$U_2$}}
\put(3,0){\makebox(0,0){$U_3$}}
\put(4,0){\makebox(0,0){$U_4$}}
\put(5,0){\makebox(0,0){$U_5$}}
\put(6,0){\makebox(0,0){$U_6$}}
\put(7,0){\makebox(0,0){$U_7$}}
\put(8,0){\makebox(0,0){$U_8$}}
\put(9,0){\makebox(0,0){$U_9$}}
\put(10,0){\makebox(0,0){$U_{10}$}}
\multiput(0.35,0)(1,0){10}{\vector(1,0){0.3}}
\put(1,-1){\makebox(0,0){$V_1$}}
\put(2,-1){\makebox(0,0){$V_2$}}
\put(3,-1){\makebox(0,0){$V_3$}}
\put(4,-1){\makebox(0,0){$V_4$}}
\put(5,-1){\makebox(0,0){$V_5$}}
\put(6,-1){\makebox(0,0){$V_6$}}
\put(7,-1){\makebox(0,0){$V_7$}}
\put(8,-1){\makebox(0,0){$V_8$}}
\put(9,-1){\makebox(0,0){$V_9$}}
\put(10,-1){\makebox(0,0){$V_{10}$}}
\multiput(1,-0.25)(1,0){10}{\vector(0,-1){0.5}}
\put(1,-2){\makebox(0,0){$W_1$}}
\put(2,-2){\makebox(0,0){$W_2$}}
\put(3,-2){\makebox(0,0){$W_3$}}
\put(4,-2){\makebox(0,0){$W_4$}}
\put(5,-2){\makebox(0,0){$W_5$}}
\put(6,-2){\makebox(0,0){$W_6$}}
\put(7,-2){\makebox(0,0){$W_7$}}
\put(8,-2){\makebox(0,0){$W_8$}}
\put(9,-2){\makebox(0,0){$W_9$}}
\put(10,-2){\makebox(0,0){$W_{10}$}}
\multiput(1,-1.25)(1,0){10}{\vector(0,-1){0.5}}
\end{picture}
\end{center}
\caption{Graphical structure of the aster model for the simulated data 
example. The top layer corresponds to survival; these random variables are 
Bernoulli. The middle layer corresponds to whether or not an individual 
reproduced; these random variables are also Bernoulli. The bottom layer 
corresponds to offspring count; these random variables are zero-truncated 
Poisson.}
\label{graph3}
\end{figure}


\subsection{Aster Calculations}

The \texttt{envlpaster} and \texttt{aster2} packages are needed to carry out 
all calculations. The \texttt{envlpaster} packages requires \texttt{aster}, 
\texttt{MASS}, \texttt{caTools} and \texttt{trust} to be installed. 


\begin{knitrout}
\definecolor{shadecolor}{rgb}{0.969, 0.969, 0.969}\color{fgcolor}\begin{kframe}
\begin{alltt}
\hlkwd{library}\hlstd{(aster2)}
\hlkwd{library}\hlstd{(envlpaster)}
\hlkwd{data}\hlstd{(simdata30nodes)}
\hlstd{data} \hlkwb{<-} \hlstd{simdata30nodes.asterdata}
\end{alltt}
\end{kframe}
\end{knitrout}


The \texttt{simdata30nodes} dataset contains almost all that is necessary to 
fit the aster model used in the paper. The \texttt{root} and \texttt{pred} 
vectors and the model matrix \texttt{modmat} necessary for aster model fitting 
are included directly in the \texttt{asterdata} object included in the 
\texttt{simdata30nodes} dataset within the \texttt{envlpaster} package. 


\begin{knitrout}
\definecolor{shadecolor}{rgb}{0.969, 0.969, 0.969}\color{fgcolor}\begin{kframe}
\begin{alltt}
\hlstd{vars} \hlkwb{<-} \hlstd{variables}
\hlstd{pred} \hlkwb{<-} \hlstd{data}\hlopt{$}\hlstd{pred}
\hlstd{families} \hlkwb{<-} \hlstd{data}\hlopt{$}\hlstd{families}
\hlstd{code} \hlkwb{<-} \hlstd{data}\hlopt{$}\hlstd{code}
\end{alltt}
\end{kframe}
\end{knitrout}



The \texttt{fam} vector, which specifies the exponential families appearing in 
the graphical structure of an aster model, is not included but is constructed 
using the \texttt{code} vector. In this example, the \texttt{code} vector has 
entries that are 1-2 valued where a 1 corresponds to Bernoulli distributions 
and a 2 corresponds to zero-truncated Poisson distributions. In the 
\texttt{aster} library, a zero-truncated Poisson distribution is coded with 
the number 3 so we change the 2's that appear in \texttt{code} to 3's when we 
construct the \texttt{fam} vector.


\begin{knitrout}
\definecolor{shadecolor}{rgb}{0.969, 0.969, 0.969}\color{fgcolor}\begin{kframe}
\begin{alltt}
\hlstd{fam} \hlkwb{<-} \hlstd{code}
\hlstd{fam[fam} \hlopt{==} \hlnum{2}\hlstd{]} \hlkwb{<-} \hlnum{3}
\end{alltt}
\end{kframe}
\end{knitrout}


The model fitting is not of class \texttt{aster.formula} and is fit using the 
aster data directly. The first \texttt{nnode} columns of the 
\texttt{simdata30nodes} dataset contain the responses necessary to fit the 
aster submodel of interest. Here \texttt{nnode} = 30. 
Recall that the full aster model is a saturated statistical model (one unknown 
parameter per observed response value) and is of little scientific use; the 
aster submodel has a reduced parameter count where the parameters of this 
model correspond to the relationship between Darwinian fitness with both 
phenotypic traits and covariates of interest. The different parameterizations 
of the aster model are seen in Figure~\ref{parm}. The relationship between the 
aster model parameterizations displayed in Figure~\ref{parm} are further 
explained in \citet*{geyer3}.


\begin{knitrout}
\definecolor{shadecolor}{rgb}{0.969, 0.969, 0.969}\color{fgcolor}\begin{kframe}
\begin{alltt}
\hlstd{nnode} \hlkwb{<-} \hlkwd{length}\hlstd{(vars)}
\hlstd{xnew} \hlkwb{<-} \hlkwd{as.matrix}\hlstd{(simdata30nodes[,}\hlkwd{c}\hlstd{(}\hlnum{1}\hlopt{:}\hlstd{nnode)])}
\hlstd{m1} \hlkwb{<-} \hlkwd{aster}\hlstd{(xnew, root, pred, fam, modmat)}
\hlstd{m1}\hlopt{$}\hlstd{formula} \hlkwb{<-} \hlstd{formula}
\hlstd{m1}\hlopt{$}\hlstd{xlevels} \hlkwb{<-} \hlstd{xlevels}
\hlstd{m1}\hlopt{$}\hlstd{terms} \hlkwb{<-} \hlstd{terms}
\hlkwd{summary}\hlstd{(m1)}
\end{alltt}
\begin{verbatim}
## 
## Call:
## NULL
## 
##              Estimate Std. Error z value Pr(>|z|)    
## (Intercept) -0.414198   0.029610 -13.989  < 2e-16 ***
## vtypev       1.240414   0.040284  30.792  < 2e-16 ***
## vtypew       0.464049   0.030958  14.990  < 2e-16 ***
## uyear        0.129734   0.004546  28.539  < 2e-16 ***
## z1           0.013889   0.003181   4.366 1.26e-05 ***
## z2          -0.010015   0.003197  -3.133  0.00173 ** 
## I(z1^2)     -0.015673   0.001915  -8.185 2.72e-16 ***
## I(z2^2)     -0.010950   0.001884  -5.813 6.14e-09 ***
## I(z1 * z2)   0.015817   0.002966   5.333 9.66e-08 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
\end{verbatim}
\end{kframe}
\end{knitrout}


There are nine components in both the submodel canonical parameter vector 
$\beta$ and the submodel mean-value parameter vector $\tau$. We write 
$$  
  \tau = \left(\begin{array}{c}
     \gamma \\
     \upsilon
  \end{array}\right) = M^T\mu
$$
where $\gamma \in \R^{p-k}$ are nuisance parameters and $\upsilon \in \R^k$ 
are the parameters of interest that are relevant to the estimation of 
expected Darwinian fitness, $M$ is the model matrix, and $\mu$ is the 
mean-value parameter vector of the unconditional aster model, see 
Figure~\ref{parm}. In our example, $p = 9$ and $k = 5$. The parameter vector 
$\upsilon$ are the five unknown slope coefficients for the full quadratic 
model in terms of the covariates $z_1$ and $z_2$. Ultimately, there are four  
parameterizations of interest that are present in the aster analyses we 
consider. These four parameterizations are: 
\begin{itemize}
\item[1)] The aster model canonical parameter vector $\beta \in \R^p$. 
\item[2)] The aster model mean-value parameter vector $\tau \in \R^p$. 
\item[3)] The unconditional aster model canonical parameter vector 
  $\varphi \in \R^m$.
\item[4)] The unconditional aster model mean-value parameter vector 
  $\mu \in \R^m$. 
%\item[5)] The best surrogate of Darwinian fitness defined in \citet*{eck} as 
%  $h(\mu) \in \R^d$ or $g(\tau) \in \R^d$ where $h$ and $g$ are a 
%  differentiable functions and $d$ is the number of hypothetical individuals 
%  that expected Darwinian fitness is estimated for. 
\end{itemize}
Consult Figure~\ref{parm} to see the relationships between these 
parameterizations. 

% Our methods are not applicable with respect to conditional aster models (needs to be included).
Darwinian fitness can be defined as a mapping from either $\mu$ or $\tau$. 
It is often easier to conceptulize the map from $\mu$ to Darwinian fitness.
In this example, true Darwinian fitness is unknown. The best available 
surrogate of Darwinian fitness is total offspring count. When Darwinian 
fitness is taken to be total offspring count, the function 
$h(\mu) = A\mu$ where $A$ is a matrix with entries $a_{i,j} = 1$ when $\mu_j$ 
corresponds to offspring for individual $i$. The matrix $A$ specifies that 
expected Darwinian fitness for every individual is the sum of offspring nodes 
for that individual. Our methods will use the map from $\tau$ to Darwinian 
fitness given by $g(\tau)$. This is because $\tau$ is of a dimension that is 
less than the sample size and is feasible for envelope methodology. 
%Aster software can be used 
%to directly program in $h(\mu)$ and our theory makes use of $g(\tau)$ where 
%$\tau$ is of a far lower dimension than $\mu$ in the applications used




\begin{figure}[t!]
\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 = a + M\beta$} ([yshift= -5pt]b.east)
([xshift= 5pt]a.south) edge node[right] {$\xi_j = c_j^{\prime}(\theta)$} ([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_{\beta} c(a + M\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}
\end{center}
\caption{A depiction of the transformations necessary to change aster model 
  parameterizations. Unlabeled arrows represent transformations that need to 
  be solved using optimization software. $M$ is a known model matrix, $a$ is a 
  known offset vector, and $c$ is the cumulant function for the unconditional 
  aster model.}
\label{parm}
\end{figure}




\begin{knitrout}
\definecolor{shadecolor}{rgb}{0.969, 0.969, 0.969}\color{fgcolor}\begin{kframe}
\begin{alltt}
\hlstd{beta} \hlkwb{<-} \hlstd{m1}\hlopt{$}\hlstd{coef; p} \hlkwb{<-} \hlkwd{length}\hlstd{(beta)}
\hlstd{target} \hlkwb{<-} \hlnum{5}\hlopt{:}\hlnum{9}\hlstd{; m} \hlkwb{<-} \hlkwd{length}\hlstd{(target)}
\hlstd{x} \hlkwb{<-} \hlstd{m1}\hlopt{$}\hlstd{x}
\hlstd{nind} \hlkwb{<-} \hlstd{n} \hlkwb{<-} \hlkwd{nrow}\hlstd{(x)}
\end{alltt}
\end{kframe}
\end{knitrout}


%<<>>=
%offset <- as.vector(m1$origin)
%@


We obtain the maximum likelihood estimator (MLE) of the submodel mean-value 
parameter, denoted $\tauhat$ using \texttt{aster} software and the relations 
in Figure~\ref{parm}. 

\begin{knitrout}
\definecolor{shadecolor}{rgb}{0.969, 0.969, 0.969}\color{fgcolor}\begin{kframe}
\begin{alltt}
\hlstd{mu} \hlkwb{<-} \hlkwd{predict}\hlstd{(m1,} \hlkwc{parm.type} \hlstd{=} \hlstr{"mean.value"}\hlstd{,}
  \hlkwc{model.type} \hlstd{=} \hlstr{"unconditional"}\hlstd{)}
\hlstd{modmat.mat} \hlkwb{<-} \hlkwd{matrix}\hlstd{(modmat,} \hlkwc{ncol} \hlstd{= p)}
\hlstd{tau} \hlkwb{<-} \hlkwd{crossprod}\hlstd{(modmat.mat, mu)}
\end{alltt}
\end{kframe}
\end{knitrout}

We now estimate expected Darwinian fitness for 144 hypothetical individuals 
possessing unique values of \texttt{z1} and \texttt{z2} suggested by the 
original data. These hypothetical individuals will be used to build a fitness 
landscape for estimated expected Darwinian fitness. The fitness landscape will 
be estimated using both maximum likelihood estimation and envelope estimation. 
Before we can build the dataframe corresponding to these 144 hypothetical 
individuals, the \texttt{vars} vector containing the names of all the nodes 
in Figure~\ref{graph3} is specified. 

\begin{knitrout}
\definecolor{shadecolor}{rgb}{0.969, 0.969, 0.969}\color{fgcolor}\begin{kframe}
\begin{alltt}
\hlstd{vars} \hlkwb{<-} \hlkwd{as.vector}\hlstd{(}\hlkwd{outer}\hlstd{(}\hlkwd{c}\hlstd{(}\hlstr{"u"}\hlstd{,} \hlstr{"v"}\hlstd{,} \hlstr{"w"}\hlstd{),} \hlnum{1}\hlopt{:}\hlnum{10}\hlstd{, paste,}
  \hlkwc{sep} \hlstd{=} \hlstr{""}\hlstd{))}
\end{alltt}
\end{kframe}
\end{knitrout}

The dataframe corresponding to the 144 hypothetical individuals is constructed 
using the code below.

\begin{knitrout}
\definecolor{shadecolor}{rgb}{0.969, 0.969, 0.969}\color{fgcolor}\begin{kframe}
\begin{alltt}
\hlstd{nx} \hlkwb{<-} \hlnum{12}\hlstd{; npop} \hlkwb{<-} \hlstd{nx}\hlopt{^}\hlnum{2}
\hlstd{z1.cand} \hlkwb{<-} \hlkwd{seq}\hlstd{(}\hlkwc{from} \hlstd{=} \hlopt{-}\hlnum{6}\hlstd{,} \hlkwc{to} \hlstd{=} \hlnum{6}\hlstd{,} \hlkwc{length} \hlstd{= nx)}
\hlstd{z2.cand} \hlkwb{<-} \hlkwd{seq}\hlstd{(}\hlkwc{from} \hlstd{=} \hlopt{-}\hlnum{6}\hlstd{,} \hlkwc{to} \hlstd{=} \hlnum{6}\hlstd{,} \hlkwc{length} \hlstd{= nx)}
\hlstd{fred} \hlkwb{<-} \hlkwd{data.frame}\hlstd{(} \hlkwd{cbind}\hlstd{(}\hlkwd{rep}\hlstd{(z1.cand,} \hlkwc{each} \hlstd{= nx),}
  \hlkwd{rep}\hlstd{(z2.cand, nx),} \hlkwd{matrix}\hlstd{(}\hlnum{1}\hlstd{,}\hlkwc{nrow} \hlstd{= nx}\hlopt{^}\hlnum{2}\hlstd{,} \hlkwc{ncol} \hlstd{=} \hlnum{30}\hlstd{)))}
\hlkwd{colnames}\hlstd{(fred)} \hlkwb{<-} \hlkwd{c}\hlstd{(}\hlstr{"z1"}\hlstd{,}\hlstr{"z2"}\hlstd{,vars)}
\hlstd{fred}\hlopt{$}\hlstd{root} \hlkwb{<-} \hlnum{1}
\end{alltt}
\end{kframe}
\end{knitrout}

Next, \texttt{amat} needs to be specified. \texttt{amat} specifies which 
nodes in the graphical structure for the hypothetical individuals are 
inputs for Darwinian fitness (usually terminal nodes in an aster graph). 
\texttt{amat} can either be a three-dimensional array or a matrix. 
We construct a three-dimensional array as in \citet*{geyer4}. 

\begin{knitrout}
\definecolor{shadecolor}{rgb}{0.969, 0.969, 0.969}\color{fgcolor}\begin{kframe}
\begin{alltt}
\hlstd{amat} \hlkwb{<-} \hlkwd{array}\hlstd{(}\hlnum{0}\hlstd{,} \hlkwd{c}\hlstd{(npop, nnode, npop))}
\hlstd{foot} \hlkwb{<-} \hlkwd{grepl}\hlstd{(}\hlstr{"w"}\hlstd{, vars)}
\hlkwa{for} \hlstd{(k} \hlkwa{in} \hlnum{1}\hlopt{:}\hlstd{npop) amat[k, foot, k]} \hlkwb{<-} \hlnum{1}
\end{alltt}
\end{kframe}
\end{knitrout}

The \texttt{newdata} argument is specified with a dataframe for the 
hypothetical individuals created to build the fitness landscape. 
This dataframe has already been constructed and has been given the name 
\texttt{fred}. The model matrix and dataframe in long format, 
corresponding to the provided dataframe specified by \texttt{newdata}, 
are constructed internally. If the construction of either the new model matrix 
or the new dataframe requires more calculation than a simple 
\texttt{model.matrix} or \texttt{renewdata} call respectively, then these 
quantities can be included using the \texttt{modmat.new} or \texttt{renewdata} 
arguments respectively. 

\begin{knitrout}
\definecolor{shadecolor}{rgb}{0.969, 0.969, 0.969}\color{fgcolor}\begin{kframe}
\begin{alltt}
\hlstd{renewdata} \hlkwb{<-} \hlkwd{reshape}\hlstd{(fred,} \hlkwc{varying} \hlstd{=} \hlkwd{list}\hlstd{(vars),}
  \hlkwc{direction} \hlstd{=} \hlstr{"long"}\hlstd{,} \hlkwc{timevar} \hlstd{=} \hlstr{"varb"}\hlstd{,}
  \hlkwc{times} \hlstd{=} \hlkwd{as.factor}\hlstd{(vars),} \hlkwc{v.names} \hlstd{=} \hlstr{"resp"}\hlstd{)}
\hlstd{layer2} \hlkwb{<-} \hlkwd{gsub}\hlstd{(}\hlstr{"[0-9]"}\hlstd{,} \hlstr{""}\hlstd{,}
  \hlkwd{as.character}\hlstd{(renewdata}\hlopt{$}\hlstd{varb))}
\hlstd{renewdata} \hlkwb{<-} \hlkwd{data.frame}\hlstd{(renewdata,} \hlkwc{layer} \hlstd{= layer2)}
\hlstd{renewdata}\hlopt{$}\hlstd{vtype} \hlkwb{<-} \hlkwd{as.factor}\hlstd{(}\hlkwd{substr}\hlstd{(}\hlkwd{as.character}\hlstd{(}
  \hlstd{renewdata}\hlopt{$}\hlstd{varb),}\hlnum{1}\hlstd{,} \hlnum{1}\hlstd{))}
\hlstd{renewdata}\hlopt{$}\hlstd{year} \hlkwb{<-} \hlkwd{as.numeric}\hlstd{(}\hlkwd{substring}\hlstd{(}\hlkwd{as.character}\hlstd{(}
  \hlstd{renewdata}\hlopt{$}\hlstd{varb),} \hlnum{2}\hlstd{))}
\hlstd{renewdata}\hlopt{$}\hlstd{uyear} \hlkwb{<-} \hlstd{renewdata}\hlopt{$}\hlstd{year} \hlopt{*}
  \hlkwd{as.numeric}\hlstd{(}\hlkwd{as.character}\hlstd{(renewdata}\hlopt{$}\hlstd{vtype)} \hlopt{==} \hlstr{"u"}\hlstd{)}
\end{alltt}
\end{kframe}
\end{knitrout}

We also need to specify the \texttt{modmat.new} argument. The dataframe in 
long format is calculated as an intermediate step to the construction of this 
model matrix.

\begin{knitrout}
\definecolor{shadecolor}{rgb}{0.969, 0.969, 0.969}\color{fgcolor}\begin{kframe}
\begin{alltt}
\hlstd{modmat.renew} \hlkwb{<-} \hlkwd{model.matrix}\hlstd{(m1}\hlopt{$}\hlstd{formula,} \hlkwc{data} \hlstd{= renewdata)}
\hlstd{modmat.renew.array} \hlkwb{<-} \hlkwd{array}\hlstd{(modmat.renew,} \hlkwc{dim} \hlstd{=} \hlkwd{c}\hlstd{(npop, nnode, p))}
\end{alltt}
\end{kframe}
\end{knitrout}

We can now use \texttt{predict.aster} to estimate expected Darwinian fitness 
for all 144 hypothetical individuals. The last \texttt{rowSums} call adds 
all of the estimated mean values corresponding to $W$ in the aster graph 
for each hypothetical individual.

\begin{knitrout}
\definecolor{shadecolor}{rgb}{0.969, 0.969, 0.969}\color{fgcolor}\begin{kframe}
\begin{alltt}
\hlstd{p1} \hlkwb{<-} \hlkwd{predict}\hlstd{(m1,} \hlkwc{modmat} \hlstd{= modmat.renew.array,} \hlkwc{varvar} \hlstd{= varb,}
  \hlkwc{idvar} \hlstd{= id,} \hlkwc{root} \hlstd{= renewdata}\hlopt{$}\hlstd{root,} \hlkwc{se.fit} \hlstd{=} \hlnum{TRUE}\hlstd{)}
\hlstd{p1.mat} \hlkwb{<-} \hlkwd{matrix}\hlstd{(p1}\hlopt{$}\hlstd{fit,} \hlkwc{ncol} \hlstd{= nnode)}
\hlstd{p1.mat.fit} \hlkwb{<-} \hlkwd{rowSums}\hlstd{(}
  \hlstd{p1.mat[,} \hlkwd{which}\hlstd{(}\hlkwd{grepl}\hlstd{(}\hlstr{"w"}\hlstd{,} \hlkwd{unique}\hlstd{(renewdata}\hlopt{$}\hlstd{varb)))]}
\hlstd{)}
\end{alltt}
\end{kframe}
\end{knitrout}



%#p1$fit
%data.fit <- asterdata(fred, vars, pred, group = rep(0,30), code,
%  families)

%offset <- as.vector(m1$origin)
%rigin <- m1$origin[1:npop,]
%origin <- as.vector(origin)
%phi.MLE <- modmat.renew %*% beta + origin
%u.MLE <- transformSaturated(phi.MLE, data.fit, from = "phi",
%  to = "mu")
%mat.fit <- matrix(amat, nrow = length(mu.MLE))
%fit.MLE <- t(amat.fit) %*% mu.MLE
%fit.MLE <- as.vector(fit.MLE)
%@

%We now have the required quantities necessary to estimate the variability of 
%estimated expected Darwinian fitness using envelope methodology. The envelope 
%estimator of expected Darwinian fitness for these hypothetical individuals has asymptotic distributions given by 
%\begin{equation}\label{asymptotics-envelope-unknown}
%  \sqrt{n}\left(g(\tauenv) - g(\tau)\right) 
%   \overset{d}{\longrightarrow} 
% N\left(0,\; \nabla g(\tau)\Delta_2 \nabla g(\tau)^T \right)
%end{equation}
%where $\Delta_2$ is the unknown asymptotic covariance matrix of $\tauenv$. The asymptotic covariance matrix of estimated expected Darwinian fitness is estimated using the parametric bootstrap adjusting for model selection as seen in Algorithm 1. The dimension of the envelope space is selected using model selection criteria and is not known in advance. \citet*[Section 4]{efron} provides a double bootstrap procedure which accounts for model selection. This procedure is used to estimate the variance of estimated expected Darwinian fitness using envelope methodology. At first, data are generated as a realization from the aster model distribution evaluated at the envelope estimator. The envelope dimension is determined for each of these generated datasets. The estimator of expected Darwinian fitness is then taken to be the average of all of the envelope estimators obtained from these datasets. To estimate the variability of this envelope estimator, we generate further datasets from the aster model distribution evaluated at each separate envelope estimator of expected Darwinian fitness used to calculate the aforementioned average. The steps of this procedure are presented in \citet*[Algorithm 1]{eck}. For completeness, the steps of this algorithm are restated here. To avoid confusion, we label this algorithm as Algorithm 1.










\subsection{Envelope methodology and calculations}

We now transition from aster models to the combination of aster models and 
envelope methodology. The goal of this combination is to incorporate envelope 
methodology with respect to parameters of interest in the $\tau$ 
parameterization to yield efficiency gains when estimating expected Darwinian 
fitness. These efficiency gains have the potential for stronger inferences to 
be made about the driving features influencing Darwinian fitness. For example, 
suppose that one is interested in trait combinations that are associated with 
high values of estimated expected Darwinian fitness. There may be one distinct 
combination of traits which maximizes estimated expected Darwinian fitness and 
many others that provide estimates close to this maximum value. Each 
estimate comes with uncertainty. When the uncertainty is large then there will 
be many trait combinations that provide uncertainty interval estimators that 
intersect with uncertainty interval estimator of the maximum estimate of
expected Darwinian fitness. Efficiency gains obtained through envelope 
methodology can decrease that overall uncertainty. Therefore there may be 
fewer combinations of traits that produce estimates of expected Darwinian 
fitness that are close to the maximized value. 

We first provide some envelope methodology basics relevant to aster 
modeling scenarios. The particular branch of envelope methodology relevant to 
aster problems is envelope methodology discussed in \citet*{cook2}. The theory 
of envelope methodology requires a $\sqrt{n}$ consistent and asymptotically 
normal estimator of an unknown parameter vector and a $\sqrt{n}$ consistent 
estimator of its asymptotic covariance matrix \citep*{cook2}. Aster models, 
being an exponential family, satisfy these conditions.

Let $\T = \text{span}(\upsilon)$. The subspace $\T$ is a line in $\R^k$ since 
the unknown parameter $\upsilon$ is a point in $\R^k$. The envelope subspace 
$\mathcal{E}_{\Sigma_{\upsilon,\upsilon}}(\T)$ is defined as the intersection 
of all reducing subspaces of $\Sigma_{\upsilon,\upsilon}$ that contains $\T$. 
Reducing subspaces of $\Sigma_{\upsilon,\upsilon}$ are sums of eigenspaces of 
$\Sigma_{\upsilon,\upsilon}$ where the eigenspaces of 
$\Sigma_{\upsilon,\upsilon}$ are spans of the eigenvectors of 
$\Sigma_{\upsilon,\upsilon}$. An eigenvector $x$ of 
$\Sigma_{\upsilon,\upsilon}$ is a solution to the equation 
$\Sigma_{\upsilon,\upsilon}x = \lambda x$ where $\lambda$ is an eigenvalue 
needed to solve the equation. There are $k$ eigenvectors and eigenvalues in 
total. Eigenspaces are defined to be the space spanned by all eigenvectors 
corresponding to one particular eigenvalue. When all of the eigenvalues are 
unique, all distinct reducing subspaces of $\Sigma_{\upsilon,\upsilon}$ can 
be specified by the $2^k$ unique combinations of spans of eigenvectors 
(the corresponding paper does not consider the zero-dimensional space so that 
there are $2^k - 1$ distinct reducing subspaces under consideration). The 
decomposition of a matrix into eigenvalues and eigenvectors is a general 
mathematical decomposition that has a statistical context when applied to a 
covariance matrix. When applied to a covariance matrix, such as 
$\Sigma_{\upsilon,\upsilon}$, this decomposition provides insight into the 
directions of variability (given by the eigenvectors) and the strength of 
variability in that direction (given by the eigenvalues). 

The envelope subspace, as defined above to be the intersection of all 
reducing subspaces of $\Sigma_{\upsilon,\upsilon}$ that contain $\T$, 
satisfies both
\begin{itemize}
\item[1.] $\T \subset \mathcal{E}_{\Sigma_{\upsilon,\upsilon}}(\T)$
\item[2.] 
  $
    \Sigma_{\upsilon,\upsilon} 
      = P_{\mathcal{E}}\Sigma_{\upsilon,\upsilon}P_{\mathcal{E}} 
      + Q_{\mathcal{E}}\Sigma_{\upsilon,\upsilon}Q_{\mathcal{E}}
  $
\end{itemize}
where $\mathcal{E}$ denotes the envelope subspace when used as a subscript, 
$P_{\mathcal{E}}$ is the projection into 
$\mathcal{E}_{\Sigma_{\upsilon,\upsilon}}(\T)$, and $Q_{\mathcal{E}}$ is the 
projection into the orthogonal complement of 
$\mathcal{E}_{\Sigma_{\upsilon,\upsilon}}(\T)$. The reducing subspace 
assumption is needed for conditional 2 to hold. It is also of use to write 
these conditions in coordinate form. They are 
\begin{itemize}
\item[1.] $\T \subset \text{span}(\Gamma)$,
\item[2.]   
  $
    \Sigma_{\upsilon,\upsilon} 
      = \Gamma\Omega\Gamma^T + \Gamma_o\Omega_o\Gamma_o^T,
  $
\end{itemize}
where $\Gamma$ is a basis matrix for the envelope space 
$\mathcal{E}_{\Sigma_{\upsilon,\upsilon}}(\T)$, $(\Gamma,\Gamma_o)$ is an 
orthogonal matrix, and 
$P_{\mathcal{E}} = \Gamma\Gamma^T$ and $Q_{\mathcal{E}} = \Gamma_o\Gamma_o^T$. 
Intuitively, the envelope estimator reduces variability in estimation at no 
cost to consistency by making use of the conditions assumed of 
$\mathcal{E}_{\Sigma_{\upsilon,\upsilon}}(\T)$. The first condition allows 
for consistent estimation to take place. The second condition states that 
only a part of the asymptotic covariance matrix $\Sigma_{\upsilon,\upsilon}$, 
seen as $P_{\mathcal{E}}\Sigma_{\upsilon,\upsilon}P_{\mathcal{E}}$, is 
variability associated with the estimation of expected Darwinian fitness. The 
other part, seen as $Q_{\mathcal{E}}\Sigma_{\upsilon,\upsilon}Q_{\mathcal{E}}$, 
is immaterial to our estimation problem and can be discarded.

When $P_{\mathcal{E}}$ is known, the envelope estimator of $\upsilon$ is then 
$P_{\mathcal{E}}\hat{\upsilon}$, where
$$
   \sqrt{n}\left(P_{\mathcal{E}}\hat{\upsilon} - \upsilon\right)
     \overset{d}{\longrightarrow}
   N\left(0, P_{\mathcal{E}}\Sigma_{\upsilon,\upsilon} P_{\mathcal{E}} \right).
$$
With condition 2 in mind, we can see that the envelope estimator of $\upsilon$ 
has less asymptotic variability than the MLE when $P_{\mathcal{E}}$ is known. 
In actual problems, $P_{\mathcal{E}}$ is unknown and requires estimation. 
The \texttt{envlpaster} software can estimate $P_{\mathcal{E}}$ using the 1D 
algorithm proposed in \citet*[Algorithm 2]{cook2} and discussed in 
\citet*[Sections 3 and 4]{eck} or by using reducing subspaces directly as in 
\citet*[Section 4]{eck}. The reducing subspace approach to envelope estimation 
is of particular interest in the current applications and we only motivate 
this approach in this technical report. See the references in this paragraph 
for details on the 1D algorithm.%, specifically the last paragraph of 
%\citet*[Section 4]{eck}. 

%In the current analysis, envelope estimators are constructed via reducing subspaces as in \citet*[Section 4]{eck}. In our context, a reducing subspace is the space spanned by a subset of the eigenvectors of $\Sigma_{\upsilon,\upsilon}$ where $\Sigma_{\upsilon,\upsilon}$ is the part of the Fisher information matrix $\Sigma$ corresponding to parameters of interest. The envelope subspace is the intersection of all reducing subspaces $G$ that satisfy both
%\begin{itemize}
%\item[1.] $\T \subset G$
%\item[2.] $\Sigma_{\upsilon,\upsilon} = P_G\Sigma_{\upsilon,\upsilon}P_G + Q_G\Sigma_{\upsilon,\upsilon}Q_G$
%\end{itemize}
%where $P_G$ is the projection into the reducing subspace $G$ and $Q_G$ is the projection into the orthogonal complement of $G$. The restriction to reducing subspaces ensures that condition 2 holds. The two conditions assure that the variability represented by $P_G\Sigma_{\upsilon,\upsilon}P_G$ is the only variability necessary for the estimation of $\upsilon$. We stand to estimate $\upsilon$ at potentially much greater efficiency when the smallest reducing subspace $G$ that satisfies the two envelope conditions is identified. In applications, this space requires estimation and this estimation poses a cost on potential efficiency gains. When the cost of estimation is small relative to the discarded variability $Q_G\Sigma_{\upsilon,\upsilon}Q_G$, then we stand to estimate $\upsilon$ at much greater efficiency and expected Darwinian fitness at potentially greater efficiency.  

In our applications, envelope estimators are constructed using eigenvectors of 
the estimated Fisher information matrix $\widehat{\Sigma}_{\upsilon,\upsilon}$ 
which is the estimate of $\Sigma_{\upsilon,\upsilon}$ obtained from 
\texttt{aster} software. Therefore $\Gamma = \Gamma_{\G}$ and 
$\Gamma_o = \Gamma_{\G o}$ for some reducing subspace $\G$ of 
$\Sigma_{\upsilon,\upsilon}$ and 
$\widehat{\Gamma} = \widehat{\Gamma}_{\widehat{\G}}$ and 
$\widehat{\Gamma}_o = \widehat{\Gamma}_{\widehat{\G} o}$ for some reducing 
subspace $\widehat{\G}$ of $\widehat{\Sigma}_{\upsilon,\upsilon}$.
In this particular analysis, 
$\widehat{\Sigma}_{\upsilon,\upsilon}$ is the lower 5 by 5 block of 
$\widehat{\Sigma}$. All possible envelope estimators constructed from reducing 
subspaces of $\widehat{\Sigma}_{\upsilon,\upsilon}$ are produced and compared 
using the \texttt{selection} function in the \texttt{envlpaster} package. 
There are 31 possible reducing subspaces of positive dimension in total 
and the \texttt{combs} function in the \texttt{caTools} package is used keep 
track of the indices. 
Reducing subspaces are constructed one dimension at a time. At a specified 
dimension, we construct all of the possible envelope estimators using reducing 
subspaces. There are $\binom{5}{1}$ possible reducing subspaces of 
dimension $1$, $\binom{5}{2}$ possible reducing subspaces of dimension $2$, 
and so on. At a particular iteration of this procedure, denote the reducing 
subspace of $\widehat{\Sigma}_{\upsilon,\upsilon}$ by $\widehat{\G}$. 
The envelope estimator of $\upsilon$ is given by 
$$
  \upenv = \widehat{P}_{\widehat{\G}} \hat{\upsilon}. 
$$   
The envelope estimator of $\tau$ is 
$$
  \tauenv = \left(\begin{array}{c}
    \hat{\gamma} \\
    \hat{\upsilon}_{\text{env}}
  \end{array}\right).
$$
Estimators of this form correspond to an aster model with a new model matrix 
that incorporates the estimated projection into the envelope space. This 
model matrix is given as
$$  
  \Menv = M \left(\begin{array}{cc}
    I & 0 \\
    0 & \widehat{P}_{\widehat{\G}}
  \end{array}\right).
$$
The model matrix $M_{\text{env}}$ satisfies 
$\tauenv = \Menv^T Y$ which says that the envelope 
estimator of $\tau$ is a maximum likelihood estimator for the aster model 
with model matrix $\Menv$. 

\begin{prop} \label{Propred}
The envelope estimator $\tauenv$ is a maximum likelihood estimator of 
$\tau$ for the aster model with model matrix $\Menv$.
\end{prop}

\begin{proof}
We have
$
   l_{\text{env}}(\beta) = \inner{Y, \Menv\beta} - c(\Menv\beta) 
     = \inner{\Menv^TY, \beta} - c(\Menv\beta). 
$
and
$
  \nabla_{\beta} l(\beta) = \Menv^TY - \nabla_{\beta} c(\Menv\beta). 
$
Setting $\nabla_{\beta} l(\beta) = 0$ and solving for $\beta$ yields 
$\nabla_{\beta} c(\Menv\beta)|_{\beta = \hat{\beta}} = \Menv^TY = \tauenv$. 
\end{proof}

Once the envelope estimator $\tauenv$ is calculated, it is transformed to 
the $\beta$ parameterization using the model matrix $\Menv$ as an argument of 
\texttt{transformUnconditional}. The parameter vector $\beta$ is the canonical 
parameter vector in the log likelihood
\begin{equation} \label{likelihood}
 l(\beta) = \langle M^TY, \beta \rangle - c(\beta). 
\end{equation}

The transformation from $\hat{\tau}$ to $\hat{\beta}$ provides a maximizer of 
\eqref{likelihood}) with model matrix given by $\Menv$. However $\hat{\beta}$ 
is not unique since $\Menv$ is not of full column rank. 
This non identifiability is not problematic because final inference is not 
conducted with respect to the $\beta$ parameterization. Final inference is 
conducted with respect to the parameterization of Darwinian fitness which does 
not posses such problems; many distinct estimates of $\beta$ map to 
$\tauenv$. Each of these distinct estimated values of $\beta$ maps to the 
same estimate of $\Menv\beta$, which in turn maps to a common estimate of 
expected Darwinian fitness. The loss of 1-1 transformations is not an issue 
in this case.

When in the $\beta$ parameterization, we compute 
AIC, BIC, and the p-value of the LRT using \texttt{aster} software and the 
likelihood \eqref{likelihood}. The \texttt{selection} function in the 
\texttt{envlpaster} package below computes the AIC, BIC, and the p-value of 
the LRT for all 31 reducing subspaces including the model corresponding to the 
MLE. This is done by specifying the \texttt{method = "eigen"} and 
\texttt{type = "mean-value"} arguments in the \texttt{selection} function. 

\begin{knitrout}
\definecolor{shadecolor}{rgb}{0.969, 0.969, 0.969}\color{fgcolor}\begin{kframe}
\begin{alltt}
\hlstd{foo2} \hlkwb{<-} \hlkwd{selection}\hlstd{(}\hlkwc{parm} \hlstd{= tau,} \hlkwc{index} \hlstd{= target,} \hlkwc{model} \hlstd{= m1,}
  \hlkwc{data} \hlstd{= data,} \hlkwc{alpha} \hlstd{=} \hlnum{0.05}\hlstd{,} \hlkwc{type} \hlstd{=} \hlstr{"mean-value"}\hlstd{,}
  \hlkwc{method} \hlstd{=} \hlstr{"eigen"}\hlstd{)}
\end{alltt}
\end{kframe}
\end{knitrout}


The calculation above is not performed in this script. This calculation is 
performed in \\ \texttt{ex1-Efronboot.R} and is sourced into this document from 
the \texttt{ex1-Efronboot.RData} file.

\begin{knitrout}
\definecolor{shadecolor}{rgb}{0.969, 0.969, 0.969}\color{fgcolor}\begin{kframe}
\begin{alltt}
\hlkwd{load}\hlstd{(}\hlstr{"ex1-Efronboot.RData"}\hlstd{)}
\hlstd{foo2}\hlopt{$}\hlstd{aic; foo2}\hlopt{$}\hlstd{bic; foo2}\hlopt{$}\hlstd{LRT}
\end{alltt}
\begin{verbatim}
## [1] 1 2 4 5
## [1] 1 4 5
## [1] 1 4 5
\end{verbatim}
\end{kframe}
\end{knitrout}

In the original sample, we see that BIC and the LRT at size $\alpha = 0.05$ 
select the reducing subspace that is the sum of the first, fourth, and fifth 
eigenspaces of $\widehat{\Sigma}_{\upsilon,\upsilon}$. On the other hand, 
AIC selects the reducing subspace that is the sum of all eigenspaces of 
$\widehat{\Sigma}_{\upsilon,\upsilon}$ with the exception of the third 
eigenspace in the original sample. %Since eigenvalues of 
%$\widehat{\Sigma}_{\upsilon,\upsilon}$ are unique, the space $\widehat{G}_1$ 
%can also be defined as the space spanned by the first, fourth, and fifth 
%eigenvector of $\widehat{\Sigma}_{\upsilon,\upsilon}$ and $\widehat{G}_2$ 
%can be defined similarly. 
Both of these reducing subspaces are smaller than the full space which is a 
sum of all of the eigenspaces of $\widehat{\Sigma}_{\upsilon,\upsilon}$ and 
both are larger than the true reducing subspace. The true reducing subspace 
of $\widehat{\Sigma}_{\upsilon,\upsilon}$ used to construct this example, 
$G_{\text{true}}$, is the space spanned by the first and fourth eigenvectors 
of $\widehat{\Sigma}_{\upsilon,\upsilon}^{1}$ where 
$\widehat{\Sigma}_{\upsilon,\upsilon}^{1}$ is the estimated asymptotic 
covariance matrix of the parameters of interest constructed from a preliminary 
simulated dataset in \citet{geyer4}. The realizations in the dataset used in 
this example were generated from the aster model distribution with submodel 
mean-value parameter 
$
  \tau_{\text{true}} = 
    (\hat{\gamma}_1^T,P_{G_{\text{true}}}\hat{\upsilon}_1^T)^T.
$ 
Here, $\hat{\gamma}_1$ and $\hat{\upsilon}_1$ are the MLE of the nuisance 
parameters and parameters of interest obtained from the preliminary dataset. 

When the 1D algorithm is used, all three model selection criterion agree that 
the dimension of the envelope space is $u = 5$. Envelope methods are not 
interesting in this case. 

\begin{knitrout}
\definecolor{shadecolor}{rgb}{0.969, 0.969, 0.969}\color{fgcolor}\begin{kframe}
\begin{alltt}
\hlstd{foo3} \hlkwb{<-} \hlkwd{selection}\hlstd{(}\hlkwc{parm} \hlstd{= tau,} \hlkwc{index} \hlstd{= target,} \hlkwc{model} \hlstd{= m1,}
  \hlkwc{data} \hlstd{= data,} \hlkwc{alpha} \hlstd{=} \hlnum{0.05}\hlstd{,} \hlkwc{type} \hlstd{=} \hlstr{"mean-value"}\hlstd{,}
  \hlkwc{method} \hlstd{=} \hlstr{"1d"}\hlstd{)}
\hlstd{foo3}\hlopt{$}\hlstd{aic; foo3}\hlopt{$}\hlstd{bic; foo3}\hlopt{$}\hlstd{LRT}
\end{alltt}
\begin{verbatim}
## [1] 5
## [1] 5
## [1] 5
\end{verbatim}
\end{kframe}
\end{knitrout}

The 1D algorithm and our reducing subspace approach can be used to 
consistently select the correct dimension of the envelope model and the 
correct reducing subspace respectively when using BIC \citep{zhangmai, eck}. 
There is considerable model selection disagreement between the reducing 
subspaces approach and the 1D algorithm in this application. Consistent model 
selection was not helpful here and variability in model selection appears to 
be problematic. We alleviate these concerns through a parametric double 
bootstrap procedure discussed in the next section.


\subsubsection{Appendix 1: Consistency of the reducing subspace approach}

There is a close connection between envelope estimation using reducing 
subspaces and envelope estimation using the 1D algorithm. In the population, 
$\upenv = \upD$ (where $\upD$ is defined in the last paragraph of page 4 in 
\citet{eck}). 
The connection between both estimation methods exists in finite samples as 
seen in Theorem~\ref{connection}. 
In preparation, let $\widehat{\Gamma}_u$ be the estimator of the basis matrix 
of the envelope space of dimension $u$ obtained from the 1D algorithm, 
let $\widehat{\Gamma}_{uo}$ be such that 
$(\widehat{\Gamma}_{u},\widehat{\Gamma}_{uo})$ is an orthogonal matrix, 
and define orthogonal matrices 
$\widehat{O}_u = \left(\widehat{\Gamma}_u,\widehat{\Gamma}_{uo}\right)$,  
%be an orthogonal matrix where 
%$\widehat{\Gamma}_u$ is the output of the 1D algorithm with inputs 
%$M = \widehat{\Sigma}_{\upsilon,\upsilon}$ and 
%$\widehat{U} = \hat{\upsilon}\hat{\upsilon}^T$ at dimension $u$. 
%Define orthogonal matrices  
$
  \widehat{O}_{\widehat{\G}} 
    = \left(\widehat{\Gamma}_{\widehat{\G}},
      \widehat{\Gamma}_{\widehat{\G}o}\right), 
$ 
%be an orthogonal matrix where $\widehat{\Gamma}_{\widehat{\G}}$ is comprised 
%of the eigenvectors of $\widehat{\Sigma}_{\upsilon,\upsilon}$ such that 
%span$(\widehat{\Gamma}_{\widehat{\G}}) = \widehat{\G}$. Define 
and $\widehat{O} = \widehat{O}_{\widehat{\G}}\widehat{O}_u^T$. 
The matrices $\widehat{\Gamma}_{\widehat{\G}}$, 
$\widehat{\Gamma}_{\G o}$, $\widehat{\Gamma}_{u}$, and $\widehat{\Gamma}_{uo}$ 
converges in probability to $\Gamma_{\G}$, $\Gamma_{\G o}$, $\Gamma_{u}$, 
and $\Gamma_{uo}$ respectively. The convergence results for 
$\Gamma_{\G}$ and $\Gamma_{\G o}$ follow from the fact that the reducing 
subspaces of $\widehat{\Sigma}_{\upsilon,\upsilon}$ are  
$\surd{n}$ consistent estimators of the reducing subspaces of 
$\Sigma_{\upsilon,\upsilon}$. Now 
$\Gamma_{\G}^T\Gamma_{u}$ and $\Gamma_{\G o}^T\Gamma_{uo}$ are both $0$-$1$ 
valued rotation matrices and 
$\Gamma_{\G}^T\Gamma_{uo} = \Gamma_{\G o}^T\Gamma_{u} = 0$. These facts 
imply that $\widehat{O}$ converges in probability to a $0$-$1$ valued rotation  
matrix. We will assume that $\widehat{O} \overset{p}{\to} I$ without loss 
of generality.

\begin{thm} \label{connection}
The basis matrix $\widehat{\Gamma}_{\widehat{\G}}$ is the output of the 1D 
algorithm with inputs 
$\widehat{M} = \widehat{O}\widehat{\Sigma}_{\upsilon,\upsilon}\widehat{O}^T$ 
and 
$\widehat{U} = \widehat{O}\hat{\upsilon}\hat{\upsilon}^T\widehat{O}^T$ 
at dimension $u$. 
\end{thm}

\begin{proof}
Let $\widehat{M}_2 = \widehat{\Sigma}_{\upsilon,\upsilon}$ and 
$\widehat{U}_2 = \hat{\upsilon}\hat{\upsilon}^T$.
Similar to the proof of \cite[Proposition 6]{cook3}, let 
\begin{align*}
  Q_n(g) &= -n/2\log(g^T\widehat{O}\widehat{M}_2\widehat{O}^Tg) 
    - n/2\log\left[g^T\left\{\widehat{O}\left(\widehat{M}_2 + \widehat{U}_2\right)\widehat{O}^T\right\}^{-1}g\right]
    + n \log(g^T g) \\
  &= -n/2\log(g^T\widehat{O}\widehat{M}_2\widehat{O}^Tg) 
    - n/2\log\left\{g^T\widehat{O}\left(\widehat{M}_2 + \widehat{U}_2\right)^{-1}\widehat{O}^Tg\right\}
    + n \log(g^T\widehat{O}\widehat{O}^T g).  
\end{align*}
Now let $\hat{v}_k$, $k = 1$, $\ldots$, $u$ be the $k$th column of 
$\widehat{\Gamma}_{\G}$ and define $A \in \R^{p\times p}$ to be a matrix of 
$0$'s with $1$'s occupying the first $k$ diagonal entries. Then 
\begin{align*}
  Q_n(\hat{v}_k) 
  &= -n/2\log(\hat{v}_k^T\widehat{O}\widehat{M}_2\widehat{O}^T\hat{v}_k) 
    - n/2\log\left\{\hat{v}_k^T\widehat{O}\left(\widehat{M}_2 
       + \widehat{U}_2\right)^{-1}\widehat{O}^T\hat{v}_k\right\} \\
    &\qquad+ n \log(\hat{v}_k^T\widehat{O}^T\widehat{O}\hat{v}_k) \\
  &= - n/2\log(A_k \widehat{O}_u\widehat{M}_2\widehat{O}_u^TA_k) 
     - n/2\log\left\{A_k\widehat{O}_u\left(\widehat{M}_2 
       + \widehat{U}_2\right)^{-1}\widehat{O}_u^TA_k\right\} \\
    &\qquad+ n \log(A_k\widehat{O}_u^T\widehat{O}_u A_k) \\
  &= - n/2\log(\hat{g}_{uk}^T\widehat{M}_2\hat{g}_{uk}) 
     - n/2\log\left\{\hat{g}_{uk}^T\left(\widehat{M}_2 
       + \widehat{U}_2\right)^{-1}\hat{g}_{uk}\right\} \\
    &\qquad+ n \log(\hat{g}_{uk}^T\hat{g}_{uk})
\end{align*}
where $\hat{g}_{uk}$ is the $k$th column of $\widehat{\Gamma}_u$, 
the output of the 1D algorithm with $\widehat{M}_2$ 
and $\widehat{U}_2$ as inputs. Therefore $\hat{v}_k$ is a maximizer 
of $Q_n(g)$ and this completes the proof.
\end{proof}

Theorem~\ref{connection} in combination with \citet[Theorem 2]{zhangmai} 
allows for us to estimate $\G$ consistently. %through the connection 
%between $\widehat{\Gamma}_{\G}$ and $\widehat{\Gamma}_u$. 
Thus the variability associated with the estimation of $\upenv$ and 
$\G$ decreases as $n\to\infty$. However in practical applications correct 
model selection cannot be guaranteed.
%The envelope estimator $\upenv$ requires that the space $\G$ is selected 
%correctly. In practical settings this is cannot be guaranteed. 
Therefore the envelope estimator of expected Darwinian fitness 
$g(\tauenv)$ has an extra source of variability due to model selection 
uncertainty.





\subsection{Double bootstrap procedure}

The parametric bootstrap procedure in \citet{eck} is used to estimate 
expected Darwinian fitness and its asymptotic variability. This procedure 
accounts for model selection variability by using a specified model selection 
criterion at every iteration of the first level of the bootstrap. Let 
$B$ be the total number of iterations for the first level of our bootstrap 
procedure. Define $g$ to be the map from $\tau$ to the parameterization of 
Darwinian fitness. Then the envelope estimator of estimated expected Darwinian 
fitness using our double bootstrap procedure is 
%the envelope estimator selected by BIC and the LRT of size $\alpha = 0.05$. 
%The steps of this procedure are outlined in Algorithm 1 of the current document. 
%The steps and quantities needed to calculate the bootstrapped estimator of the variability of expected Darwinian fitness incorporating envelope methodology are outlined in this paragraph. 
%When the top-level of our bootstrap procedure in Algorithm 1 has run for a 
%total of $B$ iterations, we obtain the envelope estimator
\begin{equation} \label{genv}
  \genv = \frac{1}{B} \sum_{b=1}^B g(\tauenv^{(b)}).
\end{equation}
The individual estimates of expected Darwinian fitness in \eqref{genv} 
involve model selection which is smoothed out by averaging \citep*{efron}. 
This envelope estimator obtained from the parametric bootstrap in Algorithm 1 
has variability analogous to that in \citet*[equation (4.15)]{efron}. 
As in \citet*{efron}, for $b = 1,...,B$ we define the matrix 
$\textbf{B}^{(b)} \in \R^{K\times p}$ which has rows 
$$
  \hat{\tau}_{\text{env}}^{(b)^{(k)}} 
    - \sum_{k=1}^K\hat{\tau}_{\text{env}}^{(b)^{(k)}}/K
$$ 
and the matrix $C^{(b)} \in \R^{K \times d}$ which has columns 
$$
  g\left(\tau_{\text{env}}^{(b)^{(k)}}\right) 
    - g\left(\tau_{\text{env}}^{(b)}\right).
$$ 
Let $\Delta$ be the asymptotic variability of \eqref{genv}. The variability 
$\Delta$ is estimated by %$g\left(\tau_{\text{env}}^{(b)^{(k)}}\right) - \hat{\mu}_g$
\begin{equation} \label{Delta.boot}
  \widehat{\Delta} = \sum_{b=1}^B\left[\widehat{\text{cov}}^{(b)^T}
    \widehat{V}^{(b)^{-1}}\widehat{\text{cov}}^{(b)}\right] / B
\end{equation}
where 
\begin{equation} \label{cov.hat}
  \widehat{\text{cov}}^{(b)} = \textbf{B}^{(b)^T}C^{(b)}/K
\end{equation}
and 
\begin{equation} \label{V.hat}
  \widehat{V}^{(b)} = \textbf{B}^{(b)^T}\textbf{B}^{(b)}/K.
\end{equation}
The estimator \eqref{Delta.boot} takes into account the 
variability of model selection when estimating the variability of estimated 
expected Darwinian fitness using envelope methodology. The method of maximum 
likelihood estimation does not have the added model selection step that 
envelope estimation has. The bootstrap procedure outlined in Algorithm 1 
efficiently estimates expected Darwinian fitness and accounts for model 
selection variability.



%We now construct the envelope estimator of $\tau$ with respect to 
%$\widehat{G}_1$ as suggested by BIC and the LRT of size $\alpha = 0.05$. 
%The envelope estimator is given by $\tauenv = \Menv^T Y$ where
%$$  
%  \Menv = M \left(\begin{array}{cc}
%   I & 0 \\
%    0 & \widehat{P}_{\widehat{G}_1}
%  \end{array}\right). 
%$$



\begin{table}[h!]
\begin{tabular}{|p{14cm}|}
\hline
\textbf{Algorithm 1: Parametric bootstrap reducing subspace envelope estimator 
  of $\upsilon$} \\
\hline
\begin{enumerate}
\item[1.] Fit the aster model to the data and obtain $\hat{\upsilon}$ and 
  $\widehat{\Sigma}_{\upsilon,\upsilon}$ from the aster model fit. 
\item[2.] Compute the envelope estimator of $\upsilon$ in the original sample, 
  given as $\upenv = P_{\widehat{G}}\hat{\upsilon}$ where $P_{\widehat{G}}$ is 
  computed using reducing subspaces and selected via a model selection 
  criterion of choice.
\item[3.] Perform a parametric bootstrap by generating samples from the 
  distribution of the aster submodel evaluated at 
  $\tauenv = (\hat{\gamma}^T,\upenv^T)^T$. For iteration $b=1,...,B$ of the 
  procedure: 
  \begin{enumerate}
  {\setlength\itemindent{10pt} \item[(3a)] Compute $\hat{\tau}^{(b)}$ and 
    $\widehat{\Sigma}_{\upsilon,\upsilon}^{(b)}$ from the aster model fit 
    of the simulated data.}
  {\setlength\itemindent{10pt} \item[(3b)] Build $P_{\widehat{G}}^{(b)}$ using 
    the indices of $\widehat{\Sigma}_{\upsilon,\upsilon}^{(b)}$ that are 
    selected using the same model selection criterion as Step 2 to construct  
    $\widehat{G}$. }
  {\setlength\itemindent{10pt} \item[(3c)] Compute 
    $\upenv^{(b)} = P_{\hat{\mathcal{E}}}^{(b)}\hat{\upsilon}^{(b)}$, 
    $\tauenv^{(b)} = \left(\hat{\gamma}^{(b)^T},\upenv^{(b)^T}\right)^T$, and 
    $g\left(\hat{\tau}_{\text{env}}^{(b)}\right)$ where $g$ maps $\tau$ to 
    Darwinian fitness.}
  \end{enumerate}
\item[4.] After $B$ steps, the bootstrap estimator of expected Darwinian 
  fitness is \eqref{genv}. This completes the first part of the bootstrap 
  procedure. 
\item[5.] We now proceed with the second level of bootstrapping. 
  At $k = 1$, $\ldots$, $K$, for each $b=1$, $\ldots$, $B$ we:
  \begin{enumerate}
  {\setlength\itemindent{10pt} \item[(5a)] Generate data from the distribution 
    of the aster submodel evaluated at $\hat{\tau}_{\text{env}}^{(b)}$. }
  {\setlength\itemindent{10pt} \item[(5b)] Perform Steps 3a through 3c with 
    respect to the dataset obtained in Step 5a.}
  {\setlength\itemindent{10pt} \item[(5c)] Calculate 
    $\hat{\tau}_{\text{env}}^{(b)^{(k)}}$ and 
    $g\left(\hat{\tau}_{\text{env}}^{(b)^{(k)}}\right)$. }
  \end{enumerate}
\item[6.] Compute the standard deviation \eqref{Delta.boot}.   
\end{enumerate} \\
\hline 
\end{tabular}
\end{table}  


There are two functions in the \texttt{envlpaster} package necessary to 
perform all of the steps in Algorithm 1. The reason for separating this chore 
into two functions is computational. This separation allows users to 
parallelize this bootstrap using \texttt{mclapply}. The first of these 
functions is \texttt{fit.boot.Efron}. This function implements the first level
of the parametric bootstrap procedure given by either Algorithm 1 here 
(reducing subspaces) or the algorithm in \citet*{eck} (1D algorithm). 
This is detailed in Steps 1 through 3c outlined in Algorithm 1. 
This parametric bootstrap generates simulated data from the distribution 
evaluated at an envelope estimator of $\tau$. This procedure accounts for 
model selection variability by using the same model selection criterion at 
every iteration $b = 1$, $\ldots$, $B$, see step 3b of Algorithm 1. 

To use the \texttt{fit.boot.Efron} function, the user must specify whether 
they want to use AIC, BIC, or the LRT of size $\alpha$ as a model selection 
criterion. The user also specifies which method, either the 1D algorithm or 
the reducing subspace approach, is to be used to calculate envelope 
estimators. When one is using a partial envelope model, then this function 
constructs envelope estimators of $\upsilon$ where we write 
$\tau = (\gamma^T,\upsilon^T)^T$ and $\upsilon$ corresponds to aster model 
parameters of interest. In the original sample, the reducing subspaces are 
indices of eigenvectors of $\widehat{\Sigma}_{\upsilon,\upsilon}$ where 
$\widehat{\Sigma}_{\upsilon,\upsilon}$ is the part of $\widehat{\Sigma}$ 
corresponding to our parameters of interest. These indices are specified by 
\texttt{vectors}. When all of the components of $\tau$ are components of 
interest, then we write 
$\widehat{\Sigma}_{\upsilon,\upsilon} = \widehat{\Sigma}$. %When data is generated via the parametric bootstrap, it is the indices of the eigenvectors used to form the reducing subspace (not the original reducing subspaces) that are used to construct envelope estimators constructed using the generated data

We now call on \texttt{fit.boot.Efron} to perform the top level of 
bootstrapping for our example. The calculation above is not performed in this 
script. It is performed in \texttt{ex1-Efronboot.R} and is sourced into this 
document from the \texttt{ex1-Efronboot.RData} file.


\begin{knitrout}
\definecolor{shadecolor}{rgb}{0.969, 0.969, 0.969}\color{fgcolor}\begin{kframe}
\begin{alltt}
\hlkwd{set.seed}\hlstd{(}\hlnum{13}\hlstd{)}
\hlstd{nboot} \hlkwb{<-} \hlnum{100}
\hlstd{blah} \hlkwb{<-} \hlkwd{fit.boot.Efron}\hlstd{(}\hlkwc{model} \hlstd{= m1,} \hlkwc{nboot} \hlstd{= nboot,} \hlkwc{index} \hlstd{= target,}
  \hlkwc{dim} \hlstd{= foo2}\hlopt{$}\hlstd{bic,} \hlkwc{data} \hlstd{= data,} \hlkwc{amat} \hlstd{= amat,} \hlkwc{newdata} \hlstd{= fred,}
  \hlkwc{modmat.new} \hlstd{= modmat.renew,} \hlkwc{renewdata} \hlstd{= renewdata,} \hlkwc{criterion} \hlstd{=} \hlstr{"BIC"}\hlstd{,}
  \hlkwc{method} \hlstd{=} \hlkwd{c}\hlstd{(}\hlstr{"eigen"}\hlstd{),} \hlkwc{quiet} \hlstd{=} \hlnum{FALSE}\hlstd{)}
\end{alltt}
\end{kframe}
\end{knitrout}

The contents stored in \texttt{blah}, and previously computed quantities, are 
all that is needed to perform the second level of bootstrapping necessary for 
the estimation of the variability of estimated expected Darwinian fitness 
incorporating envelope methodology. The function \texttt{secondboot} in the 
\texttt{envlpaster} package performs this second level of bootstrapping. This 
function implements the second level of the parametric bootstrap procedure 
given by either Algorithm 1 or Algorithm 1 in \citet*{eck} with respect to the 
mean-value parameterization. This is detailed in Steps 4 through 5c 
in Algorithm 1. At iteration $b$, this parametric bootstrap generates samples 
from the distribution evaluated at the envelope estimator 
$\hat{\tau}_{env}^{(b)}$. In this case, the selected indices producing the 
reducing subspace which was used to construct the envelope estimator 
$\hat{\tau}_{env}^{(b)}$ are used to construct envelope estimators for the 
generated data. These data are used to estimate the variability of 
$\hat{\tau}_{env}^{(b)}$. 
 
We now call on \texttt{secondboot} to perform the second level of 
bootstrapping for our example. The calculations involving \texttt{secondboot} 
are not performed in this script. These calculations are performed in 
\texttt{ex1-secondboot.R} and are sourced into this document from the 
\texttt{ex1-secondboot.RData} file. The bootstrap sample size for the second 
level of bootstrapping is 500. The script below performs this task. This 
script also calculates the bootstrapped estimator of the variability of 
expected Darwinian fitness incorporating envelope methodology. 
 
\begin{knitrout}
\definecolor{shadecolor}{rgb}{0.969, 0.969, 0.969}\color{fgcolor}\begin{kframe}
\begin{alltt}
\hlkwd{load}\hlstd{(}\hlstr{"ex1-secondboot.RData"}\hlstd{)}
\end{alltt}
\end{kframe}
\end{knitrout}

\begin{knitrout}
\definecolor{shadecolor}{rgb}{0.969, 0.969, 0.969}\color{fgcolor}\begin{kframe}
\begin{alltt}
\hlstd{nboot2} \hlkwb{<-} \hlnum{500}
\hlkwd{set.seed}\hlstd{(}\hlnum{13}\hlstd{)}
\hlstd{internal} \hlkwb{<-} \hlkwa{function}\hlstd{(}\hlkwc{k}\hlstd{)\{}
  \hlkwd{set.seed}\hlstd{(}\hlnum{13}\hlstd{)}
  \hlstd{bar} \hlkwb{<-} \hlkwd{secondboot}\hlstd{(k,} \hlkwc{out} \hlstd{= blah,} \hlkwc{model} \hlstd{= m1,} \hlkwc{data} \hlstd{= data,}
    \hlkwc{nboot2} \hlstd{= nboot2,} \hlkwc{index} \hlstd{= target,} \hlkwc{newdata} \hlstd{= fred,}
    \hlkwc{amat} \hlstd{= amat,} \hlkwc{method} \hlstd{=} \hlstr{"eigen"}\hlstd{)}
  \hlkwd{return}\hlstd{(bar}\hlopt{$}\hlstd{sd.Efron)}
\hlstd{\}}
\end{alltt}
\end{kframe}
\end{knitrout}







\subsection{Results}
\sloppy

The code below calculates \eqref{Delta.boot}. These calculations are performed 
in \texttt{ex1-secondboot.R} and are sourced into this document from the 
\texttt{ex1-secondboot.RData} file. A machine with five available 
cores was used to perform this calculation.

\begin{knitrout}
\definecolor{shadecolor}{rgb}{0.969, 0.969, 0.969}\color{fgcolor}\begin{kframe}
\begin{alltt}
\hlstd{y} \hlkwb{<-} \hlstd{nboot} \hlopt{/} \hlnum{5}
\hlstd{blah.second} \hlkwb{<-} \hlkwd{mclapply}\hlstd{(}\hlnum{1}\hlopt{:}\hlnum{5}\hlstd{,} \hlkwc{mc.cores} \hlstd{=} \hlnum{5}\hlstd{,} \hlkwc{FUN} \hlstd{=} \hlkwa{function}\hlstd{(}\hlkwc{x}\hlstd{)\{}
  \hlstd{y} \hlkwb{<-} \hlstd{nboot} \hlopt{/} \hlnum{5}
  \hlstd{int} \hlkwb{<-} \hlstd{((y}\hlopt{*}\hlstd{(x}\hlopt{-}\hlnum{1}\hlstd{)} \hlopt{+} \hlnum{1}\hlstd{)} \hlopt{:} \hlstd{(y}\hlopt{*}\hlstd{x))}
  \hlstd{out} \hlkwb{<-} \hlkwd{lapply}\hlstd{(int,} \hlkwc{FUN} \hlstd{=  internal)}
  \hlkwd{return}\hlstd{(out)}
  \hlstd{\}}
\hlstd{)}
\end{alltt}
\end{kframe}
\end{knitrout}

\begin{knitrout}
\definecolor{shadecolor}{rgb}{0.969, 0.969, 0.969}\color{fgcolor}\begin{kframe}
\begin{alltt}
\hlstd{sd.Efron.mean} \hlkwb{<-} \hlkwd{rep}\hlstd{(}\hlnum{0}\hlstd{,}\hlnum{144}\hlstd{)}
\hlkwa{for}\hlstd{(k} \hlkwa{in} \hlnum{1}\hlopt{:}\hlnum{5}\hlstd{)\{}
  \hlkwa{for}\hlstd{(j} \hlkwa{in} \hlnum{1}\hlopt{:}\hlstd{y)\{}
    \hlstd{sd.Efron.mean} \hlkwb{<-} \hlstd{sd.Efron.mean} \hlopt{+}
      \hlstd{blah.second[[k]][[j]]} \hlopt{/} \hlstd{nboot}
  \hlstd{\}}
\hlstd{\}}
\end{alltt}
\end{kframe}
\end{knitrout}

The file \texttt{techreport\_eigenboot\_fitboot.RData} contains output 
corresponding to the bootstrapped estimator of the variability of expected 
Darwinian fitness using the MLE. The bootstrap sample size for this particular 
bootstrap was 500.

\begin{knitrout}
\definecolor{shadecolor}{rgb}{0.969, 0.969, 0.969}\color{fgcolor}\begin{kframe}
\begin{alltt}
\hlstd{nboot1} \hlkwb{<-} \hlstd{nboot}
\hlkwd{load}\hlstd{(}\hlstr{"techreport_eigenboot_fitboot.RData"}\hlstd{)}
\end{alltt}
\end{kframe}
\end{knitrout}

We now compute the ratio of standard errors. Values above 1 indicate that the 
envelope estimator of expected Darwinian fitness using Algorithm 1 is less 
variable than the MLE. The code below computes the bootstrapped standard 
errors of estimated expected Darwinian fitness using the MLE. 

\begin{knitrout}
\definecolor{shadecolor}{rgb}{0.969, 0.969, 0.969}\color{fgcolor}\begin{kframe}
\begin{alltt}
\hlstd{a} \hlkwb{<-} \hlstd{bar}\hlopt{$}\hlstd{MLE.boot.out}

\hlstd{data.fit} \hlkwb{<-} \hlkwd{asterdata}\hlstd{(fred, vars, pred,}
  \hlkwc{group} \hlstd{=} \hlkwd{rep}\hlstd{(}\hlnum{0}\hlstd{,}\hlnum{30}\hlstd{), code,}
  \hlstd{families)}

\hlstd{phi.MLE} \hlkwb{<-} \hlstd{origin} \hlopt{+} \hlstd{modmat.renew} \hlopt{%*%} \hlstd{beta}
\hlstd{mu.MLE} \hlkwb{<-} \hlkwd{transformSaturated}\hlstd{(phi.MLE, data.fit,} \hlkwc{from} \hlstd{=} \hlstr{"phi"}\hlstd{,}
  \hlkwc{to} \hlstd{=} \hlstr{"mu"}\hlstd{)}
\hlstd{amat.fit} \hlkwb{<-} \hlkwd{matrix}\hlstd{(amat,} \hlkwc{nrow} \hlstd{=} \hlkwd{length}\hlstd{(mu.MLE))}
\hlstd{fit.MLE} \hlkwb{<-} \hlkwd{t}\hlstd{(amat.fit)} \hlopt{%*%} \hlstd{mu.MLE}
\hlstd{fit.MLE} \hlkwb{<-} \hlkwd{as.vector}\hlstd{(fit.MLE)}
\hlstd{a} \hlkwb{<-} \hlstd{a} \hlopt{-} \hlstd{fit.MLE}

\hlstd{boot.var.MLE} \hlkwb{<-} \hlstd{(a[,} \hlnum{1}\hlstd{]} \hlopt{%o%} \hlstd{a[,} \hlnum{1}\hlstd{])} \hlopt{/} \hlkwd{ncol}\hlstd{(a)}
\hlkwa{for}\hlstd{(j} \hlkwa{in} \hlnum{2}\hlopt{:}\hlkwd{ncol}\hlstd{(a))\{}
  \hlstd{boot.var.MLE} \hlkwb{<-} \hlstd{boot.var.MLE} \hlopt{+} \hlstd{(a[, j]} \hlopt{%o%} \hlstd{a[, j])} \hlopt{/} \hlkwd{ncol}\hlstd{(a)}
\hlstd{\}}

\hlstd{sd.MLE} \hlkwb{<-} \hlkwd{sqrt}\hlstd{(}\hlkwd{diag}\hlstd{(boot.var.MLE))}
\hlstd{ratio} \hlkwb{<-} \hlstd{sd.MLE} \hlopt{/} \hlstd{sd.Efron.mean}
\end{alltt}
\end{kframe}
\end{knitrout}

Output is shown for hypothetical individuals who have estimated expected 
Darwinian fitness greater than 8 where estimation is calculated using the MLE. 
This table has 5 columns. The first two columns correspond to estimation of 
expected Darwinian fitness using envelope methodology and its bootstrapped 
standard error respectively. The next two columns correspond to estimation of 
expected Darwinian fitness using the MLE and its bootstrapped standard error 
respectively. The fifth column displays the ratios of estimated standard 
errors where ratios greater than 1 indicate that the envelope estimator has 
lower estimated variability. We can see that efficiency gains are obtained for 
the displayed output.

\begin{knitrout}
\definecolor{shadecolor}{rgb}{0.969, 0.969, 0.969}\color{fgcolor}\begin{kframe}
\begin{alltt}
\hlkwd{load}\hlstd{(}\hlstr{"techreport_eigenboot_fitboot.RData"}\hlstd{)}
\hlstd{nboot} \hlkwb{<-} \hlstd{nboot1}
\end{alltt}
\end{kframe}
\end{knitrout}

\begin{knitrout}
\definecolor{shadecolor}{rgb}{0.969, 0.969, 0.969}\color{fgcolor}\begin{kframe}
\begin{alltt}
\hlstd{tab} \hlkwb{<-} \hlkwd{matrix}\hlstd{(}\hlnum{0}\hlstd{,} \hlkwc{nrow} \hlstd{=} \hlnum{144}\hlstd{,} \hlkwc{ncol} \hlstd{=} \hlnum{5}\hlstd{)}
\hlstd{tab[,} \hlnum{1}\hlstd{]} \hlkwb{<-} \hlkwd{rowSums}\hlstd{(blah}\hlopt{$}\hlstd{env.boot.out)} \hlopt{/} \hlstd{nboot}
\hlstd{tab[,} \hlnum{2}\hlstd{]} \hlkwb{<-} \hlstd{sd.Efron.mean}
\hlstd{tab[,} \hlnum{3}\hlstd{]} \hlkwb{<-} \hlstd{fit.MLE}
\hlstd{tab[,} \hlnum{4}\hlstd{]} \hlkwb{<-} \hlstd{sd.MLE}
\hlstd{tab[,} \hlnum{5}\hlstd{]} \hlkwb{<-} \hlstd{ratio}
\hlkwd{colnames}\hlstd{(tab)} \hlkwb{<-} \hlkwd{c}\hlstd{(}\hlstr{"env"}\hlstd{,}\hlstr{"se(env)"}\hlstd{,}\hlstr{"MLE"}\hlstd{,}\hlstr{"se(MLE)"}\hlstd{,}\hlstr{"ratio"}\hlstd{)}
\hlstd{tab[tab[,} \hlnum{3}\hlstd{]} \hlopt{>} \hlnum{8}\hlstd{, ]}
\end{alltt}
\begin{verbatim}
##           env   se(env)      MLE   se(MLE)    ratio
## [1,] 8.556409 0.1742400 8.674101 0.2597591 1.490812
## [2,] 9.013616 0.1105248 8.938662 0.1350701 1.222079
## [3,] 7.817013 0.4135980 8.030673 0.4420406 1.068769
## [4,] 9.174295 0.1627223 9.193001 0.1700020 1.044737
## [5,] 9.018329 0.1126568 9.126989 0.1276560 1.133141
## [6,] 8.611867 0.1624729 8.541299 0.2776820 1.709097
## [7,] 7.761383 0.2153252 8.109067 0.3302726 1.533831
\end{verbatim}
\end{kframe}
\end{knitrout}

We plot the ratio of estimated standard errors in 
Figures~\ref{contour-plot} and~\ref{heatmap}. We can see that efficiency 
gains are obtained when using envelope methodology for hypothetical 
individuals that overlap the majority of the original data. 
%The gains are less substantial than those found when using the indices 
%suggested by BIC and the LRT of size $\alpha = 0.05$.


\begin{knitrout}
\definecolor{shadecolor}{rgb}{0.969, 0.969, 0.969}\color{fgcolor}\begin{kframe}
\begin{alltt}
\hlstd{z1} \hlkwb{<-} \hlstd{simdata30nodes}\hlopt{$}\hlstd{z1}
\hlstd{z2} \hlkwb{<-} \hlstd{simdata30nodes}\hlopt{$}\hlstd{z2}
\hlstd{max.ind} \hlkwb{<-} \hlkwd{which}\hlstd{(tab[,}\hlnum{1}\hlstd{]} \hlopt{==} \hlkwd{max}\hlstd{(tab[,}\hlnum{1}\hlstd{]))}
\hlstd{x.max} \hlkwb{=} \hlstd{fred[max.ind,}\hlnum{1}\hlstd{]}
\hlstd{y.max} \hlkwb{=} \hlstd{fred[max.ind,}\hlnum{2}\hlstd{]}
\hlkwd{plot}\hlstd{(z1,z2,} \hlkwc{pch} \hlstd{=} \hlnum{20}\hlstd{,} \hlkwc{col} \hlstd{=} \hlstr{"grey"}\hlstd{)}
\hlkwd{points}\hlstd{(}\hlkwc{x} \hlstd{= x.max,} \hlkwc{y} \hlstd{= y.max,} \hlkwc{col} \hlstd{=} \hlstr{"red"}\hlstd{,} \hlkwc{pch} \hlstd{=} \hlnum{20}\hlstd{)}
\hlstd{max.ind.env} \hlkwb{<-} \hlkwd{which}\hlstd{(tab[,}\hlnum{3}\hlstd{]} \hlopt{==} \hlkwd{max}\hlstd{(tab[,}\hlnum{3}\hlstd{]))}
\hlstd{x.max.env} \hlkwb{=} \hlstd{fred[max.ind.env,}\hlnum{1}\hlstd{]}
\hlstd{y.max.env} \hlkwb{=} \hlstd{fred[max.ind.env,}\hlnum{2}\hlstd{]}
\hlkwd{points}\hlstd{(}\hlkwc{x} \hlstd{= x.max.env,} \hlkwc{y} \hlstd{= y.max.env,} \hlkwc{col} \hlstd{=} \hlstr{"blue"}\hlstd{,} \hlkwc{pch} \hlstd{=} \hlnum{20}\hlstd{)}
\hlkwd{contour}\hlstd{(z1.cand,z2.cand,} \hlkwd{t}\hlstd{(}\hlkwd{matrix}\hlstd{(tab[,}\hlnum{5}\hlstd{],}\hlkwc{nrow} \hlstd{= nx)),}
  \hlkwc{levels} \hlstd{=} \hlkwd{c}\hlstd{(}\hlnum{0.9}\hlstd{,}\hlnum{1}\hlstd{,}\hlnum{1.5}\hlstd{,}\hlnum{2}\hlstd{,}\hlnum{3}\hlstd{),} \hlkwc{labcex} \hlstd{=} \hlnum{1}\hlstd{,} \hlkwc{add}\hlstd{=}\hlnum{TRUE}\hlstd{)}
\end{alltt}
\end{kframe}
\end{knitrout}


\begin{figure}[h!]
\begin{center}
\begin{knitrout}
\definecolor{shadecolor}{rgb}{0.969, 0.969, 0.969}\color{fgcolor}
\includegraphics[width=\maxwidth]{figure/ex1-ratios-Efron-1-1} 

\end{knitrout}
\end{center}
\caption{Contour plot for the ratios of $se\left(h(\hat{\tau})\right)$ to $se\left(h(\hat{\tau}_{\text{env}})\right)$. Ratios greater than 1 indicate efficiency gains using envelope methodology. The point in blue corresponds to the highest estimated expected Darwinian fitness value using envelope methodology and the MLE. }
\label{contour-plot}
\end{figure}


\newpage


\begin{knitrout}
\definecolor{shadecolor}{rgb}{0.969, 0.969, 0.969}\color{fgcolor}\begin{kframe}
\begin{alltt}
\hlkwd{filled.contour}\hlstd{(z1.cand,z2.cand,} \hlkwd{t}\hlstd{(}\hlkwd{matrix}\hlstd{(tab[,}\hlnum{5}\hlstd{],}\hlkwc{nrow} \hlstd{= nx)),}
  \hlkwc{levels} \hlstd{=} \hlkwd{c}\hlstd{(}\hlnum{0}\hlstd{,}\hlnum{0.9}\hlstd{,}\hlnum{1}\hlstd{,}\hlnum{1.5}\hlstd{,}\hlnum{2}\hlstd{,}\hlnum{5}\hlstd{))}
\end{alltt}
\end{kframe}
\end{knitrout}


\begin{figure}[h!]
\begin{center}
\begin{knitrout}
\definecolor{shadecolor}{rgb}{0.969, 0.969, 0.969}\color{fgcolor}
\includegraphics[width=\maxwidth]{figure/ex1-ratios-Efron-2-1} 

\end{knitrout}
\end{center}
\caption{Heat map for the ratios of $se\left(h(\hat{\tau})\right)$ to $se\left(h(\hat{\tau}_{\text{env}})\right)$. Ratios greater than 1 indicate efficiency gains using envelope methodology. }
\label{heatmap}
\end{figure}






\newpage
\subsection{Results with AIC in place of BIC}

We now perform our methods with AIC as the model selection criterion used 
in Steps 2 and 3b of Algorithm 1. Output is shown for hypothetical individuals 
who have estimated expected Darwinian fitness greater than 8 where estimation 
is calculated using the MLE.

\begin{knitrout}
\definecolor{shadecolor}{rgb}{0.969, 0.969, 0.969}\color{fgcolor}\begin{kframe}
\begin{alltt}
\hlkwd{load}\hlstd{(}\hlstr{"ex1-AIC.RData"}\hlstd{)}
\end{alltt}
\end{kframe}
\end{knitrout}

\begin{knitrout}
\definecolor{shadecolor}{rgb}{0.969, 0.969, 0.969}\color{fgcolor}\begin{kframe}
\begin{alltt}
\hlstd{tab2} \hlkwb{<-} \hlkwd{matrix}\hlstd{(}\hlnum{0}\hlstd{,} \hlkwc{nrow} \hlstd{=} \hlnum{144}\hlstd{,} \hlkwc{ncol} \hlstd{=} \hlnum{5}\hlstd{)}
\hlkwd{colnames}\hlstd{(tab2)} \hlkwb{<-} \hlkwd{c}\hlstd{(}\hlstr{"env"}\hlstd{,}\hlstr{"se(env)"}\hlstd{,}\hlstr{"MLE"}\hlstd{,}\hlstr{"se(MLE)"}\hlstd{,}\hlstr{"ratio"}\hlstd{)}
\hlstd{tab2[,} \hlnum{1}\hlstd{]} \hlkwb{<-} \hlkwd{rowSums}\hlstd{(blah}\hlopt{$}\hlstd{env.boot.out)} \hlopt{/} \hlstd{nboot}
\hlstd{tab2[,} \hlnum{2}\hlstd{]} \hlkwb{<-} \hlstd{sd.Efron.mean}
\hlstd{tab2[,} \hlnum{3}\hlstd{]} \hlkwb{<-} \hlstd{fit.MLE}
\hlstd{tab2[,} \hlnum{4}\hlstd{]} \hlkwb{<-} \hlstd{sd.MLE}
\hlstd{tab2[,} \hlnum{5}\hlstd{]} \hlkwb{<-} \hlstd{ratio}
\hlstd{tab2[tab2[,} \hlnum{3}\hlstd{]} \hlopt{>} \hlnum{8}\hlstd{, ]}
\end{alltt}
\begin{verbatim}
##           env   se(env)      MLE   se(MLE)    ratio
## [1,] 8.749144 0.1964036 8.674101 0.2599230 1.323412
## [2,] 9.025514 0.1116299 8.938662 0.1351961 1.211109
## [3,] 7.940958 0.4204229 8.030673 0.4424621 1.052422
## [4,] 9.159294 0.1638821 9.193001 0.1701628 1.038324
## [5,] 9.047524 0.1138559 9.126989 0.1277835 1.122326
## [6,] 8.393870 0.2005390 8.541299 0.2778990 1.385760
## [7,] 7.833197 0.2152651 8.109067 0.3305165 1.535393
\end{verbatim}
\end{kframe}
\end{knitrout}

The ratios of estimated standard errors are plotted in 
Figures~\ref{contour-plot-AIC} and~\ref{heatmap-AIC}. We can see in these 
figures and in the table above that efficiency gains are obtained using 
envelope methodology. The gains are less substantial than those found when 
using BIC. 


\begin{knitrout}
\definecolor{shadecolor}{rgb}{0.969, 0.969, 0.969}\color{fgcolor}\begin{kframe}
\begin{alltt}
\hlstd{z1} \hlkwb{<-} \hlstd{simdata30nodes}\hlopt{$}\hlstd{z1}
\hlstd{z2} \hlkwb{<-} \hlstd{simdata30nodes}\hlopt{$}\hlstd{z2}
\hlstd{max.ind} \hlkwb{<-} \hlkwd{which}\hlstd{(tab[,}\hlnum{1}\hlstd{]} \hlopt{==} \hlkwd{max}\hlstd{(tab[,}\hlnum{1}\hlstd{]))}
\hlstd{x.max} \hlkwb{=} \hlstd{fred[max.ind,}\hlnum{1}\hlstd{]}
\hlstd{y.max} \hlkwb{=} \hlstd{fred[max.ind,}\hlnum{2}\hlstd{]}
\hlkwd{plot}\hlstd{(z1,z2,} \hlkwc{pch} \hlstd{=} \hlnum{20}\hlstd{,} \hlkwc{col} \hlstd{=} \hlstr{"grey"}\hlstd{)}
\hlkwd{points}\hlstd{(}\hlkwc{x} \hlstd{= x.max,} \hlkwc{y} \hlstd{= y.max,} \hlkwc{col} \hlstd{=} \hlstr{"red"}\hlstd{,} \hlkwc{pch} \hlstd{=} \hlnum{20}\hlstd{)}
\hlstd{max.ind.env} \hlkwb{<-} \hlkwd{which}\hlstd{(tab[,}\hlnum{3}\hlstd{]} \hlopt{==} \hlkwd{max}\hlstd{(tab[,}\hlnum{3}\hlstd{]))}
\hlstd{x.max.env} \hlkwb{=} \hlstd{fred[max.ind.env,}\hlnum{1}\hlstd{]}
\hlstd{y.max.env} \hlkwb{=} \hlstd{fred[max.ind.env,}\hlnum{2}\hlstd{]}
\hlkwd{points}\hlstd{(}\hlkwc{x} \hlstd{= x.max.env,} \hlkwc{y} \hlstd{= y.max.env,} \hlkwc{col} \hlstd{=} \hlstr{"blue"}\hlstd{,} \hlkwc{pch} \hlstd{=} \hlnum{20}\hlstd{)}
\hlkwd{contour}\hlstd{(z1.cand,z2.cand,} \hlkwd{t}\hlstd{(}\hlkwd{matrix}\hlstd{(tab[,}\hlnum{5}\hlstd{],}\hlkwc{nrow} \hlstd{= nx)),}
  \hlkwc{levels} \hlstd{=} \hlkwd{c}\hlstd{(}\hlnum{0.9}\hlstd{,}\hlnum{1}\hlstd{,}\hlnum{1.5}\hlstd{,}\hlnum{2}\hlstd{,}\hlnum{3}\hlstd{),} \hlkwc{labcex} \hlstd{=} \hlnum{1}\hlstd{,} \hlkwc{add}\hlstd{=}\hlnum{TRUE}\hlstd{)}
\end{alltt}
\end{kframe}
\end{knitrout}


\begin{figure}[h!]
\begin{center}
\begin{knitrout}
\definecolor{shadecolor}{rgb}{0.969, 0.969, 0.969}\color{fgcolor}
\includegraphics[width=\maxwidth]{figure/ex1-AIC-1-1} 

\end{knitrout}
\end{center}
\caption{Contour plot for the ratios of $se\left(h(\hat{\tau})\right)$ to $se\left(h(\hat{\tau}_{\text{env}})\right)$ using AIC. Ratios greater than 1 indicate efficiency gains using envelope methodology. The point in blue corresponds to the highest estimated expected Darwinian fitness value using envelope methodology and MLE. }
\label{contour-plot-AIC}
\end{figure}


\newpage


\begin{knitrout}
\definecolor{shadecolor}{rgb}{0.969, 0.969, 0.969}\color{fgcolor}\begin{kframe}
\begin{alltt}
\hlkwd{filled.contour}\hlstd{(z1.cand,z2.cand,} \hlkwd{t}\hlstd{(}\hlkwd{matrix}\hlstd{(tab[,}\hlnum{5}\hlstd{],}\hlkwc{nrow} \hlstd{= nx)),}
  \hlkwc{levels} \hlstd{=} \hlkwd{c}\hlstd{(}\hlnum{0}\hlstd{,}\hlnum{0.9}\hlstd{,}\hlnum{1}\hlstd{,}\hlnum{1.1}\hlstd{,}\hlnum{1.25}\hlstd{,}\hlnum{1.5}\hlstd{,}\hlnum{2}\hlstd{))}
\end{alltt}
\end{kframe}
\end{knitrout}



\newpage


\begin{figure}[h!]
\begin{center}
\begin{knitrout}
\definecolor{shadecolor}{rgb}{0.969, 0.969, 0.969}\color{fgcolor}
\includegraphics[width=\maxwidth]{figure/ex1-AIC-2-1} 

\end{knitrout}
\end{center}
\caption{Heat map for the ratios of $se\left(h(\hat{\tau})\right)$ to $se\left(h(\hat{\tau}_{\text{env}})\right)$ using AIC. Ratios greater than 1 indicate efficiency gains using envelope methodology. }
\label{heatmap-AIC}
\end{figure}




















\newpage















\begin{figure}
\begin{center}
\setlength{\unitlength}{0.45 in}
\thicklines
\begin{picture}(4.0,0.75)(-0.1,-0.3)
  \put(0,0){\makebox(0,0){$1_{\hphantom{0}}$}}
  \put(1,0){\makebox(0,0){$Y_1$}}
  \put(2,0){\makebox(0,0){$Y_2$}}
  \multiput(0.35,0)(1,0){2}{\vector(1,0){0.3}}
\end{picture}
\end{center}
\caption{Graphical structure of the aster model for the data in Example 2. The first arrow corresponds to survival which is a Bernoulli random variable. The second arrow corresponds to reproduction count conditional on survival which is a zero-truncated Poisson random variable.}
\label{astergraph}
\end{figure}







\section{Replicating Example 2}

We now reproduce the second example in \citet*{eck}. The dataset comes from 
\citet*{lowry} and the study in which the data is obtained investigates the 
role of chromosomal inversions in adaptation and speciation. Phenotypic traits 
and covariates are recorded for 2313 yellow monkeyflowers, 
\emph{Mimulus guttatus}. The lifecycle of the individual \emph{M. guttatus} 
flowers is depicted in Figure~\ref{astergraph}. Much of the details necessary 
to the reproduction of this example are similar to those in the first example. 
The data and corresponding \texttt{.RData} file are loaded in. The data is 
then manipulated so that \texttt{aster} software can be used and then the 
quantities needed for the bootstrap functions \texttt{fit.Efron} and 
\texttt{secondboot} are computed.

The code below obtains the relevant quantities needed to fit the aster model 
corresponding to this data and manipulates the dataset so that aster models 
can be fit.


\begin{knitrout}
\definecolor{shadecolor}{rgb}{0.969, 0.969, 0.969}\color{fgcolor}\begin{kframe}
\begin{alltt}
\hlkwd{library}\hlstd{(aster2)}
\hlkwd{library}\hlstd{(envlpaster)}
\hlkwd{data}\hlstd{(Mguttatus)}
\end{alltt}
\end{kframe}
\end{knitrout}

\begin{knitrout}
\definecolor{shadecolor}{rgb}{0.969, 0.969, 0.969}\color{fgcolor}\begin{kframe}
\begin{alltt}
\hlstd{data} \hlkwb{<-} \hlstd{Mguttatus}
\hlstd{redata} \hlkwb{<-} \hlstd{Mguttatus.redata}
\hlstd{vars} \hlkwb{<-} \hlstd{quantities}\hlopt{$}\hlstd{vars}
\hlstd{pred} \hlkwb{<-} \hlstd{quantities}\hlopt{$}\hlstd{pred}
\hlstd{group} \hlkwb{<-} \hlstd{quantities}\hlopt{$}\hlstd{group}
\hlstd{code} \hlkwb{<-} \hlstd{quantities}\hlopt{$}\hlstd{code}
\hlstd{fam} \hlkwb{<-} \hlstd{quantities}\hlopt{$}\hlstd{fam}
\hlstd{nnode} \hlkwb{<-} \hlkwd{length}\hlstd{(vars)}
\hlstd{n} \hlkwb{<-} \hlkwd{nrow}\hlstd{(redata)} \hlopt{/} \hlstd{nnode}
\hlstd{families} \hlkwb{<-} \hlstd{quantities}\hlopt{$}\hlstd{families}
\hlstd{root} \hlkwb{<-} \hlstd{redata}\hlopt{$}\hlstd{root}
\hlstd{fit} \hlkwb{<-} \hlstd{redata}\hlopt{$}\hlstd{fit}
\hlstd{varvar} \hlkwb{<-} \hlstd{quantities}\hlopt{$}\hlstd{varvar}
\hlstd{idvar} \hlkwb{<-} \hlstd{quantities}\hlopt{$}\hlstd{idvar}
\end{alltt}
\end{kframe}
\end{knitrout}


The \texttt{.RData} file containing the calculations done in this example is 
loaded in.


\begin{knitrout}
\definecolor{shadecolor}{rgb}{0.969, 0.969, 0.969}\color{fgcolor}\begin{kframe}
\begin{alltt}
\hlkwd{load}\hlstd{(}\hlstr{"main-all-terms-DavidLowrydata.RData"}\hlstd{)}
\end{alltt}
\end{kframe}
\end{knitrout}


We fit the main effects model. This model includes genetic background, site, 
inversion, and ecotype as predictors.


\begin{knitrout}
\definecolor{shadecolor}{rgb}{0.969, 0.969, 0.969}\color{fgcolor}\begin{kframe}
\begin{alltt}
\hlstd{m.main} \hlkwb{<-} \hlkwd{aster}\hlstd{(resp} \hlopt{~} \hlstd{varb} \hlopt{+} \hlstd{fit}\hlopt{:}\hlstd{(gen_bac} \hlopt{+} \hlstd{site} \hlopt{+}
  \hlstd{inversion} \hlopt{+} \hlstd{type), pred, fam,} \hlkwc{varvar} \hlstd{= varvar,}
  \hlkwc{idvar} \hlstd{= idvar,} \hlkwc{data} \hlstd{= redata,} \hlkwc{root} \hlstd{= root)}
\end{alltt}
\end{kframe}
\end{knitrout}


There are eight predictors comprising $\tau$ and we partition $\tau$ into 
$(\gamma^T, \upsilon^T)^T$ where $\gamma \in \R^2$ are nuisance parameters 
and $\upsilon \in \R^6$ are relevant to the estimation of expected Darwinian 
fitness. In this example $\Sigma_{\upsilon,\upsilon}$ is the lower 6 by 6 
block of the Fisher information matrix $\Sigma$. These quantities, as well as 
$\beta$ are extracted from the aster model fit.


\begin{knitrout}
\definecolor{shadecolor}{rgb}{0.969, 0.969, 0.969}\color{fgcolor}\begin{kframe}
\begin{alltt}
\hlstd{target} \hlkwb{<-} \hlkwd{which}\hlstd{(}\hlkwd{grepl}\hlstd{(}\hlstr{"fit"}\hlstd{,} \hlkwd{names}\hlstd{(m.main}\hlopt{$}\hlstd{coef)))}
\hlstd{avar.targ} \hlkwb{<-} \hlstd{m.main}\hlopt{$}\hlstd{fisher[target, target]}
\hlstd{modmat.mat} \hlkwb{<-} \hlkwd{matrix}\hlstd{(m.main}\hlopt{$}\hlstd{modmat,} \hlkwc{nrow} \hlstd{= n}\hlopt{*}\hlstd{nnode)}
\hlstd{tau} \hlkwb{<-} \hlkwd{crossprod}\hlstd{(modmat.mat, redata}\hlopt{$}\hlstd{resp)}
\hlstd{tau.targ} \hlkwb{<-} \hlstd{tau[target]}
\hlstd{beta} \hlkwb{<-} \hlstd{m.main}\hlopt{$}\hlstd{coef}
\hlstd{p} \hlkwb{<-} \hlkwd{length}\hlstd{(beta)}
\end{alltt}
\end{kframe}
\end{knitrout}


We now call on \texttt{selection} to see which reducing subspace is suggested 
for envelope estimation.


\begin{knitrout}
\definecolor{shadecolor}{rgb}{0.969, 0.969, 0.969}\color{fgcolor}\begin{kframe}
\begin{alltt}
\hlstd{foo} \hlkwb{<-} \hlkwd{selection}\hlstd{(tau, target, m.main, data,} \hlkwc{type} \hlstd{=} \hlstr{"mean-value"}\hlstd{,}
  \hlkwc{method} \hlstd{=} \hlstr{"eigen"}\hlstd{,} \hlkwc{alpha} \hlstd{=} \hlnum{0.01}\hlstd{)}
\end{alltt}
\end{kframe}
\end{knitrout}


All of the selection criteria are in alignment, and they suggest using the 
reducing subspace that is the sum of the first, second, third, and sixth 
eigenspaces of $\widehat{\Sigma}_{\upsilon,\upsilon}$.


\begin{knitrout}
\definecolor{shadecolor}{rgb}{0.969, 0.969, 0.969}\color{fgcolor}\begin{kframe}
\begin{alltt}
\hlstd{foo}\hlopt{$}\hlstd{bic; foo}\hlopt{$}\hlstd{aic; foo}\hlopt{$}\hlstd{LRT}
\end{alltt}
\begin{verbatim}
## [1] 1 2 3 6
## [1] 1 2 3 6
## [1] 1 2 3 6
\end{verbatim}
\end{kframe}
\end{knitrout}


The variability of estimated expected Darwinian fitness is now estimated using 
using both \texttt{fit.boot.Efron} and \texttt{secondboot} to perform the 
bootstrap procedure outlined in Algorithm 1. BIC is used in Steps 2 and 3b of 
Algorithm 1 to select the reducing subspaces used to construct envelope 
estimators. First, the hypothetical \emph{M. guttatus} flowers in which 
expected Darwinian fitness is estimated are created. All variables in our 
original aster model are factor variables, therefore each hypothetical 
\emph{M. guttatus} flower possesses a unique factor-level combination.


\begin{knitrout}
\definecolor{shadecolor}{rgb}{0.969, 0.969, 0.969}\color{fgcolor}\begin{kframe}
\begin{alltt}
\hlstd{a} \hlkwb{<-} \hlkwd{levels}\hlstd{(test}\hlopt{$}\hlstd{gen_bac)}
\hlstd{b} \hlkwb{<-} \hlkwd{levels}\hlstd{(test}\hlopt{$}\hlstd{site)}
\hlstd{c} \hlkwb{<-} \hlkwd{levels}\hlstd{(test}\hlopt{$}\hlstd{inversion)}
\hlstd{d} \hlkwb{<-} \hlkwd{levels}\hlstd{(test}\hlopt{$}\hlstd{type)}
\hlstd{fred} \hlkwb{<-} \hlkwd{expand.grid}\hlstd{(}\hlkwc{a} \hlstd{= a,} \hlkwc{b} \hlstd{= b,} \hlkwc{c} \hlstd{= c,} \hlkwc{d} \hlstd{= d)}
\hlkwd{colnames}\hlstd{(fred)} \hlkwb{<-} \hlkwd{c}\hlstd{(}\hlstr{"gen_bac"}\hlstd{,}\hlstr{"site"}\hlstd{,}\hlstr{"inversion"}\hlstd{,}\hlstr{"type"}\hlstd{)}
\hlstd{fred}\hlopt{$}\hlstd{sur_flw} \hlkwb{<-} \hlnum{1}
\hlstd{fred}\hlopt{$}\hlstd{flws} \hlkwb{<-} \hlnum{1}
\end{alltt}
\end{kframe}
\end{knitrout}


The dataset in long form and the model matrix corresponding to the 
main-effects model with these hypothetical \emph{M. guttatus} flowers is now 
created.


\begin{knitrout}
\definecolor{shadecolor}{rgb}{0.969, 0.969, 0.969}\color{fgcolor}\begin{kframe}
\begin{alltt}
\hlstd{renewdata} \hlkwb{<-} \hlkwd{reshape}\hlstd{(fred,} \hlkwc{varying} \hlstd{=} \hlkwd{list}\hlstd{(vars),}
  \hlkwc{direction} \hlstd{=} \hlstr{"long"}\hlstd{,} \hlkwc{timevar} \hlstd{=} \hlstr{"varb"}\hlstd{,}
  \hlkwc{times} \hlstd{=} \hlkwd{as.factor}\hlstd{(vars),} \hlkwc{v.names} \hlstd{=} \hlstr{"resp"}\hlstd{)}
\hlstd{fit.renew} \hlkwb{<-} \hlkwd{as.numeric}\hlstd{(}\hlkwd{grepl}\hlstd{(}\hlstr{"flws"}\hlstd{, renewdata}\hlopt{$}\hlstd{varb))}
\hlstd{renewdata}\hlopt{$}\hlstd{fit} \hlkwb{<-} \hlstd{fit.renew}
\hlstd{renewdata}\hlopt{$}\hlstd{root} \hlkwb{<-} \hlnum{1}
\hlstd{npop} \hlkwb{<-} \hlkwd{nrow}\hlstd{(fred)}

\hlstd{modmat.renew} \hlkwb{<-} \hlkwd{model.matrix}\hlstd{(m.main}\hlopt{$}\hlstd{formula,} \hlkwc{data} \hlstd{= renewdata)}
\hlstd{modmat.renew[,} \hlnum{2}\hlstd{]} \hlkwb{<-} \hlnum{1} \hlopt{-} \hlstd{modmat.renew[,} \hlnum{2}\hlstd{]}
\hlstd{index} \hlkwb{<-} \hlopt{!}\hlstd{(}\hlkwd{colnames}\hlstd{(modmat.renew)} \hlopt{%in%} \hlstd{m.main}\hlopt{$}\hlstd{dropped)}
\hlstd{modmat.renew} \hlkwb{<-} \hlstd{modmat.renew[, index]}
\end{alltt}
\end{kframe}
\end{knitrout}


Next, \texttt{amat} needs to be specified. \texttt{amat} specifies which nodes 
in the graphical structure for the hypothetical \emph{M. guttatus} flowers 
correspond to Darwinian fitness nodes. \texttt{amat} can either be a 
three-dimensional array or a matrix. We construct a three-dimensional array as 
in \citet*{geyer4}.


\begin{knitrout}
\definecolor{shadecolor}{rgb}{0.969, 0.969, 0.969}\color{fgcolor}\begin{kframe}
\begin{alltt}
\hlstd{amat} \hlkwb{<-} \hlkwd{matrix}\hlstd{(}\hlnum{0}\hlstd{,} \hlkwc{nrow} \hlstd{= npop,} \hlkwc{ncol} \hlstd{= nnode} \hlopt{*} \hlstd{npop)}
\hlkwa{for}\hlstd{(j} \hlkwa{in} \hlnum{1}\hlopt{:}\hlstd{npop) amat[j, npop} \hlopt{+} \hlstd{j]} \hlkwb{<-} \hlstd{fit.renew[npop} \hlopt{+} \hlstd{j]}
\end{alltt}
\end{kframe}
\end{knitrout}


The \texttt{fit.boot.Efron} function performs the top level of bootstrapping 
for our example. The calculation below is not performed in this script. The 
calculation is performed in \texttt{main-all-terms-DavidLowrydata.R} and is 
sourced into this document from the 
\texttt{main-all-terms-DavidLowrydata.RData} file.


\begin{knitrout}
\definecolor{shadecolor}{rgb}{0.969, 0.969, 0.969}\color{fgcolor}\begin{kframe}
\begin{alltt}
\hlkwd{set.seed}\hlstd{(}\hlnum{13}\hlstd{)}
\hlstd{nboot} \hlkwb{<-} \hlnum{100}
\hlstd{blah} \hlkwb{<-} \hlkwd{fit.boot.Efron}\hlstd{(}\hlkwc{model} \hlstd{= m.main,} \hlkwc{nboot} \hlstd{= nboot,} \hlkwc{index} \hlstd{= target,}
  \hlkwc{vectors} \hlstd{= foo}\hlopt{$}\hlstd{bic,} \hlkwc{data} \hlstd{= data,} \hlkwc{amat} \hlstd{= amat,} \hlkwc{newdata} \hlstd{= fred,}
  \hlkwc{modmat.new} \hlstd{= modmat.renew,} \hlkwc{renewdata} \hlstd{= renewdata,}
  \hlkwc{criterion} \hlstd{=} \hlstr{"BIC"}\hlstd{,} \hlkwc{method} \hlstd{=} \hlstr{"eigen"}\hlstd{,} \hlkwc{quiet} \hlstd{=} \hlnum{FALSE}\hlstd{)}
\end{alltt}
\end{kframe}
\end{knitrout}


We now call on \texttt{secondboot} to perform the second level of 
bootstrapping for our example. The calculations involving \texttt{secondboot} 
are not performed in this script. These calculations are also performed in 
\texttt{main-all-terms-DavidLowrydata.R} and are sourced into this document 
from the \texttt{main-all-terms-DavidLowrydata.RData} file. The bootstrap 
sample size for the second level of bootstrapping is 500. The script below 
performs this task. This script also calculates the bootstrapped estimator of 
the variability of expected Darwinian fitness incorporating envelope 
methodology.


\begin{knitrout}
\definecolor{shadecolor}{rgb}{0.969, 0.969, 0.969}\color{fgcolor}\begin{kframe}
\begin{alltt}
\hlstd{nboot2} \hlkwb{<-} \hlnum{500}
\hlstd{internal} \hlkwb{<-} \hlkwa{function}\hlstd{(}\hlkwc{k}\hlstd{)\{}
  \hlkwd{set.seed}\hlstd{(}\hlnum{13}\hlstd{)}
  \hlstd{bar} \hlkwb{<-} \hlkwd{secondboot}\hlstd{(k,} \hlkwc{out} \hlstd{= blah,} \hlkwc{model} \hlstd{= m.main,} \hlkwc{data} \hlstd{= data,}
    \hlkwc{nboot2} \hlstd{= nboot2,} \hlkwc{index} \hlstd{= target,} \hlkwc{newdata} \hlstd{= fred,}
    \hlkwc{amat} \hlstd{= amat,} \hlkwc{method} \hlstd{=} \hlstr{"eigen"}\hlstd{)}
  \hlkwd{return}\hlstd{(bar}\hlopt{$}\hlstd{sd.Efron)}
\hlstd{\}}

\hlcom{## 5 cores are used to perform the second level of bootstrapping}
\hlstd{ncores} \hlkwb{<-} \hlnum{5}\hlstd{; nsamp} \hlkwb{<-} \hlstd{nboot} \hlopt{/} \hlstd{ncores}
\hlstd{blah.second} \hlkwb{<-} \hlkwd{mclapply}\hlstd{(}\hlnum{1}\hlopt{:}\hlstd{ncores,} \hlkwc{mc.cores} \hlstd{= ncores,}
  \hlkwc{FUN} \hlstd{=} \hlkwa{function}\hlstd{(}\hlkwc{x}\hlstd{)\{}
    \hlstd{int} \hlkwb{<-} \hlstd{((nsamp}\hlopt{*}\hlstd{(x}\hlopt{-}\hlnum{1}\hlstd{)} \hlopt{+} \hlnum{1}\hlstd{)} \hlopt{:} \hlstd{(nsamp}\hlopt{*}\hlstd{x))}
    \hlstd{out} \hlkwb{<-} \hlkwd{lapply}\hlstd{(int,} \hlkwc{FUN} \hlstd{= internal)}
    \hlkwd{return}\hlstd{(out)}
  \hlstd{\}}
\hlstd{)}

\hlcom{## estimated standard error using }
\hlstd{sd.Efron.mean} \hlkwb{<-} \hlkwd{rep}\hlstd{(}\hlnum{0}\hlstd{, npop)}
\hlkwa{for}\hlstd{(k} \hlkwa{in} \hlnum{1}\hlopt{:}\hlstd{ncores)\{}
  \hlkwa{for}\hlstd{(j} \hlkwa{in} \hlnum{1}\hlopt{:}\hlstd{nsamp)\{}
    \hlstd{sd.Efron.mean} \hlkwb{<-} \hlstd{sd.Efron.mean} \hlopt{+}
      \hlstd{blah.second[[k]][[j]]} \hlopt{/} \hlstd{nboot}
  \hlstd{\}}
\hlstd{\}}
\end{alltt}
\end{kframe}
\end{knitrout}


Expected Darwinian fitness and its variability are now estimated for the 
hypothetical \emph{M. guttatus} flowers using the MLE. In this setting, 
the variability is estimated using the output from \texttt{fit.boot.Efron} 
to construct the bootstrap estimator of variability given by
$$ 
  \frac{1}{99} \sum_{b=1}^{100} 
    \left(g(\hat{\tau}^{(b)}) - g(\hat{\tau})\right)
    \left(g(\hat{\tau}^{(b)}) - g(\hat{\tau})\right)^T
$$ 
where $g(\hat{\tau})$ is estimated expected Darwinian fitness for the 
hypothetical \emph{M. guttatus} flowers using the MLE and 
$g(\hat{\tau}^{(b)})$ are the bootstrapped estimates.


\begin{knitrout}
\definecolor{shadecolor}{rgb}{0.969, 0.969, 0.969}\color{fgcolor}\begin{kframe}
\begin{alltt}
\hlstd{data.renew} \hlkwb{<-} \hlkwd{asterdata}\hlstd{(fred,} \hlkwc{vars} \hlstd{= vars,} \hlkwc{pred} \hlstd{= pred,}
  \hlkwc{group} \hlstd{=} \hlkwd{rep}\hlstd{(}\hlnum{0}\hlstd{,} \hlkwd{length}\hlstd{(vars)),} \hlkwc{code} \hlstd{= code,}
  \hlkwc{families} \hlstd{= families)}
\hlstd{origin.renew} \hlkwb{<-} \hlstd{m.main}\hlopt{$}\hlstd{origin[}\hlnum{1}\hlopt{:}\hlstd{npop,]}
\hlstd{origin.renew} \hlkwb{<-} \hlkwd{as.vector}\hlstd{(origin.renew)}
\hlstd{mu.renew} \hlkwb{<-} \hlkwd{transformUnconditional}\hlstd{(beta, data.renew,}
  \hlkwc{modmat} \hlstd{= modmat.renew,} \hlkwc{from} \hlstd{=} \hlstr{"beta"}\hlstd{,} \hlkwc{to} \hlstd{=} \hlstr{"mu"}\hlstd{,}
  \hlkwc{offset} \hlstd{= origin.renew)}
\hlstd{fit} \hlkwb{<-} \hlstd{amat} \hlopt{%*%} \hlstd{mu.renew}


\hlstd{dev} \hlkwb{<-} \hlstd{blah}\hlopt{$}\hlstd{MLE.boot.out} \hlopt{-} \hlkwd{as.vector}\hlstd{(fit)}
\hlstd{x} \hlkwb{<-} \hlstd{(dev[,} \hlnum{1}\hlstd{]} \hlopt{%*%} \hlkwd{t}\hlstd{(dev[,} \hlnum{1}\hlstd{]))}
\hlkwa{for}\hlstd{(j} \hlkwa{in} \hlnum{2}\hlopt{:}\hlstd{nboot)\{}
  \hlstd{x} \hlkwb{<-} \hlstd{x} \hlopt{+} \hlstd{(dev[, j]} \hlopt{%*%} \hlkwd{t}\hlstd{(dev[, j]))}
\hlstd{\}}
\end{alltt}
\end{kframe}
\end{knitrout}


The ratio of the the estimated standard error using envelope estimation 
obtained from Algorithm 1 to the estimated standard error using the MLE is 
computed.


\begin{knitrout}
\definecolor{shadecolor}{rgb}{0.969, 0.969, 0.969}\color{fgcolor}\begin{kframe}
\begin{alltt}
\hlstd{sd.MLE} \hlkwb{<-} \hlkwd{sqrt}\hlstd{(}\hlkwd{diag}\hlstd{(x)} \hlopt{/} \hlstd{(nboot} \hlopt{-} \hlnum{1}\hlstd{))}
\hlstd{ratio} \hlkwb{<-} \hlstd{sd.MLE} \hlopt{/} \hlstd{sd.Efron.mean}
\end{alltt}
\end{kframe}
\end{knitrout}



The output below is shown for hypothetical \emph{M. guttatus} flowers who 
have estimated expected Darwinian fitness greater than 7 where estimation is 
calculated using MLE. The output table table has 5 columns. The first two 
columns correspond to estimation of expected Darwinian fitness using envelope 
methodology and its bootstrapped standard error respectively. The next two 
columns correspond to estimation of expected Darwinian fitness using MLE and 
its bootstrapped standard error respectively. The fifth column displays the 
ratios of estimated standard errors where ratios greater than 1 indicate that 
the envelope estimator has lower estimated variability. We can see that 
efficiency gains are obtained for the displayed output. More importantly, 
this variance reduction implies more precise inference in this life history 
analysis. For example, the envelope estimator can statistically distinguish 
($\alpha = 0.05$, unadjusted for multiple comparisons) the second row of 
the output table from the fifth row of the output table. 
%Therefore, one could additionally 
%conclude that the trait values of the individual in the second row of 
%Table~\ref{Tab6:target.env} are not expected to maximize expected Darwinian 
%fitness. 
The combination of envelope methodology into the aster model framework 
allowed for us to consider a smaller set of traits associated with high 
Darwinian fitness. 


\begin{knitrout}
\definecolor{shadecolor}{rgb}{0.969, 0.969, 0.969}\color{fgcolor}\begin{kframe}
\begin{alltt}
\hlstd{tab3} \hlkwb{<-} \hlkwd{matrix}\hlstd{(}\hlnum{0}\hlstd{,} \hlkwc{nrow} \hlstd{= npop,} \hlkwc{ncol} \hlstd{=} \hlnum{5}\hlstd{)}
\hlkwd{colnames}\hlstd{(tab3)} \hlkwb{<-} \hlkwd{c}\hlstd{(}\hlstr{"env"}\hlstd{,}\hlstr{"se(env)"}\hlstd{,}\hlstr{"MLE"}\hlstd{,}\hlstr{"se(MLE)"}\hlstd{,}\hlstr{"ratio"}\hlstd{)}
\hlstd{tab3[,} \hlnum{1}\hlstd{]} \hlkwb{<-} \hlkwd{rowSums}\hlstd{(blah}\hlopt{$}\hlstd{env.boot.out)} \hlopt{/} \hlstd{nboot}
\hlstd{tab3[,} \hlnum{2}\hlstd{]} \hlkwb{<-} \hlstd{sd.Efron.mean}
\hlstd{tab3[,} \hlnum{3}\hlstd{]} \hlkwb{<-} \hlstd{fit}
\hlstd{tab3[,} \hlnum{4}\hlstd{]} \hlkwb{<-} \hlstd{sd.MLE}
\hlstd{tab3[,} \hlnum{5}\hlstd{]} \hlkwb{<-} \hlstd{ratio}
\hlstd{tab3[(tab3[,} \hlnum{3}\hlstd{]} \hlopt{>} \hlnum{7}\hlstd{), ]}
\end{alltt}
\begin{verbatim}
##            env   se(env)       MLE   se(MLE)    ratio
## [1,]  9.646009 0.3255334  9.170619 0.6421457 1.972595
## [2,]  8.639696 0.2999494  8.887339 0.3688551 1.229724
## [3,]  7.658725 0.3154779  7.602779 0.3608801 1.143916
## [4,]  7.517441 0.5386544  7.009952 0.6493285 1.205464
## [5,] 10.942830 0.6068946 10.474696 0.8958160 1.476065
## [6,]  7.497749 0.5208250  7.521516 0.6579239 1.263234
## [7,]  9.646009 0.3255334  9.170619 0.6421457 1.972595
## [8,]  9.646009 0.3255334  9.170619 0.6421457 1.972595
\end{verbatim}
\end{kframe}
\end{knitrout}






\begin{thebibliography}{}

%\bibitem[Cook, et al.(2010)Cook, Li, and Chiaromonte]{cook}
%Cook, R.~D., Li, B., Chiaromonte, F. (2010).
%\newblock Envelope models for parsimonious and 
%efficient multivariate linear regression
%\newblock \emph{Statistica Sinica}, \textbf{20}, 927-1010.

\bibitem[Cook and Zhang(2014)Cook and Zhang]{cook2}
Cook, R.~D., Zhang, X. (2014).
\newblock Foundations for Envelope Models and Methods.
\newblock \emph{JASA}, \textbf{110:510}: 599-611.

\bibitem[Cook and Zhang(2016)Cook and Zhang]{cook3}
Cook, R.~D., Zhang, X. (2016).
\newblock Algorithms for Envelope Estimation.
\newblock \emph{Journal of Computational and Graphical Statistics}, 
Published online. \url{DOI:10.1080/10618600.2015.1029577}.

%\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[Eck(2015)]{envlp-package}
Eck, D.~J. (2015).
\newblock R package \texttt{envlpaster}, version 0.1-1.
\newblock \url{http://cran.r-project.org/package=envlpaster}.

\bibitem[Eck, Geyer, and Cook(2018)Eck, Geyer, and Cook]{eck}
Eck, D.~J., Geyer, C.~J., and Cook, R.~D. (2018).
\newblock Combining Envelope Methodology and Aster Models for Variance 
  Reduction in Life History Analyses.
\newblock \emph{In prep}.

\bibitem[Efron(2014)Efron]{efron}
Efron, B. (2014).
\newblock Estimation and Accuracy After Model Selection.
\newblock \emph{JASA}, \textbf{109:507}, 991-1007.

%\bibitem[Etterson, J.~R., et al.]{etterson}
%Etterson, J.~R., et al. (2001).
%\newblock Constraint to Adaptive Evolution in Response to
%Global Warming
%\newblock \emph{Science}, \textbf{294}, 151-154.

%\bibitem[Fisher (1930)]{fisher}
%Fisher, R.~A. (1930).
%\newblock The Genetical Theory of Natural Selection.
%\newblock Oxford: Clarendon Press.

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

%\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, et al.(2007)Geyer, Wagenius, and Shaw]{geyer}
%Geyer, C.~J., Wagenius, Stuart, Shaw, R. (2007).
%\newblock Aster models for life history analysis.
%\newblock \emph{Biometrika}, \textbf{94}, 415-426.

%\bibitem[Geyer, et al.(2013)]{geyer2}
%Geyer, C.~J., Ridley, C.~E., Latta, R.~G., Etterson, J.~R., Shaw, R.~G. (2013). 
%\newblock Local adaptation and genetic effects on fitness: calculations 
%for exponential family models with random effects
%\newblock \emph{Annals of Applied Statistics}, \textbf{7}, 1778-1795.

\bibitem[Geyer, C.~J.(2010)Geyer]{geyer3}
Geyer, C.~J. (2010).
\newblock A Philosophical Look at Aster Models. Technical Report
No. 676. School of Statistics, University of Minnesota.
\newblock http://purl.umn.edu/57163.

\bibitem[Geyer, et al.(2009)Geyer and Shaw]{geyer4}
Geyer, C.~J. and Shaw, R. (2009).
\newblock Model Selection in Estimation of Fitness Landscapes. Technical Report
No. 671. School of Statistics, University of Minnesota.
\newblock http://conservancy.umn.edu/handle/11299/56219.

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

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

\bibitem[Lowry and Willis(2010)Lowry and Willis]{lowry}
Lowry, D.~B. and Willis, J.~H. (2010)
\newblock A Widespread Chromosomal Inversion Polymorphism Contributes 
to a Major Life-History Transition, Local Adaptation, and Reproductive 
Isolation.
\newblock \emph{PLoS Biol}, 8(9): e1000500. \url{doi:10.1371/journal.pbio.1000500}.

%\bibitem[Shaw, et al.(2010)Shaw and Geyer]{shaw}
%Shaw, R., Geyer, C.~G. (2010).
%\newblock Inferring fitness landscapes.
%\newblock \emph{Evolution}, \textbf{64}, 2510-2520.

%\bibitem[Shaw, et al.(2008)Shaw, Geyer, Wagenius, Hangelbroek, and Etterson]{shaw2}
%Shaw, R., Geyer, C.~G., Wagenius, S., Hangelbroek, H., and Etterson, J.~R. (2010).
%\newblock Unifying life-history analyses for 
%inference of fitness and population growth
%\newblock \emph{The American Naturalist}, \textbf{172} E35-E47.

%\bibitem[Su, et al.(2011)Su and Cook]{su}
%Su, Z., Cook, R.~D. (2011).
%\newblock Partial envelopes for efficient estimation in
%multivariate linear regression
%\newblock \emph{Biometrika}, \textbf{98}, 133-146.

\bibitem[Zhang and Mai(2017)Zhang and Mai]{zhangmai}
Zhang, X. and Mai, Z. (2017).
\newblock Model-free Envelope Dimension Selection.
\newblock \url{https://arxiv.org/pdf/1709.03945.pdf}. \emph{Submitted}. 

\end{thebibliography}





\section*{A. make the Mguttatus data}

In this Appendix we show how the Mguttatus dataset in the \texttt{envlpaster} 
package was initially created starting from the data for the original analysis
(Lowry, et al., 2010), which was a comma separated values (CSV) file.
This file can be found at the same location at the University of Minnesota
Digital Conservancy (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{test} \hlkwb{<-} \hlkwd{read.csv}\hlstd{(}\hlstr{"Aster_across_sites_DIV1.csv"}\hlstd{,}
  \hlkwc{sep} \hlstd{=} \hlstr{","}\hlstd{,} \hlkwc{header} \hlstd{= T)}
\end{alltt}
\end{kframe}
\end{knitrout}


Now we need to specify the rest of the graph shown in Figure~\ref{astergraph}. 
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. In this analysis, there are no groups. 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}. \texttt{code} is an index vector into the families list; 
\texttt{families[[code[j]]]} is the family for the conditional distribution 
of the dependence group containing node \texttt{j} given its predecessor node.


\begin{knitrout}
\definecolor{shadecolor}{rgb}{0.969, 0.969, 0.969}\color{fgcolor}\begin{kframe}
\begin{alltt}
\hlstd{vars} \hlkwb{<-} \hlkwd{c}\hlstd{(}\hlstr{"sur_flw"}\hlstd{,}\hlstr{"flws"}\hlstd{)}
\hlstd{pred} \hlkwb{<-} \hlkwd{c}\hlstd{(}\hlnum{0}\hlstd{,}\hlnum{1}\hlstd{)}
\hlstd{group} \hlkwb{<-} \hlkwd{rep}\hlstd{(}\hlnum{0}\hlstd{,}\hlkwd{length}\hlstd{(pred))}
\hlstd{code} \hlkwb{<-} \hlkwd{c}\hlstd{(}\hlnum{1}\hlstd{,}\hlnum{2}\hlstd{)}
\hlstd{families} \hlkwb{<-} \hlkwd{list}\hlstd{(}\hlstr{"bernoulli"}\hlstd{,} \hlkwd{fam.zero.truncated.poisson}\hlstd{())}
\end{alltt}
\end{kframe}
\end{knitrout}


Now we make the asterdata object and verify that this aster object is the 
as the asterdata provided by the \texttt{envlpaster} package.


\begin{knitrout}
\definecolor{shadecolor}{rgb}{0.969, 0.969, 0.969}\color{fgcolor}\begin{kframe}
\begin{alltt}
\hlstd{data} \hlkwb{<-} \hlkwd{asterdata}\hlstd{(}\hlkwc{data} \hlstd{= test,} \hlkwc{vars} \hlstd{= vars,} \hlkwc{pred} \hlstd{= pred,}
  \hlkwc{group} \hlstd{= group,} \hlkwc{code} \hlstd{= code,} \hlkwc{families} \hlstd{= families)}
\hlkwd{all.equal}\hlstd{(data}\hlopt{$}\hlstd{redata, Mguttatus}\hlopt{$}\hlstd{redata)}
\end{alltt}
\begin{verbatim}
## [1] TRUE
\end{verbatim}
\begin{alltt}
\hlkwd{all.equal}\hlstd{(data}\hlopt{$}\hlstd{repred, Mguttatus}\hlopt{$}\hlstd{repred)}
\end{alltt}
\begin{verbatim}
## [1] TRUE
\end{verbatim}
\begin{alltt}
\hlkwd{all.equal}\hlstd{(data}\hlopt{$}\hlstd{recode, Mguttatus}\hlopt{$}\hlstd{recode)}
\end{alltt}
\begin{verbatim}
## [1] TRUE
\end{verbatim}
\begin{alltt}
\hlkwd{all.equal}\hlstd{(data}\hlopt{$}\hlstd{regroup, Mguttatus}\hlopt{$}\hlstd{regroup)}
\end{alltt}
\begin{verbatim}
## [1] TRUE
\end{verbatim}
\begin{alltt}
\hlkwd{all.equal}\hlstd{(data}\hlopt{$}\hlstd{redelta, Mguttatus}\hlopt{$}\hlstd{redelta)}
\end{alltt}
\begin{verbatim}
## [1] TRUE
\end{verbatim}
\begin{alltt}
\hlkwd{all.equal}\hlstd{(data}\hlopt{$}\hlstd{families, Mguttatus}\hlopt{$}\hlstd{families)}
\end{alltt}
\begin{verbatim}
## [1] TRUE
\end{verbatim}
\end{kframe}
\end{knitrout}






\section*{B. make the simulated data}
The last line of code shows that the generated dataset analyzed in Example 1 is the same as the dataset included in the \texttt{envlpaster} package.


\begin{knitrout}
\definecolor{shadecolor}{rgb}{0.969, 0.969, 0.969}\color{fgcolor}\begin{kframe}
\begin{alltt}
\hlkwd{set.seed}\hlstd{(}\hlnum{13}\hlstd{)}
\hlstd{nind} \hlkwb{<-} \hlnum{3000}
\hlstd{ntime} \hlkwb{<-} \hlnum{10}
\hlstd{psurv} \hlkwb{<-} \hlnum{0.95}
\hlstd{prepr} \hlkwb{<-} \hlnum{0.7}
\hlstd{mpois} \hlkwb{<-} \hlnum{1}
\hlstd{theta.surv} \hlkwb{<-} \hlkwd{log}\hlstd{(psurv)} \hlopt{-} \hlkwd{log}\hlstd{(}\hlnum{1} \hlopt{-} \hlstd{psurv)}
\hlstd{theta.repr} \hlkwb{<-} \hlkwd{log}\hlstd{(prepr)} \hlopt{-} \hlkwd{log}\hlstd{(}\hlnum{1} \hlopt{-} \hlstd{prepr)}
\hlstd{theta.pois} \hlkwb{<-} \hlkwd{log}\hlstd{(mpois)}
\hlstd{tau.pois} \hlkwb{<-} \hlstd{mpois}\hlopt{/}\hlstd{(}\hlnum{1} \hlopt{-} \hlkwd{exp}\hlstd{(}\hlopt{-}\hlstd{mpois))}


\hlstd{vars} \hlkwb{<-} \hlkwd{as.vector}\hlstd{(}\hlkwd{outer}\hlstd{(}\hlkwd{c}\hlstd{(}\hlstr{"u"}\hlstd{,} \hlstr{"v"}\hlstd{,} \hlstr{"w"}\hlstd{),} \hlnum{1}\hlopt{:}\hlnum{10}\hlstd{, paste,}
  \hlkwc{sep} \hlstd{=} \hlstr{""}\hlstd{))}
\hlstd{vtype} \hlkwb{<-} \hlkwd{as.factor}\hlstd{(}\hlkwd{substr}\hlstd{(}\hlkwd{as.character}\hlstd{(vars),} \hlnum{1}\hlstd{,}\hlnum{1}\hlstd{))}
\hlstd{fam} \hlkwb{<-} \hlkwd{rep}\hlstd{(}\hlnum{1}\hlstd{,} \hlkwd{length}\hlstd{(vars))}
\hlstd{fam[vtype} \hlopt{==} \hlstr{"w"}\hlstd{]} \hlkwb{<-} \hlnum{3}
\hlstd{pred} \hlkwb{<-} \hlkwd{seq}\hlstd{(}\hlkwc{along} \hlstd{= vars)} \hlopt{-} \hlnum{1}
\hlstd{pred[vtype} \hlopt{==} \hlstr{"u"}\hlstd{]} \hlkwb{<-} \hlkwd{seq}\hlstd{(}\hlkwc{along} \hlstd{= vars)[vtype} \hlopt{==} \hlstr{"u"}\hlstd{]} \hlopt{-} \hlnum{3}
\hlstd{pred[}\hlnum{1}\hlstd{]} \hlkwb{<-} \hlnum{0}

\hlstd{root} \hlkwb{<-} \hlkwd{matrix}\hlstd{(}\hlnum{1}\hlstd{,} \hlkwc{nrow} \hlstd{= nind,} \hlkwc{ncol} \hlstd{=} \hlkwd{length}\hlstd{(vars))}
\hlstd{theta} \hlkwb{<-} \hlstd{root}
\hlstd{theta[, vtype} \hlopt{==} \hlstr{"u"}\hlstd{]} \hlkwb{<-} \hlstd{theta.surv}
\hlstd{theta[, vtype} \hlopt{==} \hlstr{"v"}\hlstd{]} \hlkwb{<-} \hlstd{theta.repr}
\hlstd{theta[, vtype} \hlopt{==} \hlstr{"w"}\hlstd{]} \hlkwb{<-} \hlstd{theta.pois}
\hlstd{x} \hlkwb{<-} \hlkwd{raster}\hlstd{(theta, pred, fam, root)}
\hlkwd{dimnames}\hlstd{(x)} \hlkwb{<-} \hlkwd{list}\hlstd{(}\hlkwa{NULL}\hlstd{, vars)}
\hlstd{dat} \hlkwb{<-} \hlkwd{as.data.frame}\hlstd{(x)}
\hlstd{dat} \hlkwb{<-} \hlkwd{as.list}\hlstd{(dat)}
\hlstd{dat[[}\hlstr{"root"}\hlstd{]]} \hlkwb{<-} \hlkwd{rep}\hlstd{(}\hlnum{1}\hlstd{, nind)}


\hlstd{npheno} \hlkwb{<-} \hlnum{10}
\hlstd{zbase} \hlkwb{<-} \hlkwd{rnorm}\hlstd{(nind)}
\hlkwa{for} \hlstd{(i} \hlkwa{in} \hlnum{1}\hlopt{:}\hlstd{npheno) \{}
  \hlstd{labz} \hlkwb{<-} \hlkwd{paste}\hlstd{(}\hlstr{"z"}\hlstd{, i,} \hlkwc{sep} \hlstd{=} \hlstr{""}\hlstd{)}
  \hlstd{dat[[labz]]} \hlkwb{<-} \hlstd{zbase} \hlopt{+} \hlkwd{rnorm}\hlstd{(nind)}
\hlstd{\}}


\hlstd{dat} \hlkwb{<-} \hlkwd{as.data.frame}\hlstd{(dat)}
\hlstd{redata} \hlkwb{<-} \hlkwd{reshape}\hlstd{(dat,} \hlkwc{varying} \hlstd{=} \hlkwd{list}\hlstd{(vars),} \hlkwc{direction} \hlstd{=} \hlstr{"long"}\hlstd{,}
  \hlkwc{timevar} \hlstd{=} \hlstr{"varb"}\hlstd{,} \hlkwc{times} \hlstd{=} \hlkwd{as.factor}\hlstd{(vars),} \hlkwc{v.names} \hlstd{=} \hlstr{"resp"}\hlstd{)}
\hlstd{wind} \hlkwb{<-} \hlkwd{grep}\hlstd{(}\hlstr{"w"}\hlstd{,} \hlkwd{as.character}\hlstd{(redata}\hlopt{$}\hlstd{varb))}
\hlkwa{for} \hlstd{(labz} \hlkwa{in} \hlkwd{grep}\hlstd{(}\hlstr{"z"}\hlstd{,} \hlkwd{names}\hlstd{(redata),} \hlkwc{value} \hlstd{=} \hlnum{TRUE}\hlstd{)) \{}
  \hlstd{redata[[labz]][}\hlopt{-}\hlstd{wind]} \hlkwb{<-} \hlnum{0}
\hlstd{\}}


\hlstd{redata}\hlopt{$}\hlstd{vtype} \hlkwb{<-} \hlkwd{as.factor}\hlstd{(}\hlkwd{substr}\hlstd{(}\hlkwd{as.character}\hlstd{(redata}\hlopt{$}\hlstd{varb),}
  \hlnum{1}\hlstd{,} \hlnum{1}\hlstd{))}
\hlstd{out2} \hlkwb{<-} \hlkwd{aster}\hlstd{(resp} \hlopt{~} \hlstd{vtype} \hlopt{+} \hlnum{0}\hlstd{, pred, fam, varb,}
  \hlstd{id, root,} \hlkwc{data} \hlstd{= redata,} \hlkwc{type} \hlstd{=} \hlstr{"conditional"}\hlstd{)}
\hlstd{out3} \hlkwb{<-} \hlkwd{aster}\hlstd{(resp} \hlopt{~} \hlstd{vtype} \hlopt{+} \hlnum{0}\hlstd{, pred, fam, varb,}
  \hlstd{id, root,} \hlkwc{data} \hlstd{= redata)}


\hlstd{redata}\hlopt{$}\hlstd{year} \hlkwb{<-} \hlkwd{as.numeric}\hlstd{(}\hlkwd{substring}\hlstd{(}\hlkwd{as.character}\hlstd{(redata}\hlopt{$}\hlstd{varb),} \hlnum{2}\hlstd{))}
\hlstd{redata}\hlopt{$}\hlstd{uyear} \hlkwb{<-} \hlstd{redata}\hlopt{$}\hlstd{year} \hlopt{*} \hlkwd{as.numeric}\hlstd{(}\hlkwd{as.character}\hlstd{(redata}\hlopt{$}\hlstd{vtype)} \hlopt{==}
  \hlstr{"u"}\hlstd{)}
\hlstd{out4} \hlkwb{<-} \hlkwd{aster}\hlstd{(resp} \hlopt{~} \hlstd{vtype} \hlopt{+} \hlstd{uyear, pred, fam, varb,}
  \hlstd{id, root,} \hlkwc{data} \hlstd{= redata)}


\hlcom{###############################}
\hlcom{##### Section 4.2 model 1 #####}
\hlstd{out5} \hlkwb{<-} \hlkwd{aster}\hlstd{(resp} \hlopt{~} \hlstd{vtype} \hlopt{+} \hlstd{uyear} \hlopt{+} \hlstd{z1} \hlopt{+} \hlstd{z2} \hlopt{+} \hlkwd{I}\hlstd{(z1}\hlopt{^}\hlnum{2}\hlstd{)} \hlopt{+}
  \hlkwd{I}\hlstd{(z2}\hlopt{^}\hlnum{2}\hlstd{)} \hlopt{+} \hlkwd{I}\hlstd{(z1} \hlopt{*} \hlstd{z2), pred, fam, varb, id, root,}
  \hlkwc{data} \hlstd{= redata)}

\hlcom{# change the covariate data so that it is mean 0}
\hlstd{z1} \hlkwb{<-} \hlstd{dat}\hlopt{$}\hlstd{z1}
\hlstd{z2} \hlkwb{<-} \hlstd{dat}\hlopt{$}\hlstd{z2}
\hlstd{ascal} \hlkwb{<-} \hlnum{0.013}
\hlstd{quad} \hlkwb{<-} \hlstd{ascal} \hlopt{*} \hlstd{((z1} \hlopt{+} \hlstd{z2)} \hlopt{-} \hlstd{(z1}\hlopt{^}\hlnum{2} \hlopt{+} \hlstd{z2}\hlopt{^}\hlnum{2}\hlstd{)} \hlopt{+} \hlstd{z1} \hlopt{*} \hlstd{z2)}
\hlstd{con} \hlkwb{<-} \hlkwd{mean}\hlstd{(quad)}
\hlstd{beta.new} \hlkwb{<-} \hlstd{out5}\hlopt{$}\hlstd{coefficients}
\hlstd{beta.new[}\hlnum{1}\hlopt{:}\hlnum{4}\hlstd{]} \hlkwb{<-} \hlstd{out4}\hlopt{$}\hlstd{coefficients}
\hlstd{beta.new[}\hlnum{3}\hlstd{]} \hlkwb{<-} \hlstd{beta.new[}\hlnum{3}\hlstd{]} \hlopt{-} \hlstd{con}
\hlstd{beta.new[}\hlnum{5}\hlopt{:}\hlnum{6}\hlstd{]} \hlkwb{<-} \hlstd{ascal}
\hlstd{beta.new[}\hlnum{7}\hlopt{:}\hlnum{8}\hlstd{]} \hlkwb{<-} \hlstd{(}\hlopt{-}\hlstd{ascal)}
\hlstd{beta.new[}\hlnum{9}\hlstd{]} \hlkwb{<-} \hlstd{ascal}
\hlstd{beta.new} \hlkwb{<-} \hlkwd{round}\hlstd{(beta.new,} \hlnum{3}\hlstd{)}


\hlstd{theta} \hlkwb{<-} \hlkwd{predict}\hlstd{(out5,} \hlkwc{model.type} \hlstd{=} \hlstr{"conditional"}\hlstd{,}
  \hlkwc{parm.type} \hlstd{=} \hlstr{"canonical"}\hlstd{,} \hlkwc{newcoef} \hlstd{= beta.new)}
\hlstd{theta} \hlkwb{<-} \hlkwd{matrix}\hlstd{(theta,} \hlkwc{nrow} \hlstd{=} \hlkwd{nrow}\hlstd{(out5}\hlopt{$}\hlstd{x),} \hlkwc{ncol} \hlstd{=} \hlkwd{ncol}\hlstd{(out5}\hlopt{$}\hlstd{x))}
\hlstd{xnew} \hlkwb{<-} \hlkwd{raster}\hlstd{(theta, pred, fam, root)}
\hlkwd{dimnames}\hlstd{(xnew)} \hlkwb{<-} \hlkwd{list}\hlstd{(}\hlkwa{NULL}\hlstd{, vars)}
\hlstd{dnew1} \hlkwb{<-} \hlkwd{as.data.frame}\hlstd{(xnew)}
\hlstd{renew} \hlkwb{<-} \hlkwd{reshape}\hlstd{(dnew1,} \hlkwc{varying} \hlstd{=} \hlkwd{list}\hlstd{(vars),}
  \hlkwc{direction} \hlstd{=} \hlstr{"long"}\hlstd{,} \hlkwc{timevar} \hlstd{=} \hlstr{"varb"}\hlstd{,}
  \hlkwc{times} \hlstd{=} \hlkwd{as.factor}\hlstd{(vars),} \hlkwc{v.names} \hlstd{=} \hlstr{"resp"}\hlstd{)}
\hlstd{redata}\hlopt{$}\hlstd{resp1} \hlkwb{<-} \hlstd{renew}\hlopt{$}\hlstd{resp}
\hlstd{redata}\hlopt{$}\hlstd{tau1} \hlkwb{<-} \hlkwd{predict}\hlstd{(out5,} \hlkwc{newcoef} \hlstd{= beta.new)}
\hlstd{redata}\hlopt{$}\hlstd{phi1} \hlkwb{<-} \hlkwd{predict}\hlstd{(out5,} \hlkwc{parm.type} \hlstd{=} \hlstr{"canonical"}\hlstd{,}
  \hlkwc{newcoef} \hlstd{= beta.new)}
\hlstd{beta1true} \hlkwb{<-} \hlstd{beta.new}


\hlstd{out6} \hlkwb{<-} \hlkwd{aster}\hlstd{(resp1} \hlopt{~} \hlstd{vtype} \hlopt{+} \hlstd{uyear} \hlopt{+} \hlstd{z1} \hlopt{+} \hlstd{z2} \hlopt{+} \hlkwd{I}\hlstd{(z1}\hlopt{^}\hlnum{2}\hlstd{)} \hlopt{+}
  \hlkwd{I}\hlstd{(z2}\hlopt{^}\hlnum{2}\hlstd{)} \hlopt{+} \hlkwd{I}\hlstd{(z1} \hlopt{*} \hlstd{z2), pred, fam, varb, id, root,}
  \hlkwc{data} \hlstd{= redata)}


\hlcom{######################}
\hlcom{# tau envelope model stuff}
\hlstd{target} \hlkwb{<-} \hlnum{5}\hlopt{:}\hlnum{9}
\hlstd{m} \hlkwb{<-} \hlkwd{length}\hlstd{(target)}
\hlstd{p} \hlkwb{<-} \hlkwd{length}\hlstd{(out6}\hlopt{$}\hlstd{coef)}
\hlstd{avar} \hlkwb{<-} \hlstd{out6}\hlopt{$}\hlstd{fisher}
\hlstd{beta} \hlkwb{<-} \hlstd{out6}\hlopt{$}\hlstd{coef}
\hlstd{avar.targ} \hlkwb{<-} \hlstd{avar[target,target]}
\hlstd{beta.targ} \hlkwb{<-} \hlstd{beta[target]}

\hlstd{mu} \hlkwb{<-} \hlkwd{predict}\hlstd{(out6,} \hlkwc{parm.type} \hlstd{=} \hlstr{"mean.value"}\hlstd{,}
  \hlkwc{model.type} \hlstd{=} \hlstr{"unconditional"}\hlstd{)}
\hlstd{modmat} \hlkwb{<-} \hlstd{out6}\hlopt{$}\hlstd{modmat}
\hlstd{modmat.mat} \hlkwb{<-} \hlkwd{matrix}\hlstd{(modmat,} \hlkwc{ncol} \hlstd{=} \hlkwd{length}\hlstd{(beta))}
\hlstd{offset} \hlkwb{<-} \hlkwd{as.vector}\hlstd{(out6}\hlopt{$}\hlstd{origin)}
\hlstd{x} \hlkwb{<-} \hlstd{out6}\hlopt{$}\hlstd{x}
\hlstd{tau} \hlkwb{<-} \hlkwd{crossprod}\hlstd{(modmat.mat, mu)}
\hlstd{tau.targ} \hlkwb{<-} \hlstd{tau[target]}


\hlcom{# check whether transformUnconditional and predict.aster}
\hlcom{# are returning the same quantity}
\hlstd{fam[}\hlkwd{which}\hlstd{(fam} \hlopt{==} \hlnum{3}\hlstd{)]} \hlkwb{<-} \hlnum{2}
\hlstd{code} \hlkwb{<-} \hlstd{fam}
\hlstd{families} \hlkwb{=} \hlkwd{list}\hlstd{(}\hlkwd{fam.bernoulli}\hlstd{(),}\hlkwd{fam.zero.truncated.poisson}\hlstd{())}
\hlstd{data} \hlkwb{<-} \hlkwd{asterdata}\hlstd{(}\hlkwc{data} \hlstd{= dat,} \hlkwc{vars} \hlstd{= vars,} \hlkwc{pred} \hlstd{= pred,}
  \hlkwc{group} \hlstd{=} \hlkwd{rep}\hlstd{(}\hlnum{0}\hlstd{,} \hlkwd{length}\hlstd{(vars)),} \hlkwc{code} \hlstd{= code,}
  \hlkwc{families} \hlstd{= families)}
\hlstd{tau.foo} \hlkwb{<-} \hlkwd{transformUnconditional}\hlstd{(}\hlkwc{parm} \hlstd{= beta, modmat.mat, data,}
  \hlkwc{from} \hlstd{=} \hlstr{"beta"}\hlstd{,} \hlkwc{to} \hlstd{=} \hlstr{"tau"}\hlstd{,} \hlkwc{offset} \hlstd{= offset,} \hlkwc{tol} \hlstd{=} \hlnum{1e-10}\hlstd{)}


\hlcom{# put the envelope truth into the model}
\hlstd{G} \hlkwb{<-} \hlkwd{eigen}\hlstd{(avar.targ)}\hlopt{$}\hlstd{vec[,}\hlkwd{c}\hlstd{(}\hlnum{1}\hlstd{,}\hlnum{4}\hlstd{)]}
\hlstd{P} \hlkwb{<-} \hlkwd{tcrossprod}\hlstd{(G)}
\hlstd{tau_true} \hlkwb{<-} \hlstd{tau.foo}
\hlstd{M} \hlkwb{<-} \hlkwd{t}\hlstd{(modmat.mat)}
\hlstd{M2} \hlkwb{<-} \hlstd{M[target,]; M2} \hlkwb{<-} \hlstd{P} \hlopt{%*%} \hlstd{M2}
\hlstd{M[target,]} \hlkwb{<-} \hlstd{M2}
\hlstd{modmat.mat.env} \hlkwb{<-} \hlkwd{t}\hlstd{(M)}
\hlstd{tau_true[target]} \hlkwb{<-} \hlstd{P} \hlopt{%*%} \hlstd{tau.targ}
\hlstd{beta.foo} \hlkwb{<-} \hlkwd{transformUnconditional}\hlstd{(}\hlkwc{parm} \hlstd{= tau_true,}
  \hlstd{modmat.mat, data,} \hlkwc{from} \hlstd{=} \hlstr{"tau"}\hlstd{,} \hlkwc{to} \hlstd{=} \hlstr{"beta"}\hlstd{,}
  \hlkwc{offset} \hlstd{= offset)}


\hlstd{theta} \hlkwb{<-} \hlkwd{predict}\hlstd{(out6,} \hlkwc{model.type} \hlstd{=} \hlstr{"conditional"}\hlstd{,}
  \hlkwc{parm.type} \hlstd{=} \hlstr{"canonical"}\hlstd{,} \hlkwc{newcoef} \hlstd{= beta.foo)}
\hlstd{theta} \hlkwb{<-} \hlkwd{matrix}\hlstd{(theta,} \hlkwc{nrow} \hlstd{=} \hlkwd{nrow}\hlstd{(out6}\hlopt{$}\hlstd{x),} \hlkwc{ncol} \hlstd{=} \hlkwd{ncol}\hlstd{(out6}\hlopt{$}\hlstd{x))}
\hlstd{xnew} \hlkwb{<-} \hlkwd{raster}\hlstd{(theta, out6}\hlopt{$}\hlstd{pred, out6}\hlopt{$}\hlstd{fam, out6}\hlopt{$}\hlstd{root)}
\hlkwd{dimnames}\hlstd{(xnew)} \hlkwb{<-} \hlkwd{list}\hlstd{(}\hlkwa{NULL}\hlstd{, vars)}
\hlstd{xnew} \hlkwb{<-} \hlkwd{cbind}\hlstd{(xnew, dat[,}\hlkwd{c}\hlstd{(}\hlnum{32}\hlopt{:}\hlnum{41}\hlstd{)])}
\hlstd{xnew}\hlopt{$}\hlstd{root} \hlkwb{<-} \hlnum{1}
\hlstd{data} \hlkwb{<-} \hlkwd{asterdata}\hlstd{(}\hlkwc{data} \hlstd{= xnew,} \hlkwc{vars} \hlstd{= vars,} \hlkwc{pred} \hlstd{= pred,}
  \hlkwc{group} \hlstd{=} \hlkwd{rep}\hlstd{(}\hlnum{0}\hlstd{,} \hlkwd{length}\hlstd{(vars)),} \hlkwc{code} \hlstd{= code,}
  \hlkwc{families} \hlstd{= families)}
\end{alltt}
\end{kframe}
\end{knitrout}


\begin{knitrout}
\definecolor{shadecolor}{rgb}{0.969, 0.969, 0.969}\color{fgcolor}\begin{kframe}
\begin{alltt}
\hlkwd{all.equal}\hlstd{(data}\hlopt{$}\hlstd{redata, simdata30nodes.asterdata}\hlopt{$}\hlstd{redata)}
\end{alltt}
\begin{verbatim}
## [1] TRUE
\end{verbatim}
\begin{alltt}
\hlkwd{all.equal}\hlstd{(data}\hlopt{$}\hlstd{repred, simdata30nodes.asterdata}\hlopt{$}\hlstd{repred)}
\end{alltt}
\begin{verbatim}
## [1] TRUE
\end{verbatim}
\begin{alltt}
\hlkwd{all.equal}\hlstd{(data}\hlopt{$}\hlstd{recode, simdata30nodes.asterdata}\hlopt{$}\hlstd{recode)}
\end{alltt}
\begin{verbatim}
## [1] TRUE
\end{verbatim}
\begin{alltt}
\hlkwd{all.equal}\hlstd{(data}\hlopt{$}\hlstd{regroup, simdata30nodes.asterdata}\hlopt{$}\hlstd{regroup)}
\end{alltt}
\begin{verbatim}
## [1] TRUE
\end{verbatim}
\begin{alltt}
\hlkwd{all.equal}\hlstd{(data}\hlopt{$}\hlstd{redelta, simdata30nodes.asterdata}\hlopt{$}\hlstd{redelta)}
\end{alltt}
\begin{verbatim}
## [1] TRUE
\end{verbatim}
\begin{alltt}
\hlkwd{all.equal}\hlstd{(data}\hlopt{$}\hlstd{families, simdata30nodes.asterdata}\hlopt{$}\hlstd{families)}
\end{alltt}
\begin{verbatim}
## [1] TRUE
\end{verbatim}
\end{kframe}
\end{knitrout}


\end{document}
