[R-sig-ME] Altering the intercept-vector in lme4?
Douglas Bates
bates at stat.wisc.edu
Fri Sep 19 20:08:07 CEST 2008
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
>
More information about the R-sig-mixed-models
mailing list