Validating the methods of Leibold et al. (2010) and Peres-Neto et al. (2012)

The functions phylo_vector(), trait_partition() and node_test() are translations of MATLAB code from the supplements of these articles. We checked them in three ways. Each help page lists the differences between the text of the article and the code.

Comparison with the original code

We ran the MATLAB code in Octave on the numerical example of Appendix A of Peres-Neto et al. (2012) and on simulated data. The supplement lacks the functions eigen and WeightedMeanVar, which we wrote. Weighted coordinates, correlations, Moran’s eigenvector maps and the raw and adjusted fractions agreed to within 1e-13. Correlations also agreed with ade4::fourthcorner(). The package does not include the maps, as those of adespatial::dbmem() are the same as the ones of the supplement.

Type I error and power

The script data-raw/validation_gradient.R follows the simulations of the article, with 1000 data sets and 999 permutations (100 sites and 50 species, one environmental variable). The article stacks the two null distributions and uses alpha / 2, shown here as “stacked”. Rates near 0.05 are expected for a true null hypothesis, with a standard error near 0.007.

Trait and environment (Table 4 of the article), with the correlation inside phylo_vector(). The package has no function for this test, which is available in ade4::fourthcorner():

Case	Tolerance	Row, article	Row, here	Column, article	Column, here	Stacked, article	Stacked, here	Larger p, here
Random trait	1	0.836	0.812	0.056	0.050	0.030	0.048	0.050
Random trait	20	0.916	0.923	0.052	0.042	0.046	0.042	0.042
Random trait	60	0.984	0.985	0.058	0.048	0.048	0.047	0.048
Random environment	1	0.052	0.053	0.678	0.695	0.052	0.052	0.053
Random environment	20	0.056	0.052	0.832	0.866	0.038	0.051	0.052
Random environment	60	0.058	0.046	0.954	0.962	0.036	0.046	0.046

Global test of phylo_vector(), with one Brownian trait:

Case	Tolerance	Stacked, article	Stacked, here	Larger p, here
Random environment (Table 1)	1	0.028	0.039	0.051
Random phylogeny (Table 2)	1	0.042	0.040	0.040
Power (Table 3)	1	0.880	0.875	0.877
Power (Table 3)	20	0.948	0.934	0.936
Power (Table 3)	60	0.948	0.940	0.941

With tolerances of 20 and 60, which the article does not report for Tables 1 and 2, the larger p-value gave rates from 0.037 to 0.055. The test of the final vector after a significant global test kept these rates, because the global test comes first. A test of that vector alone is not valid, as it was chosen to fit the predictor.

Weights of the eigenvectors

The correlation between the vector and the trait that generated the data depends on how the eigenvectors are weighted. The text uses correlations and the code their absolute values.

Tolerance	PEV with absolute values, here	PEV in Table 3	PEV with signs, here	PEV-B, here	PEV-B in Table 3
1	0.480	0.514	0.685	0.698	0.676
20	0.455	0.526	0.722	0.707	0.703
60	0.473	0.508	0.740	0.715	0.688

Absolute values give about what the article reports for PEV, which suggests that its Table 3 used them (this is our interpretation). With signs, PEV matches PEV-B, so the gain of the Brownian vectors reported in the article may come from this difference.

Partition

Table 5 uses the logistic generator of Appendix C. The appendix leaves out the intercepts and the scale of the slopes, so we took them from the simulation code of Ives & Helmus (2011): intercepts from N(0, 0.5), slopes with a standard deviation of 1 among species, and presences from a uniform draw. Neither the article nor the appendix says how the two vectors were obtained. The fractions of Table 5 agree with the vectors that generated the data, and its rejection rates for E and S agree with the global test of phylo_vector(), so the table below combines both (data-raw/validation_table5.R, 1000 data sets, 999 permutations). The values of the article are in parentheses:

Scenario	[E, S]	[a]	[b]	[c]	Rejection of E	Rejection of S
I	0.490 (0.478)	0.223 (0.216)	0.043 (0.069)	0.223 (0.193)	0.706 (0.687)	0.710 (0.679)
II	0.473 (0.589)	0.008 (0.061)	0.458 (0.488)	0.007 (0.036)	0.700 (0.750)	0.705 (0.751)
III	0.383 (0.349)	0.239 (0.234)	0.030 (0.036)	0.114 (0.079)	0.691 (0.701)	0.592 (0.418)
IV	0.047 (0.032)	0.023 (0.015)	0.000 (0.001)	0.024 (0.016)	0.043 (0.051)	0.047 (0.061)
V	0.524 (0.457)	0.219 (0.202)	0.049 (0.054)	0.255 (0.201)	0.683 (0.597)	0.724 (0.637)
VI	0.518 (0.440)	0.217 (0.201)	0.051 (0.063)	0.251 (0.176)	0.707 (0.632)	0.697 (0.627)
VII	0.549 (0.396)	0.255 (0.194)	0.041 (0.000)	0.253 (0.202)	0.700 (0.562)	0.708 (0.597)

The pattern among scenarios is the same, and scenarios I, III and IV are the closest. Three choices are ours. In scenario III, the N(0, 1) error of Appendix C has no effect on a gradient from -15 to 15, so we added it to y with unit variance, and only in the analysis. In scenarios V and VI, x is in random order along the transect, as otherwise it is correlated with space. In scenario II we kept the text, and the gradients seem more correlated than those of the article, as fractions a and c are smaller.

Three parts remain different. The share of the compositional component is near 0.62 in all scenarios, whereas the article gives 0.39 to 0.94, so that value must depend on something that the appendix does not describe. Scenario VII has a larger b than the zero of the article, probably because the maps here are orthogonal only without site weights. The rejection of [E, S] is 0.71 to 0.91, against 0.77 to 0.84, here as the rejection of either gradient, which gives 0.085 in scenario IV against 0.055 in the article. Rejections of a and c, with vectors estimated only after a significant global test, were 0.66 to 0.71 and 0.48 to 0.72 outside scenarios II and IV, against 0.56 to 0.70 and 0.30 to 0.60 in the article; in scenario II they were zero, against 0.29 and 0.25.

Following the demonstration script of the authors instead, with a vector fixed in advance and predictors without relation to the data, rejection rates of the larger p-value were 0.023 to 0.028, because both null hypotheses are true. With environment and space acting, the global test for environment passed in 37% of the data sets. In those, the rejection rate was 0.99 for fraction a and 0.03 for fraction c, as expected for a vector estimated for environment alone.

The daphniids of Leibold et al. (2010)

The data set daphniids rebuilds the data of the article from the public survey. With the larger p-value, no predictor is significant at any node, as species permutations have little power with 13 species. The article permutes rows and columns together, which node_test() gives as p_joint. Rates of rejection per node at the 5% level (one Yule tree of 24 species, 60 sites, 300 data sets, 99 permutations, data-raw/validation_node.R):

What acts	Sites	Species	Larger p	Joint
Neither	0.053	0.049	0.021	0.050
Environment	0.777	0.049	0.049	0.764
Phylogeny	0.046	0.050	0.021	0.050
Both	0.803	0.073	0.073	0.795

The joint test is valid only when neither acts. With the environment alone it behaves as the permutation of sites, so it does not separate an effect of the environment from an association between environment and phylogeny.

In the daphniids, the joint test with 999 permutations found fewer predictors than sites alone at the 0.05 level (PHEQ, AV_DEP, DOC, TN_TP and PTL at four nodes, with none left after the false discovery rate). It therefore did not recover Figure 3a better. With sites alone, at a level of 0.05 without correction and 999 permutations, environmental results compare with Figure 3a as follows:

Node	Figure 3a	Here
retrocurva, parvula, catawba	COND, NTL	COND, NTL
Same plus pulicaria, pulex	AV_DEP, PHEQ	AV_DEP, PHEQ
longispina, galeata, rosea	AV_DEP, ELEV	AV_DEP, ELEV
All species	TN:TP, AV_DEP, ELEV, DOC, COLOR	COLOR, DOC, AV_DEP
All Daphnia except laevis, dubia	AV_DEP, PHEQ, ELEV, DOC	AV_DEP, PHEQ, ELEV, DOC, NTL, TN_TP
Clade with ambigua	ELEV	ELEV, TN_TP
retrocurva, parvula	PHEQ, AV_DEP	AV_DEP
pulicaria, pulex	ELEV, NTL	none
galeata, rosea	COLOR, DOC	none
laevis, dubia	COND	NTL, PHEQ, PTL, TN_TP

Nodes with rare species differ most: D. laevis and D. rosea occur in one and two lakes of the rebuilt set. Biogeographic results differ more (no match at most nodes), which we expected, because the article does not give the area of each lake and we approximated rivers by longitude.

References

Leibold, M. A., Economo, E. P., & Peres-Neto, P. (2010). Metacommunity phylogenetics: separating the roles of environmental filters and historical biogeography. Ecology Letters, 13, 1290–1299. https://doi.org/10.1111/j.1461-0248.2010.01523.x

Peres-Neto, P. R., Leibold, M. A., & Dray, S. (2012). Assessing the effects of spatial contingency and environmental filtering on metacommunity phylogenetics. Ecology, 93, S14–S30. https://doi.org/10.1890/11-0494.1