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)
