[R-sig-ME] Random Intercept + random slope model yields exactly the same results as random slope
Ben Bolker
bbolker at gmail.com
Sat Jun 14 18:17:19 CEST 2014
On 14-06-14 09:51 AM, Ulf Mertens wrote:
> Dear Thierry,
>
> thanks for your answer. The same also applies for those two models:
>
> m1 <- lmer(rt ~ treatment + (1+treatment|subject),data=df,REML=F)
> m2 <- lmer(rt ~ treatment + (1|subject)+ (0+treatment|subject)
> ,data=df,REML=F)
There are two possibilities, the first seems to be the case:
(1) (1|subject) is confounded with (0+treatment|subject) -- that is,
the variation of treatments by subject includes the intercept variation
among subjects. If this is the case, you should be getting warnings
about eigenvalues, and eigen(m2 at optinfo$derivs$Hessian)$values should
include a zero. This seems to be the case (see example code below).
(2) the other would be that the model is well-behaved/well-defined, but
if the among-subject variance was zero, then you'd get the same fit to
both models.
===========
## set up data
dd <- expand.grid(treatment=LETTERS[1:4],subject=factor(1:30),rep=1:10)
## v-cov matrix
mm <- matrix(0.5,nrow=4,ncol=4)
diag(mm) <- 1:4
## transform to Cholesky-factor
theta <- t(chol(mm))[lower.tri(mm,diag=TRUE)]
## simulate response
library(lme4)
dd$s1 <- simulate(~treatment+(1+treatment|subject),
newdata=dd,
newparams=list(beta=rep(0,4),
theta=theta,
sigma=0.5),
family=gaussian, seed=101)[[1]]
## utility
library(nloptr)
nloptWrap <- function(fn, par, lower, upper, control=list(), ...) {
defaultControl <- list(xtol_rel = 1e-6, maxeval = 1e5)
for (n in names(defaultControl))
if (is.null(control[[n]])) control[[n]] <- defaultControl[[n]]
res <- nloptr(x0=par, eval_f=fn, lb=lower, ub=upper, opts=control, ...)
## ------
with(res,list(par=solution,
fval=objective,
feval=iterations,
conv=if (status>0) 0 else status,
message=message))
}
## basic fit
m1 <- lmer(s1~treatment+(1+treatment|subject),data=dd,REML=FALSE)
m1B <- lmer(s1~treatment+(0+treatment|subject),data=dd,REML=FALSE)
nlc <- lmerControl(optimizer=nloptWrap,
optCtrl=list(algorithm="NLOPT_LN_BOBYQA"))
## none of these models converge/fit well, even playing around
## with optimizers etc.
m2 <-
lmer(s1~treatment+(1|subject)+(1+treatment|subject),data=dd,REML=FALSE,
control=nlc)
m3 <-
lmer(s1~treatment+(1|subject)+(0+treatment|subject),data=dd,REML=FALSE)
library(Matrix)
tmpf <- function(x) {
h <- Matrix(x at optinfo$derivs$Hessian)
ee <- eigen(h)$values
## force pos def
if (min(ee)<0) h <- Matrix::nearPD(Matrix(h))$mat
cc <- cov2cor(solve(h))
cc[row(cc)==col(cc)] <- NA ## diag() <- NA doesn't work for Matrix?
print(image(cc))
ee
}
tmpf(m1)
tmpf(m1B) ## rearranged, but not fundamentally different
tmpf(m2)
tmpf(m3)
>
> deviance(m1)
> [1] 32769.35
>
> deviance(m2)
> [1] 32769.35
>
> anova(m1,m2)
> Data: df
> Models:
> ml3: rt ~ treatment + (1 + treatment | subject)
> ml2: rt ~ treatment + (1 | subject) + (0 + treatment | subject)
> Df AIC BIC logLik deviance Chisq Chi Df Pr(>Chisq)
> ml3 10 32789 32848 -16385 32769
> ml2 11 32791 32856 -16385 32769 0 1 1
>
> Can you tell me why this is?
>
> Thanks for your help!
>
>
> On Sat, Jun 14, 2014 at 12:59 PM, ONKELINX, Thierry <
> Thierry.ONKELINX at inbo.be> wrote:
>
>> Dear Ulf,
>>
>> Treatment is a factor. Hence the random 'slope' is a random intercept for
>> each level of treatment. The difference between 1 + treatment and 0 +
>> treatment (in case of a factor) is that 1 + treament uses the first level
>> as a reference (random 'intercept'). The other levels are coded as the
>> difference from the reference levels.
>>
>> 0 + treatment uses no reference level. All levels are coded as the direct
>> effect of the level. Hence you get the same fit but with different
>> parametrisation. Note that the variance-covariance matrix of the random
>> effects are different.
>>
>> Best regards,
>>
>> ir. Thierry Onkelinx
>> Instituut voor natuur- en bosonderzoek / Research Institute for Nature and
>> Forest
>> team Biometrie & Kwaliteitszorg / team Biometrics & Quality Assurance
>> Kliniekstraat 25
>> 1070 Anderlecht
>> Belgium
>> + 32 2 525 02 51
>> + 32 54 43 61 85
>> Thierry.Onkelinx at inbo.be
>> www.inbo.be
>> To call in the statistician after the experiment is done may be no more
>> than asking him to perform a post-mortem examination: he may be able to say
>> what the experiment died of. ~ Sir Ronald Aylmer Fisher
>> The plural of anecdote is not data. ~ Roger Brinner
>> The combination of some data and an aching desire for an answer does not
>> ensure that a reasonable answer can be extracted from a given body of data.
>> ~ John Tukey
>>
>> ________________________________________
>> Van: r-sig-mixed-models-bounces at r-project.org [
>> r-sig-mixed-models-bounces at r-project.org] namens Ulf Mertens [
>> mertens.ulf at gmail.com]
>> Verzonden: zaterdag 14 juni 2014 10:56
>> Aan: r-sig-mixed-models at r-project.org
>> Onderwerp: [R-sig-ME] Random Intercept + random slope model yields exactly
>> the same results as random slope
>>
>> Hi there,
>>
>> When I run a random intercept + random slope model, I get exactly the same
>> result as when running a model where there is no random intercept. I can't
>> figure out why this is.
>>
>> Model 1:
>>
>> *m1 <- lmer(rt ~ treatment + (1+treatment|subject),data=df)*
>> *m1*
>> *Linear mixed model fit by REML ['lmerMod']*
>> *Formula: rt ~ treatment + (1 + treatment | subject)*
>> * Data: df*
>> *REML criterion at convergence: 32558.52*
>> *Random effects:*
>> * Groups Name Std.Dev. Corr *
>> * subject (Intercept) 71.64 *
>> * treatment2 54.55 -0.31 *
>> * treatment3 54.79 -0.34 0.89*
>> * Residual 97.38 *
>> *Number of obs: 2700, groups: subject, 30*
>> *Fixed Effects:*
>> *(Intercept) treatment2 treatment3 *
>> * 598.6675 -0.9037 50.9770 *
>>
>> Model 2:
>>
>> *m2 <- lmer(rt ~ treatment + (0+treatment|subject) ,data=df)*
>> *m2*
>> *Linear mixed model fit by REML ['lmerMod']*
>> *Formula: rt ~ treatment + (0 + treatment | subject)*
>> * Data: df*
>> *REML criterion at convergence: 32558.52*
>> *Random effects:*
>> * Groups Name Std.Dev. Corr *
>> * subject treatment1 71.64 *
>> * treatment2 75.18 0.72 *
>> * treatment3 74.03 0.72 0.94*
>> * Residual 97.38 *
>> *Number of obs: 2700, groups: subject, 30*
>> *Fixed Effects:*
>> *(Intercept) treatment2 treatment3 *
>> * 598.6675 -0.9037 50.9770 *
>>
>> Compare both models:
>>
>> *anova(m1,m2)*
>> *refitting model(s) with ML (instead of REML)*
>> *Data: df*
>> *Models:*
>> *m1: rt ~ treatment + (1 + treatment | subject)*
>> *m2: rt ~ treatment + (0 + treatment | subject)*
>> * Df AIC BIC logLik deviance Chisq Chi Df Pr(>Chisq)*
>> *m1 10 32598 32657 -16289 32578 *
>> *m2 10 32598 32657 -16289 32578 0 0 1*
>>
>> Thanks in advance
>>
>> [[alternative HTML version deleted]]
>>
>> _______________________________________________
>> R-sig-mixed-models at r-project.org mailing list
>> https://stat.ethz.ch/mailman/listinfo/r-sig-mixed-models
>> * * * * * * * * * * * * * D I S C L A I M E R * * * * * * * * * * * * *
>> Dit bericht en eventuele bijlagen geven enkel de visie van de schrijver
>> weer en binden het INBO onder geen enkel beding, zolang dit bericht niet
>> bevestigd is door een geldig ondertekend document.
>> The views expressed in this message and any annex are purely those of the
>> writer and may not be regarded as stating an official position of INBO, as
>> long as the message is not confirmed by a duly signed document.
>>
>
> [[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