correlation structures for residuals

[ebmtnprof]

It would be wonderful to have some basic correlation structures available for the residuals in brms (compound symmetric, heterogeneous compound symmetric, unstructured, and auto-correlation; see gls examples below). With nlme, any of these can be combined with various random effect structures. Any chance brms could support these? Cheers, Emily

Start with a "compound symmetric" model that only has one covariance (shared across all time points) and one variance (also shared across all times). This is the model assumed by traditional repeated measures ANOVA.

cs <- gls(hr ~ time, correlation=corCompSymm(form= ~1 |subject), 
    na.action=na.omit, method="ML", data=data)

compare to a "heterogeneous compound symmetric" model that holds covariances equal but allows variances to be unique

hcs <- gls(hr ~ time, correlation=corCompSymm(form= ~1 |subject), 
  weights=varIdent(form = ~ 1|time), 
  na.action=na.omit, method="ML", data=data)

compare to an "unstructured" model that allows both the variances and covariances to be anything

unstruct <- gls(hr ~ time, correlation=corSymm(form= ~1 |subject), 
       weights=varIdent(form = ~ 1|time), 
       na.action=na.omit, method="ML", data=data)

how about an autocorrelation model?

ar <- gls(hr ~ time, correlation=corAR1(form= ~1 |subject), 
   weights=varIdent(form = ~ 1|time), 
   na.action=na.omit, method="ML", data=data)

[paul-buerkner]

Thanks for the suggestions. Just to order this a bit: Which of those correlation structures are usually used for residuals and which for "random" effects?

[ebmtnprof]

They are all used for both residuals and group effects (I agree with you that is a better term than random - I just have the bad habit due to trying to not confuse students), except for autocorrelation, which only applies to the residuals (and I think that is implemented in brms, isnt' it?). The additional wish-list I have is that in nlme you have to have the same nesting for the group effects and the residuals, so for example, if I have time in individuals in families, I can't have an unstructured group-effects matrix at the family level and then an autocorrelation residual matrix at the individual level. So, if you get into working on this, maybe that's an additional extension that could be addressed..... Cheers, Emily

On Sat, Apr 14, 2018 at 12:56 AM, Paul-Christian Bürkner < [email protected]> wrote:

Thanks for the suggestions. Just to order this a bit: Which of those correlation structures are usually used for residuals and which for "random" effects?

— You are receiving this because you authored the thread. Reply to this email directly, view it on GitHub https://github.com/paul-buerkner/brms/issues/403#issuecomment-381311644, or mute the thread https://github.com/notifications/unsubscribe-auth/AkNSOnWZFEeAoOyWXTs3d8e5VQCFZfp2ks5toavKgaJpZM4TT0iy .

--

Emily A. Butler

Associate Professor & Graduate Director Family Studies and Human Development College of Agriculture & Life Sciences University of Arizona Tucson, AZ, 85721-0033

[jgabry]

On Sat, Apr 14, 2018 at 3:56 AM Paul-Christian Bürkner < [email protected]> wrote:

Thanks for the suggestions. Just to order this a bit: Which of those correlation structures are usually used for residuals and which for "random" effects?

I would assume both. In the frequenting formulation this is based on, they’re all basically part of the error, right?

—

You are receiving this because you are subscribed to this thread. Reply to this email directly, view it on GitHub https://github.com/paul-buerkner/brms/issues/403#issuecomment-381311644, or mute the thread https://github.com/notifications/unsubscribe-auth/AHb4Q6o9lGN9dxrIU0mTKB5q13ehQYvEks5toavKgaJpZM4TT0iy .

[jgabry]

frequentist not frequenting

On Sat, Apr 14, 2018 at 1:58 PM Jonah Sol Gabry [email protected] wrote:

On Sat, Apr 14, 2018 at 3:56 AM Paul-Christian Bürkner < [email protected]> wrote:

Thanks for the suggestions. Just to order this a bit: Which of those correlation structures are usually used for residuals and which for "random" effects?

I would assume both. In the frequenting formulation this is based on, they’re all basically part of the error, right?

—

You are receiving this because you are subscribed to this thread. Reply to this email directly, view it on GitHub https://github.com/paul-buerkner/brms/issues/403#issuecomment-381311644, or mute the thread https://github.com/notifications/unsubscribe-auth/AHb4Q6o9lGN9dxrIU0mTKB5q13ehQYvEks5toavKgaJpZM4TT0iy .

[ebmtnprof]

you are correct Jonah, except that auto-correlation doesn't make sense for the group-level effects (at least not in any application I've seen or thought of)....

And just for some context, the reason I'm interested in having these implemented in brms is that I teach multilevel modeling to social science grad students (e.g., they don't have any background in math or computing, and most of them have pretty limited statistics training too) and I am very pleased to finally be able to do it from both a frequentist and a Bayesian perspective - brms has made that possible - and it would be great to be able to cover all typical models both ways. So far brms allows most of the complicated models (nonlinear, etc) but not some of the simple ones....

On Sat, Apr 14, 2018 at 10:59 AM, Jonah Gabry [email protected] wrote:

frequentist not frequenting

On Sat, Apr 14, 2018 at 1:58 PM Jonah Sol Gabry [email protected] wrote:

On Sat, Apr 14, 2018 at 3:56 AM Paul-Christian Bürkner < [email protected]> wrote:

Thanks for the suggestions. Just to order this a bit: Which of those correlation structures are usually used for residuals and which for "random" effects?

I would assume both. In the frequenting formulation this is based on, they’re all basically part of the error, right?

—

You are receiving this because you are subscribed to this thread. Reply to this email directly, view it on GitHub <https://github.com/paul-buerkner/brms/issues/403# issuecomment-381311644>, or mute the thread <https://github.com/notifications/unsubscribe-auth/ AHb4Q6o9lGN9dxrIU0mTKB5q13ehQYvEks5toavKgaJpZM4TT0iy> .

— You are receiving this because you authored the thread. Reply to this email directly, view it on GitHub https://github.com/paul-buerkner/brms/issues/403#issuecomment-381347925, or mute the thread https://github.com/notifications/unsubscribe-auth/AkNSOigqi6cbF9TluUpZk4JTDwohfHUmks5tojkBgaJpZM4TT0iy .

--

Emily A. Butler

Associate Professor & Graduate Director Family Studies and Human Development College of Agriculture & Life Sciences University of Arizona Tucson, AZ, 85721-0033

[paul-buerkner]

There are situations in meta-analysis where autocorrelation makes sense for group-level effects.

So basically here is the list of correlation structures to be implemented (with already implemented stuff checked):

Group-level effects:

• [x] unstructured
• [x] compound symmetry
• [ ] heterogenous compound symmetry
• [x] autoregressive of order 1
• [ ] Toeplitz

Residuals:

• [x] unstructured
• [x] compound symmetry
• [x] heterogenous compound symmetry
• [x] autoregressive of order 1
• [ ] Toeplitz

[ebmtnprof]

thanks! Yes, that would be great. And I just thought of one more, which is usually called "variance components" and has variances for each effect, but no covariances. I think you said that is already possible anyway?

On Sun, Apr 15, 2018 at 6:03 AM, Paul-Christian Bürkner < [email protected]> wrote:

There are situations in meta-analysis where autocorrelation makes sense for group-level effects.

So basically here is the list of correlation structures to be implemented (with already implemented stuff checked):

Group-level effects:

• unstructured
• compound symmetry
• heterogenous compound symmetry
• autoregressive of order 1 (AR1)

Residuals:

• unstructured
• compound symmetry
• heterogenous compound symmetry
• autoregressive of order 1 (AR1)

— You are receiving this because you authored the thread. Reply to this email directly, view it on GitHub https://github.com/paul-buerkner/brms/issues/403#issuecomment-381404989, or mute the thread https://github.com/notifications/unsubscribe-auth/AkNSOq0Dd0BE1BcPtCDURqRLqBsl5pmyks5to0UdgaJpZM4TT0iy .

--

Emily A. Butler

Associate Professor & Graduate Director Family Studies and Human Development College of Agriculture & Life Sciences University of Arizona Tucson, AZ, 85721-0033

[paul-buerkner]

Yes indeed. Instead of | use || in group-level effects.

[ebmtnprof]

thanks - that's what I thought.... :)

On Sun, Apr 15, 2018 at 10:47 AM, Paul-Christian Bürkner < [email protected]> wrote:

Yes indeed. Instead of | use || in group-level effects.

— You are receiving this because you authored the thread. Reply to this email directly, view it on GitHub https://github.com/paul-buerkner/brms/issues/403#issuecomment-381424501, or mute the thread https://github.com/notifications/unsubscribe-auth/AkNSOmQv1C4igapuyLdLN4CX7zXw3Ssoks5to4eZgaJpZM4TT0iy .

--

Emily A. Butler

Associate Professor & Graduate Director Family Studies and Human Development College of Agriculture & Life Sciences University of Arizona Tucson, AZ, 85721-0033

[belayb]

Hi,

Is there any update about ARMA models for family 'cumulative'?

[paul-buerkner]

Not yet. And it depends on the exact specification on whether this will ever come. As soon as we have AR(MA) for random effects this will work for cumulative models as well. However, if you think of ARMA models for residuals, we are limited by the fact that that we don't have the cumulative distribution function of the multivariate normal distribution in Stan, yet (see #323).

[paul-buerkner]

I think an appropriate syntax to specify correlation structures of group-level terms would be

y ~ x + (x | gr(group, cor = "<cor type>"))

If || instead of | is used, the cor argument will be overwritten and no correlations will be modeled.

[strengejacke]

either this, or if you don't want to "break" the formula syntax, add an additional argument for group correlation structures, which accepts a named vector or list, e.g.:

group_cor = c(group = "<cor type 1>", grouptwo = "<cor type 2>")

where names represent the name of the group-level term in the formula.

[paul-buerkner]

It won't break anything. The gr() function already exists and is added internally if not included explicitely by the user. We just add the cor argument.

Adding an additional argument to brm (or brmsformula) is not able to match the complexity of brms' multilevel syntax since the same grouping factor might have several associated group-level terms each having their own correlation structure.

[ebmtnprof]

so do these recent thoughts suggest you are working on implementing correlation structures???? I am close to needing to decide what I will cover in my class in the fall, so if these were going to be in place I would definitely include them.... Cheers, Emily

[paul-buerkner]

I can tell you more the week after next week I guess. Out of interest, do you want discuss those things primarily for residuals or for group-level effects (or both)? Of course, they are conceptually similar, but in brms they require a very different treatment (residuals via autocor and group-level effects via brms' multilevel syntax).

[ebmtnprof]

Sounds hopeful :) I am mostly interested in them for the residuals, because then it would give most of the functionality of gls, which is a real work-horse in my areas. But having the full suite of options for the group-level effects would also be nice.... Thanks for considering implementing some of this! Cheers, Emily

On Sat, Jun 9, 2018 at 2:17 AM, Paul-Christian Bürkner < [email protected]> wrote:

I can tell you more the week after next week I guess. Out of interest, do you want discuss those things primarily for residuals or for group-level effects (or both)? Of course, they are conceptually similar, but in brms they require a very different treatment (residuals via autocor and group-level effects via brms' multilevel syntax).

— You are receiving this because you authored the thread. Reply to this email directly, view it on GitHub https://github.com/paul-buerkner/brms/issues/403#issuecomment-395950454, or mute the thread https://github.com/notifications/unsubscribe-auth/AkNSOoMejZtgcg_d2j9LWQkQ56Jmp4UOks5t64SDgaJpZM4TT0iy .

--

Emily A. Butler

Associate Professor & Graduate Director Family Studies and Human Development College of Agriculture & Life Sciences University of Arizona Tucson, AZ, 85721-0033

[ebmtnprof]

hi again - Same time next year and I am again working on updating my course materials :) Any chance any further work has been done on correlation structures for residuals, other than cor_ar and cor_arma? I am still interested in the first 3 options on the checklist from last year (unstructured, compound symmetric and heterogenous compound symmetric). Of course, given my weak math abilities, it is possible those are more specific cases of cor_arma models? If so, just say so and I'll go figure it out :) Thanks for any updates! Cheers, Emily

[ghost]

Is there a link to your course content?

[paul-buerkner]

I don't think they are special cases of cor_arma. I can't promise anything, but I will look into it once more next week as to whether I can implement something soonish.

[ebmtnprof]

thanks for thinking about it! Compound symmetric and unstructured residuals would be the highest on the wish list. And many thanks again and again for your awesome contribution to the world!! brms is really fantastic. Cheers, Emily

[ebmtnprof]

oh, and sorry, I didn't see the question about the course content. I don't have a link, but I am attaching my syllabi - they are both under construction, but give the flavor of the classes. I do an "Intro" class (which isn't really very intro) and a multilevel modeling class. Cheers, Emily MLM_Syllabus2018.pdf 537A_Syllabus2019.pdf

[ebmtnprof]

PS - I am spending today wrapping my brain around the current multivariate modeling capability in brms (wow, you've added a lot of power!) and I may be able to model everything on my wish list that way, or at least be able to better articulate what is missing. So I'll post an update after I give it my best shot so you don't waste your time in the interim. Thanks, Cheers, Emily

[ebmtnprof]

Update: Ok, I did my best and am really impressed with the multivariate modeling power now available! For true multivariate models, I can't think of anything else I'd want. It made me realize that what is distinct about the situation I'm trying to model is that the observations are data over time (e.g., classic repeated-measures ANOVA type temporal data) and I would like to be able to structure the residuals either alone (e.g., no group-level parameters) or on top of group-level intercepts/slopes. The behavior of cor_ar for autocorrelation is exactly what I'm after, I just need more options (unstructured, compound symmetric, etc).

If doing it with autocor is a lot of work, would it be any easier to allow constraints on the residual structure using the multivariate approach (or maybe it's already possible?). I was able to use that approach to get unstructured residuals, but it doesn't exactly mimic gls since you get the population level-estimates for each time point, rather than a population slope over time (see attached data and syntax below). Thanks yet again for engaging with this.... Cheers, Emily

wide.xlsx

unstructured residuals

fit2 <- brm( mvbind(hr01, hr02, hr03, hr04, hr05) ~ 1 , data = wide, chains=2 )

[paul-buerkner]

I agree with you that the best place to incorporate compound symmetry and unstructured is via the autocor argument. If I find the time this week, I may mention to implement at least the compound symmetry structure. Fingers crossed :-)

[ebmtnprof]

thanks so much!! That would be great.... I'll stay tuned. Cheers, Emily

[paul-buerkner]

The compound symmetry structure should now be working via function cor_cosy in the github version of brms. I haven't tried it with real data yet so I would greatly appreciate your feedback on whether it returns sensible results also for non-artificial data.

[ebmtnprof]

awesome!! And yes, it is giving very sensible results with at least one real data example. I'll let you know if I run into trouble, but otherwise you can assume it's working.

While I have your attention, can I ask one confirmation question? It seems to me that the following syntax (see below for syntax and output) models heterogeneous residuals with respect to a between-person predictor (e.g. residuals either increasing or decreasing with scores on trait anxiety). 1) am I correct about that? and 2) is "sigma_intercept" the sigma at anxiety = 0 and "sigma_anxiety" the slope of the residuals on anxiety? Many thanks for your efforts! Cheers, Emily

fit1 <- brm(bf(hr ~ anxiety + time + (time|subject), sigma ~ anxiety), data=data, chains=2, control=list(max_treedepth=15) ) Group-Level Effects: ~subject (Number of levels: 187) Estimate Est.Error l-95% CI u-95% CI Eff.Sample Rhat sd(Intercept) 11.39 0.68 10.19 12.85 300 1.01 sd(time) 0.50 0.16 0.17 0.79 80 1.05 cor(Intercept,time) -0.61 0.16 -0.95 -0.30 211 1.03

Population-Level Effects: Estimate Est.Error l-95% CI u-95% CI Eff.Sample Rhat Intercept 73.05 4.31 65.20 82.23 147 1.03 sigma_Intercept 1.04 0.15 0.76 1.33 1659 1.00 anxiety 0.40 0.42 -0.47 1.16 147 1.03 time -1.09 0.09 -1.27 -0.92 2300 1.00 sigma_anxiety 0.02 0.01 -0.01 0.05 1989 1.00

[paul-buerkner]

This is slightly off-topic, but I am making an exception this time. Yes, I think your interpretation is correct. When combining predictions for sigma and cor_cosy we will also get heterogenous compound symmtry "for free", but this combination is not yet support. Fortunately, it just needs a small (I hope) refactor to make it possible

[ebmtnprof]

Thanks for the confirmation. Much appreciated.... And I'll stay tuned for the heterogenous compound symmetry. In the meantime, I have plenty of options now to meet my pedagogical needs. Cheers, Emily