Applying GCA to VWP Data: Example 2

Individual Differences

Dan Mirman

2020-08-22

Preliminaries

library(tidyverse)
library(lme4)
library(lmerTest)
library(gazer) # helper functions and example data sets

Mirman, D., Yee, E., Blumstein, S., and Magnuson, J.S. (2011). Theories of spoken word recognition deficits in aphasia: Evidence from eye-tracking and computational modeling. Brain and Language, 117, 53-68. pdf

These are 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 aphasia exhibited more cohort competition and less rhyme competition, whereas participants with Broca’s aphasia exhibited the opposite pattern.

The starting data frame is shown below.

summary(CohortRhyme)
##      subjID          Group          Time         timeBin             Object   
##  101    :  80   Control :960   Min.   :  50   Min.   : 1.00   Competitor:800  
##  102    :  80   Broca   :400   1st Qu.: 525   1st Qu.: 5.75   Unrelated :800  
##  103    :  80   Wernicke:240   Median :1000   Median :10.50   Target    :  0  
##  104    :  80                  Mean   :1000   Mean   :10.50                   
##  105    :  80                  3rd Qu.:1475   3rd Qu.:15.25                   
##  106    :  80                  Max.   :1950   Max.   :20.00                   
##  (Other):1120                                                                 
##      Type        FixProp      
##  Cohort:800   Min.   :0.0000  
##  Rhyme :800   1st Qu.:0.0000  
##               Median :0.0417  
##               Mean   :0.0665  
##               3rd Qu.:0.1000  
##               Max.   :0.5833  
## 

The first step is to create a fourth-order polynomial based on the timeBin values to capture the rise and fall of the competition curves:

CohortRhyme.gca <- code_poly(CohortRhyme, predictor = "timeBin", poly.order = 4, draw.poly = FALSE)

Group differences

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 ~ (poly1+poly2+poly3+poly4)*Object + 
                      (1+poly1+poly2+poly3+poly4 | subjID) + 
                      (1+poly1+poly2 | subjID:Object), 
                    data=subset(CohortRhyme.gca, Type == "Cohort"), REML=F)
## boundary (singular) fit: see ?isSingular
cohort.group <- lmer(FixProp ~ (poly1+poly2+poly3+poly4)*Object*Group +
                       (1+poly1+poly2+poly3+poly4 | subjID) + 
                       (1+poly1+poly2 | subjID:Object), 
                     data=subset(CohortRhyme.gca, Type == "Cohort"), REML=F)
## boundary (singular) fit: see ?isSingular
## Warning: Model failed to converge with 1 negative eigenvalue: -2.2e-01
anova(cohort.base, cohort.group)
## Data: subset(CohortRhyme.gca, Type == "Cohort")
## Models:
## cohort.base: FixProp ~ (poly1 + poly2 + poly3 + poly4) * Object + (1 + poly1 + 
## cohort.base:     poly2 + poly3 + poly4 | subjID) + (1 + poly1 + poly2 | subjID:Object)
## cohort.group: FixProp ~ (poly1 + poly2 + poly3 + poly4) * Object * Group + 
## cohort.group:     (1 + poly1 + poly2 + poly3 + poly4 | subjID) + (1 + poly1 + 
## cohort.group:     poly2 | subjID:Object)
##              npar     AIC     BIC logLik deviance  Chisq Df Pr(>Chisq)    
## cohort.base    32 -2277.6 -2127.7 1170.8  -2341.6                         
## cohort.group   52 -2303.6 -2060.0 1203.8  -2407.6 65.997 20  8.089e-07 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

From the critical Object:Group parameter estimates below, we can see that this overall Group difference is due to the Wernicke’s aphasia group differing from the Control (baseline reference) group. We’re just going to look at the effects up to the quadratic term because the higher order terms are just capturing differences in the tails. Note that the parameter estimates are for the Unrelated object relative to the Related object, so a negative estimate for the Intercept term (e.g., ObjectUnrelated:GroupWernicke) corresponds to a larger competition effect (i.e., for the Wernicke’s group compared to the Control group, the Unrelated object has an even lower intercept than the Competitor does).

##                                         Estimate Std. Error        df
## ObjectUnrelated:GroupBroca           0.007168833 0.01684196  22.85703
## ObjectUnrelated:GroupWernicke       -0.049674167 0.02042388  22.85703
## poly1:ObjectUnrelated:GroupBroca    -0.047473215 0.05987084  36.03035
## poly1:ObjectUnrelated:GroupWernicke  0.009277717 0.07260405  36.03035
## poly2:ObjectUnrelated:GroupBroca     0.044092440 0.05884975  41.18294
## poly2:ObjectUnrelated:GroupWernicke  0.279205510 0.07136580  41.18294
## poly3:ObjectUnrelated:GroupBroca    -0.043153806 0.03734619 702.60837
## poly3:ObjectUnrelated:GroupWernicke -0.050128667 0.04528891 702.60837
## poly4:ObjectUnrelated:GroupBroca     0.025128978 0.03734619 702.60837
## poly4:ObjectUnrelated:GroupWernicke -0.177046393 0.04528891 702.60837
##                                        t value     Pr(>|t|)
## ObjectUnrelated:GroupBroca           0.4256531 0.6743417917
## ObjectUnrelated:GroupWernicke       -2.4321613 0.0232575583
## poly1:ObjectUnrelated:GroupBroca    -0.7929272 0.4330107992
## poly1:ObjectUnrelated:GroupWernicke  0.1277851 0.8990295952
## poly2:ObjectUnrelated:GroupBroca     0.7492376 0.4579725128
## poly2:ObjectUnrelated:GroupWernicke  3.9123153 0.0003353343
## poly3:ObjectUnrelated:GroupBroca    -1.1555075 0.2482754687
## poly3:ObjectUnrelated:GroupWernicke -1.1068640 0.2687315346
## poly4:ObjectUnrelated:GroupBroca     0.6728658 0.5012538929
## poly4:ObjectUnrelated:GroupWernicke -3.9092657 0.0001015433

The same steps for the rhyme competition data reveal that only the Broca’s group differed from controls:

rhyme.base <- lmer(FixProp ~ (poly1+poly2+poly3+poly4)*Object + 
                      (1+poly1+poly2+poly3+poly4 | subjID) + 
                      (1+poly1+poly2 | subjID:Object), 
                    data=subset(CohortRhyme.gca, Type == "Rhyme"), REML=F)
## boundary (singular) fit: see ?isSingular
rhyme.group <- lmer(FixProp ~ (poly1+poly2+poly3+poly4)*Object*Group +
                       (1+poly1+poly2+poly3+poly4 | subjID) + 
                       (1+poly1+poly2 | subjID:Object), 
                     data=subset(CohortRhyme.gca, Type == "Rhyme"), REML=F)
## boundary (singular) fit: see ?isSingular
## Warning: Model failed to converge with 1 negative eigenvalue: -2.7e+01
anova(rhyme.base, rhyme.group)
## Data: subset(CohortRhyme.gca, Type == "Rhyme")
## Models:
## rhyme.base: FixProp ~ (poly1 + poly2 + poly3 + poly4) * Object + (1 + poly1 + 
## rhyme.base:     poly2 + poly3 + poly4 | subjID) + (1 + poly1 + poly2 | subjID:Object)
## rhyme.group: FixProp ~ (poly1 + poly2 + poly3 + poly4) * Object * Group + 
## rhyme.group:     (1 + poly1 + poly2 + poly3 + poly4 | subjID) + (1 + poly1 + 
## rhyme.group:     poly2 | subjID:Object)
##             npar     AIC     BIC logLik deviance  Chisq Df Pr(>Chisq)  
## rhyme.base    32 -2269.6 -2119.7 1166.8  -2333.6                       
## rhyme.group   52 -2263.2 -2019.6 1183.6  -2367.2 33.646 20    0.02862 *
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
coef(summary(rhyme.group))[21:30,]
##                                         Estimate Std. Error        df
## ObjectUnrelated:GroupBroca          -0.035564333 0.01690859  39.83764
## ObjectUnrelated:GroupWernicke        0.017755000 0.02050467  39.83764
## poly1:ObjectUnrelated:GroupBroca     0.122036460 0.06499866  37.19534
## poly1:ObjectUnrelated:GroupWernicke -0.128382147 0.07882246  37.19534
## poly2:ObjectUnrelated:GroupBroca    -0.006877810 0.05623612  35.52041
## poly2:ObjectUnrelated:GroupWernicke -0.011219074 0.06819632  35.52041
## poly3:ObjectUnrelated:GroupBroca    -0.059553382 0.03763343 682.46842
## poly3:ObjectUnrelated:GroupWernicke -0.020485110 0.04563724 682.46842
## poly4:ObjectUnrelated:GroupBroca     0.016535910 0.03763343 682.46842
## poly4:ObjectUnrelated:GroupWernicke -0.009081926 0.04563724 682.46842
##                                        t value   Pr(>|t|)
## ObjectUnrelated:GroupBroca          -2.1033296 0.04179987
## ObjectUnrelated:GroupWernicke        0.8659001 0.39173002
## poly1:ObjectUnrelated:GroupBroca     1.8775226 0.06830529
## poly1:ObjectUnrelated:GroupWernicke -1.6287508 0.11180762
## poly2:ObjectUnrelated:GroupBroca    -0.1223024 0.90334892
## poly2:ObjectUnrelated:GroupWernicke -0.1645114 0.87026117
## poly3:ObjectUnrelated:GroupBroca    -1.5824594 0.11400802
## poly3:ObjectUnrelated:GroupWernicke -0.4488683 0.65366919
## poly4:ObjectUnrelated:GroupBroca     0.4393941 0.66051507
## poly4:ObjectUnrelated:GroupWernicke -0.1990025 0.84232006

So far, we have established a double dissociation – one group has larger cohort competition, the other group has larger rhyme competition – but we can also ask whether there is an association between the effects. That is, do participants with larger cohort effects tend to show smaller rhyme effects? To test this kind of “experiment internal” individual differences, we can use the random effects to estimate each participants effect size and then test the correlation between the effect sizes.

Individual differences

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). You can extract random effects from a model using the lme4::ranef function, but it’s somewhat more convenient to use gazer::get_ranef, which allows you to specify which random effects you want and makes them a little more convenient to work with:

re.cohort <- get_ranef(cohort.base, 'subjID:Object')
re.rhyme <- get_ranef(rhyme.base, 'subjID:Object')

Both data frames will have the same structure:

head(re.cohort)
##      Intercept        poly1        poly2 subjID     Object
## 1 -0.005219780  0.001170704  0.019152320    101 Competitor
## 2  0.009802661 -0.032138493 -0.007075727    101  Unrelated
## 3 -0.005881789  0.096949384 -0.070701728    102 Competitor
## 4  0.007690752 -0.039752372  0.008477722    102  Unrelated
## 5  0.042558550 -0.046323959 -0.120663482    103 Competitor
## 6 -0.012434554  0.066828944 -0.016174013    103  Unrelated

To estimate individual participant effect sizes, we need to compute the differences between random effect estimates (Intercept, poly1, poly2) for Competitor vs. Unrelated objects for each participant:

ES.coh <- re.cohort %>% group_by(subjID) %>% 
  summarise(Coh.Intercept = Intercept[Object=="Competitor"] -
              Intercept[Object=="Unrelated"], 
            Coh.Quadratic = poly2[Object=="Competitor"] -
              poly2[Object=="Unrelated"])
ES.rhy <- re.rhyme %>% group_by(subjID) %>% 
  summarise(Rhy.Intercept = Intercept[Object=="Competitor"] -
              Intercept[Object=="Unrelated"], 
            Rhy.Quadratic = poly2[Object=="Competitor"] -
              poly2[Object=="Unrelated"])
ES <- merge(ES.coh, ES.rhy, by="subjID")
group <- unique(subset(CohortRhyme, select=c(subjID, Group))) #get group assignments from original data frame
ES <- merge(ES, group, by="subjID")
summary(ES)
##     subjID          Coh.Intercept       Coh.Quadratic      Rhy.Intercept      
##  Length:20          Min.   :-0.049232   Min.   :-0.30990   Min.   :-0.046640  
##  Class :character   1st Qu.:-0.017123   1st Qu.:-0.08317   1st Qu.:-0.017738  
##  Mode  :character   Median :-0.009845   Median : 0.02630   Median :-0.009355  
##                     Mean   : 0.000000   Mean   : 0.00000   Mean   : 0.000000  
##                     3rd Qu.: 0.021994   3rd Qu.: 0.08961   3rd Qu.: 0.016755  
##                     Max.   : 0.068690   Max.   : 0.14444   Max.   : 0.067391  
##  Rhy.Quadratic           Group   
##  Min.   :-0.12373   Control :12  
##  1st Qu.:-0.05222   Broca   : 5  
##  Median : 0.01677   Wernicke: 3  
##  Mean   : 0.00000                
##  3rd Qu.: 0.05412                
##  Max.   : 0.11931

Now we can examine a scatterplot of the individual cohort and rhyme competition effect sizes (effects on the intercept, shown on the right) and compute correlations across the complete sample and just for the participants with aphasia:

with(ES, cor.test(Coh.Intercept, Rhy.Intercept))
## 
##  Pearson's product-moment correlation
## 
## data:  Coh.Intercept and Rhy.Intercept
## t = -1.7573, df = 18, p-value = 0.09586
## alternative hypothesis: true correlation is not equal to 0
## 95 percent confidence interval:
##  -0.70569266  0.07204616
## sample estimates:
##       cor 
## -0.382675
with(subset(ES, Group != "Control"), cor.test(Coh.Intercept, Rhy.Intercept))
## 
##  Pearson's product-moment correlation
## 
## data:  Coh.Intercept and Rhy.Intercept
## t = -4.138, df = 6, p-value = 0.006092
## alternative hypothesis: true correlation is not equal to 0
## 95 percent confidence interval:
##  -0.9743601 -0.3959868
## sample estimates:
##       cor 
## -0.860535
with(subset(ES, Group == "Control"), cor.test(Coh.Intercept, Rhy.Intercept))
## 
##  Pearson's product-moment correlation
## 
## data:  Coh.Intercept and Rhy.Intercept
## t = 1.1352, df = 10, p-value = 0.2828
## alternative hypothesis: true correlation is not equal to 0
## 95 percent confidence interval:
##  -0.2928018  0.7636921
## sample estimates:
##       cor 
## 0.3378764

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.