[R-sig-ME] within subjects, 3 conditions, 3 lines

ONKELINX, Thierry Thierry.ONKELINX at inbo.be
Thu Jan 22 17:16:44 CET 2015

Dear Stan,

Please keep the mailing in cc.  You use rnorm(nsubj,mean=0,sd=.1) as the random intercept. However rnorm() will be evaluate again in a1, a2, ..., a6 and thus yielding different values of the random intercept in a1, a2, ..., a6. You want

Rf <- rnorm(nsubj,mean=0,sd=.1)
a1<- rep(1+ rf,each=5)  #inta=1
a2<- rep(1 + rf,each=5)  #intb-inta=2-1=1

Note that is more clear to define the fixed, random effect and noise separately and add those components in a separate step.

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
+ 32 2 525 02 51
+ 32 54 43 61 85
Thierry.Onkelinx op inbo.be

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

Van: Stanislav Aggerwal [mailto:stan.aggerwal op gmail.com]
Verzonden: donderdag 22 januari 2015 12:25
Aan: ONKELINX, Thierry
Onderwerp: Re: [R-sig-ME] within subjects, 3 conditions, 3 lines

Thanks Thierry for your help. Could I get some clarification please?

You get the right answer for the wrong question. Because you have a bug in the code that simulates the data. You added a second noise term instead of a random intercept. Hint: set all parameters not related to the random intercept to 0 in a1 to a6. Then a1 to a6 should be equal.

There are 30*3 (subjects * conditions) random intercepts. Also 30 random slopes. There are 6 params because there are 3 conditions: (slope and intercept)*3

a1-a6 have the population values: 1, 1, 2, 3, -1, -2
Then each subject has a1-a6 values that randomly vary about these population values.
It is the model I meant, and it describes my data.
Are you saying I implemented it wrong? I don't understand how.

a1<- rep(1+rnorm(nsubj,mean=0,sd=.1),each=5)  #inta=1
a2<- rep(1+rnorm(nsubj,mean=0,sd=.1),each=5)  #intb-inta=2-1=1
a3<- rep(2+rnorm(nsubj,mean=0,sd=.1),each=5)  #intc-inta=3-1=2
a4<- rep(3+rnorm(nsubj,mean=0,sd=.1),each=5)  #slopea=3
a5<- rep(-1+rnorm(nsubj,mean=0,sd=.1),each=5) #slopeb-slopea
a6<- rep(-2+rnorm(nsubj,mean=0,sd=.1),each=5) #slopec-slopea

y[cond=="a"]<-a1    + a4*x[cond=="a"]
y[cond=="b"]<-a1+a2 + (a4+a5)*x[cond=="b"]
y[cond=="c"]<-a1+a3 + (a4+a6)*x[cond=="c"]

Then I add errors. For each subject the data are not perfectly described by a line, due to measurement error. If I did not add this, the only source of noise in the model is due to the variation in slopes and intercepts between subjects.
Are you saying I should eliminate these measurement errors?
This is the model I intend:
For each subject in each condition:
y = int + slope*x + error
each subject has 3 lines (for a, b, c) with params a1-a6

#autocorrelated errors
for(i in 1:nsubj)
  y[subj==i]<-y[subj==i] +

If you look at the plot I think it looks correct.
Thanks very much!


More information about the R-sig-mixed-models mailing list