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

[Maciej Beręsewicz]

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

	[[alternative HTML version deleted]]