Does Pagel's lambda create a positive semi-definite covariance matrix?

In phylogenetic regression, which is a form of linear regression involving multiple species observations, Pagel’s lambda indicates the extent to which covariance between observations is necessary to maximize the data’s likelihood. This concept is known as the phylogenetic signal in this context and has links to linear mixed models, as we will explore. However, since it is a parameter that affects only the covariances (i.e., the off-diagonal values of the covariance matrix) through multiplication, it was initially unclear to me why the resulting new covariance matrix would still be a positive semi-definite (PSD) matrix. Being PSD confers great properties, such as invertibility, real eigenvalues, and spectral decomposition. With this in mind, I would like to present a brief proof of this property and follow up with some remarks on how to estimate this parameter.

Proposition. Let $\textbf{C}_{\lambda} \in \mathbb{R}^{n \times n}$ be the transformed covariance matrix $\Omega \in \mathbb{R}^{n \times n}$ with off-diagonal elements multiplied by the scalar $\lambda \in [0,1]$ such that:

Then, $\textbf{C}_{\lambda}$ is positive semi-definite matrix.

Proof. We can re-write $\mathbf{C}_{\lambda}$ as the following:

Now, $\Omega$ can be defined as the covariance matrix of some random variable $\mathbf{x} \in \mathbb{R}^{n \times 1}$, with probability density $p(\mathbf{x})$, and expectation $\mu$:

Finally, let $\mathbf{u} \in \mathbb{R}^{n \times 1}$ be any vector different than 0 such that:

Since $\mathbf{u}^\top\textbf{C}_{\lambda} \mathbf{u} \geq 0$ for any $\mathbf{u} \neq 0$, then \(\textbf{C}_{\lambda}\) is positive semidefinite matrix. $\,\,\square$

Remark 1. If we considering the following linear regression model:

where $\textbf{X} \in \mathbb{R}^{n \times p}$ is matrix with observations, $\textbf{Y} \in \mathbb{R}^{n \times 1}$ is the response variable, $\beta$ is the vector with coefficients, and $\mathbf{C}_{\lambda}$ is the covariance matrix, then likelihood of the data can be expressed as:

Here is where expressions such as $\mathbf{C}_{\lambda} = \lambda\Omega + (1 - \lambda)\mathbf{W}$ start singing and reveal to us something profound and interesting about the equation above. It implies that if we allow $\lambda$ to be $0$, the variability between species observations does not influence the likelihood, and only the variance of species is important. Conversely, if $\lambda$ is $1$, then the variances of between species are significant, but the within-species variability is not as crucial. This point is where the connection with linear mixed models becomes more apparent, as the random intercepts and slopes in mixed models also strive to find a balance between variations between individuals versus variations within individuals, which usually involve repeated measurements. Lynch (1991) proposed the phylogenetic version of mixed models, called the Phylogenetic Linear Mixed Model (PLMM). Historically, however, it was less popular than Felsenstein’s (1985) Phylogenetic Independent Contrast (PIC), despite the fact that PIC can be seen as a special case of PLMM when there is not much variation within individuals (Garamszegi 2014), or, if we attempt to make a connection with $\lambda$, when $\lambda = 1$.

Obtaining a closed form for $\lambda$ is difficult

Since we talked about maximizing, let’s first define our objective function. We can obtain it from the likelihood equation:

From above we can define our objective function $J(\beta,s,\mathbf{C}_{\lambda})$ as:

Recall that $\varepsilon$ is the error from our initial linear regression model, as defined in Remark 1, and we replaced $\sigma^2$ simply with $s$. We will use either the first or second line of the objective function definition interchangeably, depending on the variable to be differentiated.

Differentiating $J(\beta,s,\mathbf{C}_{\lambda})$ with respect to $s$ and setting it to zero we have:

Now, differentiating $J(\beta,s,\mathbf{C}_{\lambda})$ with respect to $\beta$ and setting it to zero we should have the weighted least square solution:

Finally, differentiating $J(\beta,s,\mathbf{C}_{\lambda})$ with respect to $\lambda$ and setting it to zero:

We used the multiplication rule property for the right hand side differentiation. From here we can start seeing troubles. How can we get the $\lambda$ out of the trace or of the matrix inverse? For the sake of completeness let’s finish the differentitation:

For the right-hand side differentiation, we used this property: $\partial (\mathbf{A}^{-1}) = - \mathbf{A}^{-1} (\partial \mathbf{A})\mathbf{A}^{-1}$. Considering $\mathbf{C}_{\lambda}^{-1} = (\lambda\Omega + (1 - \lambda)\mathbf{W})^{-1}$, it becomes apparent that isolating $\lambda$ from the expression above is challenging. Furthermore, we have the constraint that $\lambda \geq 0$ and $\lambda \leq 1$. While multiple inequality constraints can be addressed using Lagrange multipliers followed by the Karush–Kuhn–Tucker (KKT) conditions, which look for all the active constraints that fulfill the KKT conditions (e.g., positive multipliers, feasible inequalities), this usually requires isolating $\lambda$, a task that is not trivial.

But we can use numerical optimization

From above derivations, we can propose the following optimization problem:

where $\lambda$ is the unique decision variable.

Here is a set of functions to optimize the above problem:

from scipy.optimize import minimize
import numpy as np

def objective_function(X, y, Sigma, lambda_val):
    W = np.diag( np.diag(Sigma) )
    C = lambda_val * Sigma + (1 - lambda_val) * W

    C_inv = np.linalg.inv(C)
    beta  = np.linalg.inv(X.T @ C_inv @ X) @ (X.T @ C_inv @ y)
    e = y - X @ beta 
    s = (1/n)*e.T @ C_inv @ e

    objective = n*np.log(s) + np.log(np.linalg.det(C)) + (1/s)*e.T@ C_inv @e
    return objective.astype(np.float64)[0][0]

def optimize_problem(X, y, Sigma):
    n, p = X.shape
    X = np.hstack((np.ones((n, 1)), X))

    # Define the bounds for lambda
    bounds = [(0, 1)]

    # Initial guess for lambda 
    initial_lambda = np.random.uniform(0, 1)

    # Define the optimization problem
    fnc = lambda lambda_val: objective_function(X, y, Sigma, lambda_val)
    result = minimize(fnc, x0 = initial_lambda, bounds=bounds, method='SLSQP')

    return result

Testing in two trees:

np.random.seed(123)
n = 5
# tree 1
C1 = np.array([[7,0,0,0,0],
               [0,7,6,4,4],
               [0,6,7,4,4],
               [0,4,4,7,6],
               [0,4,4,6,7]])
X1 = np.random.multivariate_normal(mean = [0]*n, cov = C1, size = 2).T
y1 = np.random.multivariate_normal(mean = [0]*n, cov = C1, size = 1).T
print(optimize_problem(X1, y1, C1).x[0]) # 1.0

# tree 2
C2 = np.array([[7, 0, 0, 0, 0],
               [0, 7, 1,.5,.5],
               [0, 1, 7,.5,.5],
               [0,.5,.5, 7, 1],
               [0,.5,.5, 1, 7]])
X2 = np.random.multivariate_normal(mean = [0]*n, cov = C2, size = 2).T
y2 = np.random.multivariate_normal(mean = [0]*n, cov = C2, size = 1).T
print(optimize_problem(X2, y2, C2).x[0]) # 0.0

References

• Garamszegi, L. Z. (Ed.). (2014). Modern phylogenetic comparative methods and their application in evolutionary biology: concepts and practice. Springer.

• Housworth, E. A., Martins, E. P., & Lynch, M. (2004). The phylogenetic mixed model. The American Naturalist, 163(1), 84-96.

• Lynch, M. (1991). Methods for the analysis of comparative data in evolutionary biology. Evolution, 45(5), 1065-1080.

• Felsenstein, J. (1985). Phylogenies and the comparative method. The American Naturalist, 125(1), 1-15.