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 |
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 |
Free forum by Nabble | Edit this page |