Are parameter estimates for cosy constrained to be positive?

[mabarbour]

Are parameter estimates for cosy constrained to be positive? Below is a reproducible example with the attached dataset example_df.xlsx.

brms code:
library(brms)
summary(brm(Y ~ T1 + T2 + (1 | Group) + cosy(gr = Group), data = example_df, control = list(adapt_delta = 0.95)))

What I believe is the equivalent code in nlme:
library(nlme)
summary(lme(Y ~ T1 + T2, data = example_df, random = ~ 1 | Group, correlation = corCompSymm(form = ~ 1 | Group)))

The population-level estimates are virtually the same, but nlme estimates a negative value for the rho whereas brms estimates a positive value. In addition, the lower bound for the brms estimate is zero.

[paul-buerkner]

Yes, the correlation is bound by zero in brms. Is there any theoretical reason, why a negative value would make sense for a cosy structure? I understand that nlme may find negative estimates but that doesn't (necessarily) mean they are sensible.

[mabarbour]

It is my understanding that this correlation will be negative if the variance between groups is smaller than the variance within groups. This happens commonly in split-plot field experiments. One example is given by Nelder 1977:

It has been argued that one way to account for this "excess variance" is to fit a compound symmetry correlation structure (Ben Bolker on Stackoverflow). Without the compound symmetry correlation structure, these group-level effects would be estimated as effectively zero in a mixed-model. This occurs in both nlme or lme4:

summary(lme(Y ~ T1 + T2, data = example_df, random = ~ 1 | Group))
or
summary(lmer(Y ~ T1 + T2 + (1 | Group), data = example_df))

To my understanding, this suggests that the residual correlations of the model have been mis-specified, which will affect the uncertainty in the parameter estimates.

Below are some additional links to previous discussions I've seen on this topic, by people more knowledgeable than me:

[R-sig-ME] FW: Negative estimates of variance component

The above thread is continued more clearly here.

[paul-buerkner]

I have seen this argument about negative covariances/correlations before and I would be open to allow negative cosy correlations. However, this negative correlation cannot be arbitary as, at a certain point x between -1 and 0 (where x depends on the dimention of the correlation matrix), the implied correlation matrix will not be positive definite anymore, that is, not be a proper correlation matrix. We need to tell Stan about this x as otherwise it will likely cross this border during the exploration of the parameter space and fail with a lot of divergent transitions or just fail completely right away. I don't have time to do the math myself right now or search the literature for this x, so I would be grateful for any pointers in the right direction.

[hansvancalster]

I believe this article has much to say about this topic and contains the math you are looking for:
Molenberghs G. & Verbeke G. (2011). A note on a hierarchical interpretation for negative variance components. Statistical Modelling: An International Journal 11 (5): 389–408. URL: http://journals.sagepub.com/doi/10.1177/1471082X1001100501. DOI: https://doi.org/10.1177/1471082X1001100501.

[paul-buerkner]

Thanks. As noted in the paper, the eigenvalues of an equi-correlation matrices are 1 - r and 1 + (n - 1) * r where r is the correlation and n is the dimension of the correlation matrix. The first eigenvalue is trivially positive for all -1 < r < 1 and the second one is positive for all r > - 1 / (n - 1). The latter equation gives us the lower bound of r so that the correlation matrix is still positive definite.

[paul-buerkner]

Ok. Should be working now. You example_df seems to be corrupted when I download it. Would you mind trying out this new feature using the github version of brms?

[mabarbour]

I've tested it and the github version allows negative correlations now. Thanks for such quick implementation and thanks to @hansvancalster for pointing out where to find the math.

Note that I had to increase adapt_delta = 0.9999 and apply a strong prior on cosy. I arbitrarily chose a normal prior with mean zero and a low variance (normal(0, 0.1)). Both the low variance on the prior and high adapt_delta were necessary to get the model to converge with no divergent transitions. This also worked when I applied a non-normal family for the error distribution.

[paul-buerkner]

It would be interesting to know whether this just happens because of the negative correlations. Did the same happen before? Matt Barbour <notifications@github.com> schrieb am Sa., 18. Apr. 2020, 14:09:

…

I've tested it and the github version allows negative correlations now. Thanks for such quick implementation and thanks to @hansvancalster <https://github.com/hansvancalster> for pointing out where to find the math. Note that I had to increase adapt_delta = 0.9999 and apply a strong prior on cosy. I arbitrarily chose a normal prior with mean zero and a low variance (normal(0, 0.1)). Both the low variance on the prior and high adapt_delta were necessary to get the model to converge with no divergent transitions. This also worked when I applied a non-normal family for the error distribution. — You are receiving this because you modified the open/close state. Reply to this email directly, view it on GitHub <#878 (comment)>, or unsubscribe <https://github.com/notifications/unsubscribe-auth/ADCW2ADFNXYAYPCUUWKYYKTRNGKANANCNFSM4MJQBTHQ> .

[mabarbour]

This did not happen before. The model converged with no divergent transitions at adapt_delta = 0.95 and default priors.

[paul-buerkner]

Hmm this either means the Lower bound is incorrect or negative values cause undesired issues. I would tend to go back to the 0 lower bound if we cannot get rid of these issues. Matt Barbour <notifications@github.com> schrieb am Sa., 18. Apr. 2020, 15:17:

…

This did not happen before. The model converged with no divergent transitions at adapt_delta = 0.95 and default priors. — You are receiving this because you modified the open/close state. Reply to this email directly, view it on GitHub <#878 (comment)>, or unsubscribe <https://github.com/notifications/unsubscribe-auth/ADCW2AE325QS4HLRAKCGDS3RNGR7NANCNFSM4MJQBTHQ> .

[paul-buerkner]

I have reverted back to lower=0 for now and will reopen this issue. If we can find a way to avoid the divergent transitions without super high adapt delta, I am happy to change it to the negative boundary.

[paul-buerkner]

With PR #1324 being merged, you can now choose bounds of parameters yourself via the prior interface. Hence I think this issue can be closed.