[R] problem fitting linear mixed models

[Joan Valls]

Hello,

I'm fitting linear mixed models to gene-expression data from 
microarrays, in a data set where 4608 genes are studied.
For a sample of 5 subjects and for each gene we observe the expression 
level (Intensity) in four different tissues: N, Tp, Tx and M.
I want to test whether the expression level is different accross 
tissues. Between-subject variability is modeled with a random intercept, 
and the within-subject by allowing heteroscedastic and correlated errors 
accross tissues. The proposed model can then be fitted by

lme(Intensity~Tissue-1, 
weights=varIdent(form=~1|Tissue),correlation=corSymm(),random=~1|Subject)

I have fitted this model for each gene. As a consequence of balanced 
data, fixed-effects estimates are exactly the sample mean for each gene.
But I have found one particular gene for which this does not happen: the 
fixed-effects estimates are completely no-sense.

Finally, I haved found a solution for that: rounding data.  (see the 
code below)

I have no explanation. Is it a numerical problem?

Thank you for your time.

Joan Valls
Catalan Institute of Oncology
Barcelona


#data set for the gene
Subject<-c(rep("C4",4),rep("HM1",4),rep("C1",4),rep("997",4),rep("C3",4))
Tissue<-rep(c("N","Tp","Tx","M"),5)
IntensityA<-c(10.6720000000000010,10.564,10.6080000000000010,10.8673333333333350,10.9430000000000000,10.8910000000000000, 

11.1260000000000010,10.9693333333333330,10.7690000000000000,10.8110000000000000,10.8739999999999990,10.8890000000000010, 

11.6679999999999990,11.5320000000000000,11.7100000000000010,11.3519999999999990,11.0000000000000000,10.6300000000000010, 

10.8720000000000020,10.6133333333333350)
Gene<-data.frame(Subject=as.factor(as.character(Subject)),Tissue=as.factor(as.character(Tissue)),Intensity=IntensityA) 


# fitted model (does not work!)
mlmA<-lme(Intensity ~ 
Tissue-1,weights=varIdent(form=~1|Tissue),correlation=corSymm(),data=Gene,random=~1|Subject); 

summary(mlmA)

#sample means
lapply(split(Gene$Intensity,Gene$Tissue),mean)    


#fitted model with rounded data (it works)
Gene$Intensity<-round(Gene$Intensity,4)
mlmA<-lme(Intensity ~ 
Tissue-1,weights=varIdent(form=~1|Tissue),correlation=corSymm(),data=Gene,random=~1|Subject); 

summary(mlmA)