[R-sig-ME] time*treatment vs time + time:treatment in RCTs

Reinhold Kliegl re|nho|d@k||eg| @end|ng |rom gm@||@com
Mon Aug 29 22:19:39 CEST 2022


Instead of combining the groups at pretest, you could also switch to indicator variables. 

dat_long$Subj <- factor(dat_long$ID)

mm <- model.matrix(~ 1 + time*group, data=dat_long)
t <- mm[,2]
t_x_g <- mm[,4]
res3 <- lmer(vo2 ~  1 + t + t_x_g + (1 | Subj),
             data = dat_long, REML=FALSE,
             control=lmerControl(calc.derivs = FALSE))

anova(res3, res2)

#Data: dat_long
#Models:
#res3: vo2 ~ 1 + t + t_x_g + (1 | Subj)
#res2: vo2 ~ 1 + time + group:time + (1 | Subj)
#     npar    AIC    BIC  logLik deviance  Chisq Df Pr(>Chisq)  
#res3    5 600.56 614.49 -295.28   590.56                       
#res2    6 599.81 616.54 -293.91   587.81 2.7446  1    0.09758 .
#---
#Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

> On 29. Aug 2022, at 21:25, Karl Ove Hufthammer <karl using huftis.org> wrote:
> 
> BTW, this is actually a rather annoying feature of the ways model formulas work in R. For a *randomised* longitudinal study, the population means for the randomisation groups are *identical* at baseline (due to the randomisation). So to properly *adjust* for any (random) group unbalances at baseline in the samples, one should fit a model *without* a ‘group’ effect at baseline. (I know it sounds strange to fit a model without group differences at baseline to *adjust* for any baseline differences, but if you think about it for a while, it makes sense …)
> 
> So instead of fitting
> 
> time + group + time:group + ...
> 
> (where ‘group’ represents the population differences at baseline)
> 
> one would naively *think* that one should fit a
> 
> time + time:group + ...
> 
> model. But R changes the *meaning* of the ‘time:group’ term in the second model, so one ends up fitting the exact same model as the first model (though with a different parametrisation), i.e., a model *not* adjusting for any sample differences at baseline.
> 
> The only way I’ve found to easily fit the proper model, is to create a factor of all ‘interaction(time, group)’ values and manually collapse the baseline groups to have the same level. That works fine, but it makes some tests *much* more complication, as you can no longer use R’s formula handling for simplifying models.
> 
> 
> Karl Ove Hufthammer
> 
> Karl Ove Hufthammer skreiv 29.08.2022 20:59:
>> Hmm, the automatic HTML-to-plaintext conversion didn’t work too well. Here’s a plaintext version:
>> 
>> No, you actually get equivalent results for your two models (which is really one model, just parametrised differently). The likelihood for the two models should be identical, and the P-value for testing whether there is an interaction should be identical. Here’s a simple simulation for data very similar to the ones you have:
>> 
>> library(lmerTest)
>> 
>> d_temp = expand.grid(group = c("PLA", "FUT"),
>>                      time = c("Baseline", "3month", "4month"))
>> n_varcombo = nrow(d_temp)
>> d_temp$exp = c(29, 30, 24.5, 28, 24, 27)
>> n_ind = 30
>> dat_long = d_temp[rep(1:n_varcombo, each = n_ind), ]
>> dat_long$ID = rep(1:n_ind, each = n_varcombo)
>> set.seed(6)
>> dat_long$ID_effect = rep(rnorm(n_ind, sd = 3), each = n_varcombo)
>> dat_long$vo2 = dat_long$exp + rnorm(n_varcombo * n_ind, sd = 3)
>> 
>> # Model without interaction
>> res0 <- lmer(vo2 ~ group + time + (1 | ID), data = dat_long)
>> 
>> # Two (equivalent) models with interaction
>> res1 <- lmer(vo2 ~  group * time + (1  | ID), data = dat_long)
>> summary(res1)
>> #> Fixed effects:
>> #>                     Estimate Std. Error      df t value Pr(>|t|)
>> #> (Intercept)          28.6258     0.5439 24.0000  52.635 < 2e-16 ***
>> #> groupFUT              1.2876     0.7691 24.0000   1.674 0.1071
>> #> time3month           -4.8032     0.7691 24.0000  -6.245 1.87e-06 ***
>> #> time4month           -4.8628     0.7691 24.0000  -6.322 1.55e-06 ***
>> #> groupFUT:time3month   2.7573     1.0877 24.0000   2.535 0.0182 *
>> #> groupFUT:time4month   1.5326     1.0877 24.0000   1.409 0.1717
>> #>
>> 
>> res2 <- lmer(vo2 ~  time + group:time + (1 | ID), data = dat_long)
>> summary(res2)
>> #> Fixed effects:
>> #>                       Estimate Std. Error      df t value Pr(>|t|)
>> #> (Intercept)            28.6258     0.5439 24.0000  52.635 < 2e-16 ***
>> #> time3month             -4.8032     0.7691 24.0000  -6.245 1.87e-06 ***
>> #> time4month             -4.8628     0.7691 24.0000  -6.322 1.55e-06 ***
>> #> timeBaseline:groupFUT   1.2876     0.7691 24.0000   1.674 0.10710
>> #> time3month:groupFUT     4.0449     0.7691 24.0000   5.259 2.16e-05 ***
>> #> time4month:groupFUT     2.8202     0.7691 24.0000   3.667 0.00122 **
>> 
>> # The models have the same log likelihood (and degrees of freedom)
>> logLik(res1)
>> #> 'log Lik.' -442.2284 (df=8)
>> logLik(res2)
>> #> 'log Lik.' -442.2284 (df=8)
>> 
>> # And the P-values for the interaction are exactly the same
>> anova(res0, res1)
>> #> refitting model(s) with ML (instead of REML)
>> #> Data: dat_long
>> #> Models:
>> #> res0: vo2 ~ group + time + (1 | ID)
>> #> res1: vo2 ~ group * time + (1 | ID)
>> #>      npar    AIC    BIC  logLik deviance  Chisq Df Pr(>Chisq)
>> #> res0    6 906.62 925.78 -447.31 894.62
>> #> res1    8 903.79 929.33 -443.89   887.79 6.8379  2 0.03275 *
>> #> ---
>> #> Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
>> 
>> anova(res0, res2)
>> #> refitting model(s) with ML (instead of REML)
>> #> Data: dat_long
>> #> Models:
>> #> res0: vo2 ~ group + time + (1 | ID)
>> #> res2: vo2 ~ time + group:time + (1 | ID)
>> #>      npar    AIC    BIC  logLik deviance  Chisq Df Pr(>Chisq)
>> #> res0    6 906.62 925.78 -447.31 894.62
>> #> res2    8 903.79 929.33 -443.89   887.79 6.8379  2 0.03275 *
>> #> ---
>> #> Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
>> 
>> 
>> As you can see, the two models are actually equivalent, and any inference based on the models should give the same results.
>> 
>> 
>> Karl Ove Hufthammer
>> 
>> 
>> Karl Ove Hufthammer skreiv 29.08.2022 20:56:
>>> No, you actually get equivalent results for your two models (which is
>>> really one model, just parametrised differently). The likelihood for the
>>> two models should be identical, and the P-value for testing whether
>>> there is an interaction should be identical. Here’s a simple simulation
>>> for data very similar to the ones you have:
>>> 
>>> |library(lmerTest)
>>> <http://localhost:28019/session/grave-cobra_reprex_preview.html?viewer_pane=1&capabilities=1&host=http%3A%2F%2F127.0.0.1%3A55034#cb1-13>d_temp 
>>> = expand.grid(group = c("PLA", "FUT"),
>>> <http://localhost:28019/session/grave-cobra_reprex_preview.html?viewer_pane=1&capabilities=1&host=http%3A%2F%2F127.0.0.1%3A55034#cb1-14> 
>>> time = c("Baseline", "3month", "4month"))
>>> <http://localhost:28019/session/grave-cobra_reprex_preview.html?viewer_pane=1&capabilities=1&host=http%3A%2F%2F127.0.0.1%3A55034#cb1-15>n_varcombo 
>>> = nrow(d_temp)
>>> <http://localhost:28019/session/grave-cobra_reprex_preview.html?viewer_pane=1&capabilities=1&host=http%3A%2F%2F127.0.0.1%3A55034#cb1-16>d_temp$exp 
>>> = c(29, 30, 24.5, 28, 24, 27)
>>> <http://localhost:28019/session/grave-cobra_reprex_preview.html?viewer_pane=1&capabilities=1&host=http%3A%2F%2F127.0.0.1%3A55034#cb1-17>n_ind 
>>> = 30
>>> <http://localhost:28019/session/grave-cobra_reprex_preview.html?viewer_pane=1&capabilities=1&host=http%3A%2F%2F127.0.0.1%3A55034#cb1-18>dat_long 
>>> = d_temp[rep(1:n_varcombo, each = n_ind), ]
>>> <http://localhost:28019/session/grave-cobra_reprex_preview.html?viewer_pane=1&capabilities=1&host=http%3A%2F%2F127.0.0.1%3A55034#cb1-19>dat_long$ID 
>>> = rep(1:n_ind, each = n_varcombo)
>>> <http://localhost:28019/session/grave-cobra_reprex_preview.html?viewer_pane=1&capabilities=1&host=http%3A%2F%2F127.0.0.1%3A55034#cb1-20>set.seed(6) 
>>> <http://localhost:28019/session/grave-cobra_reprex_preview.html?viewer_pane=1&capabilities=1&host=http%3A%2F%2F127.0.0.1%3A55034#cb1-21>dat_long$ID_effect 
>>> = rep(rnorm(n_ind, sd = 3), each = n_varcombo)
>>> <http://localhost:28019/session/grave-cobra_reprex_preview.html?viewer_pane=1&capabilities=1&host=http%3A%2F%2F127.0.0.1%3A55034#cb1-22>dat_long$vo2 
>>> = dat_long$exp + rnorm(n_varcombo * n_ind, sd = 3)
>>> <http://localhost:28019/session/grave-cobra_reprex_preview.html?viewer_pane=1&capabilities=1&host=http%3A%2F%2F127.0.0.1%3A55034#cb1-23> 
>>> # Model without interaction
>>> <http://localhost:28019/session/grave-cobra_reprex_preview.html?viewer_pane=1&capabilities=1&host=http%3A%2F%2F127.0.0.1%3A55034#cb1-25>res0 
>>> <- lmer(vo2 ~ group + time + (1 | ID), data = dat_long)
>>> <http://localhost:28019/session/grave-cobra_reprex_preview.html?viewer_pane=1&capabilities=1&host=http%3A%2F%2F127.0.0.1%3A55034#cb1-26> 
>>> # Two (equivalent) models with interaction
>>> <http://localhost:28019/session/grave-cobra_reprex_preview.html?viewer_pane=1&capabilities=1&host=http%3A%2F%2F127.0.0.1%3A55034#cb1-28>res1 
>>> <- lmer(vo2 ~ group * time + (1 | ID), data = dat_long)
>>> <http://localhost:28019/session/grave-cobra_reprex_preview.html?viewer_pane=1&capabilities=1&host=http%3A%2F%2F127.0.0.1%3A55034#cb1-29>summary(res1) 
>>> #> Fixed effects:
>>> <http://localhost:28019/session/grave-cobra_reprex_preview.html?viewer_pane=1&capabilities=1&host=http%3A%2F%2F127.0.0.1%3A55034#cb1-48>#> 
>>> Estimate Std. Error df t value Pr(>|t|)
>>> <http://localhost:28019/session/grave-cobra_reprex_preview.html?viewer_pane=1&capabilities=1&host=http%3A%2F%2F127.0.0.1%3A55034#cb1-49>#> 
>>> (Intercept) 28.6258 0.5439 24.0000 52.635 < 2e-16 ***
>>> <http://localhost:28019/session/grave-cobra_reprex_preview.html?viewer_pane=1&capabilities=1&host=http%3A%2F%2F127.0.0.1%3A55034#cb1-50>#> 
>>> groupFUT 1.2876 0.7691 24.0000 1.674 0.1071
>>> <http://localhost:28019/session/grave-cobra_reprex_preview.html?viewer_pane=1&capabilities=1&host=http%3A%2F%2F127.0.0.1%3A55034#cb1-51>#> 
>>> time3month -4.8032 0.7691 24.0000 -6.245 1.87e-06 ***
>>> <http://localhost:28019/session/grave-cobra_reprex_preview.html?viewer_pane=1&capabilities=1&host=http%3A%2F%2F127.0.0.1%3A55034#cb1-52>#> 
>>> time4month -4.8628 0.7691 24.0000 -6.322 1.55e-06 ***
>>> <http://localhost:28019/session/grave-cobra_reprex_preview.html?viewer_pane=1&capabilities=1&host=http%3A%2F%2F127.0.0.1%3A55034#cb1-53>#> 
>>> groupFUT:time3month 2.7573 1.0877 24.0000 2.535 0.0182 *
>>> <http://localhost:28019/session/grave-cobra_reprex_preview.html?viewer_pane=1&capabilities=1&host=http%3A%2F%2F127.0.0.1%3A55034#cb1-54>#> 
>>> groupFUT:time4month 1.5326 1.0877 24.0000 1.409 0.1717 #>
>>> <http://localhost:28019/session/grave-cobra_reprex_preview.html?viewer_pane=1&capabilities=1&host=http%3A%2F%2F127.0.0.1%3A55034#cb1-65>res2 
>>> <- lmer(vo2 ~ time + group:time + (1 | ID), data = dat_long)
>>> <http://localhost:28019/session/grave-cobra_reprex_preview.html?viewer_pane=1&capabilities=1&host=http%3A%2F%2F127.0.0.1%3A55034#cb1-66>summary(res2) 
>>> #> Fixed effects:
>>> <http://localhost:28019/session/grave-cobra_reprex_preview.html?viewer_pane=1&capabilities=1&host=http%3A%2F%2F127.0.0.1%3A55034#cb1-85>#> 
>>> Estimate Std. Error df t value Pr(>|t|)
>>> <http://localhost:28019/session/grave-cobra_reprex_preview.html?viewer_pane=1&capabilities=1&host=http%3A%2F%2F127.0.0.1%3A55034#cb1-86>#> 
>>> (Intercept) 28.6258 0.5439 24.0000 52.635 < 2e-16 ***
>>> <http://localhost:28019/session/grave-cobra_reprex_preview.html?viewer_pane=1&capabilities=1&host=http%3A%2F%2F127.0.0.1%3A55034#cb1-87>#> 
>>> time3month -4.8032 0.7691 24.0000 -6.245 1.87e-06 ***
>>> <http://localhost:28019/session/grave-cobra_reprex_preview.html?viewer_pane=1&capabilities=1&host=http%3A%2F%2F127.0.0.1%3A55034#cb1-88>#> 
>>> time4month -4.8628 0.7691 24.0000 -6.322 1.55e-06 ***
>>> <http://localhost:28019/session/grave-cobra_reprex_preview.html?viewer_pane=1&capabilities=1&host=http%3A%2F%2F127.0.0.1%3A55034#cb1-89>#> 
>>> timeBaseline:groupFUT 1.2876 0.7691 24.0000 1.674 0.10710
>>> <http://localhost:28019/session/grave-cobra_reprex_preview.html?viewer_pane=1&capabilities=1&host=http%3A%2F%2F127.0.0.1%3A55034#cb1-90>#> 
>>> time3month:groupFUT 4.0449 0.7691 24.0000 5.259 2.16e-05 ***
>>> <http://localhost:28019/session/grave-cobra_reprex_preview.html?viewer_pane=1&capabilities=1&host=http%3A%2F%2F127.0.0.1%3A55034#cb1-91>#> 
>>> time4month:groupFUT 2.8202 0.7691 24.0000 3.667 0.00122 ** # The models
>>> have the same log likelihood (and degrees of freedom)
>>> <http://localhost:28019/session/grave-cobra_reprex_preview.html?viewer_pane=1&capabilities=1&host=http%3A%2F%2F127.0.0.1%3A55034#cb1-104>logLik(res1) 
>>> <http://localhost:28019/session/grave-cobra_reprex_preview.html?viewer_pane=1&capabilities=1&host=http%3A%2F%2F127.0.0.1%3A55034#cb1-105>#> 
>>> 'log Lik.' -442.2284 (df=8)
>>> <http://localhost:28019/session/grave-cobra_reprex_preview.html?viewer_pane=1&capabilities=1&host=http%3A%2F%2F127.0.0.1%3A55034#cb1-106>logLik(res2) 
>>> <http://localhost:28019/session/grave-cobra_reprex_preview.html?viewer_pane=1&capabilities=1&host=http%3A%2F%2F127.0.0.1%3A55034#cb1-107>#> 
>>> 'log Lik.' -442.2284 (df=8)
>>> <http://localhost:28019/session/grave-cobra_reprex_preview.html?viewer_pane=1&capabilities=1&host=http%3A%2F%2F127.0.0.1%3A55034#cb1-108> 
>>> # And the P-values for the interaction are exactly the same
>>> <http://localhost:28019/session/grave-cobra_reprex_preview.html?viewer_pane=1&capabilities=1&host=http%3A%2F%2F127.0.0.1%3A55034#cb1-110>anova(res0, 
>>> res1)
>>> <http://localhost:28019/session/grave-cobra_reprex_preview.html?viewer_pane=1&capabilities=1&host=http%3A%2F%2F127.0.0.1%3A55034#cb1-111>#> 
>>> refitting model(s) with ML (instead of REML)
>>> <http://localhost:28019/session/grave-cobra_reprex_preview.html?viewer_pane=1&capabilities=1&host=http%3A%2F%2F127.0.0.1%3A55034#cb1-112>#> 
>>> Data: dat_long
>>> <http://localhost:28019/session/grave-cobra_reprex_preview.html?viewer_pane=1&capabilities=1&host=http%3A%2F%2F127.0.0.1%3A55034#cb1-113>#> 
>>> Models:
>>> <http://localhost:28019/session/grave-cobra_reprex_preview.html?viewer_pane=1&capabilities=1&host=http%3A%2F%2F127.0.0.1%3A55034#cb1-114>#> 
>>> res0: vo2 ~ group + time + (1 | ID)
>>> <http://localhost:28019/session/grave-cobra_reprex_preview.html?viewer_pane=1&capabilities=1&host=http%3A%2F%2F127.0.0.1%3A55034#cb1-115>#> 
>>> res1: vo2 ~ group * time + (1 | ID)
>>> <http://localhost:28019/session/grave-cobra_reprex_preview.html?viewer_pane=1&capabilities=1&host=http%3A%2F%2F127.0.0.1%3A55034#cb1-116>#> 
>>> npar AIC BIC logLik deviance Chisq Df Pr(>Chisq)
>>> <http://localhost:28019/session/grave-cobra_reprex_preview.html?viewer_pane=1&capabilities=1&host=http%3A%2F%2F127.0.0.1%3A55034#cb1-117>#> 
>>> res0 6 906.62 925.78 -447.31 894.62 #> res1 8 903.79 929.33 -443.89
>>> 887.79 6.8379 2 0.03275 *
>>> <http://localhost:28019/session/grave-cobra_reprex_preview.html?viewer_pane=1&capabilities=1&host=http%3A%2F%2F127.0.0.1%3A55034#cb1-119>#> 
>>> ---
>>> <http://localhost:28019/session/grave-cobra_reprex_preview.html?viewer_pane=1&capabilities=1&host=http%3A%2F%2F127.0.0.1%3A55034#cb1-120>#> 
>>> Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
>>> <http://localhost:28019/session/grave-cobra_reprex_preview.html?viewer_pane=1&capabilities=1&host=http%3A%2F%2F127.0.0.1%3A55034#cb1-121>anova(res0, 
>>> res2)
>>> <http://localhost:28019/session/grave-cobra_reprex_preview.html?viewer_pane=1&capabilities=1&host=http%3A%2F%2F127.0.0.1%3A55034#cb1-122>#> 
>>> refitting model(s) with ML (instead of REML)
>>> <http://localhost:28019/session/grave-cobra_reprex_preview.html?viewer_pane=1&capabilities=1&host=http%3A%2F%2F127.0.0.1%3A55034#cb1-123>#> 
>>> Data: dat_long
>>> <http://localhost:28019/session/grave-cobra_reprex_preview.html?viewer_pane=1&capabilities=1&host=http%3A%2F%2F127.0.0.1%3A55034#cb1-124>#> 
>>> Models:
>>> <http://localhost:28019/session/grave-cobra_reprex_preview.html?viewer_pane=1&capabilities=1&host=http%3A%2F%2F127.0.0.1%3A55034#cb1-125>#> 
>>> res0: vo2 ~ group + time + (1 | ID)
>>> <http://localhost:28019/session/grave-cobra_reprex_preview.html?viewer_pane=1&capabilities=1&host=http%3A%2F%2F127.0.0.1%3A55034#cb1-126>#> 
>>> res2: vo2 ~ time + group:time + (1 | ID)
>>> <http://localhost:28019/session/grave-cobra_reprex_preview.html?viewer_pane=1&capabilities=1&host=http%3A%2F%2F127.0.0.1%3A55034#cb1-127>#> 
>>> npar AIC BIC logLik deviance Chisq Df Pr(>Chisq)
>>> <http://localhost:28019/session/grave-cobra_reprex_preview.html?viewer_pane=1&capabilities=1&host=http%3A%2F%2F127.0.0.1%3A55034#cb1-128>#> 
>>> res0 6 906.62 925.78 -447.31 894.62 #> res2 8 903.79 929.33 -443.89
>>> 887.79 6.8379 2 0.03275 *
>>> <http://localhost:28019/session/grave-cobra_reprex_preview.html?viewer_pane=1&capabilities=1&host=http%3A%2F%2F127.0.0.1%3A55034#cb1-130>#> 
>>> ---
>>> <http://localhost:28019/session/grave-cobra_reprex_preview.html?viewer_pane=1&capabilities=1&host=http%3A%2F%2F127.0.0.1%3A55034#cb1-131>#> 
>>> Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1|
>>> 
>>> 
>>> As you can see, the two models are actually equivalent, and any
>>> inference based on the models should give the same results.
>>> 
>>> 
>>> Karl Ove Hufthammer
>>> 
>>> Jorge Teixeira skreiv 29.08.2022 20:20:
>>>> Thank you all for the replies. Still processing them...
>>>> 
>>>> Indeed, Wolfgang, I was mainly thinking of time as a factor. Although, I
>>>> welcome comments as if it was numeric as well.
>>>> 
>>>> Your reply is surprising to me, because in my data I get different results.
>>>> The ES and p-values are very different regarding the interactions at 3 and
>>>> 4-months, which are the relevant data to me. My df has 3 time points.
>>>> 
>>>>    res1 <- lmer(vo2 ~  group*time + ( 1  | ID  ), data = dat_long )
>>>> 
>>>>      res2 <- lmer(vo2 ~  time + group:time + ( 1 | ID  ), data =  dat_long )
>>>> 
>>>> 
>>>> *res1:*
>>>> Fixed effects:
>>>>                          Estimate Std. Error      df t value Pr(>|t|)
>>>> (Intercept)             29.0705     0.9998 61.4510  29.076 < 2e-16 ***
>>>> groupFUT              1.0395     1.4140 61.4510   0.735 0.465036
>>>> time3month              -4.4917     1.0918 64.1740  -4.114 0.000113 ***
>>>> time4month              -5.0305     1.0622 63.8295  -4.736 1.26e-05 ***
>>>> *groupFUT:time3month *  2.5467     1.4396 61.8093   1.769 0.081822 .
>>>> *groupFUT:time4month*   1.7643     1.4424 61.8409   1.223 0.225909
>>>> 
>>>> 
>>>> *res2:*
>>>> Fixed effects:
>>>>                          Estimate Std. Error      df t value Pr(>|t|)
>>>> (Intercept)             29.0705     0.9998 61.4510  29.076 < 2e-16 ***
>>>> time3month              -4.4917     1.0918 64.1740  -4.114 0.000113 ***
>>>> time4month              -5.0305     1.0622 63.8295  -4.736 1.26e-05 ***
>>>> time0month:groupFUT   1.0395     1.4140 61.4510   0.735 0.465036
>>>> *time3month:groupFUT*   3.5862     1.5402 73.4895   2.328 0.022643 *
>>>> *time4month:groupFUT*  2.8038     1.5428 73.7427   1.817 0.073226 .
>>>> 
>>>> 
>>>> 
>>>> Viechtbauer, Wolfgang (NP)<wolfgang.viechtbauer using maastrichtuniversity.nl>
>>>> escreveu no dia segunda, 29/08/2022 à(s) 16:11:
>>>> 
>>>>> I strongly suspect that 'time' is treated as a factor in the examples
>>>>> Jorge is referring to. In this case, the two formulations are just
>>>>> different parameterizations of the same model. We can use the 'Orthodont'
>>>>> data to illustrate this. Think of 'age' as the time variable (as a
>>>>> four-level factor) and 'Sex' as the treatment variable (as a two-level
>>>>> factor). In fact, I will throw in a third parameterization, which I think
>>>>> is even more intuitive.
>>>>> 
>>>>> library(lme4)
>>>>> 
>>>>> data("Orthodont", package="nlme")
>>>>> 
>>>>> Orthodont$age <- factor(Orthodont$age)
>>>>> 
>>>>> res1 <- lmer(distance ~ age*Sex + (1 | Subject), data=Orthodont)
>>>>> summary(res1)
>>>>> 
>>>>> res2 <- lmer(distance ~ age + age:Sex + (1 | Subject), data=Orthodont)
>>>>> summary(res2)
>>>>> 
>>>>> res3 <- lmer(distance ~ 0 + age + age:Sex + (1 | Subject), data=Orthodont)
>>>>> summary(res3)
>>>>> 
>>>>> logLik(res1)
>>>>> logLik(res2)
>>>>> logLik(res3)
>>>>> 
>>>>> The fit is identical.
>>>>> 
>>>>> In 'res3', we get the estimated intercepts (means) of the reference group
>>>>> (in this case for 'Male') at all 4 timepoints and the age:Sex coefficients
>>>>> are the difference between the Female and Male groups at those 4 timepoints.
>>>>> 
>>>>> Since these are just all different parameterizations of the same model,
>>>>> there is no reasons for preferring one over the other.
>>>>> 
>>>>> One has to be careful though when using anova() on those models, esp. with
>>>>> respect to the age:Sex test. In anova(res1), the test examines if the
>>>>> difference between males and females is the same at all 4 timepoints, while
>>>>> in anova(res2) and anova(res3) the test examines if the difference is 0 at
>>>>> all 4 timepoints. However, one could get either test out of all three
>>>>> parameterizations, by forming appropriate contrasts. So again, no reason to
>>>>> prefer one over the other (except maybe convenience depending on what one
>>>>> would like to test).
>>>>> 
>>>>> Best,
>>>>> Wolfgang
>>>>> 
>>>>>> -----Original Message-----
>>>>>> From: R-sig-mixed-models [mailto:r-sig-mixed-models-bounces using r-project.org]
>>>>> On
>>>>>> Behalf Of Douglas Bates
>>>>>> Sent: Monday, 29 August, 2022 16:14
>>>>>> To: Phillip Alday
>>>>>> Cc: R-mixed models mailing list; Jorge Teixeira
>>>>>> Subject: Re: [R-sig-ME] time*treatment vs time + time:treatment in RCTs
>>>>>> 
>>>>>> M2 is an appropriate model if time corresponds to "time on treatment" or
>>>>> in
>>>>>> general if the covariate over which the measurements are repeated has a
>>>>>> scale where 0 is meaningful.  I think of it as the "zero dose" model
>>>>>> because zero dose of treatment 1 is the same as zero dose of treatment 2
>>>>> is
>>>>>> the same as zero dose of the placebo.  Similarly zero time on treatment is
>>>>>> the same for any of the treatments or the placebo.
>>>>>> 
>>>>>> In those cases we would not expect a main effect for treatment because
>>>>> that
>>>>>> corresponds to systematic differences before the study begins (or at zero
>>>>>> dose), but we would expect an interaction of time (or dose) with
>>>>> treatment.
>>>>>> On Mon, Aug 29, 2022 at 8:28 AM Phillip Alday<me using phillipalday.com>
>>>>> wrote:
>>>>>>> On 8/29/22 05:53, Jorge Teixeira wrote:
>>>>>>>> Hi. In medicine's RCTs, with 3 or more time-points, whenever LMMs are
>>>>>>> used
>>>>>>>> and the code is available, a variation of  y ~ time*treatment + (1 |
>>>>> ID)
>>>>>>>> *(M1)* is always used (from what I have seen).
>>>>>>>> 
>>>>>>>> Recently I came across the model  time + time:treatment + (1 | ID)*
>>>>> (M2)*
>>>>>>>> in Solomun Kurz's blog and in the book of Galecki (LMMs using R).
>>>>>>>> 
>>>>>>>> Questions:
>>>>>>>> *1)* Are there any modelling reasons for M2 to be less used in
>>>>> medicine's
>>>>>>>> RCTs?
>>>>>>> It depends a bit on what `y` is: change from baseline or the 'raw'
>>>>>>> measure. If it's the raw measure, then (M2) doesn't include a
>>>>>>> description of differences at baseline between the groups.
>>>>>>> 
>>>>>>> Perhaps most importantly though: (M2) violates the principle of
>>>>>>> marginality discussed e.g. in Venables' Exegeses on Linear Models
>>>>>>> (https://www.stats.ox.ac.uk/pub/MASS3/Exegeses.pdf)
>>>>>>> 
>>>>>>>> *2)* Can anyone explain, in layman terms, what is the estimand in M2?
>>>>> I
>>>>>>>> still struggle to understand what model is really measuring.
>>>>>>> Approximately the same thing as M1, except that the "overall" effect of
>>>>>>> treatment is assumed to be zero. "Overall" is a bit vague because it
>>>>>>> depends on the contrast coding used for time and treatment.
>>>>>>> 
>>>>>>> You can see this for yourself. M1 can also be written as:
>>>>>>> 
>>>>>>> y ~ time + time:treatment + treatment + (1|ID).
>>>>>>> 
>>>>>>> If you force the coefficient on treatment to be zero, then you have M2.
>>>>>>> 
>>>>>>>> *3)* On a general basis, in a RCT with 3 time points (baseline,
>>>>> 3-month
>>>>>>> and
>>>>>>>> 4-month), would you tend to gravitate more towards model 1 or 2?
>>>>>>> Definitely (1).
>>>>>>> 
>>>>>>> PS: When referencing a blog entry, please provide a link to it. :)
>>>>>>> 
>>>>>>>> Thank you
>>>>>>>> Jorge
>>>>     [[alternative HTML version deleted]]
>>>> 
>>>> _______________________________________________
>>>> R-sig-mixed-models using r-project.org  mailing list
>>>> https://stat.ethz.ch/mailman/listinfo/r-sig-mixed-models
>>>> 
>> 
> 
> -- 
> Karl Ove Hufthammer
> 
> _______________________________________________
> R-sig-mixed-models using r-project.org mailing list
> https://stat.ethz.ch/mailman/listinfo/r-sig-mixed-models


	[[alternative HTML version deleted]]



More information about the R-sig-mixed-models mailing list