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.
