[R-meta] Time as indicator vs time as meaning

[Stefanou Revesz]

Sure, but in the lme(), "correlation=" has to do with the structure of
V matrix (i.e., e_ij), not random-effects, no?

Say, I had used the categorical "time" (measurement occasions) with
struct = "HAR", then would it have made sense to add a "UN" on top of
that?

yi ~ time, random = list(~ time | study, ~ time | study), struct =
c("HAR", "UN")

Best of all,
Stefanou

On Tue, Oct 12, 2021 at 10:34 AM Viechtbauer, Wolfgang (SP)
<wolfgang.viechtbauer using maastrichtuniversity.nl> wrote:
>
> Yes, I discovered this bug while working with your data. So thank you - you brought something to my attention that requires fixing.
>
> And here we go with the modeling questions :)
>
> The two structures are doing rather different things. One is modeling differences in intercepts and slopes of the regression lines, the other accounting for autocorrelation in the data (to be precise, in the deviations around the study-specific regression lines).
>
> If you look at the literature on longitudinal data analysis, you will find that things like this are done commonly also with raw data. For example, I quite frequently work with data collected in diary studies (i.e., ecological momentary assessment / experience sampling studies). A model like:
>
> lme(outcome ~ between.x + within.x, random = ~ within.x | subject/day,
>     correlation=corCAR1(form = ~ hour | subject/day)
>
> is quite standard there. So I am applying the same idea to the meta-analytic context.
>
> Best,
> Wolfgang
>
> >-----Original Message-----
> >From: Stefanou Revesz [mailto:stefanourevesz using gmail.com]
> >Sent: Tuesday, 12 October, 2021 17:21
> >To: Viechtbauer, Wolfgang (SP)
> >Cc: R meta
> >Subject: Re: [R-meta] Time as indicator vs time as meaning
> >
> >Dear Wolfgang,
> >
> >Thank you very much! In your answer, you also demystified the problem
> >of why model2 crashes, if I don't round the "time_wthn" up to 8
> >decimal places (that had me thinking the whole day yesterday).
> >
> >I'm all good. But I wonder what made you add two
> >differently-structured ("GEN" and "CAR") set of random-effects to the
> >**same ID variable** (i.e., study)?
> >
> >In other words, when such a strategy is warranted?
> >
> >Best of all,
> >Stefanou
> >
> >On Tue, Oct 12, 2021 at 8:05 AM Viechtbauer, Wolfgang (SP)
> ><wolfgang.viechtbauer using maastrichtuniversity.nl> wrote:
> >>
> >> While Stefanou contacted me off-list and also shared the data with me, I want
> >to bring this back here.
> >>
> >> So, now that I have seen the data, I see that the "time_wthn" and "time_btw"
> >variables were computed based on the measurenent occasion counter (time), not the
> >time variable itself (time_wk). That mixes different things and of course then
> >the results are not the same. So, using your data, do the following:
> >>
> >> dat$time_btw <- ave(dat$time_wk, dat$study, FUN=mean)
> >> dat$time_wthn <- dat$time_wk - dat$time_btw
> >> dat$time_wthn <- round(dat$time_wthn, 8)
> >>
> >> and then run:
> >>
> >> model1 <- rma.mv(yi ~ time_btw + time_wthn, vi,
> >>                  random = list(~ time_wthn | study, ~time_wk | study),
> >>                  struct = c("GEN","CAR"), data=dat)
> >>
> >> model2 <- rma.mv(yi ~ time_btw + time_wthn, vi,
> >>                  random = list(~ time_wthn | study, ~time_wthn | study),
> >>                  struct = c("GEN","CAR"), data=dat)
> >>
> >> model1
> >> model2
> >>
> >> which yields:
> >>
> >> Variance Components:
> >>
> >> outer factor: study      (nlvls = 49)
> >> inner term:   ~time_wthn (nlvls = 44)
> >>
> >>             estim    sqrt  fixed  rho:  intr    tm_w
> >> intrcpt    0.4685  0.6845     no           -  1.0000
> >> time_wthn  0.0229  0.1513     no          no       -
> >>
> >> outer factor: study   (nlvls = 49)
> >> inner factor: time_wk (nlvls = 12)
> >>
> >>             estim    sqrt  fixed
> >> gamma^2    0.0603  0.2456     no
> >> phi        0.0295             no
> >>
> >> for model1 and the same for model2. So all is as expected.
> >>
> >> profile(model1) also indicates peaks at all the estimates. And yes, that
> >estimate of rho drifts towards 1. Not nice when that happens, but that's where
> >the peak of the likelihood surface is.
> >>
> >> I also see in the data now that there are at times multiple effects at the same
> >time point. So one could even extend the model above further, for example with:
> >>
> >> dat$obs <- 1:nrow(dat)
> >>
> >> model3 <- rma.mv(yi ~ time_btw + time_wthn, vi,
> >>           random = list(~ time_wthn | study, ~time_wk | study, ~ 1 | obs),
> >>           struct = c("GEN","CAR"), data=dat)
> >> model3
> >>
> >> This further improves the fit:
> >>
> >> anova(model1, model3)
> >>
> >> These are some suggestions - I have not looked more deeply into these data. And
> >my comment from before stands: You probably need a proper V matrix, not just the
> >sampling variances (vi) if these effects are measured in the same people within
> >studies.
> >>
> >> One final thing: The line
> >>
> >> dat$time_wthn <- round(dat$time_wthn, 8)
> >>
> >> is a bit weird and seems unncessary. Well, I just got screwed by floating point
> >issues. There is a part in rma.mv() where I figure out the *unique* time point
> >values for struct="CAR" and use unique() for that. In essence, this is what
> >happens:
> >>
> >> unique(c(1/2, sqrt(1/2)*sqrt(1/2)))
> >>
> >> Argh. Rounding time_wthn a bit solves that issue here. Of course, that needs to
> >be fixed eventually in rma.mv() itself.
> >>
> >> Best,
> >> Wolfgang