[R-sig-ME] lme4- lmer2: wrong fixed effect estimate for some parameterization

[Douglas Bates]

On 2/6/07, Ulrich Halekoh <Ulrich.Halekoh at agrsci.dk> wrote:
> Hej,
>
> The lmer2-version (not the lmer)  calculates wrong fixed effect
> estimates if
> a model is fitted  with a factor-effects in the mean-value
> parameterization:

Yes, you are correct.  Thank you for sending a bug report with a
reproducible example.

The CHOLMOD sparse matrix library has managed to outsmart me and
produce a permutation that does not follow the pattern that I require.
 I enclose a bit of output with the Cholesky decomposition at the
parameter estimates for the g2.withOutIntercept fit.  The permutation
chosen by CHOLMOD intermingles columns from the random effects with
those from the fixed effects and that shouldn't happen.  I thought I
had prevented that from happening but apparently I didn't.  I should
be able to upload a patched version of the lme4 package later today.


>
>
> Example
> library(lme4)
>
> set.seed(98)
>
> ### data set generation of a hiearchial model with
> #  one random components sub
> ####
> #  y_ijk= a_i + delta_j(i) +   e_ijk
>
>
> n.treat<-2
> n.sub<-3
> n.rep<-4
>
> d<-expand.grid(rep=1:n.rep,sub=1:n.sub,treat=1:n.treat)
>
> treat<-d$treat
> d$treat<-factor(d$treat)
> d$treat<-relevel(d$treat,1)
> d$sub<-factor(d$sub)
> d$sub<-relevel(d$sub,1)
> d$rep<-factor(d$rep)
> d$rep<-relevel(d$rep,1)
>
>
> e.sub<-rnorm(n.treat*n.sub,sd=3)[d$treat:d$sub]
> e.rep<-rnorm(n.treat*n.sub*n.rep)
>
>
> d$y<- treat + e.sub + e.rep
>
> #parameterization with intecept
> g.withIntercept<-  lmer(y~treat+(1|treat:sub),data=d,method='REML')
> g2.withIntercept<-lmer2(y~treat+(1|treat:sub),data=d,method='REML')
>
> #parameterization without intecept
> g.withOutIntercept<-  lmer(y~0+treat+(1|treat:sub),data=d,method='REML')
> g2.withOutIntercept<-lmer2(y~0+treat+(1|treat:sub),data=d,method='REML')
>
>
> print('with intecept -all OK')
> print(g.withIntercept)
> print(g2.withIntercept)
>
>
>
> print('without intecept- in the lmer2 fit wrong estiimate for treat1
> (and STDERR)')
> print(g.withOutIntercept)
> print(g2.withOutIntercept)
>
>
> Ulrich
>
> platform       i386-pc-mingw32
> arch           i386
> os             mingw32
> system         i386, mingw32
> status
> major          2
> minor          4.1
> year           2006
> month          12
> day            18
> svn rev        40228
> language       R
> version.string R version 2.4.1 (2006-12-18)
>
>
> lme4: [Package lme4 version 0.9975-11
>
>
>
>
> Ulrich Halekoh, Ph.D
> UNIVERSITY OF AARHUS
> Faculty of Agricultural Sciences
> Dept. of Genetics and Biotechnology
> E-mail: Ulrich.Halekoh at agrsci.dk <mailto:Ulrich.Halekoh at agrsci.dk>
>
>
>         [[alternative HTML version deleted]]
>
> _______________________________________________
> R-sig-mixed-models at r-project.org mailing list
> https://stat.ethz.ch/mailman/listinfo/r-sig-mixed-models
>
-------------- next part --------------
> as(g2.withOutIntercept at L, "sparseMatrix")
9 x 9 sparse Matrix of class "dtCMatrix"
                                                                          
 [1,]  6.9183964 .         .        .         .         .         .       
 [2,]  .         6.9183964 .        .         .         .         .       
 [3,]  .         .         6.918396 .         .         .         .       
 [4,]  1.9789973 1.9789973 1.978997 0.5007087 .         .         .       
 [5,]  .         .         .        .         6.918396  .         .       
 [6,]  .         .         .        .         .         6.918396  .       
 [7,]  .         .         .        .         .         .         6.918396
 [8,]  .         .         .        .         1.978997  1.978997  1.978997
 [9,] -0.9793224 0.3366513 2.006967 0.1150608 6.897771 -5.165302 -9.247812
                         
 [1,]  .         .       
 [2,]  .         .       
 [3,]  .         .       
 [4,]  .         .       
 [5,]  .         .       
 [6,]  .         .       
 [7,]  .         .       
 [8,]  0.5007087 .       
 [9,] -0.6338222 4.131376
> str(g2.withOutIntercept at L)
Formal class 'dCHMsimpl' [package "Matrix"] with 10 slots
  ..@ x       : num [1:23] 47.864  0.286 -0.142 47.864  0.286 ...
  ..@ p       : int [1:10] 0 3 6 9 11 14 17 20 22 23
  ..@ i       : int [1:23] 0 3 8 1 3 8 2 3 8 3 ...
  ..@ nz      : int [1:9] 3 3 3 2 3 3 3 2 1
  ..@ nxt     : int [1:11] 1 2 3 4 5 6 7 8 9 -1 ...
  ..@ prv     : int [1:11] 10 0 1 2 3 4 5 6 7 8 ...
  ..@ colcount: int [1:9] 3 3 3 2 3 3 3 2 1
  ..@ perm    : int [1:9] 0 1 2 6 3 4 5 7 8
  ..@ type    : int [1:4] 1 0 0 1
  ..@ Dim     : int [1:2] 9 9