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. ______________________________________________ R-help@r-project.org mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide http://www.R-project.org/posting-guide.html and provide commented, minimal, self-contained, reproducible code.