[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
> PLEASE do read the posting guide
> http://www.R-project.org/posting-guide.html
> and provide commented, minimal, self-contained, reproducible code.
>
[[alternative HTML version deleted]]
More information about the R-help
mailing list