issues with weights in glmer (or glmmADMB) [SEC=UNCLASSIFIED]

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