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

[ONKELINX, Thierry]

I made a Gist (https://gist.github.com/ThierryO/58ebc3c08da4fec58fb4) with a working example to simulate data. Note that the standard errors for the fixed effects of the mixed model are larger than those of the linear model.

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

Belgium

+ 32 2 525 02 51

+ 32 54 43 61 85

Thierry.Onkelinx op inbo.be


www.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 [stan.aggerwal op gmail.com]

Verzonden: zondag 25 januari 2015 14:00

Aan: ONKELINX, Thierry

CC: r-sig-mixed-models op r-project.org

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















Thierry, thanks for your help, but your model is completely different from what I want, so far as I can tell.






The model I want to fit is this. Each subject produces 3 lines, one for each condition



y = b1 + b2*x + error  #condition a


y = b3 + b4*x + error # condition b


y = b5 + b6*x + error #condition c


Within each subject, these three lines have differing slopes and intercepts. For each subject, each parameter (b1 to pick an example) is not identical to that of other subjects but is similar to that from other subjects.





The intercepts, b1, b3, and b5 vary between subjects. Each subject's b1, b3, and b5 is randomly sampled from N(p_i,sd_i) where p_i is the population value for each of the 3 parameters



The intercepts, b2, b4, and b6 vary between subjects. Each subject's b2, b4, and b6 is randomly sampled from N(p_i,sd_i) where p_i is the population value for each of the 3 parameters.





In my simulations I used params a1-a6 which are the R parametrisation used in fitting linear models, instead of b1-b6. E.g. a2 = b1+b2











So far as I can tell (and I looked hard), my code is producing the right sort of simulated data. And these simulated data look similar to my real data.





If we pretend that subjects is a fixed effect, we have


lm(y ~ x*cond*subj)


A different slope and intercept for each condition and each subject.





Thanks again,

Bill











On Sat, Jan 24, 2015 at 10:35 AM, Stanislav Aggerwal
<stan.aggerwal op gmail.com> wrote:



Thanks Thierry!





Stan






On Thu, Jan 22, 2015 at 4:16 PM, ONKELINX, Thierry
<Thierry.ONKELINX op inbo.be> wrote:


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,



Thierry



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

Belgium

+ 32 2 525 02 51

+ 32 54 43 61 85

Thierry.Onkelinx op inbo.be

www.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] +

as.numeric(filter(rnorm(5*3,mean=0,sd=.5),filter=0.5,method="recursive"))

  }



If you look at the plot I think it looks correct.

Thanks very much!

Cheers

Stan





Disclaimer<https://www.inbo.be/nl/disclaimer-mailberichten-van-het-inbo>




























Disclaimer<https://www.inbo.be/nl/disclaimer-mailberichten-van-het-inbo>