To: r-sig-mixed-models@r-project.org<mailto:r-sig-mixed-models@r-project.org>

Subject: Re: [R-sig-ME] lmer residual variance estimate with prior

      weights     [SEC=UNCLASSIFIED]

Message-ID: <loom.20130219T035106-107@post.gmane.org<mailto:loom.20130219T035106-107@post.gmane.org>>

Content-Type: text/plain; charset=us-ascii



Ben Bolker writes:



>I have a sneaking suspicion that lmer doesn't scale the sum of the weights to 1 before doing calculations, as it arguably should.

>What happens if you manually scale the weights to 1 (which seems appropriate in this case)?



Thanks for the suggestion Ben. Yes often weights are scaled so that they sum to 1, however, in this case I think its appropriate to not scale the sample sizes to give weights that sum to 1 and I give my reasons after first investigating what happens in comparison when I do scale them to sum to 1. For lm( ) scaling them to 1 made no difference to the regression coefficients and their standard errors since the scale factor of the sum of the sample sizes was simply absorbed into the estimate of the "Residual standard error" which was "Residual standard error"(scaled weights)= "Residual standard error"(unscaled weights)/(sum of unscaled weights)^0.5 (see output). However, this was not the case for the lmer( ) and asreml( ) models with slope point estimate unchanged but its SE dramatically different (see below last two lmer outputs). So which do we choose scaled or unscaled? I argue in this application it is the unscaled weights that should be used but to argue this I define a more adequate mixed model since the random cell intercept model is a bit naive.



The response is the log10 of the mean density for each cell by year combination. The data is quite unbalanced at both levels with quite a few combinations of cell x year missing and sample sizes of trawl shots within cell x year combinations that are present ranging from 1 to 233. Each mean has a sampling error with variance inversely proportional to sample size and if we assume a common variance across cell x year combinations we can use the fit of the LMM to estimate this sampling variance (i.e. we could alternatively use the actual empirical variances but these are not available for sample size of 1 and other issues can cause problems like zeros in the raw shot-level data). However, the LMM sigma^2 (or 'units' level variance) is inflated by lack-of-fit introduced by fitting only a linear YEAR trend. The more adequate LMM (initially ignoring sample sizes and therefore not weighting) is therefore



>lme.01 <- lmer(log10Abund ~ YEAR.vcentred + (1|YEAR.fs) + (YEAR.vcentred|cell.f), data=data.tab.sci)



or equivalently in asreml



>asreml.02 <-asreml(fixed = log10Abund ~  YEAR.vcentred,

+        random=~YEAR.fs + us(pol(YEAR.vcentred,1)):cell.f,

+        data=data.tab.sci, na.method.Y = "exclude",

+        na.method.X = "exclude", maxiter=50)



Where YEAR.fs is a factor version of YEAR.v (i.e. the uncentred YEAR variable) and as a random term gives the random lack-of-fit departures from linearity and as well now a full random coefficients model is fitted.



Both lmer( ) and asreml( ) give similar results when these unweighted versions of the model were fitted though both struggle to fit this complex model with the correlation between random cell intercepts and slopes reaching its bound of 1.00 for asreml( ) and lmer( ) returns an "Na" for this parameter.

However, the weighted versions gave very different results between asreml( ) and lmer( ). Reducing model complexity by dropping the random slopes (and therefore additionally the correlation parameter) the fits behaved much better and gave similar comparisons with unweighted fits close to equivalent and weighted fits giving very different SEs (see below).



The unscaled sample sizes represent when combined with sigma^2 (and after accounting for the lack-of-fit random effect variance) estimates of the actual sampling variance of the mean densities therefore unscaled sample sizes should be used as weights.





So the story is still that when the weighted LMM is fitted asreml( ) and lmer( ) are giving very different estimates of the "residual standard error" (i.e. sigma) and the asreml( ) estimate is more in line with expectation given the results from lm( ). For the simplified LMM with unscaled weights the asreml( ) YEAR slope estimate is -0.012489 (SE=0.004935) while for lmer( ) the point estimate is the same but the SE is half as big with SE=0.001989.



Since this appears to be a significant problem with lmer( ) I have included relevant output. Sorry if its a bit longwinded.



>#scale weights to sum to 1



data.tab.sci$scvT_noshots <- data.tab.sci$vT_noshots/sum(data.tab.sci$vT_noshots)



> #lm( ) with scaled weights

> lm.02 <- lm(log10Abund ~ YEAR.vcentred, weights=data.tab.sci$scvT_noshots, data=data.tab.sci)

>

> summary(lm.02)



Call:

lm(formula = log10Abund ~ YEAR.vcentred, data = data.tab.sci,

    weights = data.tab.sci$scvT_noshots)



Weighted Residuals:

     Min       1Q   Median       3Q      Max

-0.19009 -0.03387 -0.01061  0.02823  0.11931



Coefficients:

               Estimate Std. Error t value Pr(>|t|)

(Intercept)    1.083599   0.053572   20.23   <2e-16 ***

YEAR.vcentred -0.002986   0.007107   -0.42    0.675

---

Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1



Residual standard error: 0.05357 on 104 degrees of freedom

Multiple R-squared: 0.001695,   Adjusted R-squared: -0.007904

F-statistic: 0.1766 on 1 and 104 DF,  p-value: 0.6752

>

> #lm( ) with unscaled weights



> lm.02 <- lm(log10Abund ~ YEAR.vcentred, weights=data.tab.sci$vT_noshots, data=data.tab.sci)

>

> summary(lm.02)



Call:

lm(formula = log10Abund ~ YEAR.vcentred, data = data.tab.sci,

    weights = data.tab.sci$vT_noshots)



Weighted Residuals:

     Min       1Q   Median       3Q      Max

-11.0579  -1.9706  -0.6171   1.6421   6.9404



Coefficients:

               Estimate Std. Error t value Pr(>|t|)

(Intercept)    1.083599   0.053572   20.23   <2e-16 ***

YEAR.vcentred -0.002986   0.007107   -0.42    0.675

---

Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1



Residual standard error: 3.116 on 104 degrees of freedom

Multiple R-squared: 0.001695,   Adjusted R-squared: -0.007904

F-statistic: 0.1766 on 1 and 104 DF,  p-value: 0.6752



> summary(lm.02)[["sigma"]]/((sum(data.tab.sci$vT_noshots))^0.5)

[1] 0.05357206



>#lmer( ) without weights



> lme.01 <- lmer(log10Abund ~ YEAR.vcentred + (1|YEAR.fs) + (1|cell.f), data=data.tab.sci)

>

> summary(lme.01)

Linear mixed model fit by REML

Formula: log10Abund ~ YEAR.vcentred + (1 | YEAR.fs) + (1 | cell.f)

   Data: data.tab.sci

 AIC   BIC logLik deviance REMLdev

282 295.4   -136    263.2     272

Random effects:

Groups   Name        Variance Std.Dev.

YEAR.fs  (Intercept) 0.086775 0.29458

 cell.f   (Intercept) 0.124647 0.35305

 Residual             0.579444 0.76121

Number of obs: 106, groups: YEAR.fs, 23; cell.f, 10



Fixed effects:

               Estimate Std. Error t value

(Intercept)    0.823531   0.157875   5.216

YEAR.vcentred -0.009147   0.012447  -0.735



Correlation of Fixed Effects:

            (Intr)

YEAR.vcntrd 0.203



>#asreml( ) without weights



> asreml.02 <-asreml(fixed = log10Abund ~  YEAR.vcentred,

+        random=~YEAR.fs + cell.f,

+        data=data.tab.sci, na.method.Y = "exclude",

+        na.method.X = "exclude", maxiter=50)



asreml(): 3.0.1 Library: 3.01gl IA32  Run: Wed Feb 20 11:37:36 2013



     LogLik         S2      DF

    -40.9675      0.6239   104  11:37:36

    -40.7186      0.6097   104  11:37:36

    -40.5176      0.5933   104  11:37:36

    -40.4491      0.5801   104  11:37:36

    -40.4483      0.5794   104  11:37:36

    -40.4483      0.5794   104  11:37:36



Finished on: Wed Feb 20 11:37:36 2013

LogLikelihood Converged

> summary(asreml.02)

$call

asreml(fixed = log10Abund ~ YEAR.vcentred, random = ~YEAR.fs +

    cell.f, data = data.tab.sci, na.method.Y = "exclude", na.method.X = "exclude",

    maxiter = 50)



$loglik

[1] -40.4483



$nedf

[1] 104



$sigma

[1] 0.7612125



$varcomp

                        gamma  component  std.error  z.ratio constraint

YEAR.fs!YEAR.fs.var 0.1497568 0.08677574 0.07182311 1.208187   Positive

cell.f!cell.f.var   0.2151150 0.12464720 0.08957924 1.391474   Positive

R!variance          1.0000000 0.57944445 0.09336550 6.206195   Positive



attr(,"class")

[1] "summary.asreml"

> anova(asreml.02)

Wald tests for fixed effects



Response: log10Abund



Terms added sequentially; adjusted for those above



              Df Sum of Sq Wald statistic Pr(Chisq)

(Intercept)    1   17.3929         30.017 4.284e-08 ***

YEAR.vcentred  1    0.3129          0.540    0.4624

residual (MS)       0.5794

---

Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

> summary(asreml.02)$varcomp

                        gamma  component  std.error  z.ratio constraint

YEAR.fs!YEAR.fs.var 0.1497568 0.08677574 0.07182311 1.208187   Positive

cell.f!cell.f.var   0.2151150 0.12464720 0.08957924 1.391474   Positive

R!variance          1.0000000 0.57944445 0.09336550 6.206195   Positive

> asreml.02$coefficients$fixed

YEAR.vcentred   (Intercept)

 -0.009146783   0.823530369

> ((asreml.02$vcoeff$fixed)^0.5)

[1] 0.01635185 0.20740526

> summary(asreml.02)$sigma

[1] 0.7612125



>#lmer( ) with unscaled weights



lme.01 <- lmer(log10Abund ~ YEAR.vcentred + (1|YEAR.fs) + (1|cell.f), weights=data.tab.sci$vT_noshots, data=data.tab.sci)

>

> summary(lme.01)

Linear mixed model fit by REML

Formula: log10Abund ~ YEAR.vcentred + (1 | YEAR.fs) + (1 | cell.f)

   Data: data.tab.sci

   AIC   BIC logLik deviance REMLdev

536.2 549.5 -263.1      517   526.2

Random effects:

Groups   Name        Variance  Std.Dev.

YEAR.fs  (Intercept) 0.0040445 0.063596

cell.f   (Intercept) 0.0027164 0.052119

Residual             0.1623334 0.402906

Number of obs: 106, groups: YEAR.fs, 23; cell.f, 10



Fixed effects:

               Estimate Std. Error t value

(Intercept)    0.908852   0.025576   35.54

YEAR.vcentred -0.012489   0.001989   -6.28



Correlation of Fixed Effects:

            (Intr)

YEAR.vcntrd 0.175



>#asreml( ) with unscaled weights



> asreml.02 <-asreml(fixed = log10Abund ~  YEAR.vcentred,

+        random=~YEAR.fs + cell.f, weights=vT_noshots,

+        data=data.tab.sci, na.method.Y = "exclude",

+        na.method.X = "exclude", maxiter=50)



asreml(): 3.0.1 Library: 3.01gl IA32  Run: Wed Feb 20 11:44:21 2013



     LogLik         S2      DF

    -35.9673      4.1526   104  11:44:21  (1 component(s) restrained)

    -30.9356      5.1675   104  11:44:21

    -29.5909      5.3096   104  11:44:21

    -29.2832      5.2321   104  11:44:21

    -29.2672      5.1850   104  11:44:21

    -29.2672      5.1824   104  11:44:21

    -29.2672      5.1824   104  11:44:21



Finished on: Wed Feb 20 11:44:21 2013

LogLikelihood Converged

>

>

>

> summary(asreml.02)

$call

asreml(fixed = log10Abund ~ YEAR.vcentred, random = ~YEAR.fs +

    cell.f, data = data.tab.sci, weights = vT_noshots, na.method.Y = "exclude",

    na.method.X = "exclude", maxiter = 50)



$loglik

[1] -29.26721



$nedf

[1] 104



$sigma

[1] 2.276491



$varcomp

                         gamma  component  std.error  z.ratio constraint

YEAR.fs!YEAR.fs.var 0.02491499 0.12911967 0.05500330 2.347490   Positive

cell.f!cell.f.var   0.01673371 0.08672093 0.06160439 1.407707   Positive

R!variance          1.00000000 5.18240914 0.84279187 6.149097   Positive



attr(,"class")

[1] "summary.asreml"

> anova(asreml.02)

Wald tests for fixed effects



Response: log10Abund



Terms added sequentially; adjusted for those above



              Df Sum of Sq Wald statistic Pr(Chisq)

(Intercept)    1   224.757         43.369 4.532e-11 ***

YEAR.vcentred  1     6.403          1.236    0.2663

residual (MS)        5.182

---

Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

> summary(asreml.02)$varcomp

                         gamma  component  std.error  z.ratio constraint

YEAR.fs!YEAR.fs.var 0.02491499 0.12911967 0.05500330 2.347490   Positive

cell.f!cell.f.var   0.01673371 0.08672093 0.06160439 1.407707   Positive

R!variance          1.00000000 5.18240914 0.84279187 6.149097   Positive

> asreml.02$coefficients$fixed

YEAR.vcentred   (Intercept)

  -0.01248908    0.90885295

> ((asreml.02$vcoeff$fixed)^0.5)

[1] 0.004935411 0.063477032

> summary(asreml.02)$sigma

[1] 2.276491



># lmer( ) with scaled weights

> lme.01 <- lmer(log10Abund ~ YEAR.vcentred + (1|YEAR.fs) + (1|cell.f), weights=data.tab.sci$scvT_noshots, data=data.tab.sci)

>

> summary(lme.01)

Linear mixed model fit by REML

Formula: log10Abund ~ YEAR.vcentred + (1 | YEAR.fs) + (1 | cell.f)

   Data: data.tab.sci

    AIC    BIC logLik deviance REMLdev

-325.2 -311.9  167.6   -344.5  -335.2

Random effects:

Groups   Name        Variance Std.Dev.

YEAR.fs  (Intercept) 13.68668 3.69955

 cell.f   (Intercept)  9.19242 3.03190

 Residual              0.16233 0.40291

Number of obs: 106, groups: YEAR.fs, 23; cell.f, 10



Fixed effects:

              Estimate Std. Error t value

(Intercept)    0.90886    1.48771   0.611

YEAR.vcentred -0.01249    0.11568  -0.108



Correlation of Fixed Effects:

            (Intr)

YEAR.vcntrd 0.175



>

> Hi mixed-modellers

> I am fitting a simple linear mixed model to some abundance data which

> are mean densities on the log10 scale for a set of spatial cells with

> sample sizes used as prior weights. I define a random cell intercept

> model and fit a linear year trend. I get very similar estimates of the

> intercept and slope when using each of lmer(.) and

> asreml(.) but get much smaller estimates of the SEs of these

> parameters for lmer(.). This is due to a much smaller estimate of the

> residual standard error with estimate of 0.534 for lmer and

> 3.016 for asreml with corresponding estimates of the cell-level

> standard deviation of 0.0014 for lmer and 0.2151 (=0.04626^0.5) for

> asreml. Comparing these to a simple lm(.) fit gives a similar but

> slightly higher estimate of the residual standard error compared to

> the asreml estimate with the a value of 3.116. The lmer estimate

> appears to be orders of magnitude out. Am I interpreting these results

> correctly? Has it something to do with how the weighting is done!





___________________________________________________________________________

    Australian Antarctic Division - Commonwealth of Australia
IMPORTANT: This transmission is intended for the addressee only. If you are not the
intended recipient, you are notified that use or dissemination of this communication is
strictly prohibited by Commonwealth law. If you have received this transmission in error,
please notify the sender immediately by e-mail or by telephoning +61 3 6232 3209 and
DELETE the message.
        Visit our web site at http://www.antarctica.gov.au/
___________________________________________________________________________

	[[alternative HTML version deleted]]

