Calculating derived statistics as in BPA 2012 is wrong - proof by a simulation study and proposed remedy

[Tomas T.]

Dear all,

I have done a simulation study as some of you suggested, please find attached, with scripts. I believe it clearly shows that the method to calculate the posterior of the correlation coefficient (r) on posterior distribution of two correlated parameters, as presented in BPA 2012, critically underestimates the uncertainty of the correlation coefficent. I also propose alternative method that performs well, although I think it isn't perfect, and I would appreciate if some statistically experienced folks here could review it and say if it is good enough, or "as good as it can be", or if there is a better way to do it.

I was quite frustrated on Friday by how difficult is for people to understand each other, even in the field of statistics, where everything should be clear, shouldn't it? :-) But this is clearly a very confusing issue, and that's why noone noticed this problem till now, although great number of smart people were reviewing and using the method!

Clarity was obscured by several confusing arguments. First, I could not understand the reasoning behind an argument that several people brought - "it is because your sample variation is too little" - basically blaming the scenario for the discrepancy. How can you blame a parameter of a scenario? Shouldn't a good method work well in all scenarios? Especially, when the low sample variation means that you have a really good data. So, paradoxically, the method performs worse when the data is better, which is exactly what Figure 1 in the attached study shows. I hope the study also clearly shows that the scenario is not to blame.

Another confusing argument that appeared multiple times was the frequentist vs bayesian. Yes, there is a difference, but you cannot ascribe a huge discrepancy to this!

But probably the biggest confusion is that the correlation coefficient is just a single value. No, it is not. The correlation coefficient is an estimate of the (unknown) parameter that has its own confidence interval and uncertainty - see the output of cor.test()! It is no different from any other parameter estimate, like e.g. regression slope.

So, I really like what Marc & Michael taught me, that we can calculate derived quantities on the posterior samples! This is true and it really is a beauty of the bayesian world. However, here, the quantity that we calculate for each posterior sample is not a deterministic one. It is an estimate of a true unknown parameter, and this estimate goes with some uncertainty that cannot be neglected.

At times it seemed like people just assume the method is correct and when the results look weird, try to find problem everywhere else.

To achieve clarity, we must first agree on some premise: on the properties of a good posterior and a way to test it.  Posterior is a distribution of the __true value__ of the parameter given the data sample. I suggest to measure this by measuring coverage of the 95% CI, i.e. the probability (in the simulated data sets) that the 95% CI will cover the true value. You might again bring the frequentist vs bayesian argument. But there are established methods that test "bayesian models" on a similar principle: see e.g. Simulation-based Calibration - the SBC package, which Matthijs recommended in another thread (thanks for that, it is really good!). If you don't agree with this coverage premise, please bring an alternative one, along with the way how to test it!

So, our goal is to get the posterior distribution of the __true value__ of the parameter (correlation coefficient). But what the BPA 2012 method gives us, is the posterior of the __estimator of the true parameter__ ! It is basiclly posterior of the mean of the posterior we want, and thus the uncertainty is strongly underestimated, as confirmed by the Figure 1 in the attachment.

I believe the attached study brings an evidence that the BPA method doesn't work, and I attempted to explain why is that. I provided a corrected method. I would greatly appreciate your opinions and inputs on the performance of the "corrected" method, or perhaps a suggestion of different method to do this? 

I publicly apologize to Mark & Michael for this relatively abrupt way of bringing this up, but I am invested in this as I have a very urgent ongoing work, along with two other published papers that use this method. So I am also greatly responsible for spreading this method :-) and not seeing the mistake earlier.

Tomas

PS: John, thank you for your supporting voice, it was nice to have at least someone confirm that I am not completely crazy! :-DD

On Sat, 22 Feb 2025 at 00:35, 'John Clare' via hmecology: Hierarchical Modeling in Ecology <hmec...@googlegroups.com> wrote:

Oh, I suppose the distinction is between finite sample inference on the correlation between parameters (the parameters were correlated in these years with data) and otherwise (over many many years not yet sampled we expect the parameters will be correlated as…).

On Feb 21, 2025, at 4:25 PM, 'John C' via hmecology: Hierarchical Modeling in Ecology <hmec...@googlegroups.com> wrote:

Oh,  I see what Tomas means just playing around with this quickly. I've repeated the following a few times with different nData and so forth.  Yeah, seems like one should probably should use the correlation across the posterior of posterior-predicted parameters. Kind of like a true Monte Carlo assessment or parametric bootstrap or something.

John

###Start by simulating "parameters"

set.seed(8675309)
np <- 15
r <- 0.25 # pearson r
X <- rmvnorm(n, sigma = matrix(c(1, r, r, 1), 2))
X2<-plogis(X)

cor.test(X[,1], X[,2]) ###0.033, 95 %CI -0.49 - 0.54
cor.test(X2[,1], X2[,2]) ###0.042, 95 %CI -0.48 - 0.54

###Let's say X2[, 1] is survival probability,
###and X2[, 2] is...the probability of a successful nest or something.

nData<-1000
###we'll make the simplifying assumption that in every year,
###we observe an independent set of 100 birds that survive or do not
###and have a successful nest or not.

Data1<-matrix(NA, nData, np)
Data2<-matrix(NA, nData, np)

for (i in 1:np){
  Data1[, i]<-rbinom(nData, 1, X2[i, 1]) ##survived or not
  Data2[, i]<-rbinom(nData, 1, X2[i, 2]) ##nest or not
}

library(nimble)
blah<-nimbleCode({
  mu.phi~dlogis(0, 1)
  sigma.phi~dunif(0, 3)
  mu.gamma~dlogis(0, 1)
  sigma.gamma~dunif(0, 3)
 
  for (i in 1:15){
    logit(phi[i])~dnorm(mu.phi, sd=sigma.phi)
    logit(gamma[i])~dnorm(mu.gamma, sd=sigma.gamma)
    logit(simphi[i])~dnorm(mu.phi, sd=sigma.phi)
    logit(simgamma[i])~dnorm(mu.gamma, sd=sigma.gamma)
    for (j in 1:1000){
      Data1[j, i]~dbern(phi[i])
      Data2[j, i]~dbern(gamma[i])
    }
  }
})

Data=list(Data1=Data1, Data2=Data2)

Sim2<- nimbleModel(code = blah, name = "blah", constants = NULL,
                   data = Data, inits = NULL, calculate=FALSE)

SimMCMC2<-configureMCMC(Sim2, monitors=c("phi", "gamma", "simphi", "simgamma"), useConjugacy = FALSE)

SimRMCMC2<-buildMCMC(SimMCMC2)

SimCMCMMC2<-compileNimble(SimRMCMC2,Sim2)

samples<-runMCMC(mcmc=SimCMCMMC2$SimRMCMC2, niter=20000, nburnin = 10000, thin=20, nchains=1, samplesAsCodaMCMC=FALSE)

corMCMC<-rep(NA, 500)
for (i in 1:500){ ###number of stored MCMC samples
  corMCMC[i]<-cor.test(samples[i, 1:15], samples[i, 16:30])$estimate
}

quantile(corMCMC, probs=c(.025, .5, .975)) ###too narrow

corMCMC2<-rep(NA, 500)
for (i in 1:500){ ###number of stored MCMC samples
  corMCMC2[i]<-cor.test(samples[i, 31:45], samples[i, 46:60])$estimate
}

quantile(corMCMC2, probs=c(.025, .5, .975)) ###about right.

On Friday, February 21, 2025 at 12:44:33 PM UTC-7 hender...@gmail.com wrote:

Hi all,

Tomas, this is an interesting point you’ve brought up. I thought I'd throw in my 2 cents.

If i'm correct, your example reflects the extremely small measurement error (sd = 0.05) assigned to the 15 data points. If you check the standard deviation of each point’s values across your 5000 MCMC draws (column-wise SDs of ~0.048–0.052), you’ll find almost no variation from one draw to another. In other words, if the model believes each point is pinned down with minimal measurement error, it will naturally infer that the correlation implied by those 15 points is quite certain.

Frequentistically, cor.test() on 15 data points yields a wide confidence interval for the sample correlation i.e., “if we repeatedly collected new datasets of size 15, how variable would the estimated correlation be?” By contrast, a Bayesian approach (like Kéry & Schaub) asks, “Given the model, data, and pzriors, how plausible are different correlation values?” If the model is told there’s almost no uncertainty in each observation, it makes sense how it could produce narrowly distributed posterior for the correlation, even with just 15 points.

I don't think your toy example necessarily show that the BPA 2012 approach is wrong. Instead, i think it just highlights that extremely small measurement error leads to overly confident posterior inferences. Once you allow more realistic measurement error (or other relevant variance components), the distribution of correlation values will widen and better reflect the true uncertainty.

RE: your concern about p‐values and the probability that a correlation is positive or negative. A Bayesian workflow does let you determine how often the posterior distribution of the correlation falls above/below zero. Unlike a single “mean correlation,” you get a distribution of correlation values —one for each posterior draw- which directly answers whether the correlation is likely positive, negative, or near zero. There’s no need for a separate “SE of correlation” inflation step, because the Bayesian draws already capture all uncertainty (provided the model itself is realistic..).

Finally, I like Matt's comment suggesting a more thorough simulation study!

Nick

On Friday, February 21, 2025 at 10:48:24 AM UTC-8 matthijs....@gmail.com wrote:

Hey Tomas,

I think I understand what you're saying. I believe one way of framing it is that you don't think the posterior distribution over the empirical correlation coefficient (i.e., computing the coefficient for each posterior draw) gives you adequate uncertainty over the correlation. This may be correct, but I'm not sure; in some ways it is a hacky solution. But like I said, if you are interested in formally modeling the correlation, assuming your rates (or link-function transformed versions thereof) are multivariate normal distributed will let you actually estimate those correlations, albeit under that assumed parametric model.

What I would suggest is simulating bivariate normal data, then analysing those in (1) a bivariate normal model where you estimate the correlation as a parameter (i.e., the exact data-generating process) and (2) a model with two univariate normal parameters and then estimating the correlation coefficient in the method described in BPA. That would be a better way to test this issue compared to the simulation you provided. This is roughly what John has suggested as well, I believe.

If you do end up doing this, I'd be keen to see the results.

Cheers,

Matt

On Fri, Feb 21, 2025 at 10:04 PM Tomas T. <tomas.t...@gmail.com> wrote:

Hi Matthijs,

You're not interested in the frequentist 95% CI around the correlation coefficient of your sample, in this case the 2 vectors of length 15 that you simulated. The correlation coefficient is just a single value that you're computing for each MCMC sample from your posterior.

No, wrong. Most of the people (that includes the cited references, and also other studies which used this method) are not interested just in a single value of the correlation coefficient. Most people are interested in inference, that is, the probability that the correlation coefficient is positive/negative/different from zero, which is what you can do with the CI of the correlation coefficient returned by cor.test() - and which you aim to do with your resultant posterior distribution. Or, practically equivalently, in frequentistic world, people are interested in p-value - what is the probability that they mistakenly reject the null hypothesis (correlation equals to zero) when it is true.

 

Therefore, when you do this for each MCMC sample, you have draws from the posterior of the correlation between your time varying parameters.

No. Using this method, you have posterior of the mean estimate of the correlation coefficient, which is useless, because you cannot use it for inference of the probability that the correlation coefficient is positive/negative/different from zero. Which is, again, what the cited references and other studies using this method aim to do.

The example I provided should clearly demonstrate how this method could lead to wrong inference and thus wrong conclusions.

Tomas

  

Re: the correlation matrix, I don't think it would slow down too much. If you just model the varying rates as univariate normal distributions, this implies a prior distribution that there's 0 correlation between the parameters. If you model it multivariate normal, you can still express an uninformative prior over the correlation matrix, but at least your prior is that the rates may be correlated.

Hope this helps.

Matt

On Fri, Feb 21, 2025 at 4:36 AM Tomas T. <tomas.t...@gmail.com> wrote:

Dear Matthijs,

answering below:

On Thu, 20 Feb 2025 at 17:46, Matthijs Hollanders <matthijs....@gmail.com> wrote:

You're simulating very little random noise around your simulated X with your "MCMC samples", which is what's reflected in your "CI".

Of course, that was the intention and that's the point, to provide an example that would clearly demonstrate that the method is in principle wrong. Doesn't matter if it is realistic or not. It clearly demonstrates that the method ignores the variance of the correlation statistics, which leads to wrong inference.
 

Btw, if you want to directly estimate correlations, you can model the two (or more) time-varying demographic parameters as multivariate normal and estimate the correlation coefficient(s). 

True, you could do that, then, if you have e.g. 4 parameters, you would have to estimate the whole correlation matrix... I am not sure how much it would slow down the model fit though. And what would the change of the model mean - perhaps also the estimates might be different then?

Regards,

Tomas

 

On Thursday, February 20, 2025 at 11:43:50 PM UTC+11 tomas.t...@gmail.com wrote:

Hello all,

suppose you want to calculate correlation on two variables estimated by a model fitted in bayesian software. For example, a model that calculated two demographic parameters over years, and you have the posterior samples. How to do this?

Until now, I heavily used the method recommended in Bayesian Population Analysis using WinBUGS (Kery & Schaub 2012, code for Fig. 11-8 (Web appendix 1)), which calculates the pearson correlation coefficient for each posterior sample and then uses that as posterior distribution of the correlation. I thought this method is neat. Only now I realized this method is probably wrong?! It does take into account the uncertainty of the estimated parameters, but it completely ignores the uncertainty of the pearson correlation coefficient itself. So, if you have relatively precise estimates of the parameters, it can make very weak correlations appear very strong! I will illustrate it by an example below.

For me this realization is relatively serious, as I already used it in a published paper. Perhaps this is also case of the cited paper Schaub et al (2011)? I understand mistakes like this happen, I make them myself, but I think it would deserve some big red warning so that no more people make this mistake.

Have you used this method yourself? Or are you using a different one? How do you solve this? Below, I also proposed a remedy to this - instead of taking just the mean estimate of pearson, sample from its error distribution. I am not 100% sure this is statistically correct though.

Please let me know!

Best regards,

Tomas

PS: Here is a reproducible code that demonstrates the problem:

library(mvtnorm)

### generate two correlated vectors - true values
set.seed(123)
n <- 15
r <- 0.25 # pearson r
X <- rmvnorm(n, sigma = matrix(c(1, r, r, 1), 2))

cor.test(X[,1], X[,2])
# r = 0.1944498
# you see the CI is very wide: -0.3529724 to 0.6426943
# which means very weak correlation

### simulate the "MCMC samples" of these two vectors, with small uncertainty
samp <- 5000
A <- matrix(rnorm(n*samp, X[,1], 0.05), samp, byrow = TRUE)
B <- matrix(rnorm(n*samp, X[,2], 0.05), samp, byrow = TRUE)


### Calculate correlation coefficients as in BPA 2012, Code for Fig. 11-8 (Web appendix 1)
correl <- c()
for (i in 1:samp) {
correl[i] <- cor(A[i,], B[i,])
}
# Credible intervals of correlation coefficients
quantile(correl, c(0.05, 0.5, 0.95))

# the resultant CI is very narrow - i.e. the correlation suddenly seems super stong:
#        5%       50%       95%
# 0.1628700 0.1936207 0.2257177

#

# This is definitely wrong, as we know that the true correlation is very weak!

####################################################
#  Could this be a possible remedy? I am not sure:
#  to sample not just mean value, but from the probability distribution of each statistics
#  (it is also an approximation, since the distribution of Pearson is not exactly normal, see
# https://en.wikipedia.org/wiki/Pearson_correlation_coefficient#Using_the_exact_distribution

correl2 <- c()
for (i in 1:samp) {
st <- cor.test(A[i,], B[i,])
correl2[i] <- rnorm(1, st$estimate, diff(st$conf.int)/1.96/2)
}
quantile(correl2, c(0.05, 0.5, 0.95))

--
*** Three hierarchical modeling email lists ***
(1) unmarked: for questions specific to the R package unmarked
(2) SCR: for design and Bayesian or non-bayesian analysis of spatial capture-recapture
(3) HMecology (this list): for everything else, especially material covered in the books by Royle & Dorazio (2008), Kéry & Schaub (2012), Kéry & Royle (2016, 2021) and Schaub & Kéry (2022)
---
You received this message because you are subscribed to the Google Groups "hmecology: Hierarchical Modeling in Ecology" group.
To unsubscribe from this group and stop receiving emails from it, send an email to hmecology+...@googlegroups.com.
To view this discussion visit https://groups.google.com/d/msgid/hmecology/20b2720e-f495-4de3-b260-5cd6ee98b1fcn%40googlegroups.com.

--
*** Three hierarchical modeling email lists ***
(1) unmarked: for questions specific to the R package unmarked
(2) SCR: for design and Bayesian or non-bayesian analysis of spatial capture-recapture
(3) HMecology (this list): for everything else, especially material covered in the books by Royle & Dorazio (2008), Kéry & Schaub (2012), Kéry & Royle (2016, 2021) and Schaub & Kéry (2022)
---
You received this message because you are subscribed to the Google Groups "hmecology: Hierarchical Modeling in Ecology" group.
To unsubscribe from this group and stop receiving emails from it, send an email to hmecology+...@googlegroups.com.
To view this discussion visit https://groups.google.com/d/msgid/hmecology/ab522e71-b15a-471c-ac6c-35fab77d4a53n%40googlegroups.com.

--
*** Three hierarchical modeling email lists ***
(1) unmarked: for questions specific to the R package unmarked
(2) SCR: for design and Bayesian or non-bayesian analysis of spatial capture-recapture
(3) HMecology (this list): for everything else, especially material covered in the books by Royle & Dorazio (2008), Kéry & Schaub (2012), Kéry & Royle (2016, 2021) and Schaub & Kéry (2022)
---
You received this message because you are subscribed to the Google Groups "hmecology: Hierarchical Modeling in Ecology" group.
To unsubscribe from this group and stop receiving emails from it, send an email to hmecology+...@googlegroups.com.
To view this discussion visit https://groups.google.com/d/msgid/hmecology/45C15B2C-3339-498F-B148-C126A84C2474%40wisc.edu.

[Tomas T.]

Sorry, here is a version of the simulation script with better comments

Tomas

[Matthijs Hollanders]

I will have more time to read your simulations later. But I think there's confusion about parameters here. You say the correlation is not a single but "an estimate of the (unknown) parameter that has its own confidence interval and uncertainty". That's not really true. The correlation coefficient is cov(x, y) / (sd(x) * sd(y)). This is an empirical computation just like mean() and sd(). If you want to know the posterior distribution of the parameter, you should fit a model with a correlation coefficient or a correlation matrix. 

The best comparison I can think of is a random effects model. Let's say there are 5 grouping levels of your random effect, and you're trying to estimate the mean and standard deviation of the population, x_i ~ normal(mu, sigma). Because there's only 5 grouping levels, there will be a lot of uncertainty in the parameters mu and sigma. However, as generated quantities, you can still compute the means and SDs of the group-level effects, mean(x[1:5]) and sd(x[1:5]), for each posterior draw, which will have much less uncertainty. The exact same thing is happening when you compute the empirical correlation coefficient for each posterior draw. It is almost certainly not going to capture the posterior distribution of the parameter had you included it in the model as a parameter.

Idk what they were trying to do in BPA and if they were implying it's the same as putting the parameter in the model. If you want to evaluate the correlation statistically, you should model it.

To view this discussion visit https://groups.google.com/d/msgid/hmecology/CAGy9Uw8KWb7Y0X66uoEbERBSUD-49ebc0%2BXAnzPDVxhzYTtDEQ%40mail.gmail.com.

[Matthijs Hollanders]

Btw, for anyone wondering, modeling correlation matrices is now trivial in NIMBLE and Stan. Both include LKJ priors for the Cholesky factors of correlation matrices. In NIMBLE, where D is number of dimensions, it's an upper triangular:

mu[1:D] ~ dnorm()

sigma[1:D] ~ dexp()

Omega_U[1:D, 1:D] ~ dlkj_corr_cholesky(1, D)  # upper triangular of cor. matrix

Sigma_U[1:D, 1:D] <- diag(sigma[1:D]) %*% Omega_U[1:D, 1:D]  # upper triangular Cholesky factor of cov. matrix

for (i in 1:I) {

  y[i, 1:D] ~ dmnorm(mu[1:D], cholesky = Sigma_U[1:D, 1:D], prec_param = 0)

}

Omega[1:D, 1:D] <- t(Omega_U[1:D, 1:D]) %*% Omega_U[1:D, 1:D]  # recover correlation matrix

In Stan:

parameters {

  vector[D] mu;

  vector<lower=0>[D] sigma;

  cholesky_factor_corr[D] Omega_L;

}

model {

  y ~ multi_normal_cholesky(mu, diag_pre_multiply(sigma, Omega_L));

}

generated quantities {

  matrix[D, D] Omega = multiply_lower_tri_self_transpose(Omega_L);

}

And for random effects, you almost always want to use non-centered parameterisation. In NIMBLE, if theta were the multivariate realisations for the random effects:

theta[1:I, 1:D] <- z[1:I, 1:D] %*% Omega_U[1:D, 1:D] %*% diag(sigma[1:D])

for (i in 1:I) {

  for (d in 1:D) {

    z[i, d] ~ dnorm(0, 1)

  }

}

Or in Stan (only changes shown):

parameters {

  matrix[D, I] z;

}

model {

  to_vector(z) ~ std_normal();

  matrix[D, I] theta = diag_pre_multiply(sigma, Omega_L) * z;

}

In Stan, I wrote a function that returns the vector of correlations from the Cholesky factor of a correlation matrix, and a function that gets the number of correlations in a correlation matrix:

/**
 * Return vector of correlations from lower triangular Cholesky factor of correlation matrix
 * @param O_L  Lower triangular Cholesky factor of correlation matrix
 * @return     Vector of correlations
 */
vector rho(matrix O_L) {
  int R = cols(O_L);
  int n_rho = cors(R);
  vector[n_rho] rho;
  matrix[R, R] O = multiply_lower_tri_self_transpose(O_L);
  int i = 1;
  if (n_rho) {
    for (r in 1:(R - 1)) {
      for (rr in (r + 1):R) {
        rho[i] = O[r, rr];
        i = i + 1;
      }
    }
  }
  return rho;
}

// number of correlations

int n_cor(int R) {
  return (R * (R - 1)) %/% 2;
}

Cheers,

Matt

[Tomas T.]

Not true. The pearson correlation coefficient __is__ an estimate of the (unknown) parameter that has its own confidence interval and uncertainty. You must look at it that way if you want to calculate its posterior distribution, because posterior is the probability distribution of the true (unknown) parameter. Same way, you can look at mean() as an estimate of the true mean, and sd() is an estimate of true standard deviation. 

Feel free to fit a model that estimates the correlation as a parameter in multivariate normal and send us results of simulation study measured the same way as I did. I didn't feel the need to do so for two reasons: first, I am pretty sure results will be similar to the corrected method, maybe slightly better. Second, as I wrote in the study, my aim is to develop a post hoc method, which you can run post hoc on any already calculated models. I have already some models like that that ran for days and don't want to re-run them again just to calculate the correlation. Practical reasons.

Tomas

[John C]

Hi Tomas, 

I think I did not initially understand your question/issue, so apologies for the earlier responses (unhelpful in hindsight). If I understand now, this really is an issue of finite population inference (not finite sample, as I previously called it) vs. broader. I do not have any good idea about  how to transform the finite population estimate into  some broader estimate incorporating parametric/estimate uncertainty arising from the limited sample of years (or intervals), and can't give any advice beyond using an MVN prior as Matt suggested.

But I do think it might be worth considering whether finite population inference might be sufficient for your purposes.The Gelman/Hill book discusses this in the context of ANOVA (using an example somewhat similar to Matt's), where there might be a group specific random effect beta[group]. You might be interested in predicting new groups or trying to describe the variation across all groups, including the unsampled ones--the hyper-parameter sigma_beta might be useful for that. But, you might alternatively be interested in describing the variation across the groups sampled. Iteratively deriving sd(beta[1:ngroups]) is sufficient for that. (note, sd(X[,1]) is also just a calculation that has no error if X[,1] is measured without error). Say the groups are continents and data collection occurs in each. There are only 7, and so sigma_beta is likely fairly imprecise for describing the 7 continents sampled. There's also the question of whether some broader population of continents is of interest, and frankly, what broader population of continents might have been selected other than the 7 existing continents.

I think the method in BPA 11.8 is fine if interpreted as a finite population estimate that accounts for the fact that the data leading to the parameter estimates is a sample, but not that the strata that define the replication in the correlation are a sample. The correlation coefficient cor() is calculated based on the variance/covariance between two variables--if we are just interested in the correlation between X[,1] and X[, 2], there similarly is no error . Cor.test returns uncertainty under the assumption that one is trying to generalize to a broader population not fully encompassed by the rows of X and some parametric assumptions about estimation error. For demographic models, I think it's ok to use the BPA method to make statements about the correlation between (e.g.) survival and fecundity over the time period in which data-collection occurred, conditioning on the certain location considered.  On the other hand, one could employ some MVN prior for the parameters of interest and interpret correlation parameter rho as an estimate of the correlation between survival and fecundity across some broader population of time intervals. One question is whether this of interest. To be honest, for the demographic models that motivate the question, I can't see how it's important descriptively because this broader population/sampling frame of time intervals can be difficult to define, and it is probably not sampled in any defensible formal way. It might be useful from a forecasting perspective more reliant on posterior prediction, but there are other issues here. Does this make sense?

John

[John C]

"I think the method in BPA 11.8 is fine if interpreted as a finite population estimate that accounts for the fact that the data leading to the parameter estimates are a sample [leading to estimate error/uncertainty in the finite population parameters like survival/fecundity and their correlation], but not that the strata that define the replication in the correlation are a sample". Something something about an edit button or drink being useful. Fin.

[Tomas T.]

Dear all,

attached you will find an update to my simulation study. I have added:

1) The direct estimation of the correlation coefficient in the model using multivariate normal, as Matthijs proposed (method "model_est"). As I expected, the results are comparable to all the other methods except BPA, but are slightly better than the post hoc methods. However, we still need to use good post hoc method for practical reasons. 

2) Added a method "ignore", which just calculates the correlation on the mean estimates of the parameters, ignoring their uncertainty. 

Both newly added methods further confirm that the BPA method is off. Newly added Figure 4 shows the comparison of the CI width, which is comparable for all methods except  the  BPA method, where CIs are much narrower. Figure 1 confirms that the BPA method is strongly underestimating the CI width, whereas all the other methods are more or less fine.

Just a note on the email from Marc that this conference received on Wednesday Feb 26th - Marc has sent this email already on Friday 20th, I received it already then because it was CC'ed to me directly. So he was not responding to my simulation study, which I sent on Tuesday 25th. On the contrary, my study was already response to his email. Somehow the googlegroups.com messed this up.

Response to John C's email in this thread. I understand what you mean, I actually had a poster on this "finite population inference" in Riga (EOU 2011). However, I would argue that this "finite population inference" is of limited use. The inference we normally do, both frequentist and bayesian, corresponds to what you call "broader". The definition of posterior distribution itself corresponds to this "broader" sense. It is a distribution of the true value of some parameter of a process, that generated the data we have - thus, we are looking at our data as if it were a sample - one of many possible samples - from a process whose parameter we want to estimate. I believe this "one of many possible samples" corresponds to what you call the broader inference. And thus, as I already mentioned in my email and the study, the BPA method is not constructing the posterior of the `r` - P(r), but the probability distribution of the r̂ (r hat) - P(r̂). So your email kind of philosophically meditates on how normal statistics works, and argues why we should do it differently. I am not sure that this direction is practical, and it can be quite confusing. I think the way back to clarity is to take the definition of posterior, and then to test whether it's width is appropriate - for which I propose the 95% CI coverage - not only I, but it is commonly used in bayesian world, see e.g. the amazing SBC package! (https://hyunjimoon.github.io/SBC/articles/SBC.html).  (BTW, I actually feel like your original email - the one with the model - was very helpful and understanding!).

Have a nice weekend!

Tomas

--

*** Three hierarchical modeling email lists ***
(1) unmarked: for questions specific to the R package unmarked
(2) SCR: for design and Bayesian or non-bayesian analysis of spatial capture-recapture
(3) HMecology (this list): for everything else, especially material covered in the books by Royle & Dorazio (2008), Kéry & Schaub (2012), Kéry & Royle (2016, 2021) and Schaub & Kéry (2022)
---
You received this message because you are subscribed to the Google Groups "hmecology: Hierarchical Modeling in Ecology" group.
To unsubscribe from this group and stop receiving emails from it, send an email to hmecology+...@googlegroups.com.

To view this discussion visit https://groups.google.com/d/msgid/hmecology/483a9a22-e950-4d8a-aff2-8da12dc5697fn%40googlegroups.com.