[R-sig-ME] Altering the intercept-vector in lme4?

[Douglas Bates]

On Fri, Sep 19, 2008 at 10:25 AM, Rense Nieuwenhuis
<contact at rensenieuwenhuis.nl> wrote:
> Dear all,
>
> while working on a project regarding a new approach at evaluating
> measures on influential 'cases' in mixed effects models, we would
> like to alter the intercept in such a way that observations nested
> within a specified grouping factor have a value 0 on the intercept
> (and all the others the normal 1).
>
> Whereas this is relatively simple in some statistical packages that
> require the manual definition of the intercept, I have difficulties
> to understand how to do this in the lme4 package. This is what I tried:
>
> require(lme4)
>
> data(bdf, package='mlmRev')
> model.1 <- lmer(IQ.perf ~ sex + schoolSES + (1 | schoolNR), data=bdf)
> data.adapted <- bdf
> data.adapted$intercept.alt <- 1
> data.adapted$intercept.alt[data.adapted$schoolNR==1] <- 0
> data.adapted$dummy <- 0
> data.adapted$dummy[data.adapted$schoolNR==1] <- 1

> model.2 <- lmer(IQ.perf ~ -1 + intercept.alt + dummy + sex +
> schoolSES + (-1 + intercept.alt | schoolNR), data=data.adapted)

> (Please note that we also add a dummy with value 1 for the specified
> grouping factor, and a 0 for all other groups.)

I think that is the cause of the problem when you combine those
indicators with the indicators for sex.  The way that a model matrix
is constructed in R, if the intercept is removed then both indicators
are included for the first factor (sex in this case).  You can
circumvent this problem by creating a factor for schoolNR == 1 as
shown below.  In general you should avoid trying to create indicator
variables explicitly and use the model formula to create them

> bdf <- within(bdf, sch1 <- factor(ifelse(schoolNR == 1, "Y", "N")))
> (model.1 <- lmer(IQ.perf ~ sex + schoolSES + (1 | schoolNR), data=bdf))
Linear mixed model fit by REML
Formula: IQ.perf ~ sex + schoolSES + (1 | schoolNR)
   Data: bdf
   AIC   BIC logLik deviance REMLdev
 10058 10087  -5024    10034   10048
Random effects:
 Groups   Name        Variance Std.Dev.
 schoolNR (Intercept) 0.17049  0.41291
 Residual             4.58623  2.14155
Number of obs: 2287, groups: schoolNR, 131

Fixed effects:
            Estimate Std. Error t value
(Intercept) 10.16620    0.26013   39.08
sex1        -0.46084    0.09063   -5.08
schoolSES    0.05676    0.01333    4.26

Correlation of Fixed Effects:
          (Intr) sex1
sex1      -0.146
schoolSES -0.960 -0.020
> (model.2 <- lmer(IQ.perf ~ 0 + sch1 + sex + schoolSES + (1 | schoolNR), bdf))
Linear mixed model fit by REML
Formula: IQ.perf ~ 0 + sch1 + sex + schoolSES + (1 | schoolNR)
   Data: bdf
   AIC   BIC logLik deviance REMLdev
 10058 10092  -5023    10033   10046
Random effects:
 Groups   Name        Variance Std.Dev.
 schoolNR (Intercept) 0.16993  0.41222
 Residual             4.58667  2.14165
Number of obs: 2287, groups: schoolNR, 131

Fixed effects:
          Estimate Std. Error t value
sch1N     10.21585    0.26491   38.56
sch1Y      9.62569    0.61365   15.69
sex1      -0.46154    0.09064   -5.09
schoolSES  0.05445    0.01353    4.03

Correlation of Fixed Effects:
          sch1N  sch1Y  sex1
sch1Y      0.242
sex1      -0.145 -0.055
schoolSES -0.961 -0.241 -0.019

> Unfortunately, model.2 does not converge**. This is the error message
> that is given:
>
> Error in mer_finalize(ans, verbose) :
>   Downdated X'X is not positive definite, 4.
>
> Does anyone know if it is possible to manually change the intercept-
> vector, giving it the value 0 for a select number of cases? If it is
> possible, how is it done? This would help us greatly in translating a
> procedure to R-Project, making it available to a larger audience.
>
> Many thanks in advance,
>
> with kind regards,
>
> Rense Nieuwenhuis
>
>
>
> ** It does when I specify the sex-variable as.numeric(). This however
> results in really messy outcomes.
>
>
> - - -- --- ----- --------
> Rense Nieuwenhuis
> +31 6 481 05 683
>
> www.rensenieuwenhuis.nl
>
>
>        [[alternative HTML version deleted]]
>
> _______________________________________________
> R-sig-mixed-models at r-project.org mailing list
> https://stat.ethz.ch/mailman/listinfo/r-sig-mixed-models
>