From: <Bill.Venables_at_csiro.au>

Date: Sun, 06 Apr 2008 16:50:37 +1000

http://www.cmis.csiro.au/bill.venables/

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

Date: Sun, 06 Apr 2008 16:50:37 +1000

*> 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 7700mailto: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.
*