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

[Stanislav Aggerwal]

I am trying to figure out mixed models using a simulated example.

I get the essentially the same parameter values and SEs whether I do:
1. lm() ignoring the within subjects nature of the expt
2. lmer()
3. fit a line to each subject's data and find mean and SE of parameters
across subjects

At this point, my feeling is I should just do (1) or possibly (3) above.
Please tell me why I should do anything else. Maybe my example is somehow
weird? It is typical of expts I deal with.

- cond is a factor with levels a, b, and c
- x is a continuous indep var
- each condition has slope and intercept
- each subject receives all conditions; each subject has own slope
and intercept randomly varying about population value

==========================

set.seed(1234)
nsubj<-30
z<-expand.grid(x=1:5,cond=c("a","b","c"),subj=as.factor(1:nsubj))
x<-z$x
cond<-z$cond
subj<-z$subj

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

y<- (cond=="a")*(a1 + a4*x) + (cond=="b")*(a1+a2 + (a4+a5)*x) +
  (cond=="c")*(a1+a3 + (a4+a6)*x) + rnorm(5*3*nsubj,mean=0,sd=.5)

#### ignore within subjects design and treat as independent errors
fit1<-lm(y~cond*x)
summary(fit)

## mixed model using lmer()
library(lme4)
fit2<-lmer(y~x*cond + (x*cond|subj))
summary(fit2)

## lme from nlme package
library(nlme)
fit3<-lme(y~x*cond, random=~(x*cond)|subj)
#does not converge

##fit line to each subject; get mean and se of params across subjects
beta<-matrix(0,ncol=6,nrow=nsubj)
for(i in 1:nsubj)
  {
  f<-lm(y[subj==i]~cond[subj==i]*x[subj==i])
  beta[i,]<-coef(f)
  }

colMeans(beta)
apply(beta,2,sd)/sqrt(nsubj)
# 0.9730283  0.3577780  1.8709841  3.0270420 -0.8013505 -1.9840602
# 0.10121264 0.17860228 0.15895837 0.02820369 0.04765495 0.03834946

======================

This last method makes the most sense to me. Getting 95% CIs and p-values
for the params would be straightforward.

Thanks very much for your help.

Stan

	[[alternative HTML version deleted]]