Thank you for sending the data. It is very helpful
in understanding
the nature of the problem.
For example, in your original description of your
study you referred
to week as a factor, which is a completely
reasonable term, but I
mistakenly thought that you meant an object of class
"factor" and that
was why I replied that you would be estimating too
many variances and
covariances.
I can tell you why you are having problems fitting a
mixed-effects
model. Strangely it is because there is too little
variability in the
patterns across the replicates, especially at the
early times. The
leaf number is discrete with a small range (all the
leaffra
observations in this example are 4, 5, 6, 7 or 8)
and non-decreasing
over time. (I assume the nature of the experiment
is such that the
leaf number is necessarily non-decreasing.) That
doesn't allow for
many patterns. I'm sure some clever person reading
this will be able
to tell us exactly how many different such patterns
you could get but
I will simply say "not many".
Notice that the first 10 lines show that the leaffra
is 4 at week 4 in
*every* replicate. There isn't a whole lot of
variation here for the
random effects to model.
The best place to start is with a plot of the data.
I changed the
levels of the rep factor to "f1"-"f5" and "s1"-"s5"
to indicate that
each rep is at one level of the treat. (Those who
are playing along at
home should be careful of the ordering of the
original levels because
ID10, which I now call s5, occurs between ID1 and
ID2.)
With this set of labels the cross-tabulation and
treat should be
xtabs(~ rep + treat, leaf)
treat
rep pHf pHs
f1 4 0
f2 4 0
f3 4 0
f4 4 0
f5 4 0
s1 0 4
s2 0 4
s3 0 4
s4 0 4
s5 0 4
Now look at lattice plots such as
library(lattice)
xyplot(heightfra ~ week | rep, leaf, type = c("g",
"p", "r"), layout =
c(5,2), aspect = 'xy', groups = treat)
(I enclose PDF files of these plots for each of the
three responses.)
First you can see that there is very little
variation at the low end.
Strangely enough, this causes a problem in fitting
mixed-effects
models because the mle's of the variances of the
random effects for
the intercept will be zero. The lme function does
not handle this
gracefully. The lmer function from the lme4 package
does a better job
on this type of model.
Also, note that the pattern of heightfra over time
is not linear. It
is consistently concave down. Thuso a mixed-effects
model that is
linear in week will miss much of the structure in
the data.
The point of R is to encourage you to explore your
data rather than
subjecting it to a "canned" analysis. You could try
fitting a
mixed-effects model to these data in SAS PROC MIXED
or SPSS MIXED and
I have no doubt that those packages would give you
estimates (not to
mention p-values, something that the author of lmer
has been woefully
negligent in not providing :-) but you probably
won't get much of a
hint that the model doesn't make sense. I would
prefer to start with
the plot and see what the data have to say.
The technical problem with convergence in lme is
that the mle of the
variance of the intercept term is zero. You can see
that if you use
lmer from the lme4 package instead to fit the model.
the random effect for the intercept is estimate
On Dec 14, 2007 4:40 AM, Ilona Leyer
<ileyer at yahoo.de> wrote:
Here an simple example:
rep treat heightfra leaffra leafvim
ID1 pHf 1.54 4 4 4
ID2 pHf 1.49 4 4 4
ID3 pHf 1.57 4 5 4
ID4 pHf 1.48 4 4 4
ID5 pHf 1.57 4 4 4
ID6 pHs 1.29 4 5 4
ID7 pHs 0.97 4 5 4
ID8 pHs 2.06 4 4 4
ID9 pHs 0.88 4 4 4
ID10 pHs 1.47 4 4 4
ID1 pHf 3.53 5 6 6
ID2 pHf 4.08 6 6 6
ID3 pHf 3.89 6 6 6
ID4 pHf 3.78 5 6 6
ID5 pHf 3.92 6 6 6
ID6 pHs 2.76 5 5 6
ID7 pHs 3.31 6 7 6
ID8 pHs 4.46 6 7 6
ID9 pHs 2.19 5 5 6
ID10 pHs 3.83 5 5 6
ID1 pHf 5.07 7 7 9
ID2 pHf 6.42 7 8 9
ID3 pHf 5.43 6 8 9
ID4 pHf 6.83 6 8 9
ID5 pHf 6.26 6 8 9
ID6 pHs 4.57 6 9 9
ID7 pHs 5.05 6 7 9
ID8 pHs 6.27 6 8 9
ID9 pHs 3.37 5 7 9
ID10 pHs 5.38 6 8 9
ID1 pHf 5.58 7 9 12
ID2 pHf 7.43 8 9 12
ID3 pHf 6.18 8 10 12
ID4 pHf 6.91 7 10 12
ID5 pHf 6.78 7 10 12
ID6 pHs 4.99 6 13 12
ID7 pHs 5.50 7 8 12
ID8 pHs 6.56 7 10 12
ID9 pHs 3.72 6 10 12
ID10 pHs 5.94 6 10 12
I used the procedure described in Crawley?s new R
Book.
For two of the tree response variables
(heightfra,leaffra) it doesn?t work, while it
with leafvim (but in another R session, yesterday,
leaffra worked as well...).
Here the commands: