Hypotheses

Theoretical Framework


The socio-cognitive framework for acting includes three inter-dependent components: (i) cognitive perspective-taking; (ii) emotional perspective-taking; and (iii) executive functions. From a large-scale brain network perspective, cognitive and emotional perspective-taking have been tied to the hyper- and hypo-activation of the default mode network (DMN), while executive functions have been linked to the fronto-parietal network. Although not directly implicated, it is likely that the salience network also plays an important role in the acting process, given its role in social cognition and integration of external sensory stimuli with internal states (Toyomaki & Murohashi, 2013; Menon 2015; Banducci et al., 2017).

Specific Goals of the Current Study


We conducted an 8-week acting intervention to determine whether there any differences in the resting state functional connectivity in the “active-experiencing” [active-intervention] group, versus the control group, using graph metrics. We investigated both local and global topology. For the local topology, we restricted our exploration to the regions in the three networks that are hypothesized to be related to acting – i.e., the default mode network, the salience network, and the fronto-parietal network. Within this context, our analysis was at node-level in which we compared the two groups on metrics of clustering, degree, and eigen-centrality. At this level, we also hypothesized that the frontal pole, the middle and superior temporal gyrus (associated with cognitive empathy), as well as regions in the inferior and orbital frontal cortex and the anterior cingulate gyrus (associated with emotional empathy) would be implicated.
For the global topology, we expanded our exploration to all 264 functionally-distinct regions (“nodes”) identified in the Power et al. (2012) paper. In this regard, we performed a whole-brain analysis in which we compared the two groups on modularity, associativity, efficiency, and transitivity.

Modularity Analyses


1. We ran the Gen Louvain algorithm on the weighted seed-correlation matrix for each subject (at each time point, pre- and post-intervention) for gamma parameters of 0.5, 1, 1.5, 2, 2.5, and 3. For each gamma value, the Gen Louvain algorithm generates two output parameters: Ci - which refers to community assignment for node i, and Q, which is the modularity “quality” function. Although fast and accurate, this algorithm is stochastic (the output or community assignments and related Q value can vary from run to run).
2. To account for this, we ran the algorithm 100 times. It generated 100 sets of community assignments.
3. For each node, the median Q value from these 100 runs was recorded. All such “median Q” values (across all gamma levels) were combined.

References

  1. https://stats.stackexchange.com/questions/111010/interpreting-qqplot-is-there-any-rule-of-thumb-to-decide-for-non-normality/111013#111013 determining normality from qqplots
  2. https://stats.stackexchange.com/questions/101274/how-to-interpret-a-qq-plot
  3. https://stats.stackexchange.com/questions/61055/r-how-to-interpret-the-qqplots-outlier-numbers
  4. https://ase.tufts.edu/gsc/gradresources/guidetomixedmodelsinr/mixed%20model%20guide.html
  5. https://stats.stackexchange.com/questions/226946/r-lmer-vs-glmer
  6. https://www.r-bloggers.com/linear-mixed-models-in-r/
  7. https://stats.stackexchange.com/questions/283375/violated-normality-of-residuals-assumption-in-linear-mixed-model
  8. https://stats.stackexchange.com/questions/48671/what-is-restricted-maximum-likelihood-and-when-should-it-be-used
  9. https://stats.stackexchange.com/questions/77313/why-cant-i-match-glmer-family-binomial-output-with-manual-implementation-of-g
  10. https://support.sas.com/documentation/cdl/en/statug/63347/HTML/default/viewer.htm#statug_glimmix_a0000001433.htm For line about why crossed random effects can be handled with the Laplace approximation (nAHQ=1) but not with quadrature (nAHQ > 1); see code for glmer further below
  11. http://andrewgelman.com/2006/04/10/log_transformat/
  12. https://stats.stackexchange.com/questions/47840/linear-model-with-log-transformed-response-vs-generalized-linear-model-with-log to transform or not to transform, and choosing GLMM over LMM

Library

## Loading required package: carData
## Loading required package: Matrix
## 
## Attaching package: 'lme4'
## The following object is masked from 'package:nlme':
## 
##     lmList

Dataset

modularity<-read.csv(file="//Volumes//aishwar2_ACT//ACT//restingState//methods//graphTheory//Betzel//analyses//mixedModels_datasheets//combinedGammas//mixedModel_modularityMedianQ_combinedGamma_withoutsub6087.csv",
                                head=TRUE,sep=',', 
                                na.strings = "NA",
                                stringsAsFactors=TRUE)

str(modularity)
## 'data.frame':    1356 obs. of  4 variables:
##  $ subID  : Factor w/ 113 levels "sub01","sub04",..: 12 13 14 16 18 20 24 25 26 27 ...
##  $ 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 ...
##  $ medianQ: num  0.28 0.242 0.196 0.312 0.315 ...

Assumptions (adapted from Reference 4)

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.

1. Using qqnorm and Shapiro Wilkes test to see if the fitted [linear-mixed] model has normally-distributed residuals.

fit <- lmer(medianQ ~ groupID * timeID + (1|subID), data = modularity, REML = TRUE)
# REML = TRUE following the link and advice here: https://stats.stackexchange.com/questions/48671/what-is-restricted-maximum-likelihood-and-when-should-it-be-used (Reference #8)
qqpoints<- qqnorm(resid(fit))
qqline(resid(fit))

shapiro.test(resid(fit))
## 
##  Shapiro-Wilk normality test
## 
## data:  resid(fit)
## W = 0.99668, p-value = 0.00539
#non-normal residuals

2. Generate a number of data sets that are actually normal of the same sample size as the one you have (for why sample size is important when analyzing qqplots, see here: https://stats.stackexchange.com/questions/101274/how-to-interpret-a-qq-plot) - (for example, say 24 of them).

Plot your real data among a grid of such plots (5x5 in the case of 24 random sets). If it’s not especially unusual looking (the worst looking one, say), it’s reasonably consistent with normality.

z <- resid(fit)
n = length(z)
xz = cbind(matrix(rnorm(12*n),nr=n),z,matrix(rnorm(12*n),nr=n))
colnames(xz) = c(letters[1:12],"Z",letters[13:24])

opar = par()
par(mfrow=c(5,5));
par(mar=c(0.5,0.5,0.5,0.5))
par(oma=c(1,1,1,1));

ytpos = (apply(xz,2,min)+3*apply(xz,2,max))/4
cn = colnames(xz)

for(i in 1:25) {
  qqnorm(xz[,i],axes=FALSE,ylab= colnames(xz)[i],xlab="",main="")
  qqline(xz[,i],col=2,lty=2)
  box("figure", col="darkgreen")
  text(-1.5,ytpos[i],cn[i])
}

par(opar)
# Look at plot "z".

Initially, panel “Z” looked non-normal (i.e., with the original dataset). Following step3, we first removed sub6087 [as a potential outlier/extreme variable], which seemed to make panel “Z” more normal-like.However, after also removing sub6100 [during the second iteration of removing outliers/cleaning the data], the panel looks worse. Therefore, we’ll drop sub6187 but retain sub6100, and continue analyzing the data. Although Panel Z looks more-normal, an overly-critical eye may not think it’s normal-like, so we’ll treat the data as non-normal for now.

3. Is there a focused problem/pattern to the outliers? (Example, only 1 subject is causing the problem).

fit <- lmer(medianQ ~ groupID * timeID + (1|subID), data = modularity, REML = TRUE)  
# REML = TRUE following the link and advice here: https://stats.stackexchange.com/questions/48671/what-is-restricted-maximum-likelihood-and-when-should-it-be-used (Reference #8)
plot(fit, which=2) # QQ-Plot

fit.resid <- residuals(fit) # save the residuals
sorted_resid<- sort(abs(fit.resid), decreasing=TRUE)

#sorted_resid #print only onto console; commented out for cleaner knitting into html.

#The first time this was run [using original csv file titled: mixedANOVA_modularity_medianQInteractionEffects.csv], indices 822, 309, 879 were identified as potential outliers 
# (corresponding to rows 823, 310, 880...on the excel sheet) 

# cbind(modularity$subID,modularity$timeID, modularity$medianQ)[c(822, 309, 879), ]

#The second time this was run [using csv file titled: mixedANOVA_modularity_medianQInteractionEffects_withoutsub6087.csv], indices 337,28,140 were identified as potential outliers (corresponding to rows 338, 29, 141...on the excel sheet) 

# cbind(modularity$subID,modularity$timeID, modularity$medianQ)[c(337, 28, 140), ]

Sub6087 seemed to be causing the problem, so removed them and reanalyzed the data upto this point. Although the data failed while testing for normality using the Shapiro-Wilkes test (step 1), this panel of plots suggests that “panel z” is not distinctly different from the other panels.

3. Recent statistics literature suggests that transformations to normality make interpretation of model results more difficult, and they often make mischief with the variance of the transformed variable.However, let’s see if “simple” transformations (such as log or inverse) are helpful.

(see also: http://andrewgelman.com/2006/04/10/log_transformat/)

modularity$log_medianQ <- log(modularity$medianQ)
fit_log <- lmer(log_medianQ ~ groupID * timeID + (1|subID), data = modularity, REML = TRUE)
# REML = TRUE following the link and advice here: https://stats.stackexchange.com/questions/48671/what-is-restricted-maximum-likelihood-and-when-should-it-be-used (Reference #8)
qqnorm(resid(fit_log))
qqline(resid(fit_log))

hist(resid(fit_log))

shapiro.test(resid(fit_log))
## 
##  Shapiro-Wilk normality test
## 
## data:  resid(fit_log)
## W = 0.99768, p-value = 0.04975

For the original data, log doesn’t help (makes it worse). For the second round (after removing sub6087), log seems to help (but p<=0.05; at the threshold of significance).

modularity$inverse_medianQ <- 1/(modularity$medianQ)
fit_inverse <- lmer(inverse_medianQ ~ groupID * timeID + (1|subID), data = modularity, REML = TRUE)
# REML = TRUE following the link and advice here: https://stats.stackexchange.com/questions/48671/what-is-restricted-maximum-likelihood-and-when-should-it-be-used (Reference #8)
qqnorm(resid(fit_inverse))
qqline(resid(fit_inverse))

shapiro.test(resid(fit_inverse))
## 
##  Shapiro-Wilk normality test
## 
## data:  resid(fit_inverse)
## W = 0.9828, p-value = 1.271e-11

Inverse doesn’t help (makes it worse)



** Final Verdict: Data are considered to be non-normal.**

Choosing Appropriate Models

PQL


REML and maximum likelihood methods for estimating the effect sizes in the model make assumptions of normality that don’t apply to the data, so we have to use a different method for parameter estimation.

First, we need to test whether we can use penalized quasilikelihood (PQL) or not. PQL is a flexible technique that can deal with non-normal data, unbalanced design, and crossed random effects.

summary(PQL <- glmmPQL(fixed=medianQ~groupID + timeID, random= ~1 | subID, family = gaussian, data = modularity, verbose = FALSE))
## Linear mixed-effects model fit by maximum likelihood
##  Data: modularity 
##   AIC BIC logLik
##    NA  NA     NA
## 
## Random effects:
##  Formula: ~1 | subID
##         (Intercept)   Residual
## StdDev:  0.03481867 0.02652234
## 
## Variance function:
##  Structure: fixed weights
##  Formula: ~invwt 
## Fixed effects: medianQ ~ groupID + timeID 
##                   Value   Std.Error   DF  t-value p-value
## (Intercept)  0.28193163 0.004823477 1242 58.44987  0.0000
## groupIDUAA   0.01020235 0.006715132  111  1.51931  0.1315
## timeIDpre   -0.01130039 0.001442092 1242 -7.83611  0.0000
##  Correlation: 
##            (Intr) gIDUAA
## groupIDUAA -0.702       
## timeIDpre  -0.149  0.000
## 
## Standardized Within-Group Residuals:
##         Min          Q1         Med          Q3         Max 
## -2.70718012 -0.71067022  0.02366616  0.67111409  3.00952140 
## 
## Number of Observations: 1356
## Number of Groups: 113

The results suggest that modularity increases from pre to post (significantly)

However, PQL should not be used when the mean of the response variable is less than 5, which is the case for the medianQ.

GLMER


Glmer allows fitting a generalized linear mixed-effects model: these models include a link function that allows to predict response variables with non-Gaussian distributions (the error distribution is [also?] non-Gaussian). We use the Laplace approximation, a special case of a parameter estimation method called Gauss-Hermite quadrature (GHQ), with one iteration. GHQ is more accurate than Laplace due to repeated iterations, but becomes less flexible after the first iteration, so it can be used only for one random effect.
Sidenotes: 1. Crossed random effects, models without subjects, or models with non-nested subjects can be handled with the Laplace approximation but not with quadrature (reference 10; https://support.sas.com/documentation/cdl/en/statug/63347/HTML/default/viewer.htm#statug_glimmix_a0000001433.htm).
See also: https://stats.stackexchange.com/questions/77313/why-cant-i-match-glmer-family-binomial-output-with-manual-implementation-of-g 2. We consider a log link because a log-transformation of the response helped improve the model fit previously (log-transformations also transform the errors; see here: https://stats.stackexchange.com/questions/47840/linear-model-with-log-transformed-response-vs-generalized-linear-model-with-log)

summary(GHQ <- glmer(medianQ ~ groupID + timeID + groupID * timeID + (1 | subID), data = modularity,family = gaussian(link="log"), nAGQ = 1)) 
## Generalized linear mixed model fit by maximum likelihood (Laplace
##   Approximation) [glmerMod]
##  Family: gaussian  ( log )
## Formula: medianQ ~ groupID + timeID + groupID * timeID + (1 | subID)
##    Data: modularity
## 
##      AIC      BIC   logLik deviance df.resid 
##  -5814.8  -5783.5   2913.4  -5826.8     1350 
## 
## Scaled residuals: 
##      Min       1Q   Median       3Q      Max 
## -2.60050 -0.61070  0.06136  0.65096  2.74807 
## 
## Random effects:
##  Groups   Name        Variance  Std.Dev.
##  subID    (Intercept) 0.0053498 0.07314 
##  Residual             0.0008404 0.02899 
## Number of obs: 1356, groups:  subID, 113
## 
## Fixed effects:
##                       Estimate Std. Error t value Pr(>|z|)    
## (Intercept)          -1.289360   0.024985 -51.606  < 2e-16 ***
## groupIDUAA            0.015068   0.035498   0.424    0.671    
## timeIDpre            -0.054481   0.007243  -7.522 5.41e-14 ***
## groupIDUAA:timeIDpre  0.042437   0.010014   4.238 2.26e-05 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Correlation of Fixed Effects:
##             (Intr) grIDUAA tmIDpr
## groupIDUAA  -0.700               
## timeIDpre   -0.140  0.098        
## grpIDUAA:ID  0.101 -0.139  -0.723
#getME(GHQ,'b')

Interactions are significant. Note that groupIDUAA:timeIDpre > 0. That is, as groupID changes to AE and timeID changes to post, AE increases.

The Linear Model for “Simple” Transformations


Recognizing that “the mean of the log-transformed variable is not always the same as the log of the variable’s mean”, we applied a glmer model with family= gaussian(link=‘log’) first (see previous tab). However,we can also explore the linear model when the response variable follows a log-distribution.

count(modularity, 'groupID')# Check if data are balanced across groupID
##   groupID freq
## 1      AE  672
## 2     UAA  684
count(modularity, 'timeID') # Check if data are balanced across groupID
##   timeID freq
## 1   post  678
## 2    pre  678
lmm <- lmer(log_medianQ ~ groupID * timeID + (1 | subID), data = modularity,
            REML = FALSE) #data not balanced across groupID, so set REML to false
summary(lmm)
## Linear mixed model fit by maximum likelihood  ['lmerMod']
## Formula: log_medianQ ~ groupID * timeID + (1 | subID)
##    Data: modularity
## 
##      AIC      BIC   logLik deviance df.resid 
##  -2018.4  -1987.1   1015.2  -2030.4     1350 
## 
## Scaled residuals: 
##      Min       1Q   Median       3Q      Max 
## -3.01637 -0.65999  0.02129  0.66760  2.83950 
## 
## Random effects:
##  Groups   Name        Variance Std.Dev.
##  subID    (Intercept) 0.01727  0.1314  
##  Residual             0.01015  0.1007  
## Number of obs: 1356, groups:  subID, 113
## 
## Fixed effects:
##                       Estimate Std. Error t value
## (Intercept)          -1.265650   0.018400 -68.784
## groupIDUAA            0.018863   0.025908   0.728
## timeIDpre            -0.069243   0.007772  -8.910
## groupIDUAA:timeIDpre  0.037056   0.010943   3.386
## 
## Correlation of Fixed Effects:
##             (Intr) grIDUAA tmIDpr
## groupIDUAA  -0.710               
## timeIDpre   -0.211  0.150        
## grpIDUAA:ID  0.150 -0.211  -0.710
Anova(lmm) #No significant interaction effects
## Analysis of Deviance Table (Type II Wald chisquare tests)
## 
## Response: log_medianQ
##                  Chisq Df Pr(>Chisq)    
## groupID         2.1802  1  0.1397996    
## timeID         85.3712  1  < 2.2e-16 ***
## groupID:timeID 11.4678  1  0.0007081 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Interactions are significant. Note that groupIDUAA:timeIDpre > 0. That is, as groupID changes to AE and timeID changes to post, AE increases.

Non-Linear Mixed Effects Model


Linear refers to the relationship between the parameters that you are estimating (e.g., β) and the outcome (e.g., yi). Hence, y=e^xβ+ϵ is linear, but y=e^βx+ϵ is not. A linear model means that your estimate of your parameter vector can be written β̂ =∑iwiyi, where the {wi} are weights determined by your estimation procedure. Linear models can be solved algebraically in closed form, while many non-linear models need to be solved by numerical maximization using a computer. Let us consider the possibility that β and y have a non-linear relationship. Then:

non_lmm <- lme(fixed=medianQ ~ groupID * timeID, random= ~1|subID, data = modularity)
summary(non_lmm)
## Linear mixed-effects model fit by REML
##  Data: modularity 
##         AIC       BIC   logLik
##   -5611.392 -5580.136 2811.696
## 
## Random effects:
##  Formula: ~1 | subID
##         (Intercept)   Residual
## StdDev:  0.03515406 0.02639354
## 
## Fixed effects: medianQ ~ groupID * timeID 
##                            Value   Std.Error   DF  t-value p-value
## (Intercept)           0.28465282 0.004913376 1241 57.93426  0.0000
## groupIDUAA            0.00480771 0.006918020  111  0.69495  0.4885
## timeIDpre            -0.01674277 0.002036306 1241 -8.22213  0.0000
## groupIDUAA:timeIDpre  0.01078927 0.002867113 1241  3.76311  0.0002
##  Correlation: 
##                      (Intr) grIDUAA tmIDpr
## groupIDUAA           -0.710               
## timeIDpre            -0.207  0.147        
## groupIDUAA:timeIDpre  0.147 -0.207  -0.710
## 
## Standardized Within-Group Residuals:
##        Min         Q1        Med         Q3        Max 
## -2.6155667 -0.7203260  0.0315781  0.7048443  2.9214832 
## 
## Number of Observations: 1356
## Number of Groups: 113

Interactions are significant. Note that groupIDUAA:timeIDpre > 0. That is, as groupID changes to AE and timeID changes to post, AE increases.

However, such models are more difficult to interpret and there is no rationale for assuming that the relationship between y and β should be non-linear.

Final Model Consideration

To determined whether to consider the GLMM over the LMM or vice versa (the non-linear model is out of the running since it is more difficult to interpret and justify why the relationship between y and β should be non-linear), we evaluate the following statements:


*It boils down to whether or not you need to address linearity and heteroscedasticity (unequal variances) in your data, or just linearity. Transforming the data affects both the linearity and variance assumptions of a model. For example, if your residuals exhibit issues with both, you could consider transforming the data, which potentially could fix both. The transformation transforms the errors and thus their variance.

In contrast, using the link function only affects the linearity assumption, not the variance. The log is taken of the mean (expected value), and thus the variance of the residuals is not affected.

In summary, if you don’t have an issue with non-constant variance, she suggests using the link function over transformation, because you don’t want to change your variance in that case (you’re already meeting the assumption).* [Reference 12: https://stats.stackexchange.com/questions/47840/linear-model-with-log-transformed-response-vs-generalized-linear-model-with-log]


TESTING ASSUMPTION OF HOMOSCEDASTICITY USING LEVENE’S TEST

# Method 1:
# ---------
tapply(modularity$medianQ, list(modularity$groupID,modularity$timeID),sd)
##           post        pre
## AE  0.03760133 0.04806290
## UAA 0.03512723 0.05193188
# RESULT:
#          post        pre
#AE  0.03760133 0.04806290
#UAA 0.03512723 0.05193188
#Variances don't seem to be significantly different from each other.

# Method 2:
# ---------
preQ <- data.frame(subset(modularity, timeID == 'pre', select = c(groupID,medianQ)))
with(preQ, leveneTest(medianQ, as.factor(groupID), center="mean"))
## Levene's Test for Homogeneity of Variance (center = "mean")
##        Df F value Pr(>F)
## group   1  0.6528 0.4194
##       676
#RESULT: p-value = 0.4194

postQ <- data.frame(subset(modularity, timeID == 'post', select = c(groupID,medianQ)))
with(postQ, leveneTest(medianQ, as.factor(groupID), center="mean"))
## Levene's Test for Homogeneity of Variance (center = "mean")
##        Df F value Pr(>F)
## group   1  2.2812 0.1314
##       676
#RESULT: p-value = 0.1314
Homoscedasticity assumption is met. So let’s just go with GLMM.

GRAPH_SIGNIFICANT_INTERACTIONS

MedianQ_GLMM model


** Note that although medianQ is not negative, we are using the link=log function, rendering the fit as negative. The Y axis should probably read as “medianQ (link=log)”.**

GHQ <- glmer(medianQ ~ groupID + timeID + groupID * timeID + (1 | subID), data = modularity,family = gaussian(link="log"), nAGQ = 1) 

modularity$fit <- predict(GHQ)
time_order <- c('pre', 'post') 

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