Testing species trait values

[teixeirak]

This relates to issue #74 and #92.

From R1: "The authors might like to consider including Table 4 as a supplement, because as it stands it looks like an exploratory analysis. If the ‘single-variable approach’ is not just exploratory, state it clearly and why it is needed besides the full models with several parameters."

Thinking more about this, we need to balance two considerations: 1- Guard against over-fitting, which is a danger with just 12 species and 5 traits of interest. 2- test the explanatory power of each trait.

I'm going to open this issue now as a placeholder and come back with more discussion.

[mcgregorian1]

When I first read this comment I was a little confused to the meaning behind it - isn't the point of Table 4 that it's exploratory but as a vessel for getting to Table 5 and the rest of the analysis? I could see it becoming a supplement in order to give more room in the paper, i.e. "We first tested our hypotheses and determined which variables were flagged as candidate traits for inclusion in the full models for each scenario via instances where dAICc > 1 [see Figure S4]."

For testing the explanatory power of each trait - this was the original reason for our individual hypotheses as shown in Table 4, correct? And thus the dAICc is supposed to represent that explanatory power. As we've seen, though, this power varies a bit going from Rt to arima_ratio.

[teixeirak]

To guard against over-fitting, I think it's important to test each trait independently, and then use that info to inform which variables are allowed as candidates in the best full model, limiting to 2 or at most 3.

The independent tests of each variable were previously covered in Table 4, where the only other fixed variable in the model was height. This has the disadvantage of not accounting for other variables-- particularly TWI, which consistently comes out significant.

What if we did the following? 1- Include all individual-level variables (height, TWI, crown position) as candidates in all models. 2- Iteratively run through models of all species-level variables (the traits). Results of this analysis will be summarized in Table 1, details given in revised version of Table 4 (which may go in SI) and/or summarized in a figure (showing partial effect of variable in model with confidence bounds), similar in style to the example below. 3- For each drought year (or all combined), run combined model with top 2-3 trait variables (only those that come out significant on their own). The model with lowest AIC would be the top overall model for that drought. These top models could also be candidate for visualization, again similar in style to the example below.

Thus, the most fundamental change would be that we test each species trait in the context of a full model for individual-level variables (null in Table 4 expands). This meets the two objectives above.

[teixeirak]

When I first read this comment I was a little confused to the meaning behind it - isn't the point of Table 4 that it's exploratory but as a vessel for getting to Table 5 and the rest of the analysis? I could see it becoming a supplement in order to give more room in the paper, i.e. "We first tested our hypotheses and determined which variables were flagged as candidate traits for inclusion in the full models for each scenario via instances where dAICc > 1 [see Figure S4]."

For testing the explanatory power of each trait - this was the original reason for our individual hypotheses as shown in Table 4, correct? And thus the dAICc is supposed to represent that explanatory power. As we've seen, though, this power varies a bit going from Rt to arima_ratio.

Hopefully the comment I just sent addresses all this.

[mcgregorian1]

It does! I'll get on this

[mcgregorian1]

@teixeirak quick question - do we want to include TWI*height interaction as part of the base model? If we do then technically all our base models (before adding in species-level traits) would be:

• arima_ratio ~ TWI*height + crown.position as opposed to
• arima_ratio ~ TWI + height + crown.position

My hunch is we should not include the interaction since TWI alone is what has been consistently significant.

[teixeirak]

yes, I think we can include it, but expect it probably won't come out significant most of the time.

[mcgregorian1]

When we do this, we get rp and PLA as significant in 1966 and 1977/1999 respectively. Every other scenario for any trait, adding in the hydraulic traits makes things worse.

I will finish doing the workflow you suggested and make new tables.

[mcgregorian1]

yes, I think we can include it, but expect it probably won't come out significant most of the time.

Because it's a base model, it won't come out as significant at all, really, because I'm not testing the interaction against the individual constituents anyways.

However, I just tried re-running without the interaction and I can confirm it barely affects the results I got for rp and PLA. I'll leave as is.

[teixeirak]

When we do this, we get rp and PLA as significant in 1966 and 1977/1999 respectively. Every other scenario for any trait, adding in the hydraulic traits makes things worse.

Is this for Rt or ARIMA_ratio?

[mcgregorian1]

This is for arima_ratio

[teixeirak]

Okay. The one we want to present depends on the outcome of issue #92 .

[mcgregorian1]

I've redone Tables 4 and 5 using this new method of having the base model include the individual traits, and only testing the species-level traits against them. These can be found with the _reform suffix to the files (e.g. we can compare Rt and arimaratio using the tables-figures_original_reform.pdf and the figures-tables_arimaratio_reform.pdf). I've also updated the corresponding csvs.

[teixeirak]

Thank you; this is great!

It's notable that in both scenarios PLA and TLP are both consistently negative, as predicted, whereas the other traits have coefficients going in different directions for different droughts. This makes me think these two are more fundamental effects that are consistent across droughts. Of course, TLP doesn't come out significant.

It would be good to also see a version of top models including only the traits with consistent coefficients across droughts when tested themselves. From a quick look, that seems to be TLP and PLA, but needs double checking.

[mcgregorian1]

As from #92, I have moved the arima tables to the archive folder. With regard to this issue, in the main tables_figures folder I've kept the original Rt Tables 4-5 (using only height in the base model) plus the _reform.csv tables (using the height*TWI + crown.position as the base model) in case we want to compare more later.

It's notable that in both scenarios PLA and TLP are both consistently negative, as predicted, whereas the other traits have coefficients going in different directions for different droughts. This makes me think these two are more fundamental effects that are consistent across droughts. Of course, TLP doesn't come out significant.

With the original Rt, we see RP, TLP, and PLA as candidate traits for inclusion in the best models (see Table 4).

It would be good to also see a version of top models including only the traits with consistent coefficients across droughts when tested themselves. From a quick look, that seems to be TLP and PLA, but needs double checking.

Correct, this is TLP and PLA, as RP has an anomaly for the 1977 drought.

Thus, here are the top models when you only consider TLP and PLA.

And here are the top models when you consider TLP, PLA, plus RP.

[mcgregorian1]

I can get each version's full coefficients and R^2s if you want to see those.

[teixeirak]

As from #92, I have moved the arima tables to the archive folder. With regard to this issue, in the main tables_figures folder I've kept the original Rt Tables 4-5 (using only height in the base model) plus the _reform.csv tables (using the height*TWI + crown.position as the base model) in case we want to compare more later.

Actually, the two I think we should keep/ present in the pub are Rt_reform (main doc) and arimaratio_reform (SI).

[mcgregorian1]

Here are the coefficients and R^2s for the two scenarios:

All the top candidate vars (TLP, PLA, RP):

Only consistent coefficients across droughts when tested separately (only PLA and TLP):

[mcgregorian1]

Actually, the two I think we should keep/ present in the pub are Rt_reform (main doc) and arimaratio_reform (SI).

Gotcha, I'll put those back.

[teixeirak]

Thus, here are the top models when you only consider TLP and PLA.

This actually seems to be the most parsimonious. I've never liked RP as a predictor, as its just 2 (uneven) categories and goes in opposite directions. So, the criteria could be:

"Candidate variables were selected, based on the single-variable tests, as those whose addition to a corresponding null model improved fit (at $\Delta$AICc $\ge$ 1.0) in at least one drought year and whose direction of response was consistent across all drought years."

Do you agree?

[teixeirak]

I think this approach, which leads to the the inclusion of TLP, is justified by the fact that TLP comes out in the top model for 2 of the 4 scenarios.

[mcgregorian1]

"Candidate variables were selected, based on the single-variable tests, as those whose addition to a corresponding null model improved fit (at $\Delta$AICc $\ge$ 1.0) in at least one drought year and whose direction of response was consistent across all drought years."

Do you agree?

That sounds good. I wonder if we'd get more pushback from that? I see no reason we would but I'm not too knowledgeable on statistical rights.

[mcgregorian1]

Do you want the arimaratio_reform tables to reflect this as well? So, we'll have two versions of the same workflow with one using arimaratio and the other using Rt?

[teixeirak]

I don’t see why we’d get pushback. The important thing is to be clear and transparent about what we did. We could always present alternate results if reviewers want to see them, but I do think this criteria makes sense and leads to the most parsimonious results.

[teixeirak]

Do you want the arimaratio_reform tables to reflect this as well? So, we'll have two versions of the same workflow with one using arimaratio and the other using Rt?

Yes.

[mcgregorian1]

@teixeirak I remembered just now that one of the reasons we had not included position_all in the models with height earlier is because of the high correlation (see below). I spoke with Monika about this as she did this same thing for her dissertation, and she said for her, anything above a 0.7 was grounds for removal. I think Valentine would know.

It seems it would make sense to drop position instead of height since the measurement of height is more accurate.

I will wait to redo the tables from your last comment until we decide what to do here, and based on what we decide for #94.

[teixeirak]

This high correlation is definitely an issue, and limits our interpretation, but I'd rather treat it in the discussion than remove canopy position from the analysis. We currently say this in the discussion:

"However, the collinearity between the two variables (Fig. 2d) makes it impossible to confidently partition causality. Taller trees are more likely to be in dominant canopy positions (Fig. 2d) and, largely as a consequence of their position relative to others, face different microenvironments (Fig. 2a-b). "

[mcgregorian1]

Do you want the arimaratio_reform tables to reflect this as well? So, we'll have two versions of the same workflow with one using arimaratio and the other using Rt?

Yes.

I am currently doing this but I think it's worth pointing out that the candidate traits from arimaratio are only RP and PLA:

[teixeirak]

That’s fine.

[teixeirak]

@mcgregorian1 , I'm confused here-- about two things: 1- You mentioned that you were implementing the same workflow for Rt_ARIMA, but this table does not have rp and does have TLP. 2- By our criteria, rp should not be included in the Rt_ARIMA full models because its coefficient is inconsistent across years. This would leave just PLA.

That said, I think we can just stick with what you currently have. Rather than making the workflow parallel, we can test the same candidate variables in models for both variables. This is also informative, and its noteworthy that TLP comes out in two of the top models (although I'm not sure I understand that... table S4 indicates that TLP makes all models worse, whereas table S5 indicates it improves models for two drought years. I think its exactly the same set of variables, right?).

[mcgregorian1]

@mcgregorian1 , I'm confused here-- about two things: 1- You mentioned that you were implementing the same workflow for Rt_ARIMA, but this table does not have rp and does have TLP.

Yes, both Table 5 and S5 have TLP and PLA because for both arimaratio and Rt I manually constrained the candidate variables to be just those two.

2- By our criteria, rp should not be included in the Rt_ARIMA full models because its coefficient is inconsistent across years. This would leave just PLA.

That said, I think we can just stick with what you currently have. Rather than making the workflow parallel, we can test the same candidate variables in models for both variables.

Yes, I made sure to test arimaratio and Rt against the same model format (each is tested for PLA and TLP).

This is also informative, and its noteworthy that TLP comes out in two of the top models (although I'm not sure I understand that... table S4 indicates that TLP makes all models worse, whereas table S5 indicates it improves models for two drought years. I think its exactly the same set of variables, right?).

It is, so I guess including PLA with TLP makes TLP improve the models? Table S4 results was just TLP by itself (with TWI*height and crown position), and same thing for PLA. Actually, I'm looking at S4 and S5 right now for TLP and it seems to have very similar coefficients and the same sign. Am I looking at the wrong thing?

[teixeirak]

Let's stick with this: that is, Rt as the main variable constrains which traits we include as candidates in the model (it's similar either way--difference would be whether we include TLP or not).

[teixeirak]

It is, so I guess including PLA with TLP makes TLP improve the models? Table S4 results was just TLP by itself (with TWI*height and crown position), and same thing for PLA. Actually, I'm looking at S4 and S5 right now for TLP and it seems to have very similar coefficients and the same sign. Am I looking at the wrong thing?

I'm looking at the same thing (TLP coefficients in tables S4 and S5). What I don't understand is why the top model has TLP, when table S4 indicates that adding it in makes the model worse. Similarly for PLA across all droughts. 1966 makes sense. Does your full model allow there to be no species traits in the top model?

[mcgregorian1]

Ahh I see what you mean. Sorry, I was looking mainly at the coefficients not the dAIC.

Wow that's bad on my part. I changed the code to include TLP and PLA but not include instances of both or none. I've changed that to be correct. Here's the result for arimaratio:

And here's the result for Rt:

The following tables have been updated.

• tested_traits_best_coeff_reform_arimaratio
• tested_traits_best_reform_arimaratio
• top_models_dAIC_reform_arimaratio

Also the same for Rt:

• tested_traits_best_coeff_reform
• tested_traits_best_reform
• top_models_dAIC_reform

[teixeirak]

Wow, glad we caught that! New results look good, and seem to present a pretty consistent story that both come out in top models sometimes, but its a little stochastic which, if any, show up.