Combining Envelope Methodology and Aster Models for Variance Reduction in Life History Analyses

TL;DR

This paper introduces a new envelope methodology tailored for aster models to improve the precision of Darwinian fitness estimates in life history studies, demonstrated through simulations and real data.

Contribution

We develop an envelope framework and estimator specifically for aster models, enhancing variance reduction and statistical inference in life history analyses.

Findings

Significant variance reduction in simulated data

Improved precision in extit{Mimulus guttatus} case study

Stronger conclusions about fitness determinants

Abstract

Precise estimation of expected Darwinian fitness, the expected lifetime number of offspring of organism, is a central component of life history analysis. The aster model serves as a defensible statistical model for distributions of Darwinian fitness. The aster model is equipped to incorporate the major life stages an organism travels through which separately may effect Darwinian fitness. Envelope methodology reduces asymptotic variability by establishing a link between unknown parameters of interest and the asymptotic covariance matrices of their estimators. It is known both theoretically and in applications that incorporation of envelope methodology reduces asymptotic variability. We develop an envelope framework, including a new envelope estimator, that is appropriate for aster analyses. The level of precision provided from our methods allows researchers to draw stronger conclusions about the driving forces of Darwinian fitness from their life history analyses than they could with the aster model alone. Our methods are illustrated on a simulated dataset and a life history analysis of \emph{Mimulus guttatus} flowers is provided. Useful variance reduction is obtained in both analyses.

Tables

Table 1. Table 1: Comparison of the maximum likelihood estimator and the envelope estimator for seven individuals with high estimates of expected Darwinian fitness in Example 1. The fifth column is the ratio of bootstrapped standard errors for g ^ MLE subscript ^ 𝑔 MLE \hat{g}_{\text{MLE}} to g ^ env subscript ^ 𝑔 env \hat{g}_{\text{env}} .

${\hat{g}}_{\text{env}}$	se$\left({\hat{g}}_{\text{env}}\right)$	${\hat{g}}_{\text{MLE}}$	se$\left({\hat{g}}_{\text{MLE}}\right)$	ratio
8.556	0.174	8.701	0.260	1.491
9.014	0.111	8.939	0.135	1.222
7.817	0.414	8.054	0.442	1.069
9.174	0.163	9.193	0.170	1.045
9.018	0.113	9.120	0.128	1.133
8.612	0.162	8.518	0.278	1.709
7.761	0.215	8.096	0.331	1.534

Table 2. Table 2: Comparison of the maximum likelihood estimator and the envelope estimator for seven individuals with high estimates of expected Darwinian fitness in Example 2. The fifth column is the ratio of bootstrapped standard errors for g ^ MLE subscript ^ 𝑔 MLE \hat{g}_{\text{MLE}} to g ^ env subscript ^ 𝑔 env \hat{g}_{\text{env}} .

${\hat{g}}_{\text{env}}$	se$\left({\hat{g}}_{\text{env}}\right)$	${\hat{g}}_{\text{MLE}}$	se$\left({\hat{g}}_{\text{MLE}}\right)$	ratio
9.646	0.326	9.171	0.642	1.973
8.640	0.300	8.887	0.369	1.230
7.659	0.315	7.603	0.361	1.144
7.517	0.539	7.010	0.649	1.205
10.943	0.607	10.475	0.896	1.476
7.329	0.707	6.618	1.038	1.469
7.498	0.521	7.522	0.658	1.263

Equations

$$\surd n\left(\hat{\tau }-\tau \right)\stackrel{d}{⟶}N\left(0,\Sigma \right),$$

$$\surd n\left\{g(\hat{\tau })-g(\tau )\right\}\stackrel{d}{⟶}N\left\{0,\nabla g(\tau )\Sigma \nabla g(\tau {)}^T\right\}.$$

$$\mathcal{T}\subset {\mathcal{E}}_{{\Sigma }_{\upsilon ,\upsilon }}(\mathcal{T}),{\Sigma }_{\upsilon ,\upsilon }=P_{\mathcal{E}}{\Sigma }_{\upsilon ,\upsilon }{P}_{\mathcal{E}}+Q_{\mathcal{E}}{\Sigma }_{\upsilon ,\upsilon }{Q}_{\mathcal{E}};$$

$$\mathcal{T}\subset \text{span}(\Gamma ),{\Sigma }_{\upsilon ,\upsilon }=\Gamma \Omega {\Gamma }^T+{\Gamma }_o{\Omega }_o{\Gamma }_o^T;$$

$$\hat{\tau}_{\text{1D}}=\left(\begin{array}[]{c}\widehat{\gamma}\\ \hat{\upsilon}_{\text{1D}}\end{array}\right)=\left(\begin{array}[]{cc}I&0\\ 0&P_{\hat{\mathcal{E}}}\end{array}\right)M^{T}Y=M_{\text{1D}}^{T}Y,\qquad M_{\text{1D}}=M\left(\begin{array}[]{cc}I&0\\ 0&P_{\hat{\mathcal{E}}}\end{array}\right).$$

$${\hat{g}}_{\text{1D}}=\frac{1}{B}\sum _{b=1}^{B}g\left\{{\hat{\tau }}_{\text{1D}}^{(b)}\right\}.$$

$$\hat{\tau}_{\text{env}}=\left(\begin{array}[]{c}\widehat{\gamma}\\ \hat{\upsilon}_{\text{env}}\end{array}\right)=\left(\begin{array}[]{cc}I&0\\ 0&\widehat{P}_{\widehat{\mathcal{G}}}\end{array}\right)M^{T}Y=M_{\text{env}}^{T}Y,\qquad M_{\text{env}}=M\left(\begin{array}[]{cc}I&0\\ 0&\widehat{P}_{\widehat{\mathcal{G}}}\end{array}\right).$$

$$Q_n(g)$$

$$=-n/2\log (g^T\hat{O}{\hat{M}}_2{\hat{O}}^{T}g)-n/2\log \left\{g^T\hat{O}{\left({\hat{M}}_2+{\hat{U}}_2\right)}^{-1}{\hat{O}}^{T}g\right\}+n\log (g^T\hat{O}{\hat{O}}^{T}g).$$

$$Q_n({\hat{v}}_k)$$

$$+n\log ({\hat{v}}_k^T{\hat{O}}^T\hat{O}{\hat{v}}_k)$$

$$=-n/2\log (A_k{\hat{O}}_u{\hat{M}}_2{\hat{O}}_u^T{A}_k)-n/2\log \left\{A_k{\hat{O}}_u{\left({\hat{M}}_2+{\hat{U}}_2\right)}^{-1}{\hat{O}}_u^T{A}_k\right\}$$

$$+n\log (A_k{\hat{O}}_u^T{\hat{O}}_u{A}_k)$$

$$=-n/2\log ({\hat{g}}_{uk}^T{\hat{M}}_2{\hat{g}}_{uk})-n/2\log \left\{{\hat{g}}_{uk}^T{\left({\hat{M}}_2+{\hat{U}}_2\right)}^{-1}{\hat{g}}_{uk}\right\}$$

$$+n\log ({\hat{g}}_{uk}^T{\hat{g}}_{uk})$$

$${\hat{g}}_{\text{env}}=\frac{1}{B}\sum _{b=1}^{B}g\left\{{\hat{\tau }}_{\text{env}}^{(b)}\right\}.$$

Full text

Darwinian fitness; fitness landscape; envelope model; parametric bootstrap

1 Introduction

The estimation of expected Darwinian fitness is a very important procedure in both biology and genetics. The importance of this is not just limited to scientific disciplines, it is important for public policy. With genetic theory and simulation studies, Burger and Lynch (1995) shows that, under certain conditions, a changing environment leads to extinction of species. In a field study, Etterson and Shaw (2001) argued that the predicted evolutionary response to predicted rates of climate change are far too slow. In these papers, and all life history analyses of their kind, expected Darwinian fitness is the response variable. The interesting scientific conclusions are drawn from it.

In many life history analyses, values of expected Darwinian fitness are plotted using a fitness landscape (Lande and Arnold, 1983; Shaw and Geyer, 2010). A fitness landscape is the conditional expectation of Darwinian fitness given phenotypic trait values considered as a function of those values. When fitness is the response variable in a regression model and phenotypic traits are the covariates, the fitness landscape is the regression function. Estimation of the fitness landscape began with Lande and Arnold (1983). They use ordinary least squares regression of fitness on phenotypes to estimate the best linear approximation of the fitness landscape and quadratic regression to estimate the best quadratic approximation. Here “best” means minimum variance unbiased, as in the Gauss-Markov theorem. Their use of $t$ and $F$ tests and confidence intervals requires the assumption that fitness is conditionally homoscedastically normally distributed given phenotypic trait values. This assumption is almost always grossly incorrect when one uses a good surrogate for Darwinian fitness (Mitchell-Olds and Shaw, 1987; Shaw et al., 2008).

Aster models (Geyer et al., 2007) were designed to fix all of the problems of the Lande and Arnold (1983) approach and of all other approaches to life history analysis (Shaw et al., 2008). The aster model is the state-of-the-art model for all life history analyses in which the estimation of expected Darwinian fitness is the primary goal. Geyer et al. (2007); Shaw et al. (2008); Stanton-Geddes et al. (2012); Shaw et al. (2015); Eck et al. (2015) show various kinds of life history data for which aster models are necessary. Assumptions for aster models are given in Section 2 below.

In this article we combine envelope methodology (Cook et al., 2010; Cook and Zhang, 2015; Su and Cook, 2011) with aster models (Geyer et al., 2007; Shaw et al., 2008) to estimate the fitness landscape (Lande and Arnold, 1983; Shaw and Geyer, 2010) in life history analysis. The primary emphasis is that this combination of methods estimates the fitness landscape with less variability than is possible with aster models alone. We first show how existing envelope estimators constructed from the 1D algorithm (Cook and Zhang, 2015, 2016; Zhang and Mai, 2017) can reduce variability in estimation of the fitness landscape. We then develop a new envelope estimator that avoids the potential numerical pitfalls of the 1D algorithm. Variance reduction is assessed using parametric bootstrap techniques in Efron (2014, Section 4). These bootstrap algorithms account for variability in model selection. Our methodology provides the most precise estimation of expected Darwinian fitness to date. Researchers using our methods can therefore draw stronger conclusions about the driving forces of Darwinian fitness from their life history analyses.

In a life history analysis of M. guttatus flowers and a simulated example, we show that our methodology leads to variance reduction in estimation of expected Darwinian fitness when compared with analyses that use aster models alone. We show that this variance reduction leads to sharper scientific inferences about the potential causes of Darwinian fitness in the M. guttatus life history analysis. Our examples are fully reproducible, and the calculations necessary for their reproduction are included in an accompanying technical report (Eck et al., 2018).

2 The aster model

Aster models are regular full exponential families. Parameters are estimated by maximum likelihood. If $Y$ is the response vector, $\mu=E(Y)$ is the saturated model mean value parameter, and $M$ is the model matrix for an unconditional canonical affine submodel, then the maximum likelihood estimate of $\mu$ satisfies $M^{T}\hat{\mu}=M^{T}y$ (Geyer et al., 2007, Section 3.2). $\tau=M^{T}\mu$ is the submodel mean value parameter (Geyer et al., 2007, Section 2.4). Likelihood ratio tests for model comparison and confidence intervals for all parameters are based on the usual asymptotics and Fisher information are provided by R package aster (Geyer, 2017a). In particular, for this article we need to know that $\hat{\tau}=M^{T}y$ is is a minimum variance unbiased estimator of the parameter it estimates and its exact variance matrix is the Fisher information matrix for the submodel canonical parameter $\beta$ (the vector of regression coefficients). That is, in the usual asymptotics of maximum likelihood for this parameter the mean and variance are exact not approximate; only the normal distribution is approximate.

The aster model is a directed acyclic graphical model (Lauritzen, 1996, Section 3.2.3) in which the joint density is a product of conditional densities that are specified by the arrows depicted in the graph. Lines that appear in the graph specify nodes which are dependent. For example, an organism may have multiple paths in their life history and can only go down one of them (Eck et al., 2015). The aster model follows five assumptions which are: A1 The graph of arrows is acyclic. A2 In the graph of lines every connected component is a complete graph, which is called a dependence group. A3 Every node in a dependence group with more than one node has the same predecessor (there is an arrow from the predecessor to each node in the group). Every dependence group consisting of exactly one node has at most one predecessor. A4 The joint distribution is the product of conditional distributions, one conditional distribution for each dependence group. A5 Predecessor is sample size, meaning each conditional distribution is the distribution of the sum of $N$ independent and identically distributed random vectors, where $N$ is the value of the predecessor, the sum of zero terms being zero. A6 The conditional distributions are exponential families having the components of the response vector for the dependence group as their canonical statistics. Assumptions A5 and A6 mean for an arrow $y_{k}\rightarrow y_{j}$ that $y_{j}$ is the sum of independent and identically distributed random variables from the exponential family for the arrow and there are $y_{k}$ terms in the sum (the sum of zero terms is zero). These assumptions imply that the joint distribution of the aster model is an exponential family (Geyer et al., 2007, Section 2.3). Three of these assumptions have a clear biological meaning as well. Assumptions A1 through A3 restricts an individual from revisiting life stages that have come to pass. Assumption A5 implies that dead individuals remain dead and have no offspring though the course of the study.

These aster models are saturated, having one parameter per component of the response vector, and are not useful. Hence, as with linear and generalized linear models, canonical affine submodels are used with change of parameter $\varphi=a+M\beta$, where $\varphi$ is the saturated model parameter vector (linear predictor in the terminology of generalized linear models), $\beta$ is the submodel parameter (“coefficients” in R terminology), $a$ is the offset vector, and $M$ is the model matrix. The aster submodel has log likelihood $l(\beta)=\langle M^{T}Y,\beta\rangle-c(a+M\beta)$ where $Y$ is the response vector and $c(\cdot)$ is the cumulant function of the exponential family.

There are three parameters of interest that are present in the aster analyses we consider (see Geyer (2010) and Eck et al. (2018) for more detail on aster model parameterizations). These parameterizations are

1. the aster submodel canonical parameter vector $\beta\in\mathbb{R}^{p}$,
2. the aster submodel mean-value parameter vector $\tau\in\mathbb{R}^{p}$,
3. the saturated aster model mean-value parameter vector $\mu\in\mathbb{R}^{m}$, where $m$ is the number of individuals sampled multiplied by the number of nodes in the aster graph.

These three parameterizations are all linked via invertible 1-1 transformations when the model matrix $M$ is of full rank. The usual asymptotics of maximum likelihood estimation give

[TABLE]

where $\Sigma=\operatorname{Var}(M^{T}Y)$ is the Fisher information matrix associated with the canonical parameter vector $\beta$. The maximum likelihood estimator of $\beta$ is asymptotically normal with variance given by $\Sigma^{-1}$. From (1) and the delta method we can obtain the asymptotic distribution for any differentiable function of $\hat{\tau}$. The asymptotic distribution for a differentiable function $g$ of $\hat{\tau}$ is

[TABLE]

In particular, the asymptotic distribution of estimated expected Darwinian fitness is of interest. Let $h(\mu)$ be expected Darwinian fitness. Both $\beta$ being a function of $\tau$ and $\mu=\nabla c(a+M\beta)$ imply that g(\tau)=h\Bigl{[}\nabla c\bigl{\{}a+Mf(\tau)\bigr{\}}\Bigr{]} is expected Darwinian fitness as a function of $\tau$ and is differentiable if $h$ is differentiable. The estimator $g(\hat{\tau})$ has asymptotic distribution given by (2). We have the potential to do better through the incorporation of envelope methodology.

3 Incorporation of envelope methodology

Aster model estimates of expected Darwinian fitness may be too variable to be useful, and in consequence we may not be able to statistically distinguish estimates of expected Darwinian fitness over the fitness landscape. We address this problem through the incorporation of envelope methodology into the aster modelling framework.

Envelope models were developed originally as a variance reduction tool for the multivariate linear regression model. In this article we focus on envelope methodology for general vector-valued parameter estimation (Cook and Zhang, 2015). Envelope methodology has the potential to reduce the variability of any $\surd{n}$-consistent asymptotically normal distributed consistent estimator. We now define the envelope subspace. Let $\upsilon$ be a parameter of interest and suppose that $\hat{\upsilon}$ is a $\surd{n}$-consistent estimator of $\upsilon$ with asymptotic covariance matrix $\Sigma_{\upsilon,\upsilon}$. Let $\mathcal{T}=\text{span}(\upsilon)=\{a\upsilon:a\in\mathbb{R}\}.$ The envelope subspace $\mathcal{E}_{\Sigma_{\upsilon,\upsilon}}(\mathcal{T})$ is defined as the intersection of all reducing subspaces of $\Sigma_{\upsilon,\upsilon}$ that contain $T$ (a reducing subspace is a sum of eigenspaces if all eigenvalues of $\Sigma$ have multiplicity one). The envelope space satisfies both

[TABLE]

where $P_{\mathcal{E}}$ is the projection into the envelope subspace and $Q_{\mathcal{E}}$ is the projection into the orthogonal complement. In coordinate form, these two envelope conditions are

[TABLE]

where $(\Gamma,\Gamma_{o})$ is a partitioned orthogonal matrix, the columns of $\Gamma$ are a basis for $\mathcal{E}_{\Sigma_{\upsilon,\upsilon}}(\mathcal{T})$, and the dimensions of the positive definite matrices $\Omega$ and $\Omega_{o}$ are such that the matrix multiplications are defined. The quantities $P_{\mathcal{E}}\Sigma_{\upsilon,\upsilon}P_{\mathcal{E}}$ and $\Gamma\Omega\Gamma^{T}$ are often referred to as ’material information’ in the envelope literature (Cook et al., 2010; Cook and Zhang, 2015) since they represent the portion of variability that is necessary for the task of estimating $\upsilon$. Similarly, $Q_{\mathcal{E}}\Sigma_{\upsilon,\upsilon}Q_{\mathcal{E}}$ and $\Gamma_{o}\Omega_{o}\Gamma_{o}^{T}$ are referred to as ’immaterial information’ since they represent extraneous variability.

Intuitively, the envelope estimator reduces variability in estimation at no cost to consistency. An illuminating depiction and explanation of how an envelope analysis increases efficiency in multivariate linear regression problems was given by Su and Cook (2011, pgs. 134–135). The same intuition applies to envelope methodology of Cook and Zhang (2015). In applications, there is a cost to estimate $u=\text{dim}\left\{\mathcal{E}_{\Sigma_{\upsilon,\upsilon}}(\mathcal{T})\right\}$ and $\Gamma$. With the basis matrix $\Gamma$ estimated, we can then assess the variance reduction of the envelope estimator through the parametric bootstrap.

The 1D algorithm (Cook and Zhang, 2015, Algorithm 2) estimates a basis matrix $\Gamma$ for $\mathcal{E}_{\Sigma_{\upsilon,\upsilon}}(\mathcal{T})$ at a user-supplied envelope dimension $u$. The estimate of $\Gamma$ is obtained by providing $\widehat{\Sigma}_{\upsilon,\upsilon}$ and $\hat{\upsilon}\hat{\upsilon}^{T}$ as inputs into the 1D algorithm. The resulting estimator of $\Gamma$ obtained from the 1D algorithm, $\widehat{\Gamma}_{u}$, gives a $\surd{n}$ consistent estimator $P_{\hat{\mathcal{E}}}$ of the projection onto the envelope subspace $P_{\mathcal{E}}$ (Cook and Zhang, 2015). The 1D algorithm can be used to estimate $u$ consistently (Zhang and Mai, 2017, Theorem 2).

Define $\hat{\upsilon}_{\text{1D}}=P_{\hat{\mathcal{E}}}\hat{\upsilon}$ to be the envelope estimator of $\upsilon$ where $\upsilon$ are parameters that link the estimation of expected Darwinian fitness to covariates of interest. We write $\tau=(\gamma^{T}$, $\upsilon^{T})^{T}$ where $\gamma$ are aster model parameters not linking covariates to the estimation of Darwinian fitness. The envelope estimator of $\tau$ is given as

[TABLE]

The model matrix $M_{\text{1D}}$ corresponds to the aster model that incorporates the envelope structure.

Proposition 1**.**

The envelope estimator $\hat{\tau}_{\text{1D}}$ is a maximum likelihood estimator of $\tau$ for the aster model with model matrix $M_{\text{1D}}$.

Proof.

We have $l_{\text{env}}(\beta)=\langle Y,M_{\text{1D}}\beta\rangle-c(M_{\text{1D}}\beta)=\langle M_{\text{1D}}^{T}Y,\beta\rangle-c(M_{\text{1D}}\beta).$ and $\nabla_{\beta}l(\beta)=M_{\text{1D}}^{T}Y-\nabla_{\beta}c(M_{\text{1D}}\beta).$ Setting $\nabla_{\beta}l(\beta)=0$ and solving for $\beta$ yields $\nabla_{\beta}c(M_{\text{1D}}\beta)|_{\beta=\hat{\beta}}=M_{\text{1D}}^{T}Y=\hat{\tau}_{\text{1D}}$. ∎

This proposition justifies the use of the transformations to switch between maximum likelihood estimators of the different aster model parameterizations. We compare envelope dimensions $u$ by transforming envelope estimators of $\tau$ to envelope estimators of $\beta$ and then evaluate the log likelihood at the envelope estimator of $\beta$. The randomness inherent in $M_{\text{1D}}$ is non-problematic in most applications. This is because the 1D algorithm provides a $\surd{n}$-consistent envelope estimator of $\upsilon$ even when $u$ is estimated. Inferences about aster model parameters implicitly assumes that $n$ is large enough for the asymptotic normality to be a good approximation for the distribution of maximum likelihood estimators.

We then compute an envelope estimator of expected Darwinian fitness $g(\tau)$ using the aster model with model matrix $M_{\text{1D}}$. Note that the model matrix $M_{\text{1D}}$ is not of full column rank. Therefore the transformations used to switch between aster model parameterizations are not 1-1. In particular, many distinct estimates of $\beta$ map to $\hat{\tau}_{\text{1D}}$. Each of these distinct estimated values of $\beta$ maps to the same estimate of $M_{\text{1D}}\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.

Our estimator of estimated expected Darwinian fitness is given by Efron (2014, equation (4.5)), with $g(\cdot)$ replacing $t(\cdot)$. Steps for this algorithm are given in Algorithm 1. When the top-level of our bootstrap procedure (Steps 1 through 4 in Algorithm 1) has run for $B$ iterations, we obtain the envelope estimator

[TABLE]

The envelope estimator $\hat{g}_{\text{1D}}$ implicitly behaves as a weighted average with the weights reflecting the likelihood of observing a particular estimated value of $u$, $\hat{\upsilon}_{\text{1D}}$, and $M_{\text{1D}}$. The value of $u$ is estimated with either the Bayesian information criterion or Akaike information criterion at every iteration of the parametric bootstrap. The intuition is that the averaging in $\hat{g}_{\text{1D}}$ will smooth out variability due to the estimating $u$ with our chosen model selection criterion. Variability of $\hat{g}_{\text{1D}}$ is estimated using the double bootstrap technique in Efron (2014, Section 4), with steps shown below in Algorithm 1. This bootstrap technique accounts for all estimation error, including model selection error. The reason it accounts for all estimation error, is that all estimation, including model selection is done in each iteration of the bootstrap (nothing estimated is ever treated as known in bootstrap iterations).

Algorithm 1. Parametric bootstrap for assessing the variability of the envelope estimator $\hat{g}_{\text{1D}}$:

Fit the aster model to the data and obtain $\hat{\upsilon}$ and $\widehat{\Sigma}_{\upsilon,\upsilon}$ from the aster model fit. 2. 2.

Choose a model selection criterion. Compute the envelope estimator of $\upsilon$ in the original sample, given as $\hat{\upsilon}_{\text{1D}}=P_{\hat{\mathcal{E}}}\hat{\upsilon}$ where $P_{\hat{\mathcal{E}}}$ is obtained from the 1D algorithm and the chosen model selection criterion. 3. 3.

Perform a parametric bootstrap by generating samples from the distribution of the aster submodel evaluated at $\hat{\tau}_{\text{1D}}=(\hat{\gamma}^{T},\hat{\upsilon}_{\text{1D}}^{T})^{T}$. For $b=1$, $\ldots$, $B$ of the procedure:

1. (3a)

Compute $\hat{\tau}^{(b)}$ and $\widehat{\Sigma}_{\upsilon,\upsilon}^{(b)}$ from the aster model fit to the resampled data. 2. (3b)

Compute $P_{\hat{\mathcal{E}}}^{(b)}$ as done in Step 2. 3. (3c)

Compute $\hat{\tau}_{\text{1D}}^{(b)}=\left\{\hat{\gamma}^{(b)^{T}},\hat{\upsilon}_{\text{1D}}^{(b)^{T}}\right\}^{T}$ and $g\left\{\hat{\tau}_{\text{1D}}^{(b)}\right\}$. 4. 4.

The bootstrap estimator of expected Darwinian fitness is the average of the envelope estimators computed in Step 3c. This completes the first part of the bootstrap procedure. 5. 5.

At $k=1$, $\ldots$, $K$, for each $b=1$, $\ldots$, $B$ we:

1. (5a)

Generate data from the distribution of the aster submodel evaluated at $\hat{\tau}_{\text{1D}}^{(b)}$. 2. (5b)

Perform Steps 3a through 3c with respect to the dataset obtained in Step 5a to calculate both $\hat{\tau}_{\text{1D}}^{(b)^{(k)}}$ and $g\left\{\hat{\tau}_{\text{1D}}^{(b)^{(k)}}\right\}$. 6. 6.

Compute both $\hat{g}_{\text{1D}}$ and the standard deviation in Efron (2014, equation (4.15)).

4 A direct envelope estimator using reducing subspaces

We propose a new way of constructing envelope estimators provided that the eigenvalues of $\Sigma_{\upsilon,\upsilon}$ have multiplicity one. In this section, envelope estimators are constructed directly from the reducing subspaces of $\widehat{\Sigma}_{\upsilon,\upsilon}$. Uniqueness of the eigenvalues of $\Sigma_{\upsilon,\upsilon}$ implies that its reducing subspaces are sums of its eigenspaces. Let $\mathcal{G}$ be a reducing subspace of $\Sigma_{\upsilon,\upsilon}$. Define $\Gamma_{\mathcal{G}}$ and $P_{\mathcal{G}}$ as the basis matrix for $\mathcal{G}$ and the projection onto $\mathcal{G}$ respectively. Let $\widehat{\mathcal{G}}$ be the estimator of the reducing subspace $\mathcal{G}$, obtained by searching over the 1 dimensional eigenspaces of $\widehat{\Sigma}_{\upsilon,\upsilon}$. In applications, eigenvalues of $\widehat{\Sigma}_{\upsilon,\upsilon}$ are almost always unique. Define $\widehat{\Gamma}_{\widehat{\mathcal{G}}}$ and $\widehat{P}_{\widehat{\mathcal{G}}}=\widehat{\Gamma}_{\widehat{\mathcal{G}}}\widehat{\Gamma}^{T}_{\widehat{\mathcal{G}}}$ as estimators of $\Gamma_{\mathcal{G}}$, and $P_{\mathcal{G}}$ respectively. The basis matrix $\widehat{\Gamma}_{\widehat{\mathcal{G}}}$ is constructed from the eigenvectors of $\widehat{\Sigma}_{\upsilon,\upsilon}$. We now define the envelope estimator of $\upsilon$ constructed from reducing subspaces.

Definition 1**.**

*The envelope estimator of $\upsilon$ constructed from the reducing subspaces $\mathcal{G}$ is defined to be $\hat{\upsilon}_{\text{env}}=\widehat{P}_{\widehat{\mathcal{G}}}\hat{\upsilon}$. *

The reducing subspaces of $\widehat{\Sigma}_{\upsilon,\upsilon}$ are $\surd{n}$ consistent estimators of the reducing subspaces of $\Sigma_{\upsilon,\upsilon}$. Therefore $\widehat{\Gamma}_{\widehat{\mathcal{G}}}$, $\widehat{P}_{\widehat{\mathcal{G}}}$, and the corresponding estimator $\widehat{P}_{\widehat{\mathcal{G}}}\hat{\upsilon}$ are $\surd{n}$ consistent estimators of $\Gamma_{\mathcal{G}}$, $P_{\mathcal{G}}$, and $\upsilon$ respectively. The envelope estimator of $\tau$ constructed by reducing subspaces is given by

[TABLE]

The model matrix $M_{\text{env}}$ corresponds to the aster model that incorporates the envelope structure obtained from the reducing subspaces of $\widehat{\Sigma}_{\upsilon,\upsilon}$. We have a similar result as Proposition 1 for the aster model with model matrix $M_{\text{env}}$.

Proposition 2**.**

The envelope estimator $\hat{\tau}_{\text{env}}$ is a maximum likelihood estimator of $\tau$ for the aster model with model matrix $M_{\text{env}}$.

The proof of Proposition 2 follows the same steps as the proof for Proposition 1. There is a close connection between envelope estimation using reducing subspaces and envelope estimation using the 1D algorithm. In the population, $\hat{\upsilon}_{\text{env}}=\hat{\upsilon}_{\text{1D}}$. The connection between both estimation methods exists in finite samples as seen in Theorem 1. In preparation, define orthogonal matrices $\widehat{O}_{u}=\left(\widehat{\Gamma}_{u},\widehat{\Gamma}_{uo}\right)$, $\widehat{O}_{\widehat{\mathcal{G}}}=\left(\widehat{\Gamma}_{\widehat{\mathcal{G}}},\widehat{\Gamma}_{\widehat{\mathcal{G}}o}\right),$ and $\widehat{O}=\widehat{O}_{\widehat{\mathcal{G}}}\widehat{O}_{u}^{T}$. The matrices $\widehat{\Gamma}_{\widehat{\mathcal{G}}}$, $\widehat{\Gamma}_{\mathcal{G}o}$, $\widehat{\Gamma}_{u}$, and $\widehat{\Gamma}_{uo}$ converges in probability to $\Gamma_{\mathcal{G}}$, $\Gamma_{\mathcal{G}o}$, $\Gamma_{u}$, and $\Gamma_{uo}$ respectively. Now $\Gamma_{\mathcal{G}}^{T}\Gamma_{u}$ and $\Gamma_{\mathcal{G}o}^{T}\Gamma_{uo}$ are both [math]-$1$ valued rotation matrices and $\Gamma_{\mathcal{G}}^{T}\Gamma_{uo}=\Gamma_{\mathcal{G}o}^{T}\Gamma_{u}=0$. These facts imply that $\widehat{O}$ converges in probability to a [math]-$1$ valued rotation matrix. We will assume that $\widehat{O}\overset{p}{\to}I$ without loss of generality.

Theorem 1**.**

The basis matrix $\widehat{\Gamma}_{\widehat{\mathcal{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$.

Proof.

Let $\widehat{M}_{2}=\widehat{\Sigma}_{\upsilon,\upsilon}$ and $\widehat{U}_{2}=\hat{\upsilon}\hat{\upsilon}^{T}$. Similar to the proof of (Cook and Zhang, 2016, Proposition 6), let

[TABLE]

Now let $\hat{v}_{k}$, $k=1$, $\ldots$, $u$ be the $k$th column of $\widehat{\Gamma}_{\mathcal{G}}$ and define $A\in\mathbb{R}^{p\times p}$ to be a matrix of [math]’s with $1$’s occupying the first $k$ diagonal entries. Then

[TABLE]

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

Theorem 1 in combination with Zhang and Mai (2017, Theorem 2) allows for us to estimate $\mathcal{G}$ consistently. Thus the variability associated with the estimation of $\hat{\upsilon}_{\text{env}}$ and $\mathcal{G}$ decreases as $n\to\infty$. However in practical applications correct model selection cannot be guaranteed. Therefore the envelope estimator of expected Darwinian fitness $g(\hat{\tau}_{\text{env}})$ has an extra source of variability due to model selection uncertainty. We develop a double bootstrap procedure with steps similar to those in Algorithm 1 to account for variability in model selection. The first level of the bootstrap procedure provides the estimator of expected Darwinian fitness,

[TABLE]

The same model selection criteria is used to select the reducing subspace used to construct $\hat{\tau}_{\text{env}}^{(b)}$ at every iteration $b=1$, $\ldots$, $B$. The second of level of this bootstrap procedure estimates the variability of (4). The steps for this algorithm are provided in Eck et al. (2018). The utility of our double bootstrap procedure is shown in Section 5.1 where, in that example, there is considerable disagreement between model selection criteria.

When $k$ is small, $\hat{\upsilon}_{\text{env}}$ is preferable to $\hat{\upsilon}_{\text{1D}}$. At any iteration of the 1D algorithm, minimizers of the objective function stated in (Cook and Zhang, 2015, Algorithm 2) are pulled towards reducing subspaces of $\widehat{\Sigma}_{\upsilon,\upsilon}$. This objective function is non-convex and contains potentially many local minima. The optimizations conducted within the 1D algorithm are sensitive to starting values and can get stuck at these local minima. This undermines the 1D algorithm since it is required that users find global minima for its justification. Unlike the 1D algorithm, the reducing subspace approach does not involve any optimization routines. However when $k$ is moderately large, the computation of all of candidate envelope estimators at each reducing subspace is too computationally intensive, there are $2^{k}-1$ possible reducing subspaces in non trivial problems. The 1D algorithm may still be fast when $k$ is moderately large, only $k-1$ optimizations in non trivial problems. In most aster applications $k$ is small since data is obtained through expensive collection methods and fitness landscapes are low-dimensional (Shaw and Geyer, 2010; Eck et al., 2015).

5 Examples

5.1 Simulated Data

A population of 3000 organisms was simulated to form the dataset used in this aster analysis. These data were generated according to the graphical structure appearing in panel A of Figure 1. There are two covariates $(z_{1},z_{2})$ associated with Darwinian fitness and the aster model selected by the likelihood ratio test is a full quadratic model with respect to these covariates.

We partition $\tau$ into $(\gamma^{T},\upsilon^{T})^{T}$ where $\gamma\in\mathbb{R}^{4}$ are nuisance parameters and $\upsilon\in\mathbb{R}^{5}$ are relevant to the estimation of expected Darwinian fitness. Here, $\upsilon\in\mathbb{R}^{5}$ because our model is full quadratic in $z_{1}$ and $z_{2}$. In this example, the true reducing subspace is the space spanned by the first and fourth eigenvectors of the covariance matrix of the parameters of interest estimated from the original data. We begin by considering envelope estimators constructed using the 1D algorithm. The Akaike information criterion and the Bayesian information criterion both select $u=5$. Envelope methods are not interesting in this case.

We now consider envelope estimators constructed from reducing subspaces. In the original sample, the Bayesian information criterion selects the reducing subspace that is the sum of the first, fourth, and fifth eigenspaces of $\widehat{\Sigma}_{\upsilon,\upsilon}$ numbered in order of decreasing eigenvalues. This suggests that the dimension of the envelope space is $u=3$. The 1D algorithm and the reducing subspace approach are in disagreement, consistency of model selection is not helpful in this application. We turn to the double parametric bootstrap to estimate $\hat{g}_{\text{env}}$ and the asymptotic variability of $\hat{g}_{\text{env}}$. The Bayesian information criterion is used to select $\mathcal{G}$ at every iteration of the first level of the bootstrap.

The results are seen in Table 1. Table 1 shows seven individuals that had high values of estimated expected Darwinian fitness. Each individual has a unique set of traits. The first two columns display the envelope estimator of expected Darwinian fitness $\hat{g}_{\text{env}}$ and its bootstrapped standard error. The maximum likelihood estimator of expected Darwinian fitness and its bootstrapped standard error are displayed in the third and fourth columns respectively. The ratios of bootstrapped standard errors for $\hat{g}_{\text{MLE}}$ to $\hat{g}_{\text{env}}$ are displayed in the final column. We can see that all of the ratios are greater than 1 which indicates that the envelope estimator of expected Darwinian fitness is less variable than the maximum likelihood estimator.

Contour plots of the ratios of estimated standard errors are displayed in the technical report (Eck et al., 2018). These contour plots show that the envelope estimator of expected Darwinian fitness is less variable than the maximum likelihood estimator for the majority of the observed data. The region where the envelope estimator is less variable includes the values of $z_{1}$ and $z_{2}$ that maximize estimated expected Darwinian fitness. Variance reduction is also obtained when we use the reducing subspace suggested by the Akaike information criterion. This is also shown in Eck et al. (2018).

5.2 M. guttatus aster analysis

The yellow monkeyflower M. guttatus has been and is currently a well studied flower (Lowry and Willis, 2010; Hall et al., 2006; Ritland and Ritland, 1989; Allen and Sheppard, 1971; Vickery, 1959). M. guttatus is a species which comprises many morphologically variable populations growing in moist places such as stream banks, meadows and springs over a range that extends from the Aleutian Islands to Mexico and from the California coast to the Rocky Mountains (Vickery, 1959). The lifecycle of the individual M. guttatus flowers, for our life history analysis, is depicted in panel B of Figure 1. Lowry and Willis (2010) performed a life history analysis of M. guttatus using aster models. One of their interests was to determine which levels of genetic background, field site, inversion orientation, and ecotype of the flower are associated with high Darwinian fitness. We show that the set of candidate trait values thought to maximize expected Darwinian fitness is smaller when envelope methodology is incorporated.

Lowry and Willis (2010) collected measurements on 2313 M. guttatus. We fit a linear fitness landscape to this data. The parameters $\upsilon\in\mathbb{R}^{6}$ are relevant to the estimation of expected Darwinian fitness. In the original sample, the Bayesian information criterion leads to a selection of a reducing subspace that is the sum of all eigenspaces of $\widehat{\Sigma}_{\upsilon,\upsilon}$ with the exception of the fourth and fifth eigenspaces. The parametric double bootstrap procedure outlined in Sections 3 and 4 is used to estimate the variability of $\hat{g}_{\text{env}}$ where the Bayesian information criterion is used to select $\mathcal{G}$ at every iteration of the first level of the bootstrap.

Table 2 shows the results for seven individuals that have high values of estimated expected Darwinian fitness through maximum likelihood and envelope estimation. Each individual has a unique set of traits. We see that both methods agree on the trait values that are expected to maximize expected Darwinian fitness. We also see that all of the ratios are greater than 1. 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 Table 2 from the fifth row of Table 2. 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.

6 Software

This paper is accompanied by an R package envlpaster (Eck, 2016), which requires the two R packages for aster models: aster (Geyer, 2017a) and aster2 (Geyer, 2017b), and also a technical report (Eck et al., 2018) that reproduces the examples in this paper and shows how functions in the envlpaster package are used.

7 Discussion

One could think to perform envelope methodology with respect to the regression coefficients $\beta$ instead of $\tau$. However $\beta$ is not well-defined, one can shift $\beta$ with an arbitrarily chosen offset vector without changing the value of the mean-value parameters $\tau$ and $\mu$. In addition, whenever we have categorical predictors, R software automatically drops one category when it has an intercept in the formula, but which category it drops is arbitrary and changing which is dropped changes $\beta$. Envelope methodology is not invariant to this form of arbitrary shifting. To the best of our knowledge, the exponential family regression applications in the envelope model literature exclusively seek inference about $\beta$ because this well-definedness issue is not a problem in those applications (Cook et al., 2010; Su and Cook, 2011, 2012; Cook and Su, 2013; Cook and Helland, 2013; Cook and Zhang, 2015; Cook and Su, 2016; Park et al., 2017; Eck and Cook, 2017). The applications of envelope methods to aster models is therefore outside of the scope of previous applications to exponential family models. Additionally, the methods in this paper can be extended to functions of parameters in generalized linear regression models. Aster model are a generalization of generalized linear regression models (Shaw et al., 2008).

The consequences of potential model selection errors served as the motivation for the implementation of the bootstrap procedure in Eck and Cook (2017). In that article, inferences are only given for a canonical parameter vector in a multivariate linear regression model. We applied the Efron (2014) bootstrap procedures to alleviate possible model selection concerns. However, this particular choice of a bootstrap procedure is not without flaws. Hjort (2014) mentions that Efron does not derive the asymptotic distribution of the final estimator. The literature has not reached a consensus on the appropriate bootstrap procedure to be implemented when bootstrapping depends on data-driven model selection. As the literature currently stands, Efron (2014) provides a reasonable solution to the problem of potential model selection errors in the application of envelope methodology to aster models. The parametric bootstrap does not rely on asymptotic normality since it simulates the exact sampling distribution of the estimator for some parameter value, and the double bootstrap simulates the exact sampling distribution of the estimator for a long list of parameter values (Geyer, 2013).

Our new envelope estimator does not involve any non-convex optimization routines that are both sensitive to starting values and have potential problems with local minima. These computational problems can be detrimental to the performance of the 1D algorithm. The underlying theory of the 1D algorithm justifies the consistency properties of our new envelope estimator. In envelope modelling problems with a small number of parameters of interest the envelope estimator constructed directly from reducing subspaces is preferred since it possesses the same strengths as the 1D algorithm without its potential numerical pitfalls. However our estimator is currently expensive to compute in moderate $p$ problems. In aster analyses $p$ is typically small.

In many life history analyses, specific trait values which are estimated to produce the highest expected Darwinian fitness are of interest. It is common practice to only report such trait values (Shaw and Geyer, 2010; Eck et al., 2015). Such reporting ignores the variability associated with the estimation of expected Darwinian fitness. There are likely many trait values having estimated expected Darwinian fitness that is statistically indistinguishable from the reported values. Our methodology addresses this concern directly. The potential set of candidate traits associated with high values of expected Darwinian fitness is smaller when the combination of envelope methodology into the aster modelling framework is utilized as seen in Eck et al. (2018). Researchers using our methods will have the potential to make stronger inferences about expected Darwinian fitness through our variance reduction techniques.

Acknowledgements

Daniel J. Eck’s research is supported by NIH/NIHCD grant 1DP2HD091799-01. We would like to thank David B. Lowry for providing the dataset used in Example 2, Xin Zhang for the code that implements the 1D algorithm, and Forrest W. Crawford for helpful discussion that led to the strengthening of this paper. We would also like to especially thank Amber Eule-Nashoba for helpful comments on the technical report.