Estimability for multiple covariance transformations
By default, phylosem takes as input a phylogenetic tree
and constructs a covariance under a Brownian motion evolutionary model.
However, it also allows any combination of the following
transformations:
- Ornstein-Uhlenbeck (OU);
- Pagel’s lambda;
- Pagel’s kappa;
This then results in eight possible covariance models formed from including or excluding each transformation. We first show that these transformations are all simultaneously estimable for a reference data set.
# Compare using Pagel's kappa
library(phylopath)
# Run phylosem
model = "
DD -> RS, p1
BM -> LS, p2
BM -> NL, p3
NL -> DD, p4
"
psem = phylosem( sem = model,
data = rhino[,c("BM","NL","DD","RS","LS")],
estimate_ou = TRUE,
estimate_lambda = TRUE,
estimate_kappa = TRUE,
tree = rhino_tree,
control = phylosem_control(
getJointPrecision = TRUE,
quiet = TRUE) )
#
V = psem$sdrep$cov.fixed
Rsub = cov2cor(V)[c('lnalpha','logitlambda','lnkappa'),c('lnalpha','logitlambda','lnkappa')]
knitr::kable(c("minimum_eigenvalue"=min(eigen(psem$sdrep$jointPrecision)$values),
"maximum_eigenvalue"=max(eigen(psem$sdrep$jointPrecision)$values)), digits=3)| x | |
|---|---|
| minimum_eigenvalue | 0.035 |
| maximum_eigenvalue | 43439.097 |
| lnalpha | logitlambda | lnkappa | |
|---|---|---|---|
| lnalpha | 1.000 | 0.741 | 0.422 |
| logitlambda | 0.741 | 1.000 | 0.183 |
| lnkappa | 0.422 | 0.183 | 1.000 |
We see that the minimum eigenvalue for the precision matrix (which corresponds to the maximum eigenvalue of the covariance matrix) is greater than zero, such that the precision is invertible and the Hessian matrix is full rank. From this we can see that all parameters are estimable when applying all three transformations. We also see that transformation parameters are not perfectly correlated, although the logit-lambda and ln-alpha parameters have a correlation of approximately 0.75.
We can also extract the standard errors, calculated as the
square-root of the diagonal elements of the covariance matrix. We
extract the estimated values for lnalpha,
logitlambda, and lnkappa and plot plot the 95%
confidence interval as their estimates plus or minus 1.96 times their
standard errors. These values are all unbounded by construction.
However, we also convert these estimates, lower, and upper
confidence-interval bounds to the conventional parameters
alpha, lambda, and kappa, and
plot those as well. These derived quantities are bounded as expected for
Ornstein-Uhlenbeck and Pagel parameters:
library(ggplot2)
# Compile estimates and SEs
pdat = data.frame( "Estimate" = psem$sdrep$par.fixed[c('lnalpha','logitlambda','lnkappa')],
"StdErr" = sqrt(diag(V)[c('lnalpha','logitlambda','lnkappa')]) )
pdat = cbind( pdat, "Param" = rownames(pdat))
#
pdat$lower = pdat$Estimate - 1.96*pdat$StdErr
pdat$upper = pdat$Estimate + 1.96*pdat$StdErr
pdat$type = "estimated"
# Transform from log / logit-space to natural space
pdat2 = pdat
pdat2$Param = c("alpha", "lambda", "kappa")
pdat2['lnalpha',c("Estimate","lower","upper")] = exp(pdat2['lnalpha',c("Estimate","lower","upper")])
pdat2['lnkappa',c("Estimate","lower","upper")] = exp(pdat2['lnkappa',c("Estimate","lower","upper")])
pdat2['logitlambda',c("Estimate","lower","upper")] = plogis(as.numeric(pdat2['logitlambda',c("Estimate","lower","upper")]))
pdat2$type = "derived"
# Plot
ggplot( rbind(pdat,pdat2), aes( x=Param, y = Estimate, ymin = lower, ymax = upper)) +
geom_pointrange(position = position_dodge(width = 0.6)) +
theme_classic() +
facet_grid( rows=vars(type), scales="free" )
These confidence intervals suggest that parameter lnkappa overlaps with 0.0, suggesting that Pagel’s kappa overlaps with 1.0 and could be instead turned off and thereby fixed at that default value.
