Revisit special handling of `seedling_mortality` in run.meta.analysis

[infotroph]

Current lines 103-104 of run.meta.analysis.R:

prior_variances[names(trait_average), ] <- 0.001 * trait_average ^ 2
  prior_variances["seedling_mortality", 1] <- 1

@infotroph in a comment on #3278:

Would it be possible to add a few comments here explaining (1) where the constants 0.001 and 1 come from, and (2) why seedling mortality gets special handling and whether that's a hack or a principled design choice?

On scrolling down further I see this is pre-existing code moved up from run.meta.analysis.pft, so maybe these explanations are lost to time...

@ashiklom replied:

Yeah, this is pre-existing code. The oldest git blame gets me is here:

https://github.com/PecanProject/pecan/blame/bed27a2b37619c2297c3c66becf97b8c16233494/modules/meta.analysis/R/run.meta.analysis.R

...but I'm almost certain this logic predates that, which looks more like refactoring than development.

The more correct thing to do here is probably (1) remove the seedling_mortality hard-coded value altogether and replace it with a global default prior of some kind in the database; (2) pass the 0.001 as an (optional) argument called something like default_prior_variance_fraction with 0.001 as the default value.

[aditi25srivastava]

@infotroph Heyy!!!!! I would love to work on this issue. Could you please assign it to me ?

[infotroph]

@aditi25srivastava You're very welcome to take a crack at it, but be aware that the task is incompletely specified at the moment, so the bulk of the work will be in deciding how exactly to implement it in the context of PEcAn's statistical framework.

[aditi25srivastava]

Thank you @infotroph for the clarification, that makes sense
i understand that the task is still loosely specified and depends on the final React/framework decisions. I won't start implementing anything until the approach is clearly defined.. I'm keeping an eye on the issue and would be happy to help once the requirements are finalised. Please feel free to tag or assign me when it's ready to be picked up

[mdietze]

Just took a quick look and I'll be honest that I'm more concerned about the (undocumented) default logic around prior_variances than the seedling_mortality special case (which I think can probably be dropped).

      prior.variances[names(trait.average),] = 0.001*trait.average^2 
      prior.variances["seedling_mortality",1] = 1.0
      taupriors <- list(tauA = 0.01,
                        tauB = apply(prior.variances, 1, function(x) min(0.01, x)))

A few recommendations:

1. modify the code so that users can specify priors on tau, with default values generated if user provides ones are not includes
2. better document the default prior so the function man page explains the default and the code is better documented inline.
3. consider updating the default

A few things to know about gamma priors. The first parameter is n0/2 where n0 is the prior sample size (i.e. the current default is 0.02 data points, which is very uninformative). The second parameter is SS/2 where SS is the prior sum of squares, which should be SS=n0*prior.variance. Thus, as currently used, the prior.variance isn't a prior.variance but a prior SS/2. Next, if you assume the prior SD is proportional to the prior mean, a*mu, then the prior.variance will be a^2*mu^2, which is probably where the trait.average^2 is coming from. When you divide by tauA I think this gives a=sqrt(0.1) = 0.31 (i.e. I think the prior assumes a 31% uncertainty around the mean, but double check my math). I'm willing to bet there's more details around all of this in the LeBauer et al 2013 Ecological Applications paper. Finally, one reason to assume the prior SD varies with the mean is that tauB can be very sensitive to units (this I 100% do remember from when @dlebauer originally wrote this code). A common uninformative prior on tau might be something like gamma(0.01, 0.01) but if the units of your data are something like 10^-6 then a SS of 0.01 is massively informative and no practical amount of trait data will ever overwhelm it.

[shreyannandanwar]

@mdietze, A reasonable incremental approach is could be -
first remove undocumented special cases then turn hard-coded constants into named arguments with defaults, and clearly document the statistical meaning of the existing priors without changing the underlying math.

This will improves transparency and user control immediately, preserves backward compatibility, and creates a clean path toward a more principled redesign of the prior structure later.
May I create a PR and work on this issue

[mdietze]

@shreyannandanwar yes, that'd be great. Thanks!

In terms of the underlying math and what we expose to the user, the first tau parameter (tauA) is always an effective sample size (n), but we should think about whether the tauB priors should be specified in terms of the "raw" tauB (sum of squares), variance, standard deviation, or coefficient of variation (CV). The current underlying math is closer to CV and, while nonstandard, there are advantages to that formulation in terms of avoiding units errors. I'm leaning towards n and CV but then making the math/code really match that approach.

Now that we've talked through that, the special case for seedling_mortality starts to make more sense. If you happen to have a low mortality rate you might end up with an (over) confident prior. The "special case" of 1 ensures a broad prior.

If we wanted to get really fancy, we could set up the code to accept n, CV, tauA, or tauB as named arguments and then calculate tauA if n is provided and calculate tauB if CV is provided.

[Thomas-Sedhom]

Hi @mdietze, if you have some time, I have a questions here.

First, based on your excellent analysis of each variable, I understand that
We have a Gamma prior:

• Gamma(𝑎, 𝑏) = Gamma(tauA, tauB) =>
• tauA = 𝑛0 / 2, with n0 being the prior sample size, also the default value of tauA = 0.01, so 𝑛0 = 0.02
• tauB = SS / 2, where SS is the prior sum of squares
• SS = n0 * prior_variance => tauB = ( n0 * prior_variance / 2 ) = ( tauA * prior_variance )

In the existing code, prior variances are calculated as:
𝜎 ^ 2 = 0.001 × trait.average ^ 2 This corresponds to 𝜎^2 = CV^2 * 𝜇^2 , where CV is the coefficient of variation. From the code logic, we can infer: CV = 0.0316

For seedling_mortality, we ignore this calculation because the original value might be extremely small, which would produce an overly confident prior with very narrow uncertainty. To prevent that, we force 𝜎^2 = 1 for seedling_mortality.

From the previous discussion, it seems the next step is to make n, CV, tauA, tauB parameters that users can pass. If they are not provided, the code should fall back to the old default values. We also need to document this full flow of calculations clearly.

Finally, I have a small question regarding tauB, theoretically we expecttauB = ( n0 * prior_variance / 2 ), so why does the code calculate it as tauB = min(gamma_tau, prior_variance) ?

[dlebauer]

Thanks for bringing this up and coming up with a solution. #3850 is a good cleanup of early MA prior logic that was mostly pragmatic and aimed at numerical stability/convergence while automating the analysis.
These choices were not fully documented, but old commit messages from 2011 illustrate the rationale:

• “range checking variance priors against sample mean to ensure priors don't swamp data.” (f6bcb78)
• “minor hack to get seedling mortality to not have a singularity.” (c5bb9d2)
• “hack around tauprior to get meta-analysis running.” (5924e26)

So I’d treat the legacy code as useful heuristics, but moving to named prior_n / prior_cv (with optional tauA / tauB) and documenting the formulas is definitely an improvement.