Dear all,
For the data below, I would like to fit a model with common random
slope and common random intercept as shown below. I am interested in
obtaining
separate fixed effect estimates (intercept and slope and corresponding
hypothesis test)
for each method. Instead of performing the analysis method by method, can I
use the syntax as below, specifically, the part "fixed = response ~
method/time"?
I know this is legitimate specification if there is no random effects
involved.
Thanks so much in advance!
Hanna
lme(fixed= response ~ method/time, random=~ 1+time | lot, data=dat,
weights= varIdent(form=~1|method), control = lmeControl(opt =
"optim"),
na.action = na.exclude)
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]]
li li <hannah.hlx <at> gmail.com> writes:> > Dear all,> For the data below, I would like to fit a model with common > random slope and common random intercept as shown below. I am > interested in obtaining separate fixed effect estimates (intercept > and slope and corresponding hypothesis test) for each > method. Instead of performing the analysis method by method, can I > use the syntax as below, specifically, the part "fixed = response ~ > method/time"? I know this is legitimate specification if there is > no random effects involved. Thanks so much in advance! HannaThis is probably a better question for the r-sig-mixed-models mailing list ... followups there, please. (Also, please don't post in HTML) Yes, that seems as though it should work. dat <- read.table(header=TRUE, text="response individual time method 102.9 3 0 3 103.0 3 3 3 103.0 3 6 3 102.8 3 9 3 102.2 3 12 3 102.5 3 15 3 103.0 3 18 3 102.0 3 24 3 102.8 1 0 3 102.7 1 3 3 103.0 1 6 3 102.2 1 9 3 103.0 1 12 3 102.8 1 15 3 102.8 1 18 3 102.9 1 24 3 102.2 2 0 3 102.6 2 3 3 103.4 2 6 3 102.3 2 9 3 101.3 2 12 3 102.1 2 15 3 102.1 2 18 3 102.2 2 24 3 102.7 4 0 3 102.3 4 3 3 102.6 4 6 3 102.7 4 9 3 102.8 4 12 3 102.5 5 0 3 102.4 5 3 3 102.1 5 6 3 102.3 6 0 3 102.3 6 3 3 101.9 7 0 3 102.0 7 3 3 107.4 3 0 1 101.3 3 12 1 92.8 3 15 1 73.7 3 18 1 104.7 3 24 1 92.6 1 0 1 101.9 1 12 1 106.3 1 15 1 104.1 1 18 1 95.6 1 24 1 79.8 2 0 1 89.7 2 12 1 97.0 2 15 1 108.4 2 18 1 103.5 2 24 1 96.4 4 0 1 89.3 4 12 1 112.6 5 0 1 93.3 6 0 1 99.6 7 0 1 109.5 3 0 2 98.5 3 12 2 103.5 3 24 2 113.5 1 0 2 94.5 1 12 2 88.5 1 24 2 99.5 2 0 2 97.5 2 12 2 98.5 2 24 2 103.5 4 0 2 89.5 5 0 2 87.5 6 0 2 82.5 7 0 2 ") library(nlme) ## check design with(dat,table(method,time)) with(dat,table(method,individual)) with(dat,table(time,individual)) dat <- transform(dat, method=factor(method), individual=factor(individual)) library(ggplot2); theme_set(theme_bw()) ggplot(dat,aes(time,response,colour=method))+geom_point()+ stat_summary(fun.y=mean,geom="line", aes(group=individual),colour="black",alpha=0.6) library(nlme) Is 'lot' supposed to be 'individual' here? Or is there another variable we don't know about? m1 <- lme(fixed= response ~ method/time, random=~ 1+time | individual, data=dat, weights= varIdent(form=~1|method), control = lmeControl(opt = "optim"), na.action = na.exclude) summary(m1)$tTable