[R] Re: Fitting particular repeated measures model with lme()

[J.R. Lockwood]

Hello,

A couple weeks ago I posted a message about using lme() to fit a model
based on the following: Students are tested in two years, and are
linked to teachers in the second year only.  Thus students are not
nested within teachers in the traditional sense.  The model for
student j in class i is:

Y_{ij0} = a_0 + e_{ij0}
Y_{ij1} = a_1 + b_i + e_{ij1}

with Var(b_i) the teacher variance component and Cov(e_{ij0},e_{ij1})
unstructured.  Thus if the data are organized by student, the "Z"
matrix in the usual linear mixed model notation has every other row
equal to a row of zeros.  For reference I include at the end of this
message a function "generate.data()" that simulates (balanced) data
according to this model.

I think I figured out a way to fit this model in lme() but I have some
questions about it.  My strategy was essentially to brute-force create
the Z matrix with columns containing indicators of the second year
teacher interleaved with zeros for the first year scores, and then
force the covariance matrix of the random effects to be a constant
times the identity.  Here is the code I am using:

############
library(nlme)
library(MASS)
ntch<-30
d<-generate.data(k=ntch)
varnames<-paste("tchid",1:ntch,sep="")
fmla<-as.formula(paste("~",paste(varnames,collapse="+"),"-1"))

lme.u2a<-lme(fixed=Y~time,data=d,random=list(tchid=pdIdent(form=fmla)),weights=varIdent(form=~1|time),correlation=corSymm(form = ~1|tchid/studid))

lme.u2b<-lme(fixed=Y~time,data=d,random=list(dummy.group=pdIdent(form=fmla)),weights=varIdent(form=~1|time),correlation=corSymm(form = ~1|dummy.group/studid))
############

I have one set of questions and one observation:

1) The two model fits provide (essentially) the same estimates, and
these estimates seem reasonable.  However, they provide different DF
in the table of fixed effects, reflecting the fact that the nesting
structures specified in the two models are different.  In the first
case, I specify that students are nested within teachers which is not
exactly true.  In the second, I use a dummy grouping variable
"dummy.group" which is identically equal to 1 for all records. 

My major questions are why/how these specifications fit the same
model, whether it is possible for me to fit the model without
specifying names to the list components in the random statement, and
more generally, whether lme() *requires* a nesting structure at all.
I guess what it comes down to is that I am having trouble
understanding how exactly specifying the random statement as a list of
formulas works.

2) intervals() on either of these objects fails with the following
   error message:

Error in structure(exp(as.vector(object)), names = c(paste("sd(", deparse(formula(object)[[2]]),  : 
	names attribute must be the same length as the vector

Any insights would be greatly appreciated.

best regards,

J.R. Lockwood
412-683-2300 x4941
lockwood at rand.org
http://www.rand.org/methodology/stat/members/lockwood/

## Function for generating data
generate.data<-function(k=50,n=25,mu=c(0,10),tau=sqrt(2),esd=sqrt(c(2,4)),corr=0.8){
  ## k=number of teachers
  ## n=number of students per teacher
  Sigma<-diag(esd)%*%matrix(c(1,corr,corr,1),ncol=2)%*%diag(esd)
  theta<-rnorm(k,sd=tau)
  theta<-cbind(0,rep(theta,each=n))
  e<-mvrnorm(n*k,mu=mu,Sigma=Sigma)
  Y<-c(t(theta))+c(t(e))
  tchid<-gl(k,2*n)
  z<-model.matrix(~tchid-1)
  s<-seq(from=1,to=(2*k*n-1),by=2)
  z[s,]<-rep(0,k)
  studid<-gl(n*k,2)
  time<-rep(c(0,1),times=n*k)
  dummy.group<-factor(rep(1,2*n*k))
  return(data.frame(studid,Y,time,dummy.group,tchid,z))
}