Altering the intercept-vector in lme4?
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