[R-sig-ME] within subjects, 3 conditions, 3 lines
Stanislav Aggerwal
stan.aggerwal at gmail.com
Mon Jan 19 15:12:58 CET 2015
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]]
More information about the R-sig-mixed-models
mailing list