[R] [R-sig-ME] multcomp package
RuiRui
dianxiangan32 at sina.cn
Mon Jun 13 15:22:42 CEST 2016
Hi,
It seems you want to test the fixed effect. So I think the related p-value is appropriate. I'm not sure...
--------------------------------
----- Original Message -----
From: li li <hannah.hlx at gmail.com>
To: r-sig-mixed-models <r-sig-mixed-models at r-project.org>, r-help <r-help at r-project.org>, epalmery at cns.umass.edu
Subject: Re: [R-sig-ME] multcomp package
Date: 2016-06-07 11:32
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]]
_______________________________________________
R-sig-mixed-models at r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/r-sig-mixed-models
[[alternative HTML version deleted]]
More information about the R-help
mailing list