Lmer with weights

classic Classic list List threaded Threaded
2 messages Options
Reply | Threaded
Open this post in threaded view
|

Lmer with weights

Gregor Gorjanc
Hello!

I would like to use lmer() to fit data, which are some estimates and
their standard errors i.e kind of a "meta" analysis. I wonder if weights
argument is the right one to use to include uncertainty (standard
errors) of "data" into the model. I would like to use lmer(), since I
would like to have a "freedom" in modeling, if this is at all possible.

For example we can take schools data by Gelman from R2WinBUGS package.
As you can see bellow use of weights argument did not had influence on
results.

I do not know if my specification of weights i.e. 1 / sd^2 is ok. Under
least squares one minimizes sum(e^2_i) or sum(w_i * e^2_i) with weighted
LS. If I consider that \sigma_i represents uncertainty in my "data" then
e'_i = e_i / \sigma_i and we minimize sum(e'^2_i) = sum((e_i /
\sigma_i)^2) = sum(e_i * \sigma^-2_i). Therefore weights i.e. w_i are
equal to 1 / \sigma^2_i.

Can anyone help me with this issue?

Thank you very much!

 > library("R2WinBUGS")
 > data(schools)
 > schools
 > attach(schools)
 >
 > ## Fit simple model without "weights"
 > lmer(estimate ~ 1 + (1 | school))
Linear mixed-effects model fit by REML
Formula: estimate ~ 1 + (1 | school)
     AIC    BIC  logLik MLdeviance REMLdeviance
  58.882 59.041 -27.441     59.278       54.882
Random effects:
  Groups   Name        Variance Std.Dev.
  school   (Intercept) 80.4     8.97
  Residual             30.1     5.49
# of obs: 8, groups: school, 8

Fixed effects:
             Estimate Std. Error t value
(Intercept)     8.82       3.72    2.37

 > ## Fit simple model with "weights"
 > lmer(estimate ~ 1 + (1 | school), weights = ~ 1 / (sd^2))
Linear mixed-effects model fit by REML
Formula: estimate ~ 1 + (1 | school)
     AIC    BIC  logLik MLdeviance REMLdeviance
  58.882 59.041 -27.441     59.278       54.882
Random effects:
  Groups   Name        Variance Std.Dev.
  school   (Intercept) 80.4     8.97
  Residual             30.1     5.49
# of obs: 8, groups: school, 8

Fixed effects:
             Estimate Std. Error t value
(Intercept)     8.82       3.72    2.37

--
Lep pozdrav / With regards,
     Gregor Gorjanc

----------------------------------------------------------------------
University of Ljubljana     PhD student
Biotechnical Faculty
Zootechnical Department     URI: http://www.bfro.uni-lj.si/MR/ggorjan
Groblje 3                   mail: gregor.gorjanc <at> bfro.uni-lj.si

SI-1230 Domzale             tel: +386 (0)1 72 17 861
Slovenia, Europe            fax: +386 (0)1 72 17 888

----------------------------------------------------------------------
"One must learn by doing the thing; for though you think you know it,
  you have no certainty until you try." Sophocles ~ 450 B.C.

______________________________________________
[hidden email] mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide! http://www.R-project.org/posting-guide.html
Reply | Threaded
Open this post in threaded view
|

Re: Lmer with weights

Spencer Graves
          Will your data support using "lme" in the 'nlme' package?  If yes, I
suggest you switch.  This is consistent with a response given recently
by Doug Bates to a crudely related question
(http://finzi.psych.upenn.edu/R/Rhelp02a/archive/68340.html).

          You probably know that "lmer" and associated software are Doug Bates'
development platform.  There have been recent email exchanges on
problems with weights in "lmer".  RSiteSearch("lmer with weights") and
RSiteSearch("weights in lmer") both preduced the same 30 hits.  On 31
Jan. 2006, Patrick Connolly reported, "I suspect the weights argument is
not having any effect."
(http://finzi.psych.upenn.edu/R/Rhelp02a/archive/69397.html)  I couldn't
find any replies in the R-help archives, but my reply of 3 Feb. is
copied below.  If you would like further help from this list, please
submit another question.

          hope this helps.
          spencer graves

##############################
          I agree:  The lmer weights argument seems not to have any
effect.  To check this, I modified the first example in the "lmer"
documentation as follows:

Sleep <- sleepstudy
Sleep$wts <- 1:180
(fm1 <- lmer(Reaction ~ Days + (Days|Subject), Sleep))
(fm1w <- lmer(Reaction ~ Days + (Days|Subject),
                      weights=wts, Sleep))

          The numbers from both seemed to be the same.  To try to help
diagnose this, I listed "lmer", and found that it consisted of a call to
"standardGeneric".  Then 'getMethods("lmer")' listed only one "method"
for the case where the argument "formula" had class "formula".  I tried
to trace this further, e.g., by giving it a different name and using
"debug".  After being stopped a couple of time by functions hidden in
the "Matrix" namespace, I gave up.
       
############################
Gregor Gorjanc wrote:

> Hello!
>
> I would like to use lmer() to fit data, which are some estimates and
> their standard errors i.e kind of a "meta" analysis. I wonder if weights
> argument is the right one to use to include uncertainty (standard
> errors) of "data" into the model. I would like to use lmer(), since I
> would like to have a "freedom" in modeling, if this is at all possible.
>
> For example we can take schools data by Gelman from R2WinBUGS package.
> As you can see bellow use of weights argument did not had influence on
> results.
>
> I do not know if my specification of weights i.e. 1 / sd^2 is ok. Under
> least squares one minimizes sum(e^2_i) or sum(w_i * e^2_i) with weighted
> LS. If I consider that \sigma_i represents uncertainty in my "data" then
> e'_i = e_i / \sigma_i and we minimize sum(e'^2_i) = sum((e_i /
> \sigma_i)^2) = sum(e_i * \sigma^-2_i). Therefore weights i.e. w_i are
> equal to 1 / \sigma^2_i.
>
> Can anyone help me with this issue?
>
> Thank you very much!
>
>  > library("R2WinBUGS")
>  > data(schools)
>  > schools
>  > attach(schools)
>  >
>  > ## Fit simple model without "weights"
>  > lmer(estimate ~ 1 + (1 | school))
> Linear mixed-effects model fit by REML
> Formula: estimate ~ 1 + (1 | school)
>      AIC    BIC  logLik MLdeviance REMLdeviance
>   58.882 59.041 -27.441     59.278       54.882
> Random effects:
>   Groups   Name        Variance Std.Dev.
>   school   (Intercept) 80.4     8.97
>   Residual             30.1     5.49
> # of obs: 8, groups: school, 8
>
> Fixed effects:
>              Estimate Std. Error t value
> (Intercept)     8.82       3.72    2.37
>
>  > ## Fit simple model with "weights"
>  > lmer(estimate ~ 1 + (1 | school), weights = ~ 1 / (sd^2))
> Linear mixed-effects model fit by REML
> Formula: estimate ~ 1 + (1 | school)
>      AIC    BIC  logLik MLdeviance REMLdeviance
>   58.882 59.041 -27.441     59.278       54.882
> Random effects:
>   Groups   Name        Variance Std.Dev.
>   school   (Intercept) 80.4     8.97
>   Residual             30.1     5.49
> # of obs: 8, groups: school, 8
>
> Fixed effects:
>              Estimate Std. Error t value
> (Intercept)     8.82       3.72    2.37
>

______________________________________________
[hidden email] mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide! http://www.R-project.org/posting-guide.html