[R] lme cant get parameter estimated correctly

[Bill.Venables at csiro.au]

Here is a demo you may like to consider.  (I can see what you are trying
to do with your loops, but I prefer to do it this way.)

On 32 bit Windows, (which I am forced to use), your seed is not a valid
integer, so I have changed it to something which is.

> set.seed(7658943)
> 
> fph <- 0.4
> Sigh <- sqrt(0.0002)
> Sigi <- sqrt(0.04)
> 
> reH <- rnorm(90, fph, Sigh)  ## hospid effects
> dta <- within(expand.grid(hospid = 1:90, empid = 1:80),
            fpi1 <- reH[hospid] + rnorm(7200, fph, Sigi))
> 
> require(nlme)
> 
> Modl <- lme(fpi1 ~ 1, data = dta, random = ~1|hospid)
> summary(Modl)
Linear mixed-effects model fit by REML

.........

Random effects:
 Formula: ~1 | hospid
        (Intercept)  Residual
StdDev:  0.01593939 0.2003148

.........

Compare this with what you started with, viz:

> c(Sigh, Sigi)
[1] 0.01414214 0.20000000

and it doesn't look too bad to me.  You might like to see how well the
effects themselves agree:

> plot(ranef(Modl)[[1]], reH-fph)
> abline(0,1)

> cor(ranef(Modl)[[1]], reH)
[1] 0.6317546

Again, about what I would expect.

Note that the summary to lme objects gives standard deviations, not
variance components.  If you want to see variance components, there is a
function in the 'ape' package which does it:

> require(ape)
> varcomp(Modl)
      hospid       Within 
0.0002540641 0.0401260303 

Look familiar?

Bill Venables
CSIRO Laboratories
PO Box 120, Cleveland, 4163
AUSTRALIA
Office Phone (email preferred): +61 7 3826 7251
Fax (if absolutely necessary):  +61 7 3826 7304
Mobile:                         +61 4 8819 4402
Home Phone:                     +61 7 3286 7700
mailto:Bill.Venables at csiro.au
http://www.cmis.csiro.au/bill.venables/ 

-----Original Message-----
From: r-help-bounces at r-project.org [mailto:r-help-bounces at r-project.org]
On Behalf Of toby909
Sent: Sunday, 6 April 2008 3:27 PM
To: r-help at stat.math.ethz.ch
Subject: [R] lme cant get parameter estimated correctly


I am caught in a mental trap. Why isn't the between groups variance
estimated
(0.0038) to be around the value with which I generated the data
(0.0002)?

Thanks Toby






set.seed(76589437887)

fph = 0.4

Sigh = sqrt(0.0002)
Sigi = sqrt(0.04)

ci = 1
fpi = matrix(,7200,3)
for (i in 1:90) {
 fph = rnorm(1, fph, Sigh)
 for (k in 1:80) {
  fpi[ci,1:3] = matrix(c(i, k, rnorm(1, fph, Sigi)),1)
  ci = ci+1
 }
}

colnames(fpi) = c("hospid", "empid", "fpi1")
dta = as.data.frame(fpi)



lme = lme(fpi1 ~ 1, dta, ~1|hospid)
summary(lme)






lme = lme(fpi1 ~ 1, dta, ~1|hospid)
summary(lme)
Linear mixed-effects model fit by REML
 Data: dta 
        AIC       BIC   logLik
  -2555.416 -2534.771 1280.708

Random effects:
 Formula: ~1 | hospid
        (Intercept)  Residual
StdDev:  0.06173257 0.1997302

Fixed effects: fpi1 ~ 1 
               Value   Std.Error   DF  t-value p-value
(Intercept) 0.368082 0.006919828 7110 53.19236       0

Standardized Within-Group Residuals:
       Min         Q1        Med         Q3        Max 
-3.4870696 -0.6747173 -0.0048658  0.6838012  4.2633384 

Number of Observations: 7200
Number of Groups: 90 

0.0617^2
[1] 0.00380689

______________________________________________
R-help at 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.