[R] lme: NaNs in hierarchical model

[melanie.r]

Dear R-users,
I'm somehow new to the topic of mixed models. I'd liked to simulate data
with a specific correlation-structure and analyze it with lme(). I have a
fixed effect (group) and a random effect (block):

            y    group   block
 -5.2158969    -2     1
 -2.4174197    -2     1
 -2.6731572    -2     1
 -3.6310082    -2     1
 -0.3888706     0     2
  1.3422255     0     2
  0.8303369     0     2
  1.0985184     0     2
  2.4403000     2     3
  2.0890006     2     3
  2.8849674     2     3
  3.1946739     2     3

This is just a short data.frame, to illustrate my model. My mentor called it
'hierarchical' model, because every group is in one block and only the
observations of one group are correlated. 

When I now try to analyze this using lme with the following model

fit<-lme(y~group,random=~1|block)

there is an error message:
In pf(q, df1, df2, lower.tail, log.p) : NaNs produced

In fact I generate a certain number of that data and analyze it with 
anova(lme()) to get the p-values and then count them to get the power of the
test.
What I don't understand is, why are there NaNs? Or better: why is the
group-DF zero?
The output is the following:

>fit
Linear mixed-effects model fit by REML
  Data: NULL 
  Log-restricted-likelihood: -13.90713
  Fixed: y ~ group 
(Intercept)      group0      group2 
  -3.484370    4.204923    6.136606 

Random effects:
 Formula: ~1 | block
        (Intercept)  Residual
StdDev:    1.200732 0.9005491

Number of Observations: 12
Number of Groups: 3 

> anova(fit)
            numDF denDF  F-value p-value
(Intercept)     1     9 0.002524   0.961
group           2     0 5.986678     NaN
Warning message:
In pf(q, df1, df2, lower.tail, log.p) : NaNs produced


When I do the same with a sort of 'crossed' design

         y      group block
 -1.8892392    -2     1
  0.4329959     0     1
  3.2474885     2     1
 -3.7710187    -2     2
 -1.6624300     0     2
  1.2584779     2     2
 -1.6813929    -2     3
  1.6828940     0     3
  2.2041250     2     3
 -1.2136863    -2     4
  0.7032238     0     4
  0.9012617     2     4

the call to lme() works just fine:
lme(y~group,random=~1|block)
then anova() gives me the p-value and I'm happy=)

For my issue being to simulate the power I need to have the p-value. I tried
it with lmer():
anova(lmer(y~group+(1|block)))
which works, but with no degree of freedom I don't see myself getting a
critical F value that leads to my desired power. Or am I just wrong and
there's an easy way to compute it? Or is my model formula for the
hierarchical model wrong?
So that really confuses me and with all those things I read today concerning
this problem I'm not less confused.
The problem is: I need to get it solved very quickly. So, please, if there's
someone who can help me... I'd be so grateful.

mel
-- 
View this message in context: http://www.nabble.com/lme%3A-NaNs-in-hierarchical-model-tp19578256p19578256.html
Sent from the R help mailing list archive at Nabble.com.