Combining Envelope Methodology and Aster Models for Variance Reduction in Life History Analyses
Daniel J. Eck Charles J. Geyer R. Dennis Cook


























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




Daniel J. Eck 




Department of Biostatistics, Yale School of Public Health. 

daniel.eck@yale.edu




Charles J. Geyer     and     R. Dennis Cook




Department of Statistics, University of Minnesota

















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 Mimulus guttatus flowers is provided. Useful variance 
reduction is obtained in both analyses.



Darwinian fitness; fitness landscape; envelope model; parametric bootstrap 













§ 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, 
<cit.> shows that, under certain conditions, a changing environment 
leads to extinction of species. In a field study, <cit.> 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 <cit.>. 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 <cit.>. 
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 <cit.>.

Aster models <cit.> were designed to fix all of the problems of the 
<cit.> approach and of all other approaches to life history analysis 
<cit.>. 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. <cit.> 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 <cit.> with 
aster models <cit.> to estimate the fitness landscape 
<cit.> 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 <cit.> 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 <cit.>. 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 
<cit.>. 







































































































































































§ THE ASTER MODEL



Aster models are regular full exponential families.  Parameters are 
estimated by maximum likelihood.  If Y is the response vector, 
μ = 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 μ satisfies 
M^T μ̂ = M^T y <cit.>. 
τ = M^T μ is the submodel mean value parameter 
<cit.>.  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  <cit.>.  In particular, for this 
article we need to know that τ̂ = 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 β (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 
<cit.> 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 <cit.>. 


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 → 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
<cit.>. 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 φ = a + M β, where φ is the saturated model 
parameter vector (linear predictor in the terminology of generalized linear 
models), β 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(β) = ⟨ M^TY, β⟩ - c(a + Mβ)
where Y is the response vector and c(·) is the cumulant function of 
the exponential family. 







There are three parameters of interest that are present in the aster analyses
we consider (see <cit.> and <cit.> for more detail on aster model 
parameterizations). These parameterizations are 







1) the aster submodel canonical parameter vector β∈^p,
2) the aster submodel mean-value parameter vector τ∈^p,
3) the saturated aster model mean-value parameter vector μ∈^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  

    √n(τ̂ - τ) d⟶ 
        N(0, Σ),

where Σ = (M^T Y) is the Fisher information matrix associated with 
the canonical parameter vector β. The maximum likelihood estimator of 
β is asymptotically normal with variance given by Σ^-1. 
From (<ref>) and the delta method we can obtain the 
asymptotic distribution for any differentiable function of τ̂. 
The asymptotic distribution for a differentiable function g of τ̂ 
is

    √n{g(τ̂) - g(τ)}d⟶ 
        N{0, ∇ g(τ)Σ∇ g(τ)^T}.

In particular, the asymptotic distribution of estimated expected Darwinian 
fitness is of interest. Let h(μ) be expected Darwinian fitness.  
Both β being a function of τ and μ = ∇ c(a + M β) 
imply that 
g(τ) =  h[∇ c{a + M f(τ) }]
is expected Darwinian fitness as a function of τ and is differentiable 
if h is differentiable. 


The estimator g(τ̂) has asymptotic distribution given by 
(<ref>). We have the potential to do better through the 
incorporation of envelope methodology.



§ 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 <cit.>.
Envelope methodology has the potential to reduce the 
variability of any √n-consistent asymptotically normal distributed 


consistent estimator. 






We now define the envelope subspace.
Let υ be a parameter of interest and suppose that υ̂ is 
a √n-consistent estimator of υ with asymptotic covariance 
matrix Σ_υ,υ. Let 
= span(υ) ={aυ : a ∈}. 
The envelope subspace ℰ_Σ_υ,υ() is 
defined as the intersection of all reducing subspaces of 
Σ_υ,υ that contain T (a reducing subspace is a sum of 
eigenspaces if all eigenvalues of Σ have multiplicity one). 


The envelope space satisfies both

    ⊂ℰ_Σ_υ,υ(),
            Σ_υ,υ 
        = P_ℰΣ_υ,υP_ℰ 
        + Q_ℰΣ_υ,υQ_ℰ;

where P_ℰ is the projection into the envelope subspace and 
Q_ℰ is the projection into the orthogonal complement. In 
coordinate form, these two envelope conditions are 

    ⊂span(Γ),
            Σ_υ,υ 
        = ΓΩΓ^T + Γ_oΩ_oΓ_o^T;

where (Γ,Γ_o) is a partitioned orthogonal matrix, the columns of 
Γ are a basis for ℰ_Σ_υ,υ(), and 
the dimensions of the positive definite matrices Ω and Ω_o are 
such that the matrix multiplications are defined. The quantities 
P_ℰΣ_υ,υP_ℰ and 
ΓΩΓ^T are often referred to as 'material information' 
in the envelope literature <cit.> since they represent the 
portion of variability that is necessary for the task of estimating 
υ. Similarly, 
Q_ℰΣ_υ,υQ_ℰ
and Γ_oΩ_oΓ_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 <cit.>. The same intuition applies 
to envelope methodology of <cit.>. 
In applications, there is a cost to estimate 
u = dim{ℰ_Σ_υ,υ()} 
and Γ. With the basis matrix Γ estimated, we can then assess the 
variance reduction of the envelope estimator through the parametric bootstrap. 

The 1D algorithm <cit.> estimates a basis matrix 
Γ for ℰ_Σ_υ,υ() at a 
user-supplied envelope dimension u. The estimate of Γ is obtained by 
providing Σ_υ,υ and 
υ̂υ̂^T as inputs into the 1D algorithm. The 
resulting estimator of Γ obtained from the 1D algorithm, 
Γ_u, gives a √n consistent estimator 
P_ℰ̂ of the projection onto the envelope subspace 
P_ℰ <cit.>. 



The 1D algorithm can be used to estimate u consistently 
<cit.>. 


Define = P_ℰ̂υ̂ to be the envelope 
estimator of υ where υ are parameters that link the 
estimation of expected Darwinian fitness to covariates of interest. We write 
τ = (γ^T, υ^T)^T where γ are aster model 
parameters not linking covariates to the estimation of Darwinian fitness. 


The envelope estimator of τ is given as 

    = ([ γ;   ]) = ([    I    0;    0 P_ℰ̂ ]) M^T Y = ^T Y,
           = M([    I    0;    0 P_ℰ̂ ]).

The model matrix  corresponds to the aster model that incorporates 
the envelope structure. 


 
The envelope estimator  is a maximum likelihood estimator of τ for 
the aster model with model matrix .



We have
l_env(β) = Y, β - c(β)
     = ^TY, β - c(β).
and
∇_β l(β) = ^TY - ∇_β c(β).
Setting ∇_β l(β) = 0 and solving for β yields 
∇_β c(β)|_β = β̂ = ^TY =. 



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 
τ to envelope estimators of β and then evaluate the 
log likelihood at the envelope estimator of β.
The randomness inherent in  is non-problematic in most applications. 
This is because the 1D algorithm provides a √n-consistent envelope 
estimator of υ 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(τ) using the aster model with model matrix . 
Note that the model matrix  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 β map to 
. Each of these distinct estimated values of β maps to the 
same estimate of β, 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 <cit.>, with g(·) replacing 
t(·). 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 

    ĝ_1D 
        = 1/B∑_b=1^B g{^(b)}.

The envelope estimator ĝ_1D implicitly behaves 
as a weighted average with the weights reflecting the likelihood of observing 
a particular estimated value of u, , and . 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 ĝ_1D will smooth out variability due 
to the estimating u with our chosen model selection criterion. 

Variability of ĝ_1D is estimated using the double 
bootstrap technique in <cit.>, 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 ĝ_1D:


  1. Fit the aster model to the data and obtain υ̂ 
  and Σ_υ,υ from the aster model fit. 

  2. Choose a model selection criterion. Compute the envelope 
  estimator of υ in the original sample, given as 
  = P_ℰ̂υ̂ where 
  P_ℰ̂ is obtained from the 1D algorithm and the chosen 
  model selection criterion. 
  
  
  
  

  3. Perform a parametric bootstrap by generating samples from the 
  distribution of the aster submodel evaluated at 
  = (γ̂^T,^T)^T. For b=1, …, B 
  of the procedure: 
  
  
  (3a) Compute τ̂^(b) and 
    Σ_υ,υ^(b) from the aster model fit to 
    the resampled data.
  
  (3b) Compute P_ℰ̂^(b) as done in Step 2. 
    
    
    
    
  
  (3c) Compute 
    
    ^(b) = {γ̂^(b)^T,^(b)^T}^T  
    and g{^(b)}.
  
  
  
  

  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. 
  
  At k = 1, …, K, for each b=1, …, B we:
  
  
  (5a) Generate data from the distribution of the aster submodel 
    evaluated at ^(b). 
  
  (5b) Perform Steps 3a through 3c with respect to the dataset 
    obtained in Step 5a to calculate both 
  
  
  ^(b)^(k) and g{^(b)^(k)}. 
  

  6. Compute both ĝ_1D and the standard deviation 
  in <cit.>.  



























































§ A DIRECT ENVELOPE ESTIMATOR USING REDUCING SUBSPACES
 


We propose a new way of constructing envelope estimators provided that the 
eigenvalues of Σ_υ,υ have multiplicity one. In this 
section, envelope estimators are constructed directly from the reducing 
subspaces of Σ_υ,υ. 

Uniqueness of the eigenvalues of Σ_υ, υ 
implies that its reducing subspaces are sums of its eigenspaces.  
Let  be a reducing subspace of Σ_υ,υ. 








Define Γ_ and P_ as the basis matrix for  
and the projection onto  respectively. 

Let  be the estimator of the reducing subspace , 
obtained by searching over the 1 dimensional eigenspaces of 
Σ_υ,υ.  

In applications, eigenvalues of Σ_υ,υ 
are almost always unique. Define Γ_ and 
P_ = 
    Γ_Γ^T_ 
as estimators of Γ_, and P_ respectively. The basis matrix 
Γ_ is constructed from the eigenvectors of 
Σ_υ,υ.
We now define the envelope estimator of υ constructed from 
reducing subspaces.


The envelope estimator of υ constructed from the reducing 
subspaces  is defined to be 
= P_υ̂. 




The reducing subspaces of Σ_υ,υ are 
√n consistent estimators of the reducing subspaces of 
Σ_υ,υ. Therefore Γ_, 
P_, and the corresponding estimator 
P_υ̂ are √n consistent 
estimators of Γ_, P_, and υ respectively. 
The envelope estimator of τ constructed by reducing subspaces is given by 

    = ([ γ;   ]) = ([  I  0;  0 P_ ]) M^T Y = ^T Y,
           = M([  I  0;  0 P_ ]).

The model matrix  corresponds to the aster model that incorporates 
the envelope structure obtained from the reducing subspaces of 
Σ_υ,υ. We have a similar result as 
Proposition <ref> for the aster model with model matrix .


 
The envelope estimator  is a maximum likelihood estimator of 
τ for the aster model with model matrix .





The proof of Proposition <ref> follows the same steps as the proof 
for Proposition <ref>. 
There is a close connection between envelope estimation using reducing 
subspaces and envelope estimation using the 1D algorithm. In the population, 
=.
The connection between both estimation methods exists in finite samples as 
seen in Theorem <ref>. In preparation, define orthogonal matrices 
O_u = (Γ_u,Γ_uo),  





O_ 
    = (Γ_,
      Γ_o), 



and O = O_O_u^T. 
The matrices Γ_, 
Γ_ o, Γ_u, and Γ_uo 
converges in probability to Γ_, Γ_ o, Γ_u, 
and Γ_uo respectively. Now 
Γ_^TΓ_u and Γ_ o^TΓ_uo are both 0-1 
valued rotation matrices and 
Γ_^TΓ_uo = Γ_ o^TΓ_u = 0. These facts 
imply that O converges in probability to a 0-1 valued rotation  
matrix. We will assume that Op→ I without loss 
of generality.



 
The basis matrix Γ_ is the output of the 1D 
algorithm with inputs 
M = OΣ_υ,υO^T 
and 
U = Oυ̂υ̂^TO^T 
at dimension u. 



Let M_2 = Σ_υ,υ and 
U_2 = υ̂υ̂^T.
Similar to the proof of <cit.>, let 

    Q_n(g)    = -n/2log(g^TOM_2O^Tg) 
        - n/2log[g^T{O(M_2 + U_2)O^T}^-1g]
        + n log(g^T g) 
       = -n/2log(g^TOM_2O^Tg) 
        - n/2log{g^TO(M_2 + U_2)^-1O^Tg}
        + n log(g^TOO^T g).

Now let v̂_k, k = 1, …, u be the kth column of 
Γ_ and define A ∈^p× p to be a matrix of 
0's with 1's occupying the first k diagonal entries. Then 

    Q_n(v̂_k) 
         = -n/2log(v̂_k^TOM_2O^Tv̂_k) 
        - n/2log{v̂_k^TO(M_2 
           + U_2)^-1O^Tv̂_k}
           + n log(v̂_k^TO^TOv̂_k) 
       = - n/2log(A_k O_uM_2O_u^TA_k) 
         - n/2log{A_kO_u(M_2 
           + U_2)^-1O_u^TA_k}
           + n log(A_kO_u^TO_u A_k) 
       = - n/2log(ĝ_uk^TM_2ĝ_uk) 
         - n/2log{ĝ_uk^T(M_2 
           + U_2)^-1ĝ_uk}
           + n log(ĝ_uk^Tĝ_uk)

where ĝ_uk is the kth column of Γ_u, 
the output of the 1D algorithm with M_2 
and U_2 as inputs. Therefore v̂_k is a maximizer 
of Q_n(g) and this completes the proof.



















































Theorem <ref> in combination with <cit.> 
allows for us to estimate  consistently. 

Thus the variability associated with the estimation of  and 
 decreases as n→∞. However in practical applications correct 
model selection cannot be guaranteed.


Therefore the envelope estimator of expected Darwinian fitness 
g() 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,  

    ĝ_env 
        = 1/B∑_b=1^B g{^(b)}.

The same model selection criteria is used to select the reducing subspace used 
to construct ^(b) at every iteration b = 1, …, B.
The second of level of this bootstrap procedure estimates the variability 
of (<ref>). The steps for this algorithm are provided in <cit.>. 
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,  is preferable to .
At any iteration of the 1D algorithm, minimizers of the objective function 
stated in <cit.> are pulled towards reducing subspaces of 
Σ_υ,υ. 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 <cit.>.













§ EXAMPLES









 §.§ 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 <ref>. 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 τ into (γ^T, υ^T)^T where γ∈^4 
are nuisance parameters and υ∈^5 are relevant to the estimation 
of expected Darwinian fitness. Here, υ∈^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 Σ_υ,υ 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  and the asymptotic variability of 
. The Bayesian information criterion is used to select  at every 
iteration of the first level of the bootstrap. 

The results are seen in Table <ref>. 
Table <ref> 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  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 
ĝ_MLE to ĝ_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 <cit.>. 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 <cit.>.


























































 

































 §.§ M. guttatus aster analysis

The yellow monkeyflower M. guttatus has been and is currently a well 
studied flower <cit.>. 
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 <cit.>. The lifecycle of the 
individual M. guttatus flowers, for our life history analysis, is 
depicted in panel B of Figure <ref>. 


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



<cit.> collected measurements on 2313 M. guttatus.
We fit a linear fitness landscape to this data. 



The parameters υ∈^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 Σ_υ,υ 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  where the 
Bayesian information criterion is used to select  at every iteration of 
the first level of the bootstrap. 

Table <ref> 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 (α = 0.05, unadjusted for multiple 
comparisons) the second row of Table <ref> from the 
fifth row of Table <ref>. 



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. 












§ SOFTWARE
 

This paper is accompanied by an R package 
<cit.>, which requires the two R packages for aster models: 
 <cit.> and  
<cit.>, and also a technical report <cit.> that 
reproduces the examples in this paper and shows how functions in the 
 package are used.









§ DISCUSSION









One could think to perform 
envelope methodology with respect to the regression coefficients β 
instead of τ. However β is not well-defined, one can shift β 
with an arbitrarily chosen offset vector without changing the value of the 
mean-value parameters τ and μ. 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 β.  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 β 
because this well-definedness issue is not a problem in those applications 
<cit.>. 
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 <cit.>. 





















The consequences of potential model selection errors served as the motivation 
for the implementation of the bootstrap procedure in <cit.>. 
In that article, inferences are only given for 
a canonical parameter vector in a multivariate linear regression model. 
We applied the <cit.> bootstrap procedures to alleviate 
possible model selection concerns. However, this particular choice of a 
bootstrap procedure is not without flaws. <cit.> 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, <cit.> 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 <cit.>.






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







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.


chicago



