OverviewGrowth curve analysis (GCA) is a multilevel regression technique designed for analysis of time course or longitudinal data. A major advantage of this approach is that it can be used to simultaneously analyze both group-level effects (e.g., experimental manipulations) and individual-level effects (i.e., individual differences). We have been using this method for several years, particularly in the context of visual world paradigm (VWP) eye tracking data and learning curves, though it can be applied to any time course data.
The basics are summarized below. Detailed tutorial materials are available on the LCDL github page. All of the analyses are conducted using the R package
lme4 and the graphs generated using the ggplot2 package. For more information about R, check out our R Resources page.Key challenges in analyzing time series data
The big pictureGrowth curve analysis provides a way to address those challenges by explicitly modeling change over time and quantifying both group-level and individual-level differences. To specify a growth curve model, you need to decide on three key components:
A simple exampleTo see how these components come together, consider theChickWeight data set (part of the default R installation), which has data from an experiment on the effect of diet on early growth of chicks. We start out with just a "base" model of chick growth allowing for individual variability in weight (in technical terms, a random intercept for each Chick: (1 | Chick) ):> m.base <- lmer(weight ~ Time + (1 | Chick), data=ChickWeight, REML=F) To this, we can add fixed effects of diet on the intercept (i.e., a constant difference in weights among chicks randomly assigned to different diets):
> m.0 <- lmer(weight ~ Time + Diet + (1 | Chick), data=ChickWeight, REML=F) and on the slope (i.e., effects of diet on the rate of growth):
> m.1 <- lmer(weight ~ Time * Diet + (1 | Chick), data=ChickWeight, REML=F) We can visually examine the effect of adding these terms by plotting the model fits against the observed data. The observed data are on the left, the m.0 (intercept model) fit is in the middle, and the m.1 (linear model) fit is on the right:
To statistically evaluate the effects of adding these terms, we can examine the change in the goodness of fit (log likelihood) through model comparisons:
> anova(m.base, m.0, m.1)
Of course, in addition to different baseline weights, the chicks might have other (unmeasured) individual properties that affect their growth rate. To capture this, we can add a linear effect on time of Chick to the random effects, and we can use model comparisons to examine whether this effect improved model fit:
> m.t1 <- lmer(weight ~ Time * Diet + (1 + Time | Chick), data=ChickWeight, REML=F) > anova(m.1, m.t1)
> coef(summary(m.t1))
Adding the random slopes did improve model fit. As a general rule, the model should include all random effects that are licensed by the design, that is, all the ones that could vary across participants (Barr et al., 2013). Notice that adding the random slope effect caused slight changes in the fixed effect parameter estimates because some of the variance was now captured by that random effect. This is why Barr et al. recommend using the "maximal" random effect structure -- even when a random effect does not improve model fit, it can still affect the fixed effect estimates and excluding it can elevate the false positive rate.
Applying GCA to VWP dataOrthogonal polynomialsOften, our data are not perfectly straight, so we want to capture that curvature with higher-order polynomial terms (time squared, time cubed, etc.). Because our time variable is usually only positive, natural polynomials are correlated, so the estimated parameters will be interdependent. Given a particular time range and a maximum polynomial order, we can transform the polynomial time vectors to make them independent, that is, "orthogonal". This is illustrated in the figure below for linear and quadratic time in the range 1-10. Since orthogonal polynomial time terms are independent, the parameter estimates can be interpreted independently.
With orthogonal polynomials, the intercept term reflects the average overall curve height, rather than the height at the left edge of the time window, so if you are interested in differences at the very beginning of the time window, you may be better off sticking with natural polynomials.
A complete example
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.
The first step is to create a third-order polynomial in the range of
timeBin .> t <- poly((unique(TargetFix$timeBin)), 3) > TargetFix[,paste("ot", 1:3, sep="")] <- t[TargetFix$timeBin, 1:3] 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 ~ (ot1+ot2+ot3)*Condition + (ot1+ot2+ot3 | Subject) + (ot1+ot2 | Subject:Condition), control = lmerControl(optimizer="bobyqa"), data=TargetFix, REML=F)
> coef(summary(m.full))
Notice that the parameter estimates do not have p-values. There are good reasons for that (see this FAQ for more information), but this is cold comfort to most experimental psychologists, who need to report p-values. The quick and easy solution is to assume that, because we have relatively many observations, the t-distribution converges to the z-distribution, so we can use a normal approximation:
> coefs <- data.frame(coef(summary(m.full))) > coefs$p <- format.pval(2*(1-pnorm(abs(coefs$t.value))), digits=2, eps=0.0001) #also make the p-values a bit more readable > coefs Estimate Std..Error t.value p (Intercept) 0.4773227513 0.01385240 34.457775306 <1e-04 ot1 0.6385603705 0.05993519 10.654181583 <1e-04 ot2 -0.1095979256 0.03848819 -2.847573180 0.0044 ot3 -0.0932611870 0.02041640 -4.567955536 <1e-04 ConditionLow -0.0581122429 0.01901291 -3.056462582 0.0022 ot1:ConditionLow 0.0003188189 0.06330556 0.005036191 0.9960 ot2:ConditionLow 0.1635455113 0.05426498 3.013831365 0.0026 ot3:ConditionLow -0.0020869051 0.02014728 -0.103582452 0.9175 Another 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:
Chisq Df Pr(>Chisq) ConditionLow 4.1481 1 0.0417 * ot1:ConditionLow 3.4292 1 0.0641 . ot2:ConditionLow 7.4845 1 0.0062 ** ot3:ConditionLow 0.0107 1 0.9175 The results confirm that fixation proportion was lower overall for the Low condition than the High condition (Condition effect on the intercept - ConditionLow - had a statistically significant negative estimate) and there was a difference in curvature between the conditions (statistically significant effect of condition on the quadratic term:
ot2:ConditionLow )Reporting resultsThere are two rules of thumb for reporting growth curve analysis results:
Another exampleHere is another example based on data from an experiment examining the link between statistical learning and word learning (Mirman et al., 2008). In this experiment the critical manipulation was between-participants and the data shape is a little different (requiring only a second-order polynomial).
Analyzing individual differencesHere we will consider cohort and rhyme competition data from 5 participants with Broca's aphasia, 3 participants with Wernicke's aphasia, and 12 age-matched controls (Mirman et al., 2011). The group means are shown below, with cohort competition in the top row and rhyme competition in the bottom row. Visually, it seems that the participants with Wernicke's exhibited more cohort competition and less rhyme competition, whereas participants with Broca's aphasia exhibited the opposite pattern.
> summary(CohortRhyme)
> t <- poly(unique(CohortRhyme$timeBin), 4) > CohortRhyme[,paste("ot", 1:4, sep="")] <- t[CohortRhyme$timeBin, 1:4] When the individual differences are "experiment external", such as aphasia subtype (or continuous variables like age), the best option is to add the individual difference variable as a fixed effect to the model. So we start with a base model of the overall competition effect (Object fixed effect) and add a fixed effect of Group. It would be a good idea to build up the Group effects gradually, but we'll skip to the full model to keep this example more brief:
> cohort.base <- lmer(FixProp ~ (ot1+ot2+ot3+ot4)*Object + (1+ot1+ot2+ot3+ot4 | subjID) + (1+ot1+ot2 | subjID:Object), control = lmerControl(optimizer="bobyqa"), data=subset(CohortRhyme, Type == "Cohort"), REML=F) > cohort.group <- lmer(FixProp ~ (ot1+ot2+ot3+ot4)*Object*Group + (1+ot1+ot2+ot3+ot4 | subjID) + (1+ot1+ot2 | subjID:Object), control = lmerControl(optimizer="bobyqa"), data=subset(CohortRhyme, Type == "Cohort"), REML=F) > anova(cohort.base, cohort.group)
ObjectUnrelated:GroupBroca 0.007168833 0.016841402 0.4256673 0.67035
ObjectUnrelated:GroupWernicke -0.049674167 0.020423200 -2.4322421 0.01501
ot1:ObjectUnrelated:GroupBroca -0.047473215 0.059871004 -0.7929250 0.42782
ot1:ObjectUnrelated:GroupWernicke 0.009277717 0.072604257 0.1277848 0.89832
ot2:ObjectUnrelated:GroupBroca 0.044092440 0.058847796 0.7492624 0.45370
The same steps for the rhyme competition data reveal that only the Broca's group differed from controls:
> rhyme.base <- lmer(FixProp ~ (ot1+ot2+ot3+ot4)*Object + (1+ot1+ot2+ot3+ot4 | subjID) + (1+ot1+ot2 | subjID:Object), control = lmerControl(optimizer="bobyqa"), data=subset(CohortRhyme, Type == "Rhyme"), REML=F) > rhyme.group <- lmer(FixProp ~ (ot1+ot2+ot3+ot4)*Object*Group + (1+ot1+ot2+ot3+ot4 | subjID) + (1+ot1+ot2 | subjID:Object), control = lmerControl(optimizer="bobyqa"), data=subset(CohortRhyme, Type == "Rhyme"), REML=F) > anova(rhyme.base, rhyme.group) Df AIC BIC logLik deviance Chisq Chi Df Pr(>Chisq)
rhyme.base 32 -2269.6 -2119.7 1166.8 -2333.6
rhyme.group 52 -2265.2 -2021.6 1184.6 -2369.2 35.597 20 0.01715 *
> coefs.rhyme <- as.data.frame(coef(summary(rhyme.group))) > coefs.rhyme$p <- format.pval(2*(1-pnorm(abs(coefs.rhyme[,"t value"]))))
Estimate Std. Error t value p ObjectUnrelated:GroupBroca -0.035564333 0.015761359 -2.25642560 0.0240440
ObjectUnrelated:GroupWernicke 0.017755000 0.019113455 0.92892677 0.3529270
ot1:ObjectUnrelated:GroupBroca 0.122036460 0.064618493 1.88856865 0.0589497
ot1:ObjectUnrelated:GroupWernicke -0.128382147 0.078361433 -1.63833333 0.1013522
ot2:ObjectUnrelated:GroupBroca -0.006877810 0.055995864 -0.12282711 0.9022440
ot2:ObjectUnrelated:GroupWernicke -0.011219074 0.067904959 -0.16521729 0.8687730 The random effects are the systematic deviations from the "mean" pattern predicted by the fixed effects. So we'll use the random effects from the base model, which did not distinguish between participant groups so we can get individual effect sizes relative to the overall mean (otherwise we'd be looking at individual effect sizes relative to the sub-group mean). The key function for extracting random effects is
ranef and in our case it will be a list with two elements, because we had a random effect by subjID and a random effect by subjID:Object . We're interested in the latter, because it is the one that will capture individual differences in competition effect size. In particular, we'll want to subtract the random effect estimate for the Unrelated object from the estimate for the Competitor for each subject to get that subject's individual competition effect size. Extracting the random effects is relatively easy (e.g., ranef(cohort.base)$'subjID:Object' ), but there is a moderate amount of data manipulation involved in getting them into a convenient form (requiring the reshape2 and plyr packages). There are many ways to do the data manipulation, here is my way:#get cohort effect sizes > blup.cohort <- data.frame(colsplit(row.names(ranef(cohort.base)$'subjID:Object'), ":", c("Subject", "Object")), ranef(cohort.base)$'subjID:Object') > ES.coh <- ddply(blup.cohort, .(Subject), summarize, Intercept = X.Intercept.[Object=="Competitor"] - X.Intercept.[Object=="Unrelated"], Quadratic = ot2[Object=="Competitor"] - ot2[Object=="Unrelated"])
#get rhyme effect sizes > blup.rhyme <- data.frame(colsplit(row.names(ranef(rhyme.base)$'subjID:Object'), ":", c("Subject", "Object")), ranef(rhyme.base)$'subjID:Object') > ES.rhy <- ddply(blup.rhyme, .(Subject), summarize, Rhyme.Intercept = X.Intercept.[Object=="Competitor"] - X.Intercept.[Object=="Unrelated"])
#combine > ES <- merge(ES.coh, ES.rhy, by="Subject") > group <- unique(subset(CohortRhyme, select=c(subjID, Group))) #get group assignments from original data frame > ES <- merge(ES, group, by.x="Subject", by.y="subjID") ![]() > cor.test(ES$Intercept, ES$Rhyme.Intercept) > cor.test(ES$Intercept[ES$Group != "Control"], ES$Rhyme.Intercept[ES$Group != "Control"]) > cor.test(ES$Intercept[ES$Group == "Control"], ES$Rhyme.Intercept[ES$Group == "Control"]) The correlation is significant for participants with aphasia (r = -0.86, p = 0.006), but not for control participants (r = 0.34, p > 0.25). This suggests that there may be a single mechanism behind this pattern, in contrast to the standard separable components interpretation of a double dissociation.
Another exampleHere is another example based on data from an experiment that examined the time course of activation of function and thematic semantic knowledge in 17 participants with left hemisphere stroke (Kalenine, Mirman, & Buxbaum, 2012).
I highly recommend the ggplot2 package for plotting (if you are not already familiar with it, you may want to start with a tutorial and the online documentation is a useful reference). To start, plot the target fixation data discussed above:
ggplot(TargetFix, aes(Time, meanFix, color=Condition)) + stat_summary(fun.data=mean_se, geom="pointrange") + labs(y="Fixation Proportion", x="Time since word onset (ms)") + theme_bw() + scale_color_manual(values=c("red", "blue")) You can now add the model fits by adding another stat_summary with the y-variable mapped to the predicted values in the lmer model object, which can be accessed using the
fitted() function:last_plot() + stat_summary(aes(y=fitted(m.full)), fun.y=mean, geom="line")
The other aesthetics (mappings) are automatically inherited (e.g., mapping of color to Condition), so the plot will automatically remain consistent. This is particularly useful for multi-panel plots as in the cohort and rhyme competition example from above. When the model was fit to a data subset, it's not a good idea to rely on ggplot to align the data and the model-predicted values. Instead, create a new data frame that has just the data subset and add the model-predicted values to that data frame:
dat <- subset(CohortRhyme, Type == "Cohort") dat$mfit <- fitted(cohort.group) ggplot(dat, aes(Time, FixProp, color=Object)) + facet_wrap(~ Group) + theme_bw() + stat_summary(fun.y=mean, geom="point") + stat_summary(aes(y=mfit), fun.y=mean, geom="line") + labs(y="Fixation Proportion", x="Time since word onset (ms)") + scale_color_manual(values=c("red", "blue")) |