[R-sig-ME] lme function - fixed sigma - inconsistent results with sae and proc mixed results
Maciej Beresewicz
maciej.beresewicz at gmail.com
Tue Jun 28 18:16:14 CEST 2016
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
Pozdrawiam,
Maciej Beresewicz
[[alternative HTML version deleted]]
More information about the R-sig-mixed-models
mailing list