[R] Fixed Effects in lme function

[Ben Bolker]

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!  Hanna

  This 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