# [R] lme: NaNs in hierarchical model

melanie.r melanieraeder at gmx.net
Fri Sep 19 21:30:50 CEST 2008

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.