[R] Lmer with weights

Spencer Graves spencer.graves at pdf.com
Tue Feb 14 05:01:21 CET 2006


	  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
>




More information about the R-help mailing list