lme function - fixed sigma - inconsistent results with sae and proc mixed results

Tue, Jun 28, 2016 10:14 AM

Correct. Actually, I was very surprised to hear about the 'sigma' control argument being availble in R. That used to be only available in the S-Plus version of lme(). In fact, I just tried running that model in S-Plus and this is what I got:

Linear mixed-effects model fit by REML
 Data: dat 
      AIC      BIC   logLik 
  6.02487 14.34268 1.987565

Random effects:
 Formula:  ~ 1 | SmallArea
        (Intercept) Residual 
StdDev:   0.1361994        1

Variance function:
 Structure: fixed weights
 Formula:  ~ var 
Fixed effects: yi ~ as.factor(MajorArea) 
                           Value  Std.Error DF   t-value p-value 
          (Intercept)  0.9977953 0.03179340 39  31.38373  <.0001
as.factor(MajorArea)1  0.0663901 0.05150041 39   1.28912  0.2050
as.factor(MajorArea)2  0.0535187 0.02659572 39   2.01230  0.0511
as.factor(MajorArea)3 -0.0903025 0.01466646 39  -6.15707  <.0001
 Correlation: 
                      (Intr) a(MA)1 a(MA)2 
as.factor(MajorArea)1  0.075              
as.factor(MajorArea)2 -0.157 -0.060       
as.factor(MajorArea)3 -0.392 -0.054  0.113

Standardized Within-Group Residuals:
       Min         Q1       Med        Q3     Max 
 -1.702152 -0.2304439 0.2200099 0.3981538 1.22415

Number of Observations: 43
Number of Groups: 43

[1] 0.01855028

So that matches what we get from the other packages.

I am curious -- why do you want to use lme() for fitting F-H models anyway? Why not stick to 'sae'?

Best,
Wolfgang

-----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 18:16
To: r-sig-mixed-models at r-project.org
Subject: Re: [R-sig-ME] lme function - fixed sigma - inconsistent results
with sae and proc mixed results

Dear Viechtbauer,

Thanks! So it means that there is difference in terms of REML estimation.

On 28 Jun 2016, at 17:29, Viechtbauer Wolfgang (STAT)

<wolfgang.viechtbauer at maastrichtuniversity.nl> wrote:

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