Unfortunately, I don't know how to do this. Could you please give me a hint
on how to generate simple data with sum-to-zero contrasts so that I can
take this as a starting point?

My goal is to be able to generate data so that the assumption of compound
symmetry/sphericity is either met or violated. So far,  I sampled my random
effects from a multivariate normal distribution with either equal
correlations (met) or different correlations (violated). But I do not know
how to adjust this to sum-to-zero and I couldn't find the reason why random
intercept and random slopes are confounded (why my constructed model is
unidentifiable).


Thank you in advance


On Sat, Jun 14, 2014 at 8:01 PM, Ben Bolker <bbolker@gmail.com> wrote:

>   [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@optinfo$derivs$Hessian)$values:
> >
> > eigen(m2@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@gmail.com
> > <mailto:bbolker@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@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@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@inbo.be <mailto:Thierry.ONKELINX@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@inbo.be <mailto:Thierry.Onkelinx@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@r-project.org
> >     <mailto:r-sig-mixed-models-bounces@r-project.org> [
> >     >> r-sig-mixed-models-bounces@r-project.org
> >     <mailto:r-sig-mixed-models-bounces@r-project.org>] namens Ulf
> Mertens [
> >     >> mertens.ulf@gmail.com <mailto:mertens.ulf@gmail.com>]
> >     >> Verzonden: zaterdag 14 juni 2014 10:56
> >     >> Aan: r-sig-mixed-models@r-project.org
> >     <mailto:r-sig-mixed-models@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@r-project.org
> >     <mailto:R-sig-mixed-models@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@r-project.org
> >     <mailto:R-sig-mixed-models@r-project.org> mailing list
> >     > https://stat.ethz.ch/mailman/listinfo/r-sig-mixed-models
> >     >
> >
> >     _______________________________________________
> >     R-sig-mixed-models@r-project.org
> >     <mailto:R-sig-mixed-models@r-project.org> mailing list
> >     https://stat.ethz.ch/mailman/listinfo/r-sig-mixed-models
> >
> >
>
>

	[[alternative HTML version deleted]]

