[R] [R-sig-ME] multcomp package

li li hannah.hlx at gmail.com
Tue Jun 7 05:29:10 CEST 2016


Thanks Evan and Gabriel for the reply. I think it might help me make the
question clearer if I show the data and the model here (I actually asked
questions related to this data before but I still need some help). The data
looks like the following:

   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



I used the following random intercept and random slope model for this data.

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 common for all three methods;
          T_k is the time corresponding to y_ijk;
          e_ijk is the residual.

Here I used the following specification for the lme function

mod1 <- lme(fixed= reponse ~ method*time, random=~ 1 +time | individual,
data=one, weights= varIdent(form=~1|method),
            control = lmeControl(opt = "optim"))

I did not add the method in random effects  because here I assumed common
random slope for all three methods.

This model is used to initially check whether there fixed slopes are equal.
I wanted to further evaluate whether each fixed slope (beta_1, beta_2 and
beta 3) is significantly different from zero. I was hoping to evaluate this
based on the same model.

The output is as follows. Does the highlighted part below already gives the
result for testing beta_1=0; beta_2=0 and beta_3=0?

Thanks very much.
   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










2016-06-06 11:21 GMT-04:00 Gabriel Baud-Bovy <baud-bovy.gabriel at hsr.it>:

> On 06/06/2016 4:57 PM, li li wrote:
>
> Hi all,
>   After fitting a random slope and random intercept model using lme
> function, I want
> to test whether each of the fixed slopes is equal to zero (The output of
> model is below).
> Can this be done (testing each individual slope) using multcomp package?
>
> I don't understand what you mean by testing individual slope ?
>
> For the fixed effects, you might test whether there is a method, time or
> interaction
> effect using one of the methods described below
>
> For the randome effects, according to your model specification, the time
> dependency might vary for each individual. the sd for the time
> (0.0001841179)
> is small.  You might want to test whether to include randome slope  by
> doing a LRT
> between a model with it and another model without.
>
> Whay not include method in the random effects ?
>
> Gabriel
>
>   Thanks much for the help.
>    Hanna
>
> To get p values:
>
>
> http://stats.stackexchange.com/questions/118416/getting-p-value-with-mixed-effect-with-lme4-package
>
>
> http://mindingthebrain.blogspot.it/2014/02/three-ways-to-get-parameter-specific-p.html
>
> Using lmerTest package
> https://cran.r-project.org/web/packages/lmerTest/lmerTest.pdf
>
> or using mixed in afex package
> http://rpackages.ianhowson.com/cran/afex/man/mixed.html
>
> both use pbkrtest packages
> https://cran.r-project.org/web/packages/pbkrtest/pbkrtest.pdf
>
> a faq
> http://glmm.wikidot.com/faq
>
>
>
>
> 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 >
>
> 	[[alternative HTML version deleted]]
>
> _______________________________________________R-sig-mixed-models at r-project.org mailing listhttps://stat.ethz.ch/mailman/listinfo/r-sig-mixed-models
> .
>
>
>
>
> --
> ---------------------------------------------------------------------
> Gabriel Baud-Bovy               tel.: (+39) 02 2643 4839 (office)
> UHSR University                       (+39) 02 2643 3429 (laboratory)
> via Olgettina, 58                     (+39) 02 2643 4891 (secretary)
> 20132 Milan, Italy               fax: (+39) 02 2643 4892
> ---------------------------------------------------------------------
>
>

	[[alternative HTML version deleted]]



More information about the R-help mailing list