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

[Viechtbauer Wolfgang (STAT)]

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
[1] 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
> [1] "REML"
> 
> $convergence
> [1] TRUE
> 
> $iterations
> [1] 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
> [1] 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