Hi,

Sorry for emailing about multivariate mixed models again, but I've  
read old postings and I'm still non the wiser.  I've simulated some  
data (see below) with 3 response variables each measured ONCE for a  
thousand individuals. Individuals are grouped into 100 families (10  
per family). The aim is to estimate the (co)variances within and  
between the three traits at the family level and the individual  
level.  By concatenating the three traits into a single vector, and  
making a dummy variable indexing traits, I can get the appropriate  
random effect structures using random interactions (e.g. (trait-1| 
family)).  However, the problem arises with the individual term.  
Because each individual is only measured once per trait, the residual  
variance is not identifiable, and I cannot find a way of fixing the  
residual variance to zero.  I can fit the appropraite model simply in  
ASReml, but I would like to fit an equivalent Poisson model using  
lmer which is less straightforward in ASReml.  My guess is I would  
need to fix the scale parameter to zero, but I'm not sure if this can  
be done, or if it is safe to do so.

Thanks for any help,

Jarrod


library(MASS)
library(lme4)
library(bayesm)

G<-rwishart(4, diag(3))$IW  # random covariance matrix for family
E<-diag(3)                  # identity covariance matrix for indiviual
n<-1000                     # number of individuals

r<-mvrnorm(n/10, c(0,0,0), G) # sample random effects for family
e<-mvrnorm(n, c(0,0,0), E)    # sample random effects for individuals
fam<-rep(1:(n/10),10)         # create family ids

y<-c(r[fam,]+e)                         # create trait values on  
underlying score
yP<-rpois(n*3, exp(y))                   # sample count data

lmedata<-data.frame(Rfac=as.factor(rep(fam, 3)), y=y, yP=yP,  
tr=as.factor(rep(1:3, each=n)), id=as.factor(rep(1:n, each=3)))

model1a<-lmer(y~tr+(tr-1|fam), data=lmedata)

# gets good estimates but not very flexible because E is assumed to  
be an
# identity matrix multiplied by a scalar: which of course it is in  
this example.
# a more relaistic model would be something like:

model1b<-lmer(y~tr+(tr-1|fam)+(tr-1|id), data=lmedata)

# but the residual variance is not identifiable so the model gives  
odd results for
# the id term.  Ideally I woud fix the residual variance to zero, but  
a) I don't seem
# to be able to do this and b) I guess there may be comutational  
problems in doing this.

# the poisson model also gets good G estimates, and the scale  
parameter is greater
# then 1 as expected:

model2a<-lmer(yP~tr+(tr-1|fam), data=lmedata, family=poisson)

# but again the full multivaraite model is hard to interpret

model2b<-lmer(y~tr+(tr-1|fam)+(tr-1|id), data=lmedata, family=poisson)
	[[alternative HTML version deleted]]


