[R-meta] Time as indicator vs time as meaning
Stefanou Revesz
@te|@noureve@z @end|ng |rom gm@||@com
Tue Oct 12 17:20:33 CEST 2021
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
>
> >-----Original Message-----
> >From: Stefanou Revesz [mailto:stefanourevesz using gmail.com]
> >Sent: Sunday, 10 October, 2021 23:29
> >To: Viechtbauer, Wolfgang (SP)
> >Cc: R meta
> >Subject: Re: [R-meta] Time as indicator vs time as meaning
> >
> >Resending Model 2:
> >
> >#***********************
> >MODEL 2:
> >#***********************
> >rma.mv(yi ~ time_btw + time_wthn, vi,
> > random =
> > list(~ time_wthn | study, ~time_wthn | study), struct
> >= c("GEN","CAR"),data=data)
> >
> >Variance Components:
> >
> >outer factor: study (nlvls = 49)
> >inner term: ~time_wthn (nlvls = 9)
> >
> > estim sqrt fixed rho: intr tm_w
> >intrcpt 0.4034 0.6351 no - no
> >time_wthn 0.1141 0.3377 no 1.0000 -
> >
> >outer factor: study (nlvls = 49)
> >inner factor: time_wthn (nlvls = 9)
> >
> > estim sqrt fixed
> >gamma^2 0.0770 0.2775 no
> >phi 0.0000 no
> >
> >Model Results:
> >
> > estimate se zval pval ci.lb ci.ub
> >intrcpt 0.1155 0.1861 0.6204 0.5350 -0.2493 0.4802
> >time_btw 0.3184 0.1815 1.7540 0.0794 -0.0374 0.6741 .
> >time_wthn 0.3247 0.0658 4.9344 <.0001 0.1957 0.4536 ***
> >
> >On Sun, Oct 10, 2021 at 4:24 PM Viechtbauer, Wolfgang (SP)
> ><wolfgang.viechtbauer using maastrichtuniversity.nl> wrote:
> >>
> >> >-----Original Message-----
> >> >From: Stefanou Revesz [mailto:stefanourevesz using gmail.com]
> >> >Sent: Sunday, 10 October, 2021 23:03
> >> >To: Viechtbauer, Wolfgang (SP)
> >> >Cc: R meta
> >> >Subject: Re: [R-meta] Time as indicator vs time as meaning
> >> >
> >> >ATTACHMENT(S) REMOVED: image.png
> >> >
> >> >Thanks! Reporting back the results for two models (note: "time_wk" =
> >> >"time_meaning_wks").
> >> >In model 1, I used "time_wk" to the left of | for "CAR".
> >> >In model 2, I used "time_whtn" to the left of | for "CAR".
> >> >In both models, rho is estimated to be 1. But the likelihood profile
> >> >for rho seems ok for both models (attached).
> >> >
> >> >However, in model 1, phi is estimated to be 0.1170. In model 2, phi is
> >> >estimated to be 0. In both models, the likelihood profiles for phi
> >> >seem ok.
> >> >Which phi, then, seems appropriate?
> >> >
> >> >Thank you,
> >> >Stefanou
> >> >
> >> >#******************************
> >> >MODEL 1:
> >> >#******************************
> >> >rma.mv(yi ~ time_btw + time_wthn, vi,
> >> > random =
> >> > list(~ time_wthn | study, ~time_wk | study), struct =
> >> >c("GEN","CAR"),data=data)
> >> >
> >> >Variance Components:
> >> >
> >> >outer factor: study (nlvls = 49)
> >> >inner term: ~time_wk (nlvls = 12)
> >> >
> >> > estim sqrt fixed rho: intr tm_w
> >> >intrcpt 0.1306 0.3614 no - no
> >> >time_wk 0.0236 0.1536 no 1.0000 -
> >> >
> >> >outer factor: study (nlvls = 49)
> >> >inner factor: time_wk (nlvls = 12)
> >> >
> >> > estim sqrt fixed
> >> >gamma^2 0.0703 0.2652 no
> >> >phi 0.0544 no
> >> >
> >> >Model Results:
> >> >
> >> > estimate se zval pval ci.lb ci.ub
> >> >intrcpt 0.0049 0.1816 0.0272 0.9783 -0.3510 0.3608
> >> >time_btw 0.3399 0.2164 1.5708 0.1162 -0.0842 0.7641
> >> >time_wthn 0.2837 0.0708 4.0064 <.0001 0.1449 0.4224 ***
> >>
> >> Something doesn't match up here. The first part of the 'random' formula says '~
> >time_wthn | study' but the output says that the inner term is ~time_wk.
> >>
> >> >#******************************
> >> >MODEL 2:
> >> >#******************************
> >> >rma.mv(yi ~ time_btw + time_wthn, vi,
> >> > random =
> >> > list(~ time_wthn | study, ~time_wthn | study), struct
> >> >= c("GEN","CAR"),data=data)
> >> >Variance Components:
> >> >
> >> >outer factor: study (nlvls = 49)
> >> >inner term: ~time_wthn (nlvls = 9)
> >> >
> >> > estim sqrt fixed rho: intr tm_w
> >> >intrcpt 0.4034 0.6351 no - no
> >> >time_wthn 0.1141 0.3377 no 1.0000 -
> >> >
> >> >outer factor: study (nlvls = 49)
> >> >inner factor: time_wthn (nlvls = 9)
> >> >
> >> > estim sqrt fixed
> >> >gamma^2 0.0770 0.2775 no
> >> >phi 0.0000 no
> >> >
> >> >Test for Residual Heterogeneity:
> >> >QE(df = 402) = 1450.3879, p-val < .0001
> >> >
> >> >Test of Moderators (coefficients 2:3):
> >> >QM(df = 2) = 24.4255, p-val < .0001
> >> >
> >> >Model Results:
> >> >
> >> > estimate se zval pval ci.lb ci.ub
> >> >intrcpt 0.1155 0.1861 0.6204 0.5350 -0.2493 0.4802
> >> >time_btw 0.3184 0.1815 1.7540 0.0794 -0.0374 0.6741 .
> >> >time_wthn 0.3247 0.0658 4.9344 <.0001 0.1957 0.4536 ***
More information about the R-sig-meta-analysis
mailing list