[R-sig-ME] issues with weights in glmer (or glmmADMB) [SEC=UNCLASSIFIED]
Steve Candy
Steve.Candy at aad.gov.au
Tue Jul 9 03:27:34 CEST 2013
Ben Bolker <bbolker at ...> writes:
>
> Ben Bolker <bbolker <at> ...> writes:
>
> >
> > Steve Candy <Steve.Candy <at> ...> writes:
> >
> > >
> > > I am yet to see any response from the developers/maintainers of
> > > lme4 with respect to my posting from 25
> > > Feburary on the R-sig-mixed-models Digest, Vol 74, Issue 40
> > >
> > > 5. Re: lmer residual variance estimate with prior weights (lmer
>
> [snip]
> >
> > Sorry this slipped through the cracks. We will take a look;
> > if you can (re)post the query at
> > https://github.com/lme4/lme4/issues?state=open
> > (the new site for bug/feature requests) that would help us out ...
>
> PS I looked back at your previous post, and you don't provide a
> reproducible example. If you can either provide your data or (very
> much preferably) a minimal reproducible example
> <http://tinyurl.com/reproducible-000>, that will save us time and
> greatly increase the probability that we will be able to tackle this
> in the near future ...
>
> cheers
> Ben Bolker
>
>
I have set up this small simulation to compare lmer{lme4} with lme{nlme}
(with results compared also to asreml{asreml})
Note that lme and asreml give identical results for SEs of regression
parameters while lmer SEs are an order of magnitude smaller.
I also give an asreml model that teases out the lev1 and lev2 variance
components given a known lev1 variance.
## sample size weighted LMM
#
#
library(nlme)
library(lme4)
#generate some lognormal random effects (i.e. level-2)
Nsamp <- rep(0,50)
REs <- rnorm(n=50, mean=0, sd=0.3)
#generate some within residual errors (i.e. level-1)
RES <- NULL
Xf <- NULL
X <- NULL
for (i in 1:50) {
Nsamp[i] <- rpois(n=1,lambda=20)
RES <- c(RES,rep(REs[i],each=Nsamp[i]))
Xi <- as.integer(i/10)
X <- c(X,rep(Xi,each=Nsamp[i]))
Xf <- c(Xf,rep(i,each=Nsamp[i]))}
RES <- RES + rnorm(n=length(RES), mean=0, sd=0.1)
RE.Xf <- as.factor(Xf)
#set up regression model
Y <- 1 + 5*X + RES
# calculate means and fit weighted LMM with random lack-of-fit term
Ym <- as.vector(tapply(X=Y, INDEX=RE.Xf, FUN=mean))
Xm <- as.vector(tapply(X=X, INDEX=RE.Xf, FUN=mean))
Xm.f <- as.factor(Xm)
wts <- 1/Nsamp
nlme.av <- lme(fixed=Ym ~ Xm, random=~1|Xm.f, weights=varFixed(~wts))
summary(nlme.av)
lmer.av <- lmer(formula=Ym ~ Xm + (1|Xm.f), weights=Nsamp)
summary(lmer.av)
> ## sample size weighted LMM
> #
> #
> library(nlme)
> library(lme4)
> #generate some lognormal random effects (i.e. level-2)
> Nsamp <- rep(0,50)
> REs <- rnorm(n=50, mean=0, sd=0.3)
> #generate some within residual errors (i.e. level-1)
> RES <- NULL
> Xf <- NULL
> X <- NULL
> for (i in 1:50) {
+ Nsamp[i] <- rpois(n=1,lambda=20)
+ RES <- c(RES,rep(REs[i],each=Nsamp[i]))
+ Xi <- as.integer(i/10)
+ X <- c(X,rep(Xi,each=Nsamp[i]))
+ Xf <- c(Xf,rep(i,each=Nsamp[i]))}
>
> RES <- RES + rnorm(n=length(RES), mean=0, sd=0.1)
> RE.Xf <- as.factor(Xf)
> Y <- 1 + 5*X + RES
>
> # calculate means and fit sample size weighted LMM with random lack of
fit term
> Ym <- as.vector(tapply(X=Y, INDEX=RE.Xf, FUN=mean))
> Xm <- as.vector(tapply(X=X, INDEX=RE.Xf, FUN=mean))
>
> Xm.f <- as.factor(Xm)
> wts <- 1/Nsamp
>
> nlme.av <- lme(fixed=Ym ~ Xm, random=~1|Xm.f, weights=varFixed(~wts))
> summary(nlme.av)
Linear mixed-effects model fit by REML
Data: NULL
AIC BIC logLik
48.867 56.35181 -20.4335
Random effects:
Formula: ~1 | Xm.f
(Intercept) Residual
StdDev: 6.521819e-06 1.476964
Variance function:
Structure: fixed weights
Formula: ~wts
Fixed effects: Ym ~ Xm
Value Std.Error DF t-value p-value
(Intercept) 0.876380 0.08709383 44 10.06248 0
Xm 5.036656 0.03275373 4 153.77353 0
Correlation:
(Intr)
Xm -0.84
Standardized Within-Group Residuals:
Min Q1 Med Q3 Max
-1.7507938 -0.7926107 0.1400055 0.6415305 2.0951252
Number of Observations: 50
Number of Groups: 6
>
> lmer.av <- lmer(formula=Ym ~ Xm + (1|Xm.f), weights=Nsamp)
> summary(lmer.av)
Linear mixed model fit by REML
Formula: Ym ~ Xm + (1 | Xm.f)
AIC BIC logLik deviance REMLdev
196.2 203.8 -94.08 178.9 188.2
Random effects:
Groups Name Variance Std.Dev.
Xm.f (Intercept) 0.00000 0.00000
Residual 0.11152 0.33395
Number of obs: 50, groups: Xm.f, 6
Fixed effects:
Estimate Std. Error t value
(Intercept) 0.876380 0.019693 44.5
Xm 5.036656 0.007406 680.1
Correlation of Fixed Effects:
(Intr)
Xm -0.840
>
> library(asreml)
> asreml.av <- asreml(fixed = Ym ~ Xm, random=~Xm.f, weights=Nsamp)
asreml(): 3.0.1 Library: 3.01gl X86_64 Run: Tue Jul 09 10:34:36 2013
LogLik S2 DF
18.8188 2.1432 48 10:34:36 (1 component(s) restrained)
22.5765 2.1637 48 10:34:36 (1 component(s) restrained)
23.5813 2.1793 48 10:34:36 (1 component(s) restrained)
23.6694 2.1813 48 10:34:36 (1 component(s) restrained)
23.6752 2.1814 48 10:34:36 (1 component(s) restrained)
23.6755 2.1814 48 10:34:36
23.6755 2.1814 48 10:34:36
23.6755 2.1814 48 10:34:36
Finished on: Tue Jul 09 10:34:36 2013
LogLikelihood Converged
>
> summary(asreml.av)
$call
asreml(fixed = Ym ~ Xm, random = ~Xm.f, weights = Nsamp)
$loglik
[1] 23.67552
$nedf
[1] 48
$sigma
[1] 1.476964
$varcomp
gamma component std.error z.ratio constraint
Xm.f!Xm.f.var 1.011929e-07 2.207444e-07 4.505927e-08 4.898979 Boundary
R!variance 1.000000e+00 2.181422e+00 4.452809e-01 4.898979 Positive
attr(,"class")
[1] "summary.asreml"
>
> asreml.av$coefficients$fixed
Xm (Intercept)
5.0366561 0.8763796
> (summary(asreml.av))$sigma*((asreml.av$vcoeff$fixed)^0.5)
[1] 0.03275401 0.08709455
> (summary(asreml.av)$varcomp[["component"]][1])^0.5
[1] 0.0004698345
> (summary(asreml.av)$varcomp[["component"]][2])^0.5
[1] 1.476964
>
> # tease apart lev1 and lev2 variances using a fixed "units" level known
variance
>
> RE.Xfs <- factor(x=seq(1,50), levels=c(1:50))
> asreml.av <- asreml(fixed = Ym ~ Xm, random=~Xm.f+RE.Xfs,
family=asreml.gaussian(dispersion=0.01), weights=Nsamp)
asreml(): 3.0.1 Library: 3.01gl X86_64 Run: Tue Jul 09 10:45:14 2013
LogLik S2 DF
19.5841 0.0100 48 10:45:14 (1 component(s) restrained)
22.5914 0.0100 48 10:45:14 (1 component(s) restrained)
23.1395 0.0100 48 10:45:14 (1 component(s) restrained)
23.1885 0.0100 48 10:45:14 (1 component(s) restrained)
23.1916 0.0100 48 10:45:14 (1 component(s) restrained)
23.1918 0.0100 48 10:45:14
23.1918 0.0100 48 10:45:14
23.1918 0.0100 48 10:45:14
Finished on: Tue Jul 09 10:45:14 2013
LogLikelihood Converged
>
> summary(asreml.av)
$call
asreml(fixed = Ym ~ Xm, random = ~Xm.f + RE.Xfs, family = asreml.gaussian
(dispersion = 0.01),
weights = Nsamp)
$loglik
[1] 23.19179
$nedf
[1] 48
$sigma
[1] 0.1
$varcomp
gamma component std.error z.ratio constraint
Xm.f!Xm.f.var 1.011929e-07 1.011929e-07 NA NA Boundary
RE.Xfs!RE.Xfs.var 1.165194e-01 1.165194e-01 0.02389939 4.875413 Positive
R!variance 1.000000e-02 1.000000e-02 NA NA Fixed
attr(,"class")
[1] "summary.asreml"
> (summary(asreml.av)$varcomp[["component"]][1])^0.5
[1] 0.0003181083
> (summary(asreml.av)$varcomp[["component"]][2])^0.5
[1] 0.3413494
> (summary(asreml.av)$varcomp[["component"]][3])^0.5
[1] 0.1
> ((asreml.av$vcoeff$fixed)^0.5)
[1] 0.03346914 0.08534932
> asreml.av$coefficients$fixed
Xm (Intercept)
5.0336711 0.8798403
>
Dr Steven G. Candy
Senior Applied Statistician
Wildlife Conservation and Fisheries
Australian Antarctic Division
203 Channel Hwy, Kingston, Tasmania, 7050, Australia
ph: (+61) 3 62 323135
More information about the R-sig-mixed-models
mailing list