Preliminaries
library(tidyverse)
library(lme4)
library(lmerTest)
library(gazer) # helper functions and example data sets
For this example, we will consider the simple target fixation data shown on the right from a VWP experiment with a word frequency manipulation (High frequency words recognized faster than Low frequency words). Note that since Condition (High vs. Low frequency) was manipulated within-participants, we will want to include nested random effects of Condition with Subject.
The starting data frame is shown below. Subject is the the subject ID, Time is time (in ms) from word onset, timeBin identifies the 100ms time bin corresponding to this time point (convenient for aligning orthogonal time), Condition is the word frequency condition, and meanFix is the mean proportion of fixations to the target object.
summary(TargetFix)## Subject Time timeBin Condition meanFix
## 708 : 30 Min. : 300 Min. : 1 High:150 Min. :0.02857
## 712 : 30 1st Qu.: 450 1st Qu.: 4 Low :150 1st Qu.:0.27778
## 715 : 30 Median : 650 Median : 8 Median :0.45584
## 720 : 30 Mean : 650 Mean : 8 Mean :0.44827
## 722 : 30 3rd Qu.: 850 3rd Qu.:12 3rd Qu.:0.61111
## 725 : 30 Max. :1000 Max. :15 Max. :0.82857
## (Other):120
## sumFix N
## Min. : 1.00 Min. :33.00
## 1st Qu.:10.00 1st Qu.:35.75
## Median :16.00 Median :36.00
## Mean :15.94 Mean :35.55
## 3rd Qu.:21.25 3rd Qu.:36.00
## Max. :29.00 Max. :36.00
## The first step is to create a third-order polynomial based on the timeBin values.
TargetFix.gca <- code_poly(TargetFix, predictor="timeBin", poly.order = 3, draw.poly = FALSE)
summary(TargetFix.gca)## Subject Time timeBin Condition meanFix
## 708 : 30 Min. : 300 Min. : 1 High:150 Min. :0.02857
## 712 : 30 1st Qu.: 450 1st Qu.: 4 Low :150 1st Qu.:0.27778
## 715 : 30 Median : 650 Median : 8 Median :0.45584
## 720 : 30 Mean : 650 Mean : 8 Mean :0.44827
## 722 : 30 3rd Qu.: 850 3rd Qu.:12 3rd Qu.:0.61111
## 725 : 30 Max. :1000 Max. :15 Max. :0.82857
## (Other):120
## sumFix N timeBin.Index poly1
## Min. : 1.00 Min. :33.00 Min. : 1 Min. :-0.4183
## 1st Qu.:10.00 1st Qu.:35.75 1st Qu.: 4 1st Qu.:-0.2390
## Median :16.00 Median :36.00 Median : 8 Median : 0.0000
## Mean :15.94 Mean :35.55 Mean : 8 Mean : 0.0000
## 3rd Qu.:21.25 3rd Qu.:36.00 3rd Qu.:12 3rd Qu.: 0.2390
## Max. :29.00 Max. :36.00 Max. :15 Max. : 0.4183
##
## poly2 poly3
## Min. :-0.29063 Min. :-0.4563
## 1st Qu.:-0.22835 1st Qu.:-0.2457
## Median :-0.04152 Median : 0.0000
## Mean : 0.00000 Mean : 0.0000
## 3rd Qu.: 0.26987 3rd Qu.: 0.2457
## Max. : 0.47227 Max. : 0.4563
## Since this is a simple case with just one within-subjects fixed effect that has only two levels, we can skip to the full model and examine its parameter estimates:
m.full <- lmer(meanFix ~ (poly1+poly2+poly3)*Condition +
(poly1+poly2+poly3 | Subject) +
(poly1+poly2 | Subject:Condition),
data=TargetFix.gca, REML=FALSE)## boundary (singular) fit: see ?isSingularsummary(m.full)## Linear mixed model fit by maximum likelihood . t-tests use Satterthwaite's
## method [lmerModLmerTest]
## Formula: meanFix ~ (poly1 + poly2 + poly3) * Condition + (poly1 + poly2 +
## poly3 | Subject) + (poly1 + poly2 | Subject:Condition)
## Data: TargetFix.gca
##
## AIC BIC logLik deviance df.resid
## -815.6 -723.0 432.8 -865.6 275
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -3.0347 -0.5745 0.0072 0.5989 3.3854
##
## Random effects:
## Groups Name Variance Std.Dev. Corr
## Subject:Condition (Intercept) 1.672e-03 0.040885
## poly1 1.801e-02 0.134200 -0.43
## poly2 1.269e-02 0.112668 -0.33 0.72
## Subject (Intercept) 1.120e-04 0.010581
## poly1 1.588e-02 0.126030 0.92
## poly2 9.002e-05 0.009488 -0.69 -0.92
## poly3 2.139e-03 0.046247 -0.89 -1.00 0.94
## Residual 2.030e-03 0.045050
## Number of obs: 300, groups: Subject:Condition, 20; Subject, 10
##
## Fixed effects:
## Estimate Std. Error df t value Pr(>|t|)
## (Intercept) 4.773e-01 1.385e-02 1.993e+01 34.458 < 2e-16 ***
## poly1 6.386e-01 5.994e-02 1.673e+01 10.654 7.16e-09 ***
## poly2 -1.096e-01 3.849e-02 2.000e+01 -2.848 0.00995 **
## poly3 -9.326e-02 2.042e-02 1.827e+01 -4.568 0.00023 ***
## ConditionLow -5.811e-02 1.901e-02 1.511e+01 -3.057 0.00793 **
## poly1:ConditionLow 3.188e-04 6.331e-02 1.726e+01 0.005 0.99604
## poly2:ConditionLow 1.635e-01 5.427e-02 1.987e+01 3.014 0.00689 **
## poly3:ConditionLow -2.087e-03 2.015e-02 2.352e+02 -0.104 0.91759
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Correlation of Fixed Effects:
## (Intr) poly1 poly2 poly3 CndtnL pl1:CL pl2:CL
## poly1 -0.135
## poly2 -0.297 0.424
## poly3 -0.154 -0.475 0.053
## ConditionLw -0.686 0.206 0.207 0.000
## ply1:CndtnL 0.268 -0.528 -0.447 0.000 -0.390
## ply2:CndtnL 0.201 -0.335 -0.705 0.000 -0.294 0.634
## ply3:CndtnL 0.000 0.000 0.000 -0.493 0.000 0.000 0.000
## convergence code: 0
## boundary (singular) fit: see ?isSingularAnother approach is to selectively remove individual effects of Condition and use model comparisons to evaluate each effect. This method would have produced very similar results:
## Data: TargetFix.gca
## Models:
## m.base: meanFix ~ (poly1 + poly2 + poly3) + (poly1 + poly2 + poly3 |
## m.base: Subject) + (poly1 + poly2 | Subject:Condition)
## m.0: meanFix ~ (poly1 + poly2 + poly3) + Condition + (poly1 + poly2 +
## m.0: poly3 | Subject) + (poly1 + poly2 | Subject:Condition)
## m.1: meanFix ~ (poly1 * Condition + poly2 + poly3) + (poly1 + poly2 +
## m.1: poly3 | Subject) + (poly1 + poly2 | Subject:Condition)
## m.2: meanFix ~ (poly1 + poly2) * Condition + poly3 + (poly1 + poly2 +
## m.2: poly3 | Subject) + (poly1 + poly2 | Subject:Condition)
## m.full: meanFix ~ (poly1 + poly2 + poly3) * Condition + (poly1 + poly2 +
## m.full: poly3 | Subject) + (poly1 + poly2 | Subject:Condition)
## npar AIC BIC logLik deviance Chisq Df Pr(>Chisq)
## m.base 21 -808.57 -730.79 425.28 -850.57
## m.0 22 -810.72 -729.23 427.36 -854.72 4.1481 1 0.041681 *
## m.1 23 -812.14 -726.96 429.07 -858.14 3.4293 1 0.064050 .
## m.2 24 -817.63 -728.74 432.81 -865.63 7.4845 1 0.006223 **
## m.full 25 -815.64 -723.05 432.82 -865.64 0.0107 1 0.917471
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1Here is how we might report these results:
Mirman, D. (2014). Growth Curve Analysis and Visualization Using R. Boca Raton, FL: Chapman and Hall / CRC Press.
Growth curve analysis (Mirman, 2014) was used to analyze the target gaze data from 300ms to 1000ms after word onset. The overall time course of target fixations was modeled with a third-order (cubic) orthogonal polynomial and fixed effects of Condition (Low vs. High frequency; within-participants) on all time terms. The model also included participant random effects on all time terms and participant-by-condition random effects on all time terms except the cubic (estimating random effects is “expensive” in terms of the number of observation required, so this cubic term was excluded because it tends to capture less-relevant effects in the tails). Parameter estimate degrees of freedom and corresponding p-values were estimated using Satterthwaite’s method. There was a significant effect of Condition on the intercept term, indicating lower overall target fixation proportions for the Low condition relative to the High condition (Estimate = -0.0581, SE = 0.0196, p < 0.01). There was also a significant effect on the quadratic term, indicating shallower curvature - slower word recognition - in the Low condition relative to the High condition (Estimate = 0.164, SE = 0.0544, p < 0.01). All other effects of Condition were not significant (see Table 1 for full results).