[R] model specification using lme

Thierry Onkelinx thierry.onkelinx at inbo.be
Mon May 30 10:40:14 CEST 2016

```Dear Hanna,

None of the models are correct is you want the same random intercept for
the different methods but different random slope per method.

You can random = ~ 1 + time:method | individual

The easiest way to get alpha_0 and tau_i is to apply post-hoc contrasts.
That is fairly easy to do with the multcomp package.

alpha_0 = (m1 + m2 + m3) / 3
m1 = intercept
m2 = intercept + method2
m3 = intercept + method3
hence alpha_0 = intercept + method2/3 + method3/3

m1 = alpha_0 + tau_1
tau_1 = intercept - method2/3 - method3/3

Best regards,

ir. Thierry Onkelinx
Instituut voor natuur- en bosonderzoek / Research Institute for Nature and
Forest
team Biometrie & Kwaliteitszorg / team Biometrics & Quality Assurance
Kliniekstraat 25
1070 Anderlecht
Belgium

To call in the statistician after the experiment is done may be no more
than asking him to perform a post-mortem examination: he may be able to say
what the experiment died of. ~ Sir Ronald Aylmer Fisher
The plural of anecdote is not data. ~ Roger Brinner
The combination of some data and an aching desire for an answer does not
ensure that a reasonable answer can be extracted from a given body of data.
~ John Tukey

2016-05-29 21:23 GMT+02:00 li li <hannah.hlx op gmail.com>:

> Hi all,
>   For the following data, I consider the following random intercept and
> random slope model. Denote as y_ijk the response value from *j*th
> individual within *i*th method at time point *k*. Assume the following
> model for y_ijk:
>
>       y_ijk= (alpha_0+ tau_i +a_j(i))+(beta_i+b_j(i)) T_k + e_ijk
>
>
> Here alpha_0 is the grand mean;
>           tau_i is the fixed effect for ith method;
>           a_j(i) is random intercept corresponding to the *j*th individual
> within *i*th method, assumed to be common for all three methods;
>           beta_i is the fixed slope corresponding to the ith method;
>           b_j(i) is the random slope corresponding to jth individual for
> the ith method, assumed to be different for different methods;
>           T_k is the time corresponding to y_ijk;
>           e_ijk is the residual.
>
> For this model, I consider the three specification using  the lme function
> as follows:
>
>
> mod1 <- lme(fixed= reponse ~ method*time, random=~ 1 +time | individual,
> data=one, weights= varIdent(form=~1|method),
>             control = lmeControl(opt = "optim"))
>
> mod2 <- lme(fixed= reponse ~ method*time, random=~ 0 +time | individual,
> data=one, weights= varIdent(form=~1|method),
>             control = lmeControl(opt = "optim"))
>
> mod3 <- lme(fixed= reponse ~ method*time, random=~ method +time |
> individual, data=one, weights= varIdent(form=~1|method),
>             control = lmeControl(opt = "optim"))
>
> I think mod1 is the correct one. However, I am kind of confused with the
> right usage of lme function. Can someone familiar with this give some help
> here?
>
> Another question is regarding the fixed effect   tau_1, tau_2 and tau_3
> (corresponding to the three methods). One main question I am interested in
> is whether each of them are statistically different from zero. In the
> summary results below (shaded part), it looks that the result for method 2
> and 3 are given with reference to method 1). Is there a way to obtain
> specific result separately for alpha_0 (the overall mean) and also tau_1,
> tau_2 and tau3?
>
> Thanks very much for the help!
>    Hanna
>
> > summary(mod1)Linear mixed-effects model fit by REML
>  Data: one
>        AIC      BIC    logLik
>   304.4703 330.1879 -140.2352
>
> Random effects:
>  Formula: ~1 + time | individual
>  Structure: General positive-definite, Log-Cholesky parametrization
>             StdDev       Corr
> (Intercept) 0.2487869075 (Intr)
> time        0.0001841179 -0.056
> Residual    0.3718305953
>
> Variance function:
>  Structure: Different standard deviations per stratum
>  Formula: ~1 | method
>  Parameter estimates:
>        3        1        2
>  1.00000 26.59750 24.74476
> Fixed effects: reponse ~ method * time
>                 Value Std.Error DF   t-value p-value(Intercept)
> 96.65395  3.528586 57 27.391694  0.0000
> method2       1.17851  4.856026 57  0.242689  0.8091
> method3       5.87505  3.528617 57  1.664973  0.1014time
> 0.07010  0.250983 57  0.279301  0.7810
> method2:time -0.12616  0.360585 57 -0.349877  0.7277
> method3:time -0.08010  0.251105 57 -0.318999  0.7509
>  Correlation:
>              (Intr) methd2 methd3 time   mthd2:
> method2      -0.726
> method3      -0.999  0.726
> time         -0.779  0.566  0.779
> method2:time  0.542 -0.712 -0.542 -0.696
> method3:time  0.778 -0.566 -0.779 -0.999  0.696
>
> Standardized Within-Group Residuals:
>         Min          Q1         Med          Q3         Max
> -2.67575293 -0.51633192  0.06742723  0.59706762  2.81061874
>
> Number of Observations: 69
> Number of Groups: 7 >
>
>
>
>
>
> > one   response individual time method
> 1    102.9          3    0      3
> 2    103.0          3    3      3
> 3    103.0          3    6      3
> 4    102.8          3    9      3
> 5    102.2          3   12      3
> 6    102.5          3   15      3
> 7    103.0          3   18      3
> 8    102.0          3   24      3
> 9    102.8          1    0      3
> 10   102.7          1    3      3
> 11   103.0          1    6      3
> 12   102.2          1    9      3
> 13   103.0          1   12      3
> 14   102.8          1   15      3
> 15   102.8          1   18      3
> 16   102.9          1   24      3
> 17   102.2          2    0      3
> 18   102.6          2    3      3
> 19   103.4          2    6      3
> 20   102.3          2    9      3
> 21   101.3          2   12      3
> 22   102.1          2   15      3
> 23   102.1          2   18      3
> 24   102.2          2   24      3
> 25   102.7          4    0      3
> 26   102.3          4    3      3
> 27   102.6          4    6      3
> 28   102.7          4    9      3
> 29   102.8          4   12      3
> 30   102.5          5    0      3
> 31   102.4          5    3      3
> 32   102.1          5    6      3
> 33   102.3          6    0      3
> 34   102.3          6    3      3
> 35   101.9          7    0      3
> 36   102.0          7    3      3
> 37   107.4          3    0      1
> 38   101.3          3   12      1
> 39    92.8          3   15      1
> 40    73.7          3   18      1
> 41   104.7          3   24      1
> 42    92.6          1    0      1
> 43   101.9          1   12      1
> 44   106.3          1   15      1
> 45   104.1          1   18      1
> 46    95.6          1   24      1
> 47    79.8          2    0      1
> 48    89.7          2   12      1
> 49    97.0          2   15      1
> 50   108.4          2   18      1
> 51   103.5          2   24      1
> 52    96.4          4    0      1
> 53    89.3          4   12      1
> 54   112.6          5    0      1
> 55    93.3          6    0      1
> 56    99.6          7    0      1
> 57   109.5          3    0      2
> 58    98.5          3   12      2
> 59   103.5          3   24      2
> 60   113.5          1    0      2
> 61    94.5          1   12      2
> 62    88.5          1   24      2
> 63    99.5          2    0      2
> 64    97.5          2   12      2
> 65    98.5          2   24      2
> 66   103.5          4    0      2
> 67    89.5          5    0      2
> 68    87.5          6    0      2
> 69    82.5          7    0      2
>
>         [[alternative HTML version deleted]]
>
> ______________________________________________
> R-help op r-project.org mailing list -- To UNSUBSCRIBE and more, see
> https://stat.ethz.ch/mailman/listinfo/r-help