[ADMB Users] example

[H. Skaug]

Dave,


To check out your claim I compared to

lme(y ~ x, random=~1|fac,method="ML")

(the old mixed model function in R, tailor made for linear mixed models),
and got 3.829014 for the standard deviation of the random effect,
which is the same
as what you got. Note however, that these are ML estimates,
while what lmer and lme calculate by default are REML estimates.
I modified your ADMB program to calculate REML estimates and got
sigma_u=4.42137953753, which agrees with what I get from

lme(y ~ x, random=~1|fac,method="REML")

Hence, it seems that the optimization method in lmer() is not accurate enough
in this case, which is also indicated by the warning
"In mer_finalize(ans) : false convergence (8)".


I think this is a good example of why it is important to have a
good optimizer built into the package, and I will include it
in the ADMB example collection.


Hans







On Fri, Apr 23, 2010 at 6:03 PM, dave fournier <otter at otter-rsch.com> wrote:
> I see that at the ADMB meeting there was discussion about examples etc.
> Here is an example from the R list.  It is supposed to show that a
> linear mixed model
> with only 4 levels of the random effect is a bad idea, ostensibly
> because  lmer
> does not perform well.  You can find the discussion at.
>
> https://stat.ethz.ch/pipermail/r-sig-mixed-models/2010q2/003713.html
>
> I think the main point is at
> Random effects:
>  Groups   Name        Variance Std.Dev.
>  fac      (Intercept) 1.1162   1.0565
>  Residual             1.0298   1.0148
>
> where the estimate for the std dev of the random effects is 1.0565 while
> the value used for simulating the data was 4.0.
>
> I suspect that the real problem is with the  mediocre estimation methods
> used in lmer. I built the model in ADMBRE
> and estimated the parameters.  The estimate for sigma_u  was 3.8290
>  4   sigma_u  3.8290e+00 1.3538e+00
>
> Here is the TPL file.  Note that I sorted the R ouput data by factor
> value i.e by random effect index.
>
>
> DATA_SECTION
>  init_int n
>  init_int m
>  int nm
>  !! nm=n*m;
>  vector y
>  vector x
>  init_matrix data(1,nm,1,3)
>  !! y=column(data,1);
>  !! x=column(data,2);
> PARAMETER_SECTION
>  init_number a
>  init_number b
>  init_bounded_number sigma(.1,10.)
>  init_bounded_number sigma_u(.1,10.)
>  random_effects_vector u(1,m)
>  objective_function_value f
> PROCEDURE_SECTION
>  int ii=0;
>  for (int i=1;i<=m;i++)
>  {
>    f1(a,b,u(i),sigma,sigma_u,ii);
>  }
>
> SEPARABLE_FUNCTION void f1(const prevariable& a,const prevariable&
> b,const prevariable& ui,const prevariable& sigma,const prevariable&
> sigma_u,int& ii)
>
>  dvariable v=square(sigma);
>  f+=0.5*square(ui);
>  for (int j=1;j<=n;j++)
>  {
>    ii++;
>    dvariable pred=a+b*x(ii)+sigma_u*ui;
>    dvariable diff=y(ii)-pred;
>    f+=log(sigma)+0.5*square(diff)/v;
>  }
>
> I wonder if the R experts are interested?
>
>
>
>
>
>
> _______________________________________________
> Users mailing list
> Users at admb-project.org
> http://lists.admb-project.org/mailman/listinfo/users
>