Re: [R] lme cant get parameter estimated correctly

From: <Bill.Venables_at_csiro.au>
Date: Sun, 06 Apr 2008 16:50:37 +1000

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.

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. Received on Sun 06 Apr 2008 - 06:53:46 GMT

Archive maintained by Robert King, hosted by the discipline of statistics at the University of Newcastle, Australia.
Archive generated by hypermail 2.2.0, at Sun 06 Apr 2008 - 19:30:27 GMT.

Mailing list information is available at https://stat.ethz.ch/mailman/listinfo/r-help. Please read the posting guide before posting to the list.

list of date sections of archive