[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 20:01:15 CEST 2014
[cc'ing back to r-sig-mixed-models]
My point (which I probably should have said more directly) is that
you've constructed an unidentifiable model. It's not really clear what
you're trying to do, but as currently formulated it's just not possible
to tell the difference between "among-subject variation" and
"among-subject variation in response to treatment". I think if you want
to do this, you might have to work on constructing dummy variables
yourself: I think it might work better to use sum-to-zero contrasts ...
On 14-06-14 01:16 PM, Ulf Mertens wrote:
> Thanks for your detailed answer. Actually, I am simulating my data, too,
> so it shouldn't be very difficult to figure out what I am doing wrong.
> First, here's the output from eigen(m2 at optinfo$derivs$Hessian)$values:
>
> eigen(m2 at optinfo$derivs$Hessian)$values
> [1] 482.082704275 441.411371639 219.759912576 193.156411039
> 155.499675763 69.193189621 -0.000237569
>
> See my script to simulate the data attached. I am sampling the random
> intercept and the three random slopes from a multivariate normal
> distribution. The data, in the end, looks like this (here with 6 trials
> per subject):
>
> subject trial beta0 beta1 beta2 beta3 treatEff1 treatEff2 treatEff3
> subjectEffect error treatment rt
> 1 1 1 600 0 0 64.55883 40.01480 -14.71201 -27.38827
> 51.57127 -103.181056 1 588.4050
> 2 1 2 600 0 0 64.55883 40.01480 -14.71201 -27.38827
> 51.57127 34.425169 1 726.0112
> 3 1 3 600 0 0 64.55883 40.01480 -14.71201 -27.38827
> 51.57127 -24.821108 2 612.0382
> 4 1 4 600 0 0 64.55883 40.01480 -14.71201 -27.38827
> 51.57127 -63.495611 2 573.3637
> 5 1 5 600 0 0 64.55883 40.01480 -14.71201 -27.38827
> 51.57127 -67.078324 3 621.6635
> 6 1 6 600 0 0 64.55883 40.01480 -14.71201 -27.38827
> 51.57127 -5.059233 3 683.6826
> 7 2 1 600 0 0 64.55883 15.39385 23.04594 -78.88895
> -62.95682 -5.987356 1 546.4497
>
> and the model which created the response times (rt) is this:
>
> df <- within(df, rt <- beta0 +
> subjectEffect +
> (beta1+treatmentEffect1)*x1 +
> (beta2 + treatmentEffect2)*x2 +
> (beta3 + treatmentEffect3)*x3 +
> error)
>
> Thanks in advance
>
> Regards,
> Ulf
>
>
>
>
>
> On Sat, Jun 14, 2014 at 6:17 PM, Ben Bolker <bbolker at gmail.com
> <mailto:bbolker at gmail.com>> wrote:
>
>
> 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 <mailto: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 <tel:%2B%2032%202%20525%2002%2051>
> >> + 32 54 43 61 85 <tel:%2B%2032%2054%2043%2061%2085>
> >> Thierry.Onkelinx at inbo.be <mailto:Thierry.Onkelinx at inbo.be>
> >> www.inbo.be <http://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
> <mailto:r-sig-mixed-models-bounces at r-project.org> [
> >> r-sig-mixed-models-bounces at r-project.org
> <mailto:r-sig-mixed-models-bounces at r-project.org>] namens Ulf Mertens [
> >> mertens.ulf at gmail.com <mailto:mertens.ulf at gmail.com>]
> >> Verzonden: zaterdag 14 juni 2014 10:56
> >> Aan: r-sig-mixed-models at r-project.org
> <mailto: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
> <mailto: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
> <mailto:R-sig-mixed-models at r-project.org> mailing list
> > https://stat.ethz.ch/mailman/listinfo/r-sig-mixed-models
> >
>
> _______________________________________________
> R-sig-mixed-models at r-project.org
> <mailto: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