[R] Trouble phrasing an R command that will run the model i need (ANOVA, nested)

[David Elliott]

Hi,

I have been trying to find the appropriate R command to analyse my datasets 
according to a particular model. Unfortunately, my best attempts at doing so 
have so far not worked.

I am wondering if anybody can help me to figure out what i've been doing 
wrong, and what i need to do in order to use R correctly?

The model is an ANOVA with some crossed factors, interactions, and one 
nested factor. I have listed them here:

(Factors "foresttype" through to "experiment" are all at the same level; 
there is no nesting there, but the last factor listed is nested).

foresttype
region
foresttype:region
crosstype
foresttype:crosstype
experiment
linelabel (nested within every cell)

...and these are the appropriate F-test error terms:

foresttype -> tested against foresttype:crosstype interaction,
everything else -> tested against linelabel(nested)

I expect unequal sample sizes, so i have been trying to do it using Type III 
sums-of-squares. Now, that's the model that i'm trying to run in R. My 
attempts so far have looked like this:

>prod.anova <- aov(yproductivity ~ 
>foresttype+region+foresttype:region+crosstype+foresttype:crosstype+experiment+Error(linelabel), 
>data=Productivity)

After looking at the results of this model, my plan was to redo the model 
with the foresttype:crosstype interaction as the Error() term, to get the 
right F-test for "foresttype". Like so:

>prod2.anova <- aov(yproductivity ~ 
>foresttype+region+foresttype:region+crosstype+foresttype:crosstype+experiment+Error(foresttype:region), 
>data=Productivity)

To produce an output, i have used the Anova() command, and specified Type 
III SS:

>library(car)
>Anova(prod.anova, type=c("III"))

Unfortunately, i have run into a couple of problems with this approach. The 
first is that using the "Anova()" command produces this error message:

*
Error in Anova(prod.anova, type = c("III")) :
        no applicable method for "Anova"
*

I have looked through the manual (and various guides to using R for ANOVA), 
but haven't found anything that refers to it, let alone explains it. It 
looks (to me) like an extremely generic error message, which makes me think 
i am simply not using the "Anova()" command in the right way. But if the 
"Anova()" command is not appropriate, then i don't know how else to get the 
output that i want (and what else would allow me to specify Type III SS).

That's my first problem; here's a second one: As soon as i enter that second 
formula:

>prod2.anova <- aov(yproductivity ~ 
>foresttype+region+foresttype:region+crosstype+foresttype:crosstype+experiment+Error(foresttype:region), 
>data=Productivity)

, i get this error message:

*
Warning message:
Error() model is singular in: aov(yproductivity ~ foresttype + region + 
foresttype:region +
*

It looks to me as if R is saying that i have constructed the formula 
incorrectly. I noticied that i had "foresttype:crosstype" listed twice in 
the one formula, so i ran it again with "foresttype:crosstype" ommitted from 
the formula itself, but still in the "Error()" term. That produces the same 
error message as before.

I played around with it a bit, and found that i seemed to get this Error 
message when an interaction term is present in the "Error()" part of the 
formula. But i don't know why this should be, as i all the tutorials i have 
read have led me to beleive that R has no problem with an interaction term 
being used with the "Error()" command. And, aside from that, if i can't use 
an interaction term in "Error()", then how else can i obtain my F-test for 
"foresttype"?

I sought help for this, and it was suggested to me that i give up trying to 
use "aov()" at all, and do it all using the "lme()" command in the lme4 
package.

I downloaded "lme4", and read over this page: 
http://www.stat.ucl.ac.be/ISdidactique/Rhelp/library/lme4/html/lme.html , to 
try and get an idea of how to construct my "lme()" command.

My first try was:

>library(lme4)
>lme(yproductivity ~ foresttype+crosstype, data=productivity, yproductivity 
>~ region+crosstype+experiment)

(I excluded interactions and "Error()", as i wanted to just get a feel for 
what general formula worked before i made things more complex. My tactic was 
to write out the formula i'd used earlier, separating it into a formula 
containing only fixed factors, and a formula containing only random factors, 
in accordance with the page referenced above.) Unfortunately, doing so 
results in the error message:

*
Error: couldn't find function "lme"
*

So i couldn't even get far enough to see wether the formula i had 
constructed was laid out in a valid way.

I am really lost at this stage. I suspect that i am using all the wrong 
commands, in the wrong way. But i don't know how else to run this model; the 
commands i've tried have been my best attempts.

If anyone can provide any guidance, i would really appreciate it.

- David Elliott