Basic Design:

This code conducts a mixed levels effects analysis:

Group-level design (2x2 cell):

                             AE                              UAA
            ----------------------------------------------------------------  
   Pre (1) | Parameter Estimate @ time1,   | Parameter Estimate @ time1,   |  
           | Parameter Estimate @ time2... | Parameter Estimate @ time2... |  
           | Parameter Estimate @ time 11  | Parameter Estimate @ time 11  |  
           |---------------------------------------------------------------|  
  Post (2) | Parameter Estimate @ time1,   | Parameter Estimate @ time1,   |   
           | Parameter Estimate @ time2... | Parameter Estimate @ time2... |  
           | Parameter Estimate @ time 11  | Parameter Estimate @ time 11  |  
            ----------------------------------------------------------------  

With:
Time (Pre =1 or Post =2) = within subject or repeated measures factor (fixed factor)
Intervention (AE or UAA group) = between subject factor (fixed factor)

@ Subject-level:

Subject_ID = blocking factor (random factor)

Parameter Estimate at each time point = measured (dependent) variable

##References

  1. https://ase.tufts.edu/gsc/gradresources/guidetomixedmodelsinr/mixed%20model%20guide.html
  2. https://www.theanalysisfactor.com/the-difference-between-crossed-and-nested-factors/
  3. https://www.r-bloggers.com/how-to-get-the-frequency-table-of-a-categorical-variable-as-a-data-frame-in-r/
  4. https://stats.stackexchange.com/questions/98958/plots-to-illustrate-results-of-linear-mixed-effect-model
  5. https://stackoverflow.com/questions/31075407/plot-mixed-effects-model-in-ggplot
  6. https://stats.stackexchange.com/questions/115090/why-do-i-get-zero-variance-of-a-random-effect-in-my-mixed-model-despite-some-va

##Library

## Loading required package: Matrix
## Loading required package: carData
## Registered S3 methods overwritten by 'car':
##   method                          from
##   influence.merMod                lme4
##   cooks.distance.influence.merMod lme4
##   dfbeta.influence.merMod         lme4
##   dfbetas.influence.merMod        lme4

##Assumptions (adapted from Reference 1)

STEP 1: What distribution best fits the data?

We want to pick the distribution for which the largest number of observations falls between the dashed lines.

Standard distributions tested:

  1. Normal: Syntax: qqp(dataset$columnName, “norm”)

  2. Lognormal: Syntax: qqp(dataset$columnName, “lnorm”)

  3. Neg binomial: Syntax: fitdistr(dataset$columnName, “Negative Binomial”)

  4. Poisson: Syntax: fitdistr(dataset$columnName, “Poisson”)

  5. Gamma: Syntax: fitdistr(dataset$columnName, “gamma”)


SUMMARY CHECKLIST FOR STEP 1 (current study):

Draft Checklist: What distribution best fits the data?

  1. qqp(dataset$columnName, “norm”)

  2. qqp(dataset$columnName, “lnorm”)

  3. fitdistr(dataset$columnName, “Negative Binomial”)
    A neg binomial model doesn’t make sense here, since a negative binomial is a kind of count model. The response is supposed to be counts. A count, by definition, cannot be a fractional value (i.e.,count=non-negative whole numbers).

  4. fitdistr(dataset$columnName, “Poisson”)
    Also, doesn’t make sense here, since a poisson distribution can only handle positive whole numbers.

  5. [add 1 to dataset$columnName for gamma dist. to work; gama dist. requires values >=0]
  • fitdistr(dataset$columnName+1, “gamma”)
  • qqp(dataset$columnName+1, “gamma”, shape = gamma$estimate[[1]], rate=gamma$estimate[[2]])


FINAL CHECKLIST [For the current study]: What distribution best fits the data?

  1. qqp(dataset$columnName, “norm”)

  2. qqp(dataset$columnName, “lnorm”)

  3. gamma func:
  • fitdistr(dataset$Add1toColumnName, “gamma”)
  • qqp(dataset$Add1toColumnName, “gamma”, shape = gamma$estimate[[1]], rate = gamma$estimate[[2]])




STEP 2: How should the mixed model be fit to the data?

First, a note: if your data best fit the lognormal distribution, do not transform them. This is true for any type of transformation you might apply to your data to make them normal. If you can transform your data to normality, common wisdom says you should use the transformed data.More recent statistics literature has entirely changed stance on this matter, however, because transformation makes interpretation of model results more difficult, and it makes mischief with the variance of the transformed variable. Even if your data are transformable to normality, they are still not normal, and you should move on to tnon-linear models.

If data are normally distributed:

  • We can use a linear mixed model (LMM). We will load the lme4 package and make a call to the function lmer.

  • Syntax breakdown:

  1. The first argument to the function is a formula that takes the form y ~ x1 + x2 … etc., where y is the response variable and x1, x2, etc. are explanatory variables.

  2. Random effects are added in with the explanatory variables.

  3. Crossed random effects take the form (1 | r1) + (1 | r2), while nested random effects take the form (1 | r1 / r2).

  4. The next argument is where we designate the data frame your variables come from.

  5. The next argument allows us to designate whether the mixed model will estimate the parameters using maximum likelihood or restricted maximum likelihood.

  • If the random effects are nested, or we have only one random effect, and if the data are balanced (i.e., similar sample sizes in each factor group) we can set REML to FALSE, because you can use maximum likelihood.

  • If the random effects are crossed, we don’t set the REML argument because it defaults to TRUE anyway.


SUMMARY CHECKLIST FOR STEP 2:

Draft Checklist: How should the mixed model be fit to the data?

  1. Fitting the data: are the data normally distributed? -> If no, consider nlme. If yes, use lmer.

  2. Should we include crossed factors and/or nested factors? For the current study, because subjects are crossed across time (-although they are nested within group-), we treat them as crossed factors. (Reference 2)

  3. Should we use REML or ML?


Subchecklist:

  1. Nested random effects? -> If no (for current study: No), then go to step 2

  2. Only one random effect? -> If yes (for current study: Yes, subject) then go to step 3 [If no, do nothing]

  3. Similar sample sizes? (or are data balanced?) -> If yes, then set REML to false (else do nothing).


FINAL Checklist [For the current study]: How should the mixed model be fit to the data?

  1. Fitting the data: are the data normally distributed? -> If no, consider nlme. If yes, use lmer.

  2. Should we use REML or ML? Are the data balanced?/similar sample sizes for each fixed factor? -> If yes, then set REML to false (else do nothing).




STEP 3: Interpreting Results

  1. First we get some measures of model fit, including AIC, BIC, log likelihood, and deviance.

  2. Then we get an estimate of the variance explained by the random effect. This number is important, because if it’s indistinguishable from zero, then your random effect probably doesn’t matter and you can go ahead and do a regular linear model instead. The lesson from this model output is that although there is “obviously” variation in subject performance, the extent of this subject variation can be fully or virtually-fully explained by just the residual variance term alone. There is not enough additional subject-level variation to warrant adding an additional subject-level random effect to explain all the observed variation.
  • Gedankenexperiment: Imagine we are simulating experimental data under this same paradigm. We set up the parameters so that there is residual variation on a trial-by-trial basis, but 0 subject-level variation (i.e., all subjects have the same “true mean,” plus error). Now each time we simulate data from this set of parameters, we will of course find that subjects do not have exactly equal performance. Some end up with low scores, some with high scores. But this is all just because of the residual trial-level variation. We “know” (by virtue of having determined the simulation parameters) that there is not really any subject-level variation.
  • Because we perform a repeated measures analysis, subjects are tested twice, accounting for “not-completely-zero” variance/std.deviations or random effects. (based on smillig’s response to the question posted here: https://stats.stackexchange.com/questions/115090/why-do-i-get-zero-variance-of-a-random-effect-in-my-mixed-model-despite-some-va)
  1. Next we have estimates of the fixed effects, with standard errors. Some journals like you to report the results of these models as effect sizes with confidence intervals. Some journals want you to report p-values.The creators of the lme4 package are philosophically opposed to p-values, so we use anova from library(car). The Anova function does a Wald test, which tells us how confident we are of our estimate of the effect of groupID and time on time1, and the p-value tells us that we should not be confident at all.: Syntax: Anova(lmm)


SUMMARY CHECKLIST FOR STEP 3 (current study):

FINAL Checklist [For the current study]: How do we interpret the results?

  1. Is the variance related to the random effect close to zero? If not, consider a non-linear model.

  2. For p-value estimates of the fixed effects, use Anova(lmm). Note that for significant effects, if groupIDUAA:timeIDpre > 0, then as groupID changes to AE and timeID changes to post, AE increases (symbolized as AE^ in the Summary Table Tab). Put differently, AE increases from Pre to Post, while UAA decreases from Pre to Post.


Dataset

str(MixedANOVA_BraverSeed)
## 'data.frame':    264 obs. of  16 variables:
##  $ subID  : Factor w/ 142 levels "sub01","sub04",..: 16 17 18 21 23 25 26 27 29 33 ...
##  $ groupID: Factor w/ 2 levels "AE","UAA": 1 1 1 1 1 1 1 1 1 1 ...
##  $ timeID : Factor w/ 2 levels "post","pre": 2 2 2 2 2 2 2 2 2 2 ...
##  $ seeds  : Factor w/ 1 level "Braver": 1 1 1 1 1 1 1 1 1 1 ...
##  $ ROI    : Factor w/ 1 level "BA6_supFrontal_L_MNIxneg47yneg9z44": 1 1 1 1 1 1 1 1 1 1 ...
##  $ time1  : num  -0.134 -0.0827 0.0946 -0.1433 0.0318 ...
##  $ time2  : num  0.22457 -0.067 0.41948 0.16636 0.00399 ...
##  $ time3  : num  -0.488 -1.264 2.058 0.801 -0.478 ...
##  $ time4  : num  -0.493 -1.861 3.116 1.127 -0.972 ...
##  $ time5  : num  -0.639 -1.569 3.154 1.173 -0.837 ...
##  $ time6  : num  -0.129 -1.011 1.372 0.77 -0.536 ...
##  $ time7  : num  0.1257 -0.2327 0.5063 0.3124 -0.0327 ...
##  $ time8  : num  0.2332 -0.1194 0.148 0.3265 -0.0516 ...
##  $ time9  : num  -0.0214 -0.0975 0.0168 0.0796 -0.2316 ...
##  $ time10 : num  0.2211 0.0723 0.084 0.1224 0.2351 ...
##  $ time11 : num  0.1688 -0.0532 -0.164 -0.0547 -0.1136 ...

Timepoints

Time3

  • STEP1 Checklist:
MixedANOVA_BraverSeed$t3 <- MixedANOVA_BraverSeed$time3 + 1
qqp(MixedANOVA_BraverSeed$t3, "norm")

## [1] 209 142
qqp(MixedANOVA_BraverSeed$t3, "lnorm")

## [1] 209 142
# For the gamma distr, adding 1 still gives error, so removed this.
#gamma <- fitdistr(MixedANOVA_BraverSeed$t3, "gamma")
#qqp(MixedANOVA_BraverSeed$t3, "gamma", shape = gamma$estimate[[1]], rate = gamma$estimate[[2]])

#####Verdict: Normal distribution is best suited.


  • STEP2 Checklist:
# 1. Use lmer 
# 2a. Check if data are balanced across groups:
count(MixedANOVA_BraverSeed, 'groupID')
##   groupID freq
## 1      AE  130
## 2     UAA  134
# 2b. Check if data are balanced across time:
count(MixedANOVA_BraverSeed, 'timeID') 
##   timeID freq
## 1   post  131
## 2    pre  133

#####Verdict: Use lmer, set REML=false.


  • Model:
lmm <- lmer(time3 ~ groupID * timeID + (1 | subID), data = MixedANOVA_BraverSeed,
            REML = FALSE)
summary(lmm)
## Linear mixed model fit by maximum likelihood  ['lmerMod']
## Formula: time3 ~ groupID * timeID + (1 | subID)
##    Data: MixedANOVA_BraverSeed
## 
##      AIC      BIC   logLik deviance df.resid 
##    990.6   1012.1   -489.3    978.6      258 
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -3.4029 -0.5701  0.0261  0.5723  3.8613 
## 
## Random effects:
##  Groups   Name        Variance Std.Dev.
##  subID    (Intercept) 0.1087   0.3296  
##  Residual             2.2780   1.5093  
## Number of obs: 264, groups:  subID, 142
## 
## Fixed effects:
##                      Estimate Std. Error t value
## (Intercept)            0.5698     0.1901   2.996
## groupIDUAA            -0.5917     0.2699  -2.192
## timeIDpre             -0.6337     0.2653  -2.388
## groupIDUAA:timeIDpre   0.8430     0.3724   2.264
## 
## Correlation of Fixed Effects:
##             (Intr) grIDUAA tmIDpr
## groupIDUAA  -0.704               
## timeIDpre   -0.687  0.484        
## grpIDUAA:ID  0.489 -0.695  -0.713


#####Verdict: Use lmer, set REML=false.


  • STEP3 Checklist:
#1. Close to 0, linear model okay.
# 2.
Anova(lmm)
## Analysis of Deviance Table (Type II Wald chisquare tests)
## 
## Response: time3
##                 Chisq Df Pr(>Chisq)  
## groupID        0.7415  1    0.38917  
## timeID         1.2204  1    0.26927  
## groupID:timeID 5.1253  1    0.02358 *
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1


#####Verdict: linear-model okay, significant interactions (p=0.02). Note that groupIDUAA:timeIDpre > 0. That is, as groupID changes to AE and timeID changes to post, AE increases (symbolized as AE^ in the Summary Table Tab). Put differently, AE increases from Pre to Post, while UAA decreases from Pre to Post.

Time4

  • STEP1 Checklist:
MixedANOVA_BraverSeed$t4 <- MixedANOVA_BraverSeed$time4 + 1
qqp(MixedANOVA_BraverSeed$t4, "norm")

## [1] 209 142
qqp(MixedANOVA_BraverSeed$t4, "lnorm")

## [1] 209 142
# For the gamma distr, adding 1 still gives error, so removed this.
#gamma <- fitdistr(MixedANOVA_BraverSeed$t4, "gamma")
#qqp(MixedANOVA_BraverSeed$t4, "gamma", shape = gamma$estimate[[1]], rate = gamma$estimate[[2]])

#####Verdict: Normal distribution is best suited.


  • STEP2 Checklist:
# 1. Use lmer 
# 2a. Check if data are balanced across groups:
count(MixedANOVA_BraverSeed, 'groupID')
##   groupID freq
## 1      AE  130
## 2     UAA  134
# 2b. Check if data are balanced across time:
count(MixedANOVA_BraverSeed, 'timeID') 
##   timeID freq
## 1   post  131
## 2    pre  133

#####Verdict: Use lmer, set REML=false.


  • Model:
lmm <- lmer(time4 ~ groupID * timeID + (1 | subID), data = MixedANOVA_BraverSeed,
            REML = FALSE)
summary(lmm)
## Linear mixed model fit by maximum likelihood  ['lmerMod']
## Formula: time4 ~ groupID * timeID + (1 | subID)
##    Data: MixedANOVA_BraverSeed
## 
##      AIC      BIC   logLik deviance df.resid 
##   1240.5   1262.0   -614.3   1228.5      258 
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -2.9544 -0.6187  0.0379  0.6115  3.6910 
## 
## Random effects:
##  Groups   Name        Variance Std.Dev.
##  subID    (Intercept) 0.1358   0.3685  
##  Residual             6.0101   2.4516  
## Number of obs: 264, groups:  subID, 142
## 
## Fixed effects:
##                      Estimate Std. Error t value
## (Intercept)            1.1367     0.3052   3.725
## groupIDUAA            -1.0087     0.4332  -2.329
## timeIDpre             -1.2791     0.4305  -2.971
## groupIDUAA:timeIDpre   1.4231     0.6043   2.355
## 
## Correlation of Fixed Effects:
##             (Intr) grIDUAA tmIDpr
## groupIDUAA  -0.704               
## timeIDpre   -0.694  0.489        
## grpIDUAA:ID  0.495 -0.702  -0.712


#####Verdict: Use lmer, set REML=false.


  • STEP3 Checklist:
#1. Close to 0, linear model okay.
# 2.
Anova(lmm)
## Analysis of Deviance Table (Type II Wald chisquare tests)
## 
## Response: time4
##                 Chisq Df Pr(>Chisq)  
## groupID        0.8978  1    0.34336  
## timeID         3.3965  1    0.06533 .
## groupID:timeID 5.5461  1    0.01852 *
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1


#####Verdict: linear-model okay, significant interactions (p=0.02). Note that groupIDUAA:timeIDpre > 0. That is, as groupID changes to AE and timeID changes to post, AE increases (symbolized as AE^ in the Summary Table Tab). Put differently, AE increases from Pre to Post, while UAA decreases from Pre to Post.

Time5

  • STEP1 Checklist:
MixedANOVA_BraverSeed$t5 <- MixedANOVA_BraverSeed$time5 + 1
qqp(MixedANOVA_BraverSeed$t5, "norm")

## [1] 209 142
qqp(MixedANOVA_BraverSeed$t5, "lnorm")

## [1] 209 142
# For the gamma distr, adding 1 still gives error, so removed this.
# gamma <- fitdistr(MixedANOVA_BraverSeed$t5, "gamma")
# qqp(MixedANOVA_BraverSeed$t5, "gamma", shape = gamma$estimate[[1]], rate = gamma$estimate[[2]])

#####Verdict: Normal distribution is best suited.


  • STEP2 Checklist:
# 1. Use lmer 
# 2a. Check if data are balanced across groups:
count(MixedANOVA_BraverSeed, 'groupID')
##   groupID freq
## 1      AE  130
## 2     UAA  134
# 2b. Check if data are balanced across time:
count(MixedANOVA_BraverSeed, 'timeID') 
##   timeID freq
## 1   post  131
## 2    pre  133

#####Verdict: Use lmer, set REML=false.


  • Model:
lmm <- lmer(time5 ~ groupID * timeID + (1 | subID), data = MixedANOVA_BraverSeed,
            REML = FALSE)
summary(lmm)
## Linear mixed model fit by maximum likelihood  ['lmerMod']
## Formula: time5 ~ groupID * timeID + (1 | subID)
##    Data: MixedANOVA_BraverSeed
## 
##      AIC      BIC   logLik deviance df.resid 
##   1239.5   1261.0   -613.8   1227.5      258 
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -3.2983 -0.5580  0.0161  0.5893  3.8666 
## 
## Random effects:
##  Groups   Name        Variance Std.Dev.
##  subID    (Intercept) 0.1651   0.4063  
##  Residual             5.9592   2.4412  
## Number of obs: 264, groups:  subID, 142
## 
## Fixed effects:
##                      Estimate Std. Error t value
## (Intercept)            0.9913     0.3046   3.254
## groupIDUAA            -0.9591     0.4324  -2.218
## timeIDpre             -1.1426     0.4288  -2.665
## groupIDUAA:timeIDpre   1.3600     0.6018   2.260
## 
## Correlation of Fixed Effects:
##             (Intr) grIDUAA tmIDpr
## groupIDUAA  -0.704               
## timeIDpre   -0.693  0.488        
## grpIDUAA:ID  0.494 -0.701  -0.712


#####Verdict: Use lmer, set REML=false.


  • STEP3 Checklist:
#1. Close to 0, linear model okay.
# 2.
Anova(lmm)
## Analysis of Deviance Table (Type II Wald chisquare tests)
## 
## Response: time5
##                 Chisq Df Pr(>Chisq)  
## groupID        0.7906  1    0.37391  
## timeID         2.2591  1    0.13284  
## groupID:timeID 5.1063  1    0.02384 *
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1


#####Verdict: linear-model okay, significant interactions (p=0.02). Note that groupIDUAA:timeIDpre > 0. That is, as groupID changes to AE and timeID changes to post, AE increases (symbolized as AE^ in the Summary Table Tab). Put differently, AE increases from Pre to Post, while UAA decreases from Pre to Post.

Time6

  • STEP1 Checklist:
MixedANOVA_BraverSeed$t6 <- MixedANOVA_BraverSeed$time6 + 1
qqp(MixedANOVA_BraverSeed$t6, "norm")

## [1] 209  11
qqp(MixedANOVA_BraverSeed$t6, "lnorm")

## [1] 209  11
# For the gamma distr, adding 1 still gives error, so removed this.
# gamma <- fitdistr(MixedANOVA_BraverSeed$t6, "gamma")
# qqp(MixedANOVA_BraverSeed$t6, "gamma", shape = gamma$estimate[[1]], rate = gamma$estimate[[2]])

#####Verdict: Normal distribution is best suited.


  • STEP2 Checklist:
# 1. Use lmer 
# 2a. Check if data are balanced across groups:
count(MixedANOVA_BraverSeed, 'groupID')
##   groupID freq
## 1      AE  130
## 2     UAA  134
# 2b. Check if data are balanced across time:
count(MixedANOVA_BraverSeed, 'timeID') 
##   timeID freq
## 1   post  131
## 2    pre  133

#####Verdict: Use lmer, set REML=false.


  • Model:
lmm <- lmer(time6 ~ groupID * timeID + (1 | subID), data = MixedANOVA_BraverSeed,
            REML = FALSE)
summary(lmm)
## Linear mixed model fit by maximum likelihood  ['lmerMod']
## Formula: time6 ~ groupID * timeID + (1 | subID)
##    Data: MixedANOVA_BraverSeed
## 
##      AIC      BIC   logLik deviance df.resid 
##    840.4    861.8   -414.2    828.4      258 
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -3.0117 -0.5803  0.0202  0.5731  3.3756 
## 
## Random effects:
##  Groups   Name        Variance Std.Dev.
##  subID    (Intercept) 0.08092  0.2845  
##  Residual             1.27098  1.1274  
## Number of obs: 264, groups:  subID, 142
## 
## Fixed effects:
##                      Estimate Std. Error t value
## (Intercept)            0.3868     0.1431   2.703
## groupIDUAA            -0.4281     0.2032  -2.108
## timeIDpre             -0.4811     0.1983  -2.426
## groupIDUAA:timeIDpre   0.6324     0.2783   2.273
## 
## Correlation of Fixed Effects:
##             (Intr) grIDUAA tmIDpr
## groupIDUAA  -0.704               
## timeIDpre   -0.682  0.480        
## grpIDUAA:ID  0.486 -0.690  -0.713


#####Verdict: Use lmer, set REML=false.


  • STEP3 Checklist:
#1. Close to 0, linear model okay.
# 2.
Anova(lmm)
## Analysis of Deviance Table (Type II Wald chisquare tests)
## 
## Response: time6
##                 Chisq Df Pr(>Chisq)  
## groupID        0.5558  1    0.45596  
## timeID         1.3208  1    0.25045  
## groupID:timeID 5.1645  1    0.02305 *
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1


#####Verdict: linear-model okay, significant interactions (p=0.02). Note that groupIDUAA:timeIDpre > 0. That is, as groupID changes to AE and timeID changes to post, AE increases (symbolized as AE^ in the Summary Table Tab). Put differently, AE increases from Pre to Post, while UAA decreases from Pre to Post.


Graph_SignificantInteractions

Time3
lmm <- lmer(time3 ~ groupID * timeID + (1 | subID), data = MixedANOVA_BraverSeed,
            REML = FALSE)

MixedANOVA_BraverSeed$fit3 <- predict(lmm)
time_order <- c('pre', 'post') 

ggplot(MixedANOVA_BraverSeed, aes(x=factor(timeID, level = time_order), y=time3, group=subID, col=groupID)) + 
      scale_color_manual(values=wes_palette(n=3, name="GrandBudapest1")) + 
      facet_grid(~groupID) +
      geom_line(aes(y=time3, group=subID)) + 
      geom_line(aes(y=fit3), size=5.5) +
      geom_point(alpha = 0.3) + 
      geom_hline(yintercept=0, linetype="dashed") +
      theme_bw()

Time4
lmm <- lmer(time4 ~ groupID * timeID + (1 | subID), data = MixedANOVA_BraverSeed,
            REML = FALSE)

MixedANOVA_BraverSeed$fit4 <- predict(lmm)
time_order <- c('pre', 'post') 

ggplot(MixedANOVA_BraverSeed, aes(x=factor(timeID, level = time_order), y=time4, group=subID, col=groupID)) + 
      scale_color_manual(values=wes_palette(n=3, name="GrandBudapest1")) + 
      facet_grid(~groupID) +
      geom_line(aes(y=time4, group=subID)) + 
      geom_line(aes(y=fit4), size=5.5) +
      geom_point(alpha = 0.3) + 
      geom_hline(yintercept=0, linetype="dashed") +
      theme_bw()

Time5
lmm <- lmer(time5 ~ groupID * timeID + (1 | subID), data = MixedANOVA_BraverSeed,
            REML = FALSE)

MixedANOVA_BraverSeed$fit5 <- predict(lmm)
time_order <- c('pre', 'post') 

ggplot(MixedANOVA_BraverSeed, aes(x=factor(timeID, level = time_order), y=time5, group=subID, col=groupID)) + 
      scale_color_manual(values=wes_palette(n=3, name="GrandBudapest1")) + 
      facet_grid(~groupID) +
      geom_line(aes(y=time5, group=subID)) + 
      geom_line(aes(y=fit5), size=5.5) +
      geom_point(alpha = 0.3) + 
      geom_hline(yintercept=0, linetype="dashed") +
      theme_bw()

Time6
lmm <- lmer(time6 ~ groupID * timeID + (1 | subID), data = MixedANOVA_BraverSeed,
            REML = FALSE)

MixedANOVA_BraverSeed$fit6 <- predict(lmm)
time_order <- c('pre', 'post') 

ggplot(MixedANOVA_BraverSeed, aes(x=factor(timeID, level = time_order), y=time6, group=subID, col=groupID)) + 
      scale_color_manual(values=wes_palette(n=3, name="GrandBudapest1")) + 
      facet_grid(~groupID) +
      geom_line(aes(y=time6, group=subID)) + 
      geom_line(aes(y=fit6), size=5.5) +
      geom_point(alpha = 0.3) + 
      geom_hline(yintercept=0, linetype="dashed") +
      theme_bw()

SummaryTable_SignificantInteractions

Region Time Point Cue(C) or Probe(P) Post-Pre Trend: AE>UAA (AE^) or UAA>AE (UAA^) Expected Direction?: Y/N
BA6_subFrontal_L 3 C AE^ Y
BA6_subFrontal_L 4 C AE^ Y
BA6_subFrontal_L 5 P AE^ N
BA6_subFrontal_L 6 P AE^ N