Hi, mixed modelers,
I'm coming up with kind of an old question, I guess. I've looked through
several old mails and discussions in order to find something to help me but
I was not much good at doing this search or the problem I'm facing actually
is not really a problem and needs not to be resolved.

I'm fitting a mixed effects model using lmer, and was trying to obtain some
intervals around variance components.
The problems is that I'm evaluating reproducibility of microRNA platform, so
I have kind of a large dataset (12 chips with 46227 values on it). The
effects I'm interested into are the Probe (which comes up with 7815 levels),
the Sample (4 levels) and the replicate within sample (3 levels). Probe and
Sample are nested (the probe is the same across all chips) and Replicate is
nested within sample (replicate 1 for sample 1 can not be directly compared
with replicate 1 for sample 2).

This is the code for the model:

lme.data = lmer(log.exprs ~ (1|Probe.ID) + (1|Sample/Replicate) + (1|
Probe.ID:Sample) ,data=miRNA.exprs.all)

summary(lme.data)


Linear mixed model fit by REML

Formula: log.exprs ~ (1 | Probe.ID) + (1 | Sample/Replicate) + (1 |
Probe.ID:Sample)

   Data: miRNA.exprs.all

    AIC    BIC  logLik deviance REMLdev

 728536 728603 -364262   728521  728524

Random effects:

 Groups           Name        Variance Std.Dev.

 Probe.ID:Sample  (Intercept) 0.367684 0.60637

 Probe.ID         (Intercept) 1.245118 1.11585

 Replicate:Sample (Intercept) 0.011426 0.10689

 Sample           (Intercept) 0.027912 0.16707

 Residual                     0.172734 0.41561

Number of obs: 554724, groups: Probe.ID:Sample, 31260; Probe.ID, 7815;
Replicate:Sample, 12; Sample, 4


Fixed effects:

            Estimate Std. Error t value

(Intercept)  5.71375    0.08962   63.75


That is: Probe.ID and Sample are crossed random effects, and Replicate is
nested within Sample. I specified the model this way after some (intensive)
reading of Bates book on lme 4 and long surfing on this ML, but I suppose I
could have done crucial mistakes here...is this the case? However, after
this I extracted variance components like this


var.comp.lme = c(attributes(VarCorr(lme.data)[[1]])$stddev,attributes(
VarCorr(lme.data)[[2]])$stddev,attributes(VarCorr(lme.data)[[3]])$stddev,
attributes(VarCorr(lme.data)[[4]])$stddev,attributes(VarCorr(lme.data))$sc)

names(var.comp.lme) = c(names(VarCorr(lme.data)),"Residual")

var.comp.lme


 Probe.ID:Sample         Probe.ID Replicate:Sample           Sample
Residual

       0.6063697        1.1158488        0.1068903        0.1670674
 0.4156124



 After trying to use profile function (and not being able to install package
lme4a) I resorted to a different approach (discovered through mailing lists
also) that has been thouroughly discussed here in the past I think.

I used mcmcsamp with 30000 samples and HPDinterval. However, I found that
all intervals on variance components do not contain the estimate from the
model. Only the fixed effect estimate was contained in these intervals.
Actually, xyplot and densityplot of MCMC samples found some strangeness:
Replicate within Sample and Sample chains were actually unstable and the
samples were kind of bimodal distributed, whereas remaining chains were
quite stable after a few samples, but distributions were very skewed. The
HPD intervals were as follows:

mcmc.miRNA = mcmcsamp(lme.data,n = 30000)
HPDinterval(mcmc.miRNA)

$fixef
               lower    upper
(Intercept) 5.519983 5.915997
attr(,"Probability")
[1] 0.95

$ST
         lower     upper
[1,] 0.7575154 0.7666013
[2,] 0.9003053 0.9134174
[3,] 0.1684593 0.3649006
[4,] 0.2032486 0.6025208
attr(,"Probability")
[1] 0.95

$sigma
         lower     upper
[1,] 0.4535642 0.4554482
attr(,"Probability")
[1] 0.95

Model estimated are quite close to limits of the HPD, but inevitably they
fall outside it. Is this a problem with the model (mis-specification
possibly?) or it has to do with troubles in the mcmcsamp procedure? Before
increasing number of mcmc samples I was looking for some advice (with n =
30000 the mcmcsamp procedure takes quite a long time to run, a few hours...)


Thanks for any help!!
Niccolò

	[[alternative HTML version deleted]]


