# [R-sig-ME] lme function - fixed sigma - inconsistent results with sae and proc mixed results

Viechtbauer Wolfgang (STAT) wolfgang.viechtbauer at maastrichtuniversity.nl
Tue Jun 28 17:29:01 CEST 2016

```Since when does lme() in R have the 'sigma' control argument? Ah, since 2015-11-25 (https://cran.r-project.org/web/packages/nlme/ChangeLog). Very interesting!

But apparently this is not giving the right results here. Another check:

library(sae)
library(metafor)
data(milk)
milk\$var <- milk\$SD^2
res <- rma(yi ~ as.factor(MajorArea), var, data=milk)
res
res\$tau2

Yields:

Mixed-Effects Model (k = 43; tau^2 estimator: REML)

tau^2 (estimated amount of residual heterogeneity):     0.0185 (SE = 0.0079)
tau (square root of estimated tau^2 value):             0.1362
I^2 (residual heterogeneity / unaccounted variability): 55.21%
H^2 (unaccounted variability / sampling variability):   2.23
R^2 (amount of heterogeneity accounted for):            65.85%

Test for Residual Heterogeneity:
QE(df = 39) = 86.1840, p-val < .0001

Test of Moderators (coefficient(s) 2,3,4):
QM(df = 3) = 46.5699, p-val < .0001

Model Results:

estimate      se     zval    pval    ci.lb    ci.ub
intrcpt                  0.9682  0.0694  13.9585  <.0001   0.8322   1.1041  ***
as.factor(MajorArea)2    0.1328  0.1030   1.2891  0.1974  -0.0691   0.3347
as.factor(MajorArea)3    0.2269  0.0923   2.4580  0.0140   0.0460   0.4079    *
as.factor(MajorArea)4   -0.2413  0.0816  -2.9565  0.0031  -0.4013  -0.0813   **

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

> res\$tau2
 0.01854996

Same as in 'sae' (rounded to 6 decimals) and SAS.

Not a huge difference to lme(), but larger than one would expect due to numerical differences.

If you switch to method="ML" (for both lme() and rma()), then you get 0.1245693^2 = 0.01551751 for lme() and 0.0155174 for rma(), so that's basically the same.

Best,
Wolfgang

--
Wolfgang Viechtbauer, Ph.D., Statistician | Department of Psychiatry and
Neuropsychology | Maastricht University | P.O. Box 616 (VIJV1) | 6200 MD
Maastricht, The Netherlands | +31 (43) 388-4170 | http://www.wvbauer.com

> -----Original Message-----
> From: R-sig-mixed-models [mailto:r-sig-mixed-models-bounces at r-
> project.org] On Behalf Of Maciej Beresewicz
> Sent: Tuesday, June 28, 2016 16:02
> To: r-sig-mixed-models at r-project.org
> Subject: [R-sig-ME] lme function - fixed sigma - inconsistent results
> with sae and proc mixed results
>
> I would like to estimate Fay-Herriot class models in nlme (small area
> models). Basically, these models have fixed random error which is assumed
> to be known from sample survey (sampling error). Hence, the model I would
> like to specify should have sigma = 1 (it is not estimated).
>
> I have checked new version nlme package (3.1-128) however results are
> different from those from sae package and proc mixed when it comes to
> fitting a small area model (in particular a Fay-Herriot area model).
>
> I am not sure why these results differ. They should be the same because
> sae::eblupFH fits s mixed model with one random effect and fixed residual
> variance).
>
> I would be grateful for any help on this matter. Please find the codes to
> highlight the problem below.
>
> ##############################
> ## preparation
> ##############################
>
> library(sae)
> library(nlme)
> data(milk)
> milk\$var <- milk\$SD^2
>
>
> ##############################
> ### nlme results
> ##############################
>
> > m2 <- lme(fixed = yi ~ as.factor(MajorArea),
>            random = ~ 1 | SmallArea,
>            control = lmeControl(sigma = 1,
>                                 apVar = T),
>            weights = varFixed(~var),
>            data = milk)
>
>
> ## variance (not the same as in sae and proc mixed)
> 0.1332918^2 = 0.0177667
>
> > summary(m2)
> Linear mixed-effects model fit by REML
>  Data: milk
>         AIC       BIC   logLik
>   -10.69175 -2.373943 10.34588
>
> Random effects:
>  Formula: ~1 | SmallArea
>         (Intercept) Residual
> StdDev:   0.1332918        1
>
> Variance function:
>  Structure: fixed weights
>  Formula: ~var
> Fixed effects: yi ~ as.factor(MajorArea)
>                            Value  Std.Error DF   t-value p-value
> (Intercept)            0.9680768 0.06849017 39 14.134537  0.0000
> as.factor(MajorArea)2  0.1316132 0.10183884 39  1.292367  0.2038
> as.factor(MajorArea)3  0.2269008 0.09126952 39  2.486053  0.0173
> as.factor(MajorArea)4 -0.2415905 0.08058755 39 -2.997863  0.0047
>
> ##############################
> ### results based on sae package
> ##############################
>
> library(sae)
> va <- eblupFH(formula = yi ~ as.factor(MajorArea), vardir = var, data =
> milk, method = "REML")
>
> > va\$fit
> \$method
>  "REML"
>
> \$convergence
>  TRUE
>
> \$iterations
>  4
>
> \$estcoef
>                             beta  std.error    tvalue       pvalue
> (Intercept)            0.9681890 0.06936208 13.958476 2.793443e-44
> as.factor(MajorArea)2  0.1327801 0.10300072  1.289119 1.973569e-01
> as.factor(MajorArea)3  0.2269462 0.09232981  2.457995 1.397151e-02
> as.factor(MajorArea)4 -0.2413011 0.08161707 -2.956503 3.111496e-03
>
> \$refvar
>  0.01855022
>
> \$goodness
>    loglike        AIC        BIC        KIC       AICc      AICb1
> AICb2       KICc
>  12.677478 -15.354956  -6.548956 -10.354956         NA         NA
> NA         NA
>      KICb1      KICb2 nBootstrap
>         NA         NA   0.000000
>
> ##############################
> ### Proc mixed results are consistent with sae
> ##############################
>
> proc SmallArea data= milk order=data method=reml;
> class SmallArea;
> weight var;
> model y= MajorArea2 MajorArea3  MajorArea4 / cl solution outp=predicted;
> random SmallArea;
> parms (1) (1) / hold=2;
> run;
>
> ## Covariance Parameter Estimates
> Cov Parm Estimate
> SmallArea 0.01855
> Residual 1.0000
>
> ### fixed effects
> Intercept 0.9682
> majorarea2 0.1328
> majorarea3 0.2269
> majorarea3 -0.2413
>
> Best regards,
> Maciej

```