Slides & my Presentation style

this adds ?showNotes=true to the file

click this link for pdf printable + notes

HTML is best viewed in chrome

Overview

• Repeated measures
• An example to guide you through the steps
• The mixed model
• Inference
• Assumption checking
• Shrinkage
• Terminology
• Convergence Problems
• Multiple random variables

Repeated measures It looks like there is no difference between groups?

But what if we have data from within the same subjects?

Repeated measures Individual differences Between-Subject variability can be removed by repeated measures.

Never analyze repeated measures data without taking repeated measures into account (Disaggregation, Pseudoreplication)

Accouting for repeated measures? Black = mean of the subject.

$$n_S=5$$ subjects

Sampling distributions Multiple trials per subject with $$n_T =15$$ trials

Sampling distributions This sampling distribution reflects these exact 5 subjects. But does not generalize!

Sampling distributions (pink for aggregate, purple for disaggregate). There is a strong allure to choose the purple one, easier to get significant results!

$n_{subject}(aggregate) = df+1 = 5$ $n_{"subject"}(pseudoreplication) = df+1= 5*5 = 25$

DONT DO THIS, YOU KNOW (now) BETTER

Always necessary to aggregate? $\sigma_{within} >> \sigma_{between}$

Subjects are very similar

Not so necessary to aggregate first

A single parameter summarizes subject-means well

Should we trust all subjects equally? Question: Estimate the mean of the population

Should we trust all subjects equally? Question: Estimate the mean of the population

Should we trust all subjects equally? We should take into account sample-size and variance into our estimate of the population mean

We should trust those subjects more, where we most sure about their central tendency

Summary

• Why repeated measures
• More power to find effects that exist!
• Don’t forget about repeated measures!
• Else your type I error will be very high
• Should we trust all subjects equally?
• No, we should punish for small within-subject sample-size and punish high variance

Overview

• Repeated measures
• An example to guide you through the steps
• The mixed model
• Inference
• Assumption checking
• Shrinkage
• Terminology
• Convergence Problems
• Multiple random variables

The example EEG-Experiment: Contrast vs. P100 size

=> Very noisy subjects (babys?), between 4 and 29 trials

The data one model for all subjects $\hat{y} = \beta_0 + \beta_1 * contrast$

one model for each subject $\hat{y_j} = \beta_{0,j} + \beta_{1,j} * contrast$

Excourse: two-stage statistics One solution Linear Mixed Model aka Hierarchical Model aka Multilevel Model

We estimate the population distribution directly. Subject-means with high uncertainty are less influential for the model fit

A compromise Small second level variance (limit $$0$$)

Parameters/DF assumed to be n-trials Overview

• Repeated measures
• An example to guide you through the steps
• The mixed model
• Inference
• Assumption checking
• Shrinkage
• Terminology
• Convergence Problems
• Multiple random variables

Formal definition of the linear mixed model A word of warning

Encouraging psycholinguists to use linear mixed-effects models ‘was like giving shotguns to toddlers’ (Barr)

Because mixed models are more complex and more flexible than the general linear model, the potential for confusion and errors is higher (Hamer)

Many ongoing riddles:

Book:Richly Parameterized Linear Models

• Chp. 9: Puzzles from Analyzing Real Data sets
• Chp.10: A Random Effect Competing with a Fixed Effect
• Chp.12: Competition between Random Effects
• Chp.14: Mysterious, Inconvenient or Wrong Results from Real Datasets

The mixed model in R

$$\hat{y} = \beta_0 + \beta_1 * contrast + u_{0,j} + u_{1,j} * contrast$$

$$\hat{y}_j = \beta X_j + b_jZ$$

Fitting this in R using the package lme4

response ~ beta + (beta|group)

In our example:

response ~ 1 + contrast + (1|subject) + (0+contrast|subject)

Mixed models in R

summary(lmer(p100~1 + contrast + (1 + contrast||subject),data=d))
## Linear mixed model fit by REML ['lmerMod']
## Formula: p100 ~ 1 + contrast + ((1 | subject) + (0 + contrast | subject))
##    Data: d
##
## REML criterion at convergence: 2271.4
##
## Scaled residuals:
##      Min       1Q   Median       3Q      Max
## -2.71978 -0.62658 -0.01064  0.65739  2.91110
##
## Random effects:
##  Groups    Name        Variance Std.Dev.
##  subject   (Intercept) 10.16741 3.1886
##  subject.1 contrast     0.01353 0.1163
##  Residual              54.45510 7.3794
## Number of obs: 323, groups:  subject, 20
##
## Fixed effects:
##             Estimate Std. Error t value
## (Intercept)  4.03422    1.10536   3.650
## contrast     0.09536    0.03045   3.131
##
## Correlation of Fixed Effects:
##          (Intr)
## contrast -0.324 gray: one model, color: individual model, black: mixed model

Overview

• Repeated measures
• An example to guide you through the steps
• The mixed model
• Inference
• REML vs ML
• Assumption checking
• Shrinkage
• Terminology
• Convergence Problems
• Multiple random variables

Model-Comparisons aka p-values

Unresolved problem!

Read (three links):how to get p values, problems with stepwise regression, different ways to get p-values

One way (good and common, but not the best*): Model comparison for fixed effects:

Be sure to keep random structure equal

anova(mres0,mres1)

“refitting model(s) with ML (instead of REML)”

## Data: d
## Models:
## mres0: p100 ~ 1 + (1 + contrast | subject)
## mres1: p100 ~ 1 + contrast + (1 + contrast | subject)
##       Df    AIC    BIC  logLik deviance  Chisq Chi Df Pr(>Chisq)
## mres0  5 2286.5 2305.3 -1138.2   2276.5
## mres1  6 2280.1 2302.7 -1134.0   2268.1 8.4095      1   0.003733 **
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

* currently bootstrapping or mcmc recommended

REML VS ML

Intuition:

When estimating the variance parameter using maximum likelihood (ML) $$\sigma$$ is biased because $$\mu$$ is not known at the same time

Solution:

Fit all fixed effects $$\mu$$, then refit the variance parameters on the residing estimates.

This is called restricted maximum likelihood (REML). If $$n$$ is large, no difference between ML and REML exist. If $$n$$ is small REML gives better predictive results.

If you do model comparison and fixed effects (should) differ (usual case) you have to use ML

Bootstrapping

• resample from the data with replacement (or parametrically from $$H_0$$-Model) 1000 times
• calculate statistic of interest
b1 <- bootMer(mres1, FUN = function(x) as.numeric(logLik(x)), nsim = nboot,seed = 1)
b2 <- bootMer(mres0, FUN = function(x) as.numeric(logLik(x)), nsim = nboot,seed = 1)
bootLL =  2 * b1$t - 2 * b2$t # grab the values & plot
quantile(bootLL, probs = c(0.025, 0.975)) pbkrtest

Better & easier: Likelihood-Ratio using R-Package pbkrtest

summary(pbkrtest::PBmodcomp(mres1,mres0,nsim=nboot))
## Parametric bootstrap test; time: 30.27 sec; samples: 200 extremes: 3;
## large : p100 ~ 1 + contrast + (1 + contrast | subject)
## small : p100 ~ 1 + (1 + contrast | subject)
##            stat     df    ddf  p.value
## PBtest   8.3917               0.019900 *
## Gamma    8.3917               0.011404 *
## Bartlett 7.0365 1.0000        0.007986 **
## F        8.3917 1.0000 12.385 0.013029 *
## LRT      8.3917 1.0000        0.003769 **
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

The main con against bootstrapping is the high computational time (use at least 1000 not 200 as in this example)

Overview

• Repeated measures
• An example to guide you through the steps
• The mixed model
• Inference
• Assumption checking
• Level 1
• Level 2
• Shrinkage
• Terminology
• Convergence Problems
• Multiple random variables

Checking Assumptions Level I No strong non-linearity visible in residuals

Checking assumptions Level II

In addition: test whether normality for random-variable is true:  QQ-Plots indicate normality of both random effects.

Overview

• Repeated measures
• The mixed model
• An example to guide you through the steps
• Inference
• Assumption checking
• Shrinkage
• Terminology
• Convergence Problems
• Multiple random variables

Shrinkage The individual estimates with high variance / far off the population are shrinked inwards

Shrinkage-plot  gray: one model, color: individual model, black: mixed model (with shrinkage)

Shrinkage over multiple parameters Simulated data: Subjects with faster reaction time (intercept) show larger effects for difficulty (slope)

In reverse: If I know somebody is generally slow, he should also have a smaller slope!

The multivariate normal II $$\rho = 0$$, knowing one dimension does not tell you about the other

$\sum = \left[ {\begin{array}{*{20}{c}} 1&0 \\ 0 &1 \end{array}} \right]$

Back to mixed models $\sum = \left[ {\begin{array}{*{20}{c}} 1&\rho \\ \rho &1 \end{array}} \right]$

Overview

• Repeated measures
• The mixed model
• An example to guide you through the steps
• Inference
• Assumption checking
• Shrinkage
• Terminology
• Convergence Problems
• Multiple random variables

Naming is a bit messy

Random/Fixed Coefficients

• Random intercept
• (1|subject): Different offset for each subject
• Random slope
• (predictor|subject): Different predictor-slope for each subject
• Random correlation
• (1+predictor|subject) or (X+Y|subject): Knowing about X tells you something about Y
• Fixed coefficient
• exactly the same for all subjects

Overview

• Repeated measures
• The mixed model
• An example to guide you through the steps
• Inference
• Assumption checking
• Shrinkage
• Terminology
• Convergence Problems
• Multiple random variables

A more complex model

d$contrastSquare = d$contrast^2
lmer(p100~1+contrast+contrastSquare+(1+contrast+contrastSquare|subject),data = d)


Warning messages:
1: Some predictor variables are on very different scales: consider rescaling
2: In checkConv(attr(opt, "derivs"), opt$par, ctrl = control$checkConv,  :
3: In checkConv(attr(opt, "derivs"), opt$par, ctrl = control$checkConv,  :
Model failed to converge: degenerate  Hessian with 2 negative eigenvalues


Convergence problems

Indicator

• Convergence warning (always a problem)
• random effect correlation is either 0 or 1 (very seldom correct)
• model takes very long time to fit (could be correct though)

Reasons

• misspecified or overspecified models (specified fixed variable as random variable?)
• too many parameters and not enough data,
• data have big outliers/are noisy

• clean data
• center predictors (remove mean) and scale (divide by SD of independent variable)
• use Bayesian fit and put a soft prior on the covariance structure (MCMCglmm,rstanarm,brms)
• simplify the model (parsimoneous models; Bates 2015)

The ideal case

$\sum = \left[ {\begin{array}{*{20}{c}} 1&{{\rho_{ab}}}&{{\rho_{ac}}}&{{\rho_{ad}}}\\ {{\rho_{ab}}}&1&{{\rho_{bc}}}&{{\rho_{bd}}}\\ {{\rho_{ac}}}&{{\rho_{bc}}}&1&{{\rho_{cd}}}\\ {{\rho_{ad}}}&{{\rho_{bd}}}&{{\rho_{cd}}}&1 \end{array}} \right]$

y ~ 1 + A + B + C + (1 + A + B + C | group)

If you have enough data, try to estimate all correlations.

The simplest structure

$\Sigma = \left[ {\begin{array}{*{20}{c}} 1&0&0&0\\ 0&1&0&0\\ 0&0&1&0\\ 0&0&0&1 \end{array}} \right]$

y ~ 1 + X1 + X2 + X3 + (1|group) + (0+X1|group) + (0+X2|group) + (0+X3|group)

shorthand:

y ~ 1 + X1 + X2 + X3 + (1+X1+X2+X3||group)

note the || instead if |

Which one to choose?

• 1: You can estimate the full (enough data or put a lkj-prior if you are Bayes)
• 2: Prior assumption on shape and you don’t have to change it (non-contingent on data)
• 3: Simplify in a data-driven way

Barr promotes full models. Bates promotes a step wise reduction if effects do not significantly explain variance

Overview

• Repeated measures
• The mixed model
• An example to guide you through the steps
• Inference
• Assumption checking
• Shrinkage
• Terminology
• Convergence Problems
• Multiple random variables

Multiple random variables What a nice effect! Write paper and end of story?

But very high stimulus variance! But very high stimulus variance You should include it in your model Often it is necessary to add multiple random effects! Examples are: subjects, stimuli (images, words, etc)

y ~ 1 + condA (1+condA|subject) + (1+condA|stimulus)

When to include what

Some guidance:

• Predictor is a subset of a population?
• => random/grouping variable e.g. subjects, images, words
• Predictor is within subject/grouping?
• => random coefficients
• Predictor is between subject/grouping or aggregated?
• => only fixed coefficient

Summary

• We learned how to analyze repeated measures data
• We learned how to structure random effects
• We learned how to check assumptions of Mixed Models