lme nesting/interaction advice

The main point of my question is, having a 3 way anova (or ancova, if  
you prefer), with *no* nesting, 2 fixed effects and 1 random effect,  
why is it so boneheaded difficult to specify a bog standard fully  
crossed model? I'm not talking about some rarified esoteric model  
here, we're talking about stuff tought in a first year Biology Stats  
course here[1].

That may be so, but I've never needed to use one.  

If it's bog-standard and yet boneheaded difficult, then presumably
someone else would have had this problem before you.  Perhaps a search
of the archives will help?  If you try, you will find many qualifiers
to the effect that "lme isn't very well set up for crossed random
effects".

Now, to avoid any chances of being misunderstood in my use of the  
words 'fully crossed model', what I mean is a simple

y ~ effect1 * effect2 * effect3

with effect3 being random (all all the jazz that comes from this  
fact). I fully apprecciate that the only reasonable F-tests would be  
for effect1, effect2 and effect1:effect2, but there is no way I can  
use lme to specify such simple thing without getting the *wrong*  
denDF. I need light on this topic and I'd say it's a general enough  
question not to need much more handholding than this.

Perhaps there are some circumstances unique to your situation.

I fully apprecciate that R is developed for love, not money, 

... as is the R-help community ...

and if I  
knew how to write an user friendly frontend for nlme and lme4 (and I  
knew how to actually get the model I want) I'd be pretty happy to do  
so and submit it as a library. In any case, I feel my complaint is  
pefectly valid, because specifying such basic model should ideally  
not such a chore, and I think the powers that be might actually find  
some use from user feedback.

This is not feedback.  It is a compliant.  But, the complaint boils
down to the fact that you don't know what you're doing, and you show
no evidence of having searched the R-help archives.  How is that
helpful?

Once I have sorted how to specify such trivial model I'll face the  
horror of the nesting, in any case I attach a toy dataset I created  
especially to test how to specify the correct model (silly me).

Well, these data seem to differ.  Is replica block?  If not, then how
can we reproduce your results?  And, if I assume that it is, then the
output df differ from what you sent in your original mail.  So, I find
this confusing.

Then, from your original mail,

The easiest model ignores the nested random effects and uses just
selection, males and replica and the relative interactions. The
model

lme(y ~ selection * males, random = ~1|replica/selection/males, mydata)

forgive me, but I seem to see nesting in the random statement.  That is
what happens when we separate factors with a '/'; they are nested.  We
would expect that statement to not provide the correct df for the
bog-standard fully crossed design.

Perhaps if you were to comply with the request at the bottom of each
R-help email, and provide commented, minimal, self-contained,
reproducible code, that actually ran, ideally with fewer value
judgements, you might get more attention from the people who are
smarter than you and me, but have less time than either of us.

Andrew

Andrew Robinson  
Department of Mathematics and Statistics            Tel: +61-3-8344-6410
University of Melbourne, VIC 3010 Australia         Fax: +61-3-8344-4599
http://www.ms.unimelb.edu.au/~andrewpr
http://blogs.mbs.edu/fishing-in-the-bay/

On Sun, May 11, 2008 at 07:52:50PM +0100, Federico Calboli wrote:

The main point of my question is, having a 3 way anova (or ancova, if
you prefer), with *no* nesting, 2 fixed effects and 1 random effect,
why is it so boneheaded difficult to specify a bog standard fully
crossed model? I'm not talking about some rarified esoteric model
here, we're talking about stuff tought in a first year Biology Stats
course here[1].

That may be so, but I've never needed to use one.

So what?  This is still a standard, common, garden-variety
	model that you will encounter in exercises in many (if not
	all!) textbooks on experimental design and anova.

If it's bog-standard and yet boneheaded difficult, then presumably
someone else would have had this problem before you.  Perhaps a search
of the archives will help?  If you try, you will find many qualifiers
to the effect that "lme isn't very well set up for crossed random
effects".

But that avoids the question as to *why* it isn't very well
	set up for crossed random effects?  What's the problem?
	What are the issues?  The model is indeed bog-standard.
	It would seem not unreasonable to expect that it could be
	fitted in a straightforward manner, and it is irritating to
	find that it cannot be.  If SAS and Minitab can do it at
	the touch of a button, why can't R do it?

Now, to avoid any chances of being misunderstood in my use of the
words 'fully crossed model', what I mean is a simple

y ~ effect1 * effect2 * effect3

with effect3 being random (all all the jazz that comes from this
fact). I fully apprecciate that the only reasonable F-tests would be
for effect1, effect2 and effect1:effect2, but there is no way I can
use lme to specify such simple thing without getting the *wrong*
denDF. I need light on this topic and I'd say it's a general enough
question not to need much more handholding than this.

Perhaps there are some circumstances unique to your situation.

Huh?

I fully apprecciate that R is developed for love, not money,

... as is the R-help community ...

and if I
knew how to write an user friendly frontend for nlme and lme4 (and I
knew how to actually get the model I want) I'd be pretty happy to do
so and submit it as a library. In any case, I feel my complaint is
pefectly valid, because specifying such basic model should ideally
not such a chore, and I think the powers that be might actually find
some use from user feedback.

This is not feedback.  It is a compliant.  But, the complaint boils
down to the fact that you don't know what you're doing

That's rubbish. I think it's fairly clear that Federico does
	have a pretty good idea of what he's doing, but is flummoxed
	by the arcana of lme().  As am I.

and you show
no evidence of having searched the R-help archives.  How is that
helpful?

It doesn't seem to me to be a complaint as such.  It is a
	request for insight.  I too would like some insight as to
	what on earth is going on.  And why do you say Federico
	shows no evidence of having searched the archives?  One can
	search till one is blue in the face and come away no wiser
	on this issue.

		cheers,

			Rolf Turner

######################################################################
Attention:\ This e-mail message is privileged and confid...{{dropped:9}}

On 12/05/2008, at 9:45 AM, Andrew Robinson wrote:

On Sun, May 11, 2008 at 07:52:50PM +0100, Federico Calboli wrote:

The main point of my question is, having a 3 way anova (or ancova, if
you prefer), with *no* nesting, 2 fixed effects and 1 random effect,
why is it so boneheaded difficult to specify a bog standard fully
crossed model? I'm not talking about some rarified esoteric model
here, we're talking about stuff tought in a first year Biology Stats
course here[1].

That may be so, but I've never needed to use one.

	So what?  This is still a standard, common, garden-variety
	model that you will encounter in exercises in many (if not
	all!) textbooks on experimental design and anova.

To reply in similar vein, so what?  Why should R-core or the R
community feel it necessary to reproduce every textbook example?  How
many times have *you* used such a model in real statistical work,
Rolf?

If it's bog-standard and yet boneheaded difficult, then presumably
someone else would have had this problem before you.  Perhaps a search
of the archives will help?  If you try, you will find many qualifiers
to the effect that "lme isn't very well set up for crossed random
effects".

	But that avoids the question as to *why* it isn't very well
	set up for crossed random effects?  What's the problem?
	What are the issues?  The model is indeed bog-standard.
	It would seem not unreasonable to expect that it could be
	fitted in a straightforward manner, and it is irritating to
	find that it cannot be.  If SAS and Minitab can do it at
	the touch of a button, why can't R do it?

Bates has made no secret of the fact that lme was intended first and
foremost for nested designs, and that support for crossed designs is
not promised.  He has said so on many occasions, as a search would
find.  He is now working on lme4, which will support crossed designs.
It's not done yet.

and you show
no evidence of having searched the R-help archives.  How is that
helpful?

	It doesn't seem to me to be a complaint as such.  It is a
	request for insight.  I too would like some insight as to
	what on earth is going on.  And why do you say Federico
	shows no evidence of having searched the archives?  One can
	search till one is blue in the face and come away no wiser
	on this issue.

At least one can know that there is an issue, which apparently
Federico previously did not.

Warm wishes

Andrew

Andrew Robinson  
Department of Mathematics and Statistics            Tel: +61-3-8344-6410
University of Melbourne, VIC 3010 Australia         Fax: +61-3-8344-4599
http://www.ms.unimelb.edu.au/~andrewpr
http://blogs.mbs.edu/fishing-in-the-bay/

This leads to the notion of a fully cross design.  A fully cross
design has the characteristics that:
(a) its lineaments are not clear.
(b) it leads to heated discussion.

John Maindonald             email: john.maindonald at anu.edu.au
phone : +61 2 (6125)3473    fax  : +61 2(6125)5549
Centre for Mathematics & Its Applications, Room 1194,
John Dedman Mathematical Sciences Building (Building 27)
Australian National University, Canberra ACT 0200.

On Sun, May 11, 2008 at 07:52:50PM +0100, Federico Calboli wrote:

The main point of my question is, having a 3 way anova (or ancova, if
you prefer), with *no* nesting, 2 fixed effects and 1 random effect,
why is it so boneheaded difficult to specify a bog standard fully
crossed model? I'm not talking about some rarified esoteric model
here, we're talking about stuff tought in a first year Biology Stats
course here[1].

That may be so, but I've never needed to use one.

If it's bog-standard and yet boneheaded difficult, then presumably
someone else would have had this problem before you.  Perhaps a search
of the archives will help?  If you try, you will find many qualifiers
to the effect that "lme isn't very well set up for crossed random
effects".

Now, to avoid any chances of being misunderstood in my use of the
words 'fully crossed model', what I mean is a simple

y ~ effect1 * effect2 * effect3

with effect3 being random (all all the jazz that comes from this
fact). I fully apprecciate that the only reasonable F-tests would be
for effect1, effect2 and effect1:effect2, but there is no way I can
use lme to specify such simple thing without getting the *wrong*
denDF. I need light on this topic and I'd say it's a general enough
question not to need much more handholding than this.

Perhaps there are some circumstances unique to your situation.

I fully apprecciate that R is developed for love, not money,

... as is the R-help community ...

and if I
knew how to write an user friendly frontend for nlme and lme4 (and I
knew how to actually get the model I want) I'd be pretty happy to do
so and submit it as a library. In any case, I feel my complaint is
pefectly valid, because specifying such basic model should ideally
not such a chore, and I think the powers that be might actually find
some use from user feedback.

This is not feedback.  It is a compliant.  But, the complaint boils
down to the fact that you don't know what you're doing, and you show
no evidence of having searched the R-help archives.  How is that
helpful?

Once I have sorted how to specify such trivial model I'll face the
horror of the nesting, in any case I attach a toy dataset I created
especially to test how to specify the correct model (silly me).

Well, these data seem to differ.  Is replica block?  If not, then how
can we reproduce your results?  And, if I assume that it is, then the
output df differ from what you sent in your original mail.  So, I find
this confusing.

Then, from your original mail,

The easiest model ignores the nested random effects and uses just
selection, males and replica and the relative interactions. The
model

lme(y ~ selection * males, random = ~1|replica/selection/males,  
mydata)

forgive me, but I seem to see nesting in the random statement.  That  
is
what happens when we separate factors with a '/'; they are nested.  We
would expect that statement to not provide the correct df for the
bog-standard fully crossed design.

Perhaps if you were to comply with the request at the bottom of each

R-help email, and provide commented, minimal, self-contained,
reproducible code, that actually ran, ideally with fewer value
judgements, you might get more attention from the people who are
smarter than you and me, but have less time than either of us.

Andrew

-- 
Andrew Robinson
Department of Mathematics and Statistics            Tel:  
+61-3-8344-6410
University of Melbourne, VIC 3010 Australia         Fax:  
+61-3-8344-4599
http://www.ms.unimelb.edu.au/~andrewpr
http://blogs.mbs.edu/fishing-in-the-bay/

_______________________________________________
R-sig-mixed-models at r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/r-sig-mixed-models

The main point of my question is, having a 3 way anova (or ancova, if
you prefer), with *no* nesting, 2 fixed effects and 1 random effect,
why is it so boneheaded difficult to specify a bog standard fully
crossed model? I'm not talking about some rarified esoteric model
here, we're talking about stuff tought in a first year Biology Stats
course here[1].

y ~ effect1 * effect2 * effect3

with effect3 being random (all all the jazz that comes from this
fact). I fully apprecciate that the only reasonable F-tests would be
for effect1, effect2 and effect1:effect2, but there is no way I can
use lme to specify such simple thing without getting the *wrong*
denDF. I need light on this topic and I'd say it's a general enough

You _might_ be able to use the kinship package (lmekin, which in turn uses
lme) for this.  I'm just not sure about the "random regression" 
(fixed:random) interaction terms.

David Duffy,

| David Duffy (MBBS PhD)                                         ,-_|\
| email: davidD at qimr.edu.au  ph: INT+61+7+3362-0217 fax: -0101  /     *
| Epidemiology Unit, Queensland Institute of Medical Research   \_,-._/
| 300 Herston Rd, Brisbane, Queensland 4029, Australia  GPG 4D0B994A v

On 10 May 2008, at 07:36, Kingsford Jones wrote:

Federico,

I think you'll be more likely to receive the type of response you're
looking for if you formulate your question more clearly.  The
inclusion of "commented, minimal, self-contained, reproducible code"
(as is requested at the bottom of every email sent by r-help) is an
effective way to clarify the issues.  Also, when asking a question
about fitting a model it's helpful to describe the specific research
questions you want the model to answer.

<snip>

I apprecciate that my description of the *full* model is not 100%  
clear, but my main beef was another.

The main point of my question is, having a 3 way anova (or ancova,  
if you prefer), with *no* nesting, 2 fixed effects and 1 random  
effect, why is it so boneheaded difficult to specify a bog standard  
fully crossed model? I'm not talking about some rarified esoteric  
model here, we're talking about stuff tought in a first year Biology  
Stats course here[1].

Now, to avoid any chances of being misunderstood in my use of the  
words 'fully crossed model', what I mean is a simple

y ~ effect1 * effect2 * effect3

with effect3 being random (all all the jazz that comes from this  
fact). I fully apprecciate that the only reasonable F-tests would be  
for effect1, effect2 and effect1:effect2, but there is no way I can  
use lme to specify such simple thing without getting the *wrong*  
denDF. I need light on this topic and I'd say it's a general enough  
question not to need much more handholding than this.

There is only one random effect, so where does the crossing come  
from ? The fixed effects vary across blocks, but they are fixed so are  
just covariates. For this type of data the usual model in lme4 is  
y~fixed1+fixed2+1|group and for lme split into fixed and random parts.

Having said that, I did look at the mixed-effects mailing list  
before posting here, and it looks like it was *not* the right place  
to post anyway:

'This mailing list is primarily for useRs and programmeRs interested  
in *development* and beta-testing of the lme4 package.'

although the R-Me is now CC'd in this.

I fully apprecciate that R is developed for love, not money, and if  
I knew how to write an user friendly frontend for nlme and lme4 (and  
I knew how to actually get the model I want) I'd be pretty happy to  
do so and submit it as a library. In any case, I feel my complaint  
is pefectly valid, because specifying such basic model should  
ideally not such a chore, and I think the powers that be might  
actually find some use from user feedback.

The problems seems to be that you want lme to work in the same way as  
an ANOVA table and it doesn't. The secret with lme and lme4 is to  
think about the structure of the data and describe with an equation.  
Then each term in the equation corresponds to part of the model  
definition in R.

Once I have sorted how to specify such trivial model I'll face the  
horror of the nesting, in any case I attach a toy dataset I created  
especially to test how to specify the correct model (silly me).

I'm a bit lost with your data file, it has 4 covariates, which is more  
than enough for 2 fixed effects, assuming block is the grouping and y  
the outcome.

Ken

Best,

Federico Calboli

[1] So much bog standard that the Zar, IV ed, gives a nice table of  
how to compute the F-tests correctly, taking into account that one  
of the 3  effects is randon (I'll send the exact page and table  
number tomorrow, I don't have the book at home).

<testdat.txt>
--
Federico C. F. Calboli
Department of Epidemiology and Public Health
Imperial College, St. Mary's Campus
Norfolk Place, London W2 1PG

Tel +44 (0)20 75941602   Fax +44 (0)20 75943193

f.calboli [.a.t] imperial.ac.uk
f.calboli [.a.t] gmail.com

_______________________________________________
R-sig-mixed-models at r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/r-sig-mixed-models

lme(y ~ selection * males, random = ~1|replica/selection/males,  
mydata)

forgive me, but I seem to see nesting in the random statement.   
That is
what happens when we separate factors with a '/'; they are nested.  We
would expect that statement to not provide the correct df for the
bog-standard fully crossed design.

Please read page 23/24 of the Pinheiro and Bates book, "Mixed-Effects  
Models in S and S-PLUS". It might prove enlightening.

F

--
Federico C. F. Calboli
Department of Epidemiology and Public Health
Imperial College, St. Mary's Campus
Norfolk Place, London W2 1PG

Tel +44 (0)20 75941602   Fax +44 (0)20 75943193

f.calboli [.a.t] imperial.ac.uk
f.calboli [.a.t] gmail.com

	It doesn't seem to me to be a complaint as such.  It is a
	request for insight.  I too would like some insight as to
	what on earth is going on.  And why do you say Federico
	shows no evidence of having searched the archives?  One can
	search till one is blue in the face and come away no wiser
	on this issue.

		cheers,

			Rolf Turner

Cheers for the support Rolf. I have searched the archives, and have  
the Pinheiro and Bates book in front of my nose (plus MASS4 and many  
others).

The bottom line here is, I'm pretty cool with RTFM and all that, my  
problem is that I do not have a clear FM to read about my issue, and  
hence I have to pester (because that how people seem to feel) the  
list. I apologise for asking an inconvenient question.

Fede

--
Federico C. F. Calboli
Department of Epidemiology and Public Health
Imperial College, St. Mary's Campus
Norfolk Place, London W2 1PG

Tel +44 (0)20 75941602   Fax +44 (0)20 75943193

f.calboli [.a.t] imperial.ac.uk
f.calboli [.a.t] gmail.com

On Mon, May 12, 2008 at 10:34:40AM +1200, Rolf Turner wrote:

On 12/05/2008, at 9:45 AM, Andrew Robinson wrote:

On Sun, May 11, 2008 at 07:52:50PM +0100, Federico Calboli wrote:

The main point of my question is, having a 3 way anova (or  
ancova, if
you prefer), with *no* nesting, 2 fixed effects and 1 random  
effect,
why is it so boneheaded difficult to specify a bog standard fully
crossed model? I'm not talking about some rarified esoteric model
here, we're talking about stuff tought in a first year Biology  
Stats
course here[1].

That may be so, but I've never needed to use one.

	So what?  This is still a standard, common, garden-variety
	model that you will encounter in exercises in many (if not
	all!) textbooks on experimental design and anova.

To reply in similar vein, so what?  Why should R-core or the R
community feel it necessary to reproduce every textbook example?  How
many times have *you* used such a model in real statistical work,
Rolf?

There is a very important reason why R (or any other stats package)  
should *easily* face the challenge of bog standard models: because it  
is a *tool* for an end (i.e. the analysis of data to figure out what  
the heck they tell us) rather than a end in itself.

Bog standard models are *likely* to be used over and over again  
because they are *bog standard*, and they became such by being used  
*lots*.

If someone with a relatively easy model cannot use R for his job s/he  
will use something else, and the R community will *not* increase in  
numbers. Since R is a *community driven project*, you do the math on  
what that would mean in the long run.

Regards,

Federico

--
Federico C. F. Calboli
Department of Epidemiology and Public Health
Imperial College, St. Mary's Campus
Norfolk Place, London W2 1PG

Tel +44 (0)20 75941602   Fax +44 (0)20 75943193

f.calboli [.a.t] imperial.ac.uk
f.calboli [.a.t] gmail.com

There is only one random effect, so where does the crossing come  
from ? The fixed effects vary across blocks, but they are fixed so  
are just covariates. For this type of data the usual model in lme4  
is y~fixed1+fixed2+1|group and for lme split into fixed and random  
parts.

First off, whoa, an helpful reply! thanks for that, I hope I won't  
sound sarcastic or aggressive because I do not mean to be either.

Regarding your comment, the experiment was replicated three times, in  
3 different months. I would argue that for the fixed effects to be  
meaningful, they must have an effect over an above the effect:month  
interaction (because each fixed effect, and their interaction, might  
vary between each replicate). I would then argue I need to calculate

1) fixed.effect1:random.effect
2) fixed.effect2:random.effect
3) fixed.effect1:fixed.effect2:random.effect

to test if fixed.effect1 is meaningful (and use 1) as the error); if  
fixed.effect2 has is meaningful (and use 2) as the error);  
fixed.effect1:fixed.effect2 is meaningful (and use 3) as the error).

I'm happy to be correct if I am wrong here.

The problems seems to be that you want lme to work in the same way  
as an ANOVA table and it doesn't. The secret with lme and lme4 is  
to think about the structure of the data and describe with an  
equation. Then each term in the equation corresponds to part of  
the model definition in R.

I'll try to do that.

Once I have sorted how to specify such trivial model I'll face the  
horror of the nesting, in any case I attach a toy dataset I  
created especially to test how to specify the correct model (silly  
me).

I'm a bit lost with your data file, it has 4 covariates, which is  
more than enough for 2 fixed effects, assuming block is the  
grouping and y the outcome.

In the data file, 'selection' and 'males' are fixed effects, and  
'month' is the effect I am using for the model we are discussing  
here. The y was generatde with runif() just to have something, I'm  
not expecting any intersting result, just to understand how to fit  
the right model.

In the dataset 'line' is nested within 'selection' and 'block' is  
nested within 'month'. That's the nesting I will have to take into  
account once I get the more straightforward (sic!) model we're  
discussing right.

Best,

Federico


--
Federico C. F. Calboli
Department of Epidemiology and Public Health
Imperial College, St. Mary's Campus
Norfolk Place, London W2 1PG

Tel +44 (0)20 75941602   Fax +44 (0)20 75943193

f.calboli [.a.t] imperial.ac.uk
f.calboli [.a.t] gmail.com

On 12 May 2008, at 01:05, Andrew Robinson wrote:

On Mon, May 12, 2008 at 10:34:40AM +1200, Rolf Turner wrote:

On 12/05/2008, at 9:45 AM, Andrew Robinson wrote:

On Sun, May 11, 2008 at 07:52:50PM +0100, Federico Calboli wrote:

The main point of my question is, having a 3 way anova (or  
ancova, if
you prefer), with *no* nesting, 2 fixed effects and 1 random  
effect,
why is it so boneheaded difficult to specify a bog standard fully
crossed model? I'm not talking about some rarified esoteric model
here, we're talking about stuff tought in a first year Biology  
Stats
course here[1].

That may be so, but I've never needed to use one.

So what?  This is still a standard, common, garden-variety
model that you will encounter in exercises in many (if not
all!) textbooks on experimental design and anova.

To reply in similar vein, so what?  Why should R-core or the R
community feel it necessary to reproduce every textbook example?  How
many times have *you* used such a model in real statistical work,
Rolf?

There is a very important reason why R (or any other stats package)  
should *easily* face the challenge of bog standard models: because it  
is a *tool* for an end (i.e. the analysis of data to figure out what  
the heck they tell us) rather than a end in itself.

But a tool that mostly (entirely?) appears in textbooks.

Bog standard models are *likely* to be used over and over again  
because they are *bog standard*, and they became such by being used  
*lots*.

Well.  I have documentation relevant to nlme that goes back about 10
years.  I don't know when it was first added to S-plus, but I assume
that it was about then.  Now, do you think that if the thing that you
want to do was really bog standard, that noone would have raised a
fuss or solved it within 10 years?

If someone with a relatively easy model cannot use R for his job s/he  
will use something else, and the R community will *not* increase in  
numbers. Since R is a *community driven project*, you do the math on  
what that would mean in the long run.

Fewer pestering questions?  ;)

Andrew

Andrew Robinson  
Department of Mathematics and Statistics            Tel: +61-3-8344-6410
University of Melbourne, VIC 3010 Australia         Fax: +61-3-8344-4599
http://www.ms.unimelb.edu.au/~andrewpr
http://blogs.mbs.edu/fishing-in-the-bay/

Well.  I have documentation relevant to nlme that goes back about 10
years.  I don't know when it was first added to S-plus, but I assume
that it was about then.  Now, do you think that if the thing that you
want to do was really bog standard, that noone would have raised a
fuss or solved it within 10 years?

I'm pretty unpleasant, more so in person, so I'll tell you this. If  
people raised the issue and got the answer I got, I would not be  
surprised if they'd migrated to 'any other stats software' in droves.  
I have no doubt that, given the cryptic and sparse nature of the  
documentation of the issue, most people migrated well before asking -- 
on the grounds most people have a job to do, papers to publish,  
grants to write, kids to pick up from school and pretty little time  
for RTFM and all that sanctimonious attitude.

Once people stop nagging about 'whatever', the reason is because they  
finally got the message things ain't gonna improve, so cut your  
losses and look elsewhere.

Being unpleasant, thick skinned and cheap I will keep nagging and use  
R (the fact I do like it very much might be a factor). But given the  
selection users go through, it will be Vogon time sooner or later ;).

/F

--
Federico C. F. Calboli
Department of Epidemiology and Public Health
Imperial College, St. Mary's Campus
Norfolk Place, London W2 1PG

Tel +44 (0)20 75941602   Fax +44 (0)20 75943193

f.calboli [.a.t] imperial.ac.uk
f.calboli [.a.t] gmail.com

An embedded and charset-unspecified text was scrubbed...
Name: not available
URL: <https://stat.ethz.ch/pipermail/r-sig-mixed-models/attachments/20080512/89505c87/attachment.pl>

this is not an easy and not a standard model.

if your 3-way anova is fully factorial, with one factor being random,
then you have a LOT of random effects. you have a random main effect,
you will have a number of random 2-way interction terms, and also a
random 3-way interaction effect, their covariance matrix most likely to
be nonpositive definite.

The nested model is actually the simpler one, and I have the hint that
the "basics" book does not present such a complicated model in its
introductory chapter.

also, you don't make it easy for people to help you. I could not easily
read in your data and I had to generate own data. :-(

Maybe you are looking for something similar to:

set.seed(78987)
a = rep(1:3,each=120)
b = rep(1:3,3, each=40)
c = rep(1:10,9,each=4)
y = rnorm(360,0,10)
x = cbind(a,b,c,y)
colnames(x) = c("a", "b", "c", "y")
x = as.data.frame(x)

lme = lme(y ~ factor(a)*factor(b)*factor(c)-1, x, ~factor(c))
anova(lme)

but again, I don't know if that is what you are looking for, and it may
not be correct.

T

this is not an easy and not a standard model.

I know it is not easy, but I would be able to calculate the model by  
hand, were it not for the fact that I would like to understand how to  
do it using R

if your 3-way anova is fully factorial, with one factor being random,
then you have a LOT of random effects. you have a random main effect,
you will have a number of random 2-way interction terms, and also a
random 3-way interaction effect, their covariance matrix most  
likely to
be nonpositive definite.

The nested model is actually the simpler one, and I have the hint that
the "basics" book does not present such a complicated model in its
introductory chapter.

also, you don't make it easy for people to help you. I could not  
easily
read in your data and I had to generate own data. :-(

I am sorry about this, it might be because I wrote the data to a text  
file with write.table() in OSX. The data was nothing more than a toy  
model

Maybe you are looking for something similar to:

set.seed(78987)
a = rep(1:3,each=120)
b = rep(1:3,3, each=40)
c = rep(1:10,9,each=4)
y = rnorm(360,0,10)
x = cbind(a,b,c,y)
colnames(x) = c("a", "b", "c", "y")
x = as.data.frame(x)

lme = lme(y ~ factor(a)*factor(b)*factor(c)-1, x, ~factor(c))
anova(lme)

With the very same model you're using I get:

 > lme.mod = lme(y ~ selection * males * month-1, random = ~month,  
testdat)
Error in getGroups.data.frame(dataMix, groups) :
   Invalid formula for groups
 > lme.mod = lme(y ~ selection * males * (month-1), random = ~month,  
testdat)
Error in getGroups.data.frame(dataMix, groups) :
   Invalid formula for groups

I copy/paste the datset below, in case it makes things clearer.

Thanks,

Federico


selection	line	males	month	block	y
L	L1	1	a	1	13.8156357121188
L	L1	1	a	1	12.5678496952169
L	L1	1	a	1	17.1313698710874
L	L1	1	a	1	3.87016302696429
L	L1	1	a	1	13.2627072110772
L	L2	1	a	1	17.835768135963
L	L2	1	a	1	19.3615794742946
L	L2	1	a	1	1.73416316602379
L	L2	1	a	1	12.9440758333076
L	L2	1	a	1	2.09191741654649
S	S1	1	a	1	1.56137526640669
S	S1	1	a	1	17.6580698778853
S	S1	1	a	1	18.1417595115490
S	S1	1	a	1	15.5621050691698
S	S1	1	a	1	17.0240987658035
S	S2	1	a	1	12.4378062419128
S	S2	1	a	1	6.63962595071644
S	S2	1	a	1	16.6060689473525
S	S2	1	a	1	7.1222553497646
S	S2	1	a	1	18.0590278783347
L	L1	2	a	1	1.24710303940810
L	L1	2	a	1	4.62720696791075
L	L1	2	a	1	16.0327167815994
L	L1	2	a	1	6.12926463945769
L	L1	2	a	1	7.65810538828373
L	L2	2	a	1	7.44077128893696
L	L2	2	a	1	14.9197938004509
L	L2	2	a	1	13.4244954204187
L	L2	2	a	1	11.5361888066400
L	L2	2	a	1	2.60056478204206
S	S1	2	a	1	14.8965472229756
S	S1	2	a	1	18.777876078384
S	S1	2	a	1	6.80722737265751
S	S1	2	a	1	13.1697203880176
S	S1	2	a	1	3.74557441123761
S	S2	2	a	1	5.41025308240205
S	S2	2	a	1	19.8277674221899
S	S2	2	a	1	4.76206006342545
S	S2	2	a	1	3.08200096315704
S	S2	2	a	1	11.7220768791158
L	L1	1	a	2	17.8684629611671
L	L1	1	a	2	18.5609889889602
L	L1	1	a	2	1.33335256157443
L	L1	1	a	2	12.2590920312796
L	L1	1	a	2	10.3133576815017
L	L2	1	a	2	9.08117202203721
L	L2	1	a	2	11.8387454338372
L	L2	1	a	2	2.17258459446020
L	L2	1	a	2	7.64467771397904
L	L2	1	a	2	10.1472946784925
S	S1	1	a	2	6.33078282815404
S	S1	1	a	2	14.2109518861398
S	S1	1	a	2	2.86901426501572
S	S1	1	a	2	1.33705932833254
S	S1	1	a	2	3.62769498769194
S	S2	1	a	2	10.6116549053695
S	S2	1	a	2	19.2579759012442
S	S2	1	a	2	4.93543729488738
S	S2	1	a	2	14.0185110287275
S	S2	1	a	2	13.0477287801914
L	L1	2	a	2	7.81632485729642
L	L1	2	a	2	15.8365131700411
L	L1	2	a	2	13.6505087725818
L	L1	2	a	2	4.30545190884732
L	L1	2	a	2	5.62008981755935
L	L2	2	a	2	11.6415019945707
L	L2	2	a	2	17.4424436504487
L	L2	2	a	2	5.51907726703212
L	L2	2	a	2	3.24006642703898
L	L2	2	a	2	2.96195078082383
S	S1	2	a	2	16.2962908495683
S	S1	2	a	2	11.919979267288
S	S1	2	a	2	2.93063819734380
S	S1	2	a	2	15.6936508698855
S	S1	2	a	2	2.295168728102
S	S2	2	a	2	11.7390888168011
S	S2	2	a	2	10.4391786744818
S	S2	2	a	2	11.1727120177820
S	S2	2	a	2	11.4406719871331
S	S2	2	a	2	5.20650001661852
L	L1	1	b	3	9.53833453846164
L	L1	1	b	3	1.95979442284442
L	L1	1	b	3	8.31046994985081
L	L1	1	b	3	4.39192276610993
L	L1	1	b	3	16.0887561799027
L	L2	1	b	3	5.20481191016734
L	L2	1	b	3	8.32888073613867
L	L2	1	b	3	3.47083900799043
L	L2	1	b	3	12.2260039008688
L	L2	1	b	3	12.2011876017787
S	S1	1	b	3	10.6229240614921
S	S1	1	b	3	9.49411644623615
S	S1	1	b	3	18.4179708964657
S	S1	1	b	3	4.12228938611224
S	S1	1	b	3	5.55325478035957
S	S2	1	b	3	10.4130052563269
S	S2	1	b	3	12.6257456133608
S	S2	1	b	3	19.8305484240409
S	S2	1	b	3	2.6487210446503
S	S2	1	b	3	15.9869729799684
L	L1	2	b	3	14.6902481422294
L	L1	2	b	3	16.9090498948935
L	L1	2	b	3	3.90413530776277
L	L1	2	b	3	1.89122105599381
L	L1	2	b	3	14.2131580882706
L	L2	2	b	3	16.1121652075090
L	L2	2	b	3	16.8070402105805
L	L2	2	b	3	2.19568496150896
L	L2	2	b	3	2.96183063089848
L	L2	2	b	3	12.6824307644274
S	S1	2	b	3	12.8504637307487
S	S1	2	b	3	6.1809710978996
S	S1	2	b	3	7.47571811126545
S	S1	2	b	3	8.81466605700552
S	S1	2	b	3	12.3395618333016
S	S2	2	b	3	19.9728567516431
S	S2	2	b	3	11.3028548071161
S	S2	2	b	3	6.50141697842628
S	S2	2	b	3	13.0698232366703
S	S2	2	b	3	11.474319845438
L	L1	1	b	4	1.03187379334122
L	L1	1	b	4	9.6801089819055
L	L1	1	b	4	4.00592510355636
L	L1	1	b	4	4.63362230593339
L	L1	1	b	4	15.4505966042634
L	L2	1	b	4	14.9346919001546
L	L2	1	b	4	18.3153922208585
L	L2	1	b	4	8.57050257665105
L	L2	1	b	4	1.31918784533627
L	L2	1	b	4	19.5731912786141
S	S1	1	b	4	4.34415089036338
S	S1	1	b	4	9.66746214497834
S	S1	1	b	4	12.4202494181227
S	S1	1	b	4	17.0269973296672
S	S1	1	b	4	15.1591323411558
S	S2	1	b	4	7.3169305799529
S	S2	1	b	4	14.9331590640359
S	S2	1	b	4	13.4191052850801
S	S2	1	b	4	3.8695620871149
S	S2	1	b	4	14.8654785933904
L	L1	2	b	4	18.5652933202218
L	L1	2	b	4	1.33474090369418
L	L1	2	b	4	12.9836367180105
L	L1	2	b	4	8.14379613753408
L	L1	2	b	4	17.5245339989197
L	L2	2	b	4	19.4962712256238
L	L2	2	b	4	8.36269160197116
L	L2	2	b	4	16.4331527492031
L	L2	2	b	4	19.9028004400898
L	L2	2	b	4	8.9440936667379
S	S1	2	b	4	5.1467830883339
S	S1	2	b	4	13.0068403361365
S	S1	2	b	4	17.9727395507507
S	S1	2	b	4	7.75389035535045
S	S1	2	b	4	13.1991689843126
S	S2	2	b	4	1.86596546205692
S	S2	2	b	4	13.6726454407908
S	S2	2	b	4	5.46787238144316
S	S2	2	b	4	14.9116268747021
S	S2	2	b	4	17.4470446316991
L	L1	1	c	5	6.0449733575806
L	L1	1	c	5	9.23758197389543
L	L1	1	c	5	11.8805425714236
L	L1	1	c	5	11.097381018335
L	L1	1	c	5	6.2186355558224
L	L2	1	c	5	15.4524628254585
L	L2	1	c	5	10.8791305767372
L	L2	1	c	5	8.53259862400591
L	L2	1	c	5	17.7329147965647
L	L2	1	c	5	7.95771266217344
S	S1	1	c	5	14.2891584723257
S	S1	1	c	5	18.7194010075182
S	S1	1	c	5	5.8396331758704
S	S1	1	c	5	12.2250114120543
S	S1	1	c	5	19.9965940725524
S	S2	1	c	5	8.086333316518
S	S2	1	c	5	1.02154314704239
S	S2	1	c	5	2.71978635899723
S	S2	1	c	5	11.7734509080183
S	S2	1	c	5	10.7842262475751
L	L1	2	c	5	15.3013315075077
L	L1	2	c	5	6.12277858285233
L	L1	2	c	5	6.58965252665803
L	L1	2	c	5	13.0980647155084
L	L1	2	c	5	11.3858233471401
L	L2	2	c	5	13.8056075817440
L	L2	2	c	5	11.4665438272059
L	L2	2	c	5	6.37498314119875
L	L2	2	c	5	3.85318743507378
L	L2	2	c	5	7.959969348507
S	S1	2	c	5	13.8249646106269
S	S1	2	c	5	13.3112728327978
S	S1	2	c	5	9.67586784437299
S	S1	2	c	5	17.8595201659482
S	S1	2	c	5	3.14554875413887
S	S2	2	c	5	10.4300375301391
S	S2	2	c	5	15.4386686717626
S	S2	2	c	5	7.93477151589468
S	S2	2	c	5	3.66100515658036
S	S2	2	c	5	13.5765222932678
L	L1	1	c	6	11.3484169594012
L	L1	1	c	6	12.7286028855015
L	L1	1	c	6	17.2616694702301
L	L1	1	c	6	19.8529832861386
L	L1	1	c	6	13.6375724875834
L	L2	1	c	6	1.47241126280278
L	L2	1	c	6	1.95706973574124
L	L2	1	c	6	5.88739864598028
L	L2	1	c	6	8.16921139019541
L	L2	1	c	6	15.4799758975860
S	S1	1	c	6	2.47244948730804
S	S1	1	c	6	4.87303283042274
S	S1	1	c	6	6.15663269115612
S	S1	1	c	6	1.31718108500354
S	S1	1	c	6	16.4325702653732
S	S2	1	c	6	8.77171974466182
S	S2	1	c	6	13.4219867731445
S	S2	1	c	6	15.79551491444
S	S2	1	c	6	8.52955200499855
S	S2	1	c	6	14.4046092615463
L	L1	2	c	6	17.4255252447911
L	L1	2	c	6	5.80786765902303
L	L1	2	c	6	3.27889802469872
L	L1	2	c	6	7.55257236934267
L	L1	2	c	6	10.3537469459698
L	L2	2	c	6	10.1472219382413
L	L2	2	c	6	15.3184245729353
L	L2	2	c	6	12.6165662519634
L	L2	2	c	6	11.2075877587777
L	L2	2	c	6	17.0927850408480
S	S1	2	c	6	19.5778706537094
S	S1	2	c	6	5.56091665173881
S	S1	2	c	6	6.32830694783479
S	S1	2	c	6	7.46368952002376
S	S1	2	c	6	14.5648785342928
S	S2	2	c	6	14.7789344959892
S	S2	2	c	6	3.21725518512540
S	S2	2	c	6	2.26359746395610
S	S2	2	c	6	9.14707987429574
S	S2	2	c	6	8.6291270784568

--
Federico C. F. Calboli
Department of Epidemiology and Public Health
Imperial College, St. Mary's Campus
Norfolk Place, London W2 1PG

Tel +44 (0)20 75941602   Fax +44 (0)20 75943193

f.calboli [.a.t] imperial.ac.uk
f.calboli [.a.t] gmail.com

	But that avoids the question as to *why* it isn't very well
	set up for crossed random effects?  What's the problem?
	What are the issues?  The model is indeed bog-standard.
	It would seem not unreasonable to expect that it could be
	fitted in a straightforward manner, and it is irritating to
	find that it cannot be.  If SAS and Minitab can do it at
	the touch of a button, why can't R do it?

I haven't followed this thread carefully, so apologies if I'm too off
base. But, in response to Rolf's questions/issues. First, SAS cannot
handle models with crossed random effects (at least well at all). SAS is
horribly incapable of handling even the simplest of models (especially
generalized linear mixed models). I can cite numerous (recent) examples
of SAS coming to a complete halt (proc nlmixed) for an analyses we were
recently working on. R (and Ubuntu) was the only solution to our
problem.

Now, lme is not optimized for crossed random effects, but lmer is. That
is why lmer is supported and lme is not really supported much. lmer is
optimized for models with nested random effects and crossed random
effects. 

When working with models with nested random effects, and software
optimized for those problems (e.g., HLM, SAS, mlWin) the
variance/covariance matrix forms a special, and simple structure that
can be easily worked with. This is not the case for models with crossed
random effects.

Software packages designed for nested random effects can be tricked into
handling models with crossed random effects, but this kludge is slow and
really inefficient.

If you want complete transparency into the why and how, here is a
citation for your review.

Best
Harold

@article{Doran:Bates:Bliese:Dowling:2007:JSSOBK:v20i02,
  author =	"Harold  Doran and Douglas  Bates and Paul  Bliese and
Maritza   Dowling",
  title =	"Estimating the Multilevel Rasch Model: With the lme4
Package",
  journal =	"Journal of Statistical Software",
  volume =	"20",
  number =	"2",
  pages =	"1--18",
  day =  	"22",
  month =	"2",
  year = 	"2007",
  CODEN =	"JSSOBK",
  ISSN = 	"1548-7660",
  bibdate =	"2007-02-22",
  URL =  	"http://www.jstatsoft.org/v20/i02",
  accepted =	"2007-02-22",
  acknowledgement = "",
  keywords =	"",
  submitted =	"2006-10-01",
}

Maybe you are looking for something similar to:

set.seed(78987)
a = rep(1:3,each=120)
b = rep(1:3,3, each=40)
c = rep(1:10,9,each=4)
y = rnorm(360,0,10)
x = cbind(a,b,c,y)
colnames(x) = c("a", "b", "c", "y")
x = as.data.frame(x)

lme = lme(y ~ factor(a)*factor(b)*factor(c)-1, x, ~factor(c))
anova(lme)

With the very same model you're using I get:

lme.mod = lme(y ~ selection * males * month-1, random = ~month, 

testdat)
Error in getGroups.data.frame(dataMix, groups) :
  Invalid formula for groups

lme.mod = lme(y ~ selection * males * (month-1), random = ~month, 

testdat)
Error in getGroups.data.frame(dataMix, groups) :
  Invalid formula for groups

Try

lme = lme(y ~ factor(a)*factor(b)*factor(c)-1, x, ~1|factor(c))

Good luck,
Christian

christian forkstam____________________________________________
|________|_NETHERLANDS________________________________________|
|   work | max planck institute, PO box 310, 6500AH nijmegen  |
|        | f.c. donders centre, PO box 9101, 6500HB nijmegen  |
|________|_SWEDEN_____________________________________________|
| mobile | +46 (0)73 0242643                                  |
|___work | mrcenter, N8, karolinska hospital, 17176 stockholm |

I *think* the syntax for the model Federico wants is this:

lmer(y~selection*males+ (selection|month) + (males|month))

I'll try and check against some back of the envelope calculations -- 
as I said, the model is, per se, nothing really new, and my data is  
fully balanced.

My lme syntax is a bit rusty, so I'm not confident how to recode  
with nested random effects, as in P&B p24.


Two quick points:
1. I think Federico has caused some confusion on account of the way  
he used the term 'crossing'

Sorry about that, I'll try to avoid causing such confusion in the  
future.

Cheers,

/F

--
Federico C. F. Calboli
Department of Epidemiology and Public Health
Imperial College, St. Mary's Campus
Norfolk Place, London W2 1PG

Tel +44 (0)20 75941602   Fax +44 (0)20 75943193

f.calboli [.a.t] imperial.ac.uk
f.calboli [.a.t] gmail.com

I haven't followed this thread carefully, so apologies if I'm too off
base. But, in response to Rolf's questions/issues. First, SAS cannot
handle models with crossed random effects (at least well at all).  
SAS is
horribly incapable of handling even the simplest of models (especially
generalized linear mixed models). I can cite numerous (recent)  
examples
of SAS coming to a complete halt (proc nlmixed) for an analyses we  
were
recently working on. R (and Ubuntu) was the only solution to our
problem

First off, let's keep SAS out of this. I never used it, never wanted  
to use it and did not mention anywhere I wanted to get SAS-like  
results! Although, seeing how easily it creeps up, I can sympathise  
with those who have strog feelings about it! [for those with strong  
feelings about me, this is meant to be something joke-like]

Now, lme is not optimized for crossed random effects, but lmer is.  
That
is why lmer is supported and lme is not really supported much. lmer is
optimized for models with nested random effects and crossed random
effects.

When working with models with nested random effects, and software
optimized for those problems (e.g., HLM, SAS, mlWin) the
variance/covariance matrix forms a special, and simple structure that
can be easily worked with. This is not the case for models with  
crossed
random effects.

Software packages designed for nested random effects can be tricked  
into
handling models with crossed random effects, but this kludge is  
slow and
really inefficient.

If you want complete transparency into the why and how, here is a
citation for your review.

Thank you very much. I'll read the paper and hopefully get the  
answers I was looking for.

Best,

Federico

Best
Harold

@article{Doran:Bates:Bliese:Dowling:2007:JSSOBK:v20i02,
  author =	"Harold  Doran and Douglas  Bates and Paul  Bliese and
Maritza   Dowling",
  title =	"Estimating the Multilevel Rasch Model: With the lme4
Package",
  journal =	"Journal of Statistical Software",
  volume =	"20",
  number =	"2",
  pages =	"1--18",
  day =  	"22",
  month =	"2",
  year = 	"2007",
  CODEN =	"JSSOBK",
  ISSN = 	"1548-7660",
  bibdate =	"2007-02-22",
  URL =  	"http://www.jstatsoft.org/v20/i02",
  accepted =	"2007-02-22",
  acknowledgement = "",
  keywords =	"",
  submitted =	"2006-10-01",
}

_______________________________________________
R-sig-mixed-models at r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/r-sig-mixed-models

--
Federico C. F. Calboli
Department of Epidemiology and Public Health
Imperial College, St. Mary's Campus
Norfolk Place, London W2 1PG

Tel +44 (0)20 75941602   Fax +44 (0)20 75943193

f.calboli [.a.t] imperial.ac.uk
f.calboli [.a.t] gmail.com

Maybe you are looking for something similar to:

set.seed(78987)
a = rep(1:3,each=120)
b = rep(1:3,3, each=40)
c = rep(1:10,9,each=4)
y = rnorm(360,0,10)
x = cbind(a,b,c,y)
colnames(x) = c("a", "b", "c", "y")
x = as.data.frame(x)

lme = lme(y ~ factor(a)*factor(b)*factor(c)-1, x, ~factor(c))
anova(lme)

 With the very same model you're using I get:

 > lme.mod = lme(y ~ selection * males * month-1, random = ~month, testdat)

 Error in getGroups.data.frame(dataMix, groups) :
  Invalid formula for groups

 > lme.mod = lme(y ~ selection * males * (month-1), random = ~month,

testdat)
 Error in getGroups.data.frame(dataMix, groups) :
  Invalid formula for groups

 I copy/paste the datset below, in case it makes things clearer.

 Thanks,

 Federico

I tried the code at the top of this message--it did not work for me.
The example uses some unfortunate code.  For example, it is not a good
idea to have a data frame 'lme' and a function name 'lme'.  Also,
there are vector objects a, b, and c.  These vectors are then put into
a data frame x and then the lme call refers to 'factor(a)', but is
this the vector a or the column a in the data frame x?  It is a
confusing example.

I have several times found that lme works a lot better with
groupedData objects than data frames.  (Adds confusion in my opinion,
but has some nice features such as groupedData plot methods).

For your data, try this:

d=read.table("c:/data.txt",header=TRUE)
d2=groupedData(y~selection|block, data=d)
mod=lme(y ~ selection * males * month-1, random = ~month, d2)
anova(mod)


Kevin Wright

I'm entering this discussion late so I may be discussing issues that
have already been addressed.

As I understand it, Federico, you began by describing a model for data
in which two factors have a fixed set of levels and one factor has an
extensible, or "random", set of levels and you wanted to fit a model
that you described as

y ~ effect1 * effect2 * effect3

The problem is that this specification is not complete.  An
interaction of factors with fixed levels and a factor with random
levels can mean, in the lmer specification,

lmer(y ~ effect1 * effect2 + (1| effect3) + (1|effect1:effect2:effect3), ...)

or

lmer(y ~ effect1 * effect2 + (effect1*effect2 | effect3), ...)

or other variations.  When you specify a random effect or an random
interaction term you must, either explicitly or implicitly, specify
the form of the variance-covariance matrix associated with those
random effects.

The "advantage" that other software may provide for you is that it
chooses the model for you but that, of course, means that you only
have the one choice.

If you can describe how many variance components you think should be
estimated in your model and what they would represent then I think it
will be easier to describe how to fit the model.

I'm entering this discussion late so I may be discussing issues that
have already been addressed.

As I understand it, Federico, you began by describing a model for data
in which two factors have a fixed set of levels and one factor has an
extensible, or "random", set of levels and you wanted to fit a model
that you described as

y ~ effect1 * effect2 * effect3

The problem is that this specification is not complete.

My apologies for that, I thought that the above formula was the  
shorthand for what I would call the 'full' model, i.e. the single  
factors and the 2 and 3 ways interactions.

An
interaction of factors with fixed levels and a factor with random
levels can mean, in the lmer specification,

lmer(y ~ effect1 * effect2 + (1| effect3) + (1| 
effect1:effect2:effect3), ...)

or

lmer(y ~ effect1 * effect2 + (effect1*effect2 | effect3), ...)

or other variations.  When you specify a random effect or an random
interaction term you must, either explicitly or implicitly, specify
the form of the variance-covariance matrix associated with those
random effects.

I'll play around with this and see what I can get.

The "advantage" that other software may provide for you is that it
chooses the model for you but that, of course, means that you only
have the one choice.

I'm more than happy to stick to R, and to put more legwork into my  
models

If you can describe how many variance components you think should be
estimated in your model and what they would represent then I think it
will be easier to describe how to fit the model.

I'll work on that. Incidentally, what/where is the most comprehensive  
and up to date documentation for lme4? the pdfs coming with the  
package? I suspect knowing which are the right docs will help a lot  
in keeping me within the boundaries of civility and prevent me from  
annoying anyone (which is not something I sent forth to do on purpose).

Best regards,

Federico

--
Federico C. F. Calboli
Department of Epidemiology and Public Health
Imperial College, St. Mary's Campus
Norfolk Place, London W2 1PG

Tel +44 (0)20 75941602   Fax +44 (0)20 75943193

f.calboli [.a.t] imperial.ac.uk
f.calboli [.a.t] gmail.com

On Mon, May 12, 2008 at 11:22 AM, Federico Calboli

On 12 May 2008, at 17:09, Douglas Bates wrote:

I'm entering this discussion late so I may be discussing issues that
have already been addressed.

As I understand it, Federico, you began by describing a model for data
in which two factors have a fixed set of levels and one factor has an
extensible, or "random", set of levels and you wanted to fit a model
that you described as

y ~ effect1 * effect2 * effect3

The problem is that this specification is not complete.

My apologies for that, I thought that the above formula was the shorthand
for what I would call the 'full' model, i.e. the single factors and the 2
and 3 ways interactions.

As I indicated, the trick is that the interaction of a fixed factor
and a random factor can be defined in more than one way.

It sounds as if what you want is

lmer(y ~ factor1 * factor2 + (1|factor3) + (1|factor1:factor3) +
(1|factor2:factor3) + (1|factor1:factor2:factor3), ...)

but I'm not sure.

An interaction of factors with fixed levels and a factor with random
levels can mean, in the lmer specification,

lmer(y ~ effect1 * effect2 + (1| effect3) + (1|effect1:effect2:effect3),
...)

or

lmer(y ~ effect1 * effect2 + (effect1*effect2 | effect3), ...)

or other variations.  When you specify a random effect or an random
interaction term you must, either explicitly or implicitly, specify
the form of the variance-covariance matrix associated with those
random effects.

I'll play around with this and see what I can get.

The "advantage" that other software may provide for you is that it
chooses the model for you but that, of course, means that you only
have the one choice.

I'm more than happy to stick to R, and to put more legwork into my models

If you can describe how many variance components you think should be
estimated in your model and what they would represent then I think it
will be easier to describe how to fit the model.

I'll work on that. Incidentally, what/where is the most comprehensive and up
to date documentation for lme4? the pdfs coming with the package? I suspect
knowing which are the right docs will help a lot in keeping me within the
boundaries of civility and prevent me from annoying anyone (which is not
something I sent forth to do on purpose).

Documentation for lme4 is pretty sketchy at present.  I hope to remedy
that during our summer break.

I'm entering this discussion late so I may be discussing issues that
have already been addressed.

As I understand it, Federico, you began by describing a model for data
in which two factors have a fixed set of levels and one factor has an
extensible, or "random", set of levels and you wanted to fit a model
that you described as

y ~ effect1 * effect2 * effect3

The problem is that this specification is not complete.

At *last* (as Owl said to Rabbit) we're getting somewhere!!!

	I always knew that there was some basic fundamental point
	about this business that I (and I believe many others) were
	simply missing.  But I could not for the life of me get anyone
	to explain to me what that point was.  Or to put it another
	way, I was never able to frame a question that would illuminate
	just what it was that I wasn't getting.

	I now may be at a stage where I can start asking the right
	questions.

An interaction of factors with fixed levels and a factor with random
levels can mean, in the lmer specification,

lmer(y ~ effect1 * effect2 + (1| effect3) + (1| 
effect1:effect2:effect3), ...)

or

lmer(y ~ effect1 * effect2 + (effect1*effect2 | effect3), ...)

or other variations.  When you specify a random effect or an random
interaction term you must, either explicitly or implicitly, specify
the form of the variance-covariance matrix associated with those
random effects.

The "advantage" that other software may provide for you is that it
chooses the model for you but that, of course, means that you only
have the one choice.

Now may I start asking what I hope are questions that will lift
	the fog a bit?

	Let us for specificity consider a three-way model with two
	fixed effects and one random effect from the good old Rothamstead style
	agricultural experiment context:  Suppose we have a number of
	species/breeds of wheat (say) and a number of fertilizers.
	These are fixed effects.  And we have a number of fields (blocks)
	--- a random effect.  Each breed-fertilizer combination is
	applied a number of times in each field.  We ***assume*** that
	that the field or block effect is homogeneous throughout.  This
	may or may not be a ``good'' assumption, but it's not completely
	ridiculous and would often be made in practice.  And probably
	*was* made at Rothamstead.  The response would be something like
	yield in bushels per acre.

	The way that I would write the ``full'' model for this setting,
	in mathematical style is:

	Y_ijkl = mu + alpha_i + beta_j + (alpha.beta)_ij + C_k + (alpha.C)_ik
                     + (beta.C)_jk + (alpha.beta.C)_ijk + E_ijkl

	The alpha_i and beta_j are parameters corresponding to breed and  
fertilizer
	respectively; the C_k are random effects corresponding to fields or  
blocks.
	Any effect ``involving'' C is also random.

	The assumptions made by the Package-Which-Must-Not-Be-Named are (I  
think)
	that

		C_k ~ N(0,sigma_C^2)
		(alpha.C)_ik ~ N(0,sigma_aC^2)
		(beta.C)jk ~ N(0,sigma_bC^2)
		(alpha.beta.C)_ijk ~ N(0,sigma_abC^2)
		E_ijkl ~ N(0,sigma^2)

	and these random variables are *all independent*.

	Ahhhhhhhh ... perhaps I'm on the way to answering my own question.  Is
	it this assumption of ``all independent'' which is questionable?  It
	seemed innocent enough when I first learned about this stuff, lo these
	many years ago.  But .... mmmmmaybe not!

	To start with:  What would be the lmer syntax to fit the foregoing
	(possibly naive) model?  I am sorry, but I really cannot get my head
	around the syntax of lmer model specification, and I've tried.  I
	really have.  Hard.  I know I must be starting from the wrong place,
	but I haven't a clue as to what the *right* place to start from is.
	And if I'm in that boat, I will wager Euros to pretzels that there
	are others in it.  I know that I'm not the brightest bulb in the
	chandelier, but I'm not the dullest either.

	Having got there:  Presuming that I'm more-or-less on the right track
	in my foregoing conjecture that it's the over-simple dependence  
structure
	that is the problem with what's delivered by the Package-Which-Must- 
Not-Be-Named,
	how might one go about being less simple-minded?  I.e. what might be  
some
	more realistic dependence structures, and how would one specify  
these in lmer?
	And how would one assess whether the assumed dependence structure gives
	a reasonable fit to the data?

If you can describe how many variance components you think should be
estimated in your model and what they would represent then I think it
will be easier to describe how to fit the model.

How does this fit in with my conjecture (above) about what I've been
	missing all these years?  Does it fit?  How many variance components
	are there in the ``naive'' model?  It looks like 5 to me ... but maybe
	I'm totally out to lunch in what I think I'm understanding at this  
stage.
	(And besides --- there are three sorts of statistician; those who  
can count,
	and those who can't.)

	Thank you for your indulgence.

		 cheers,

			Rolf Turner

######################################################################
Attention:\ This e-mail message is privileged and confid...{{dropped:9}}