CCA Sensitivity

Date written: 2020-07-22
Last ran: 2020-07-25

---

Overview

Here, we perform canonical correlation analysis (CCA) as reported in https://rpubs.com/navona/CCA, but without the n=43 participants who, in the outlier analysis, demonstrated multivariate non-normal data. The data have undergone Yeo-Johnson transformation, as in the full-sample analysis.

Thus, the present includes a total of n=131 participants, comprised of n=94 SSD and 37 HC; n=91 men and n=40 women; and 67 from the SPINS study and 64 from the PNS study. As before, we include 35 variables in the CCA analysis. Our \(X\) set is comprised of FA estimates in 23 white matter tracts, representing association, projection, and commissural fibers. Our \(Y\) set is comprised of 12 neurocognition and social cognition variables.

---

Data structure

[1] Review correlation matrices. Below, we summarise correlations within the X and Y sets, and the cross correlation matrix. We see some differences with the full-sample CCA. In the present analysis, our X set displays some negative correlations, of which there were none in the full sample (min >0), tohugh the median correlation value is comparable. The Y set in the present analysis has a higher median correlation value, though similar patterns of correlation. The cross-correlation matrix in the present analysis appears similar, but has a slightly wider span of values in both directions.

#run correlations (note: identical to pairwise complete, `cor(X, use='p')`)
corr <- matcor(X, Y)

#also run correlations for separate HC and SSD
corr_HC <- matcor(X_HC, Y_HC)
corr_SSD <- matcor(X_SSD, Y_SSD)

#function to describe correlations -- excluding the 1s from the diagonal
tab_corr_fn <- function(corr){
  data.frame(as.list(summary(c(corr[upper.tri(corr)], corr [lower.tri(corr)])))) 
}

#function to plot correlation matrices
plotCorr_fn <- function(df){
  ggcorrplot(df) + 
  theme_classic(base_size = 6) +
  theme(legend.position = 'top',
        axis.text.x = element_text(angle = 90, vjust=.5, hjust=1),
        axis.title= element_blank())
}

X correlation matrix

Min.	X1st.Qu.	Median	Mean	X3rd.Qu.	Max.
-0.173	0.225	0.332	0.33	0.436	0.816

Y correlation matrix

Min.	X1st.Qu.	Median	Mean	X3rd.Qu.	Max.
0.111	0.317	0.393	0.41	0.503	0.698

XY correlation matrix

Min.	X1st.Qu.	Median	Mean	X3rd.Qu.	Max.
-0.241	-0.003	0.149	0.181	0.369	0.816

---

[2] Test of multicollinearity. We attempt to identify collinearity via the variance inflation factors (VIF) metric, which is an indicator of how much the variance of the coefficient is inflated compared to the non-significant correlations with any other variables in the model. Typically, VIF = 1 indicates no correlation; > 1 < 5 indicates moderate correlation, >5 high correlation. The full-sample CCA analysis showed no evidence of multicollinearity in our variable sets. The present analysis shows higher VIF values in both the X and Y sets, though all fall short of the 'full of thumb' threshold of 10 indicating problematic structure. The highest value, of a VIF of >5, belong to CC4.

X set

usdm::vif(X)

##          Variables      VIF
## 1          AF_left 3.967281
## 2          CB_left 4.363373
## 3         ILF_left 3.081015
## 4        IOFF_left 3.902069
## 5          UF_left 2.396230
## 6         AF_right 1.678029
## 7         CB_right 4.351766
## 8        ILF_right 3.909041
## 9       IOFF_right 4.137070
## 10        UF_right 3.669493
## 11       CR_F_left 2.596059
## 12       CR_P_left 2.389890
## 13         TF_left 2.465649
## 14      CR_F_right 2.727513
## 15      CR_P_right 1.707062
## 16        TF_right 2.719120
## 17 CC1_commissural 2.297744
## 18 CC2_commissural 3.439523
## 19 CC3_commissural 3.014969
## 20 CC4_commissural 5.819784
## 21 CC5_commissural 4.045530
## 22 CC6_commissural 2.363416
## 23 CC7_commissural 1.682168

Y set

usdm::vif(Y)

##              Variables      VIF
## 1  Attention_vigilance 2.044649
## 2       Working_memory 2.547930
## 3      Verbal_learning 2.335476
## 4     Processing_speed 2.749012
## 5      Visual_learning 1.549828
## 6      Problem_solving 1.421873
## 7              TASIT_1 1.578635
## 8              TASIT_2 2.566499
## 9              TASIT_3 2.795722
## 10                RMET 1.926499
## 11                 RAD 2.136974
## 12               ER_40 1.408829

---

[3] Test of homogeneity between groups (SSD and HC). Covariance matrices tell us how much two random variables vary together (as opposed to correlation, which tells how a change in one variable is correspondance to a change in another). In the present analysis, as well as in the full-sample CCA analysis, HC and SSD populations appear to be similar:

#make covariance matrices
cov_HC <- cov(df_HC)
cov_SSD <- cov(df_SSD)

#plot
cov_HC_plot <- plotCorr_fn(cov_HC)
cov_SSD_plot <-plotCorr_fn(cov_SSD)

#arrange side by side
grid.arrange(cov_HC_plot, cov_SSD_plot, ncol=2)

We can use Box's M-test to test the null hypothesis that the covariance matrices are equal; i.e., that the observed covariance matrices of the variable sets are equal across groups. As in the full-sample analysis, the present analysis finds that the M test is significant for X (brain) variables, but not Y (behaviour) variables, meaning that patients and controls do not have equal covariance of brain variables.

X set

rstatix::box_m(X, df_demo$group)

## # A tibble: 1 x 4
##   statistic  p.value parameter method                                           
##       <dbl>    <dbl>     <dbl> <chr>                                            
## 1      360. 0.000521       276 Box's M-test for Homogeneity of Covariance Matri…

Y set

rstatix::box_m(Y, df_demo$group)

## # A tibble: 1 x 4
##   statistic p.value parameter method                                            
##       <dbl>   <dbl>     <dbl> <chr>                                             
## 1      92.0   0.133        78 Box's M-test for Homogeneity of Covariance Matric…

---

CCA

[1] Canonical correlations. We have quite high canonical correlations. As seen in the plot below, the first variate has a correlation values of 0.6111328. This is slightly higher than the full-sample analysis, in which the first canonical correlation was 0.584854.

#manipulate canonical correlation values
correlations$variate <- seq(1:ncol(Y)) #
names(correlations)[names(correlations) == 'unlist(correlations)'] <- 'canonical correlation'
correlations$`canonical correlation` <- round(correlations$`canonical correlation`, 3)

#bar plot/scree plot
ggplot(correlations, aes(x=variate, y=`canonical correlation`)) + 
  geom_bar(stat='identity') +
  scale_x_continuous(breaks=seq(1, ncol(Y), 1)) +
  theme_minimal() +
  theme(axis.text.x=element_blank()) +
  geom_text(aes(label=`canonical correlation`), hjust = 1.2, color= 'white', size=3.5) +
  coord_flip()

---

[2] Visualization. The structure of the canonical correlations differs slightly between the full-sample and present analysis. For the most part, it appears that variables are projected in the same direction from the origin.

#visualize 'target' plot and 'user' plot
plt.cc(cca_cca, var.label = T)

---

Significance testing

#the null hypothesis is that the canonical correlations in all rows are 0

#define n, p, q
n <- nrow(df)
p <- ncol(X)
q <- ncol(Y)

#function to calculate p-values using the F-approximations of different statistical tests
p_fn <- function(tstat){
  df_p <- as.data.frame(p.asym(cca_cca$cor, n, p, q, tstat = tstat)) #test
  df_p <- setDT(df_p, keep.rownames = 'variate')[] #add row names
  df_p$variate <- paste(df_p$variate, '-', ncol(Y), sep='')
  return(df_p)
}

[1] Statistical testing, F distribution. Next, we formally test for significance. In the full-sample analysis, only the final test, Roy's Largest Root, was highly significant, p<0.001. That is again the case here.

Wilks' Lambda

p_wilks <- p_fn('Wilks')

## Wilks' Lambda, using F-approximation (Rao's F):
##                  stat    approx df1       df2   p.value
## 1 to 12:   0.08678193 0.9869421 276 1058.9898 0.5471999
## 2 to 12:   0.13851494 0.9029818 242  986.4784 0.8351084
## 3 to 12:   0.20964184 0.8162318 210  911.8632 0.9650432
## 4 to 12:   0.29945230 0.7300611 180  834.9750 0.9951887
## 5 to 12:   0.41160153 0.6321096 152  755.6142 0.9997031
## 6 to 12:   0.50117240 0.5929278 126  673.5462 0.9998074
## 7 to 12:   0.60059541 0.5397545 102  588.4982 0.9999078
## 8 to 12:   0.71342048 0.4543092  80  500.1601 0.9999844
## 9 to 12:   0.79225432 0.4182122  60  408.1917 0.9999607
## 10 to 12:  0.87140691 0.3530896  42  312.2457 0.9999441
## 11 to 12:  0.93020071 0.3003862  26  212.0000 0.9996933
## 12 to 12:  0.97009285 0.2748934  12  107.0000 0.9920141

Hotelling-Lawley Trace

p_hotelling <- p_fn('Hotelling')

##  Hotelling-Lawley Trace, using F-approximation:
##                  stat    approx df1  df2   p.value
## 1 to 12:   2.86798244 0.9785085 276 1130 0.5827629
## 2 to 12:   2.27185587 0.9027967 242 1154 0.8386828
## 3 to 12:   1.75835964 0.8219634 210 1178 0.9627329
## 4 to 12:   1.32996015 0.7400982 180 1202 0.9943996
## 5 to 12:   0.95544565 0.6422020 152 1226 0.9996694
## 6 to 12:   0.73783017 0.6099786 126 1250 0.9997307
## 7 to 12:   0.53944931 0.5614856 102 1274 0.9998506
## 8 to 12:   0.35159393 0.4753843  80 1298 0.9999750
## 9 to 12:   0.24109271 0.4426730  60 1322 0.9999326
## 10 to 12:  0.14118465 0.3770527  42 1346 0.9999079
## 11 to 12:  0.07371468 0.3236830  26 1370 0.9995386
## 12 to 12:  0.03082916 0.2984434  12 1394 0.9897574

Pillai-Bartlett Trace

p_pillai <- p_fn('Pillai')

##  Pillai-Bartlett Trace, using F-approximation:
##                  stat    approx df1  df2   p.value
## 1 to 12:   2.11213110 0.9937431 276 1284 0.5185859
## 2 to 12:   1.73864783 0.9157974 242 1308 0.8042200
## 3 to 12:   1.39936966 0.8373089 210 1332 0.9483061
## 4 to 12:   1.09945392 0.7598292 180 1356 0.9902043
## 5 to 12:   0.82698353 0.6719886 152 1380 0.9989857
## 6 to 12:   0.64826087 0.6363323 126 1404 0.9992922
## 7 to 12:   0.48272012 0.5867776 102 1428 0.9996199
## 8 to 12:   0.32457345 0.5045647  80 1452 0.9999167
## 9 to 12:   0.22506773 0.4702079  60 1476 0.9998185
## 10 to 12:  0.13423463 0.4040274  42 1500 0.9997730
## 11 to 12:  0.07102913 0.3490158  26 1524 0.9990952
## 12 to 12:  0.02990715 0.3223051  12 1548 0.9855828

Roy's Largest Root

plt.asym(p.asym(rho, n, p, q, tstat = "Roy"), rhostart=1)

##  Roy's Largest Root, using F-approximation:
##               stat   approx df1 df2          p.value
## 1 to 1:  0.3734833 5.861911  12 118 0.00000006694766
## 
##  F statistic for Roy's Greatest Root is an upper bound.

---

[2] Permuted statistical test. The full-sample analysis found that Roy's was not significant under permutation, p=.2162162. We draw the same conclusion here, and the p value is now further from significance.

#for some reason, this doesn't work in rstudio (must be packages loaded - run in R instead - takes a bit of time)
#CCP::p.perm(X, Y, nboot=500, rhostart=1, type='Roy') 

---

Model exploration

#calculate correlations
rxx <- cor(X, use = "p") #pairwise complete observations
ryy <- cor(Y, use = "p")
rxy <- cor(X,Y, use = "p")

#calculate omega
omega = t(solve(chol(rxx))) %*% rxy %*% solve(chol(ryy)) 

#decompose matrix
usv <- svd(omega)

#standardized weights (beta)
x.weights.std <- solve(chol(rxx)) %*% usv$u
y.weights.std <- solve(chol(ryy)) %*% usv$v

#structure coefficients
x.structures <- rxx %*% x.weights.std
y.structures <- ryy %*% y.weights.std

#canonical cross-loadings/redundancy
x.scores <- as.matrix(X) %*% x.weights.std
y.scores <- as.matrix(Y) %*% y.weights.std

#redundancy (?)
red.x.yscores <- cor(X, y.scores) 
red.y.xscores <- cor(Y, x.scores)

[1] Participant scores (CV1). Here, we show the association between scores in the X set and Y set, for each participant, for the first canonical variate, Rc=0.611.

#bind together data
scores.plot <- as.data.frame(cbind(x.scores[, 1], y.scores[, 1], df_demo$group))

#make sure numeric
scores.plot[1:2] <- sapply(scores.plot[1:2],as.numeric)

#change names for clarity
names(scores.plot) <- c('Xset', 'Yset', 'group')

#make a scatter plot
ggplot(scores.plot, aes(x=Xset, y=Yset)) +
  geom_point(size = 4, aes(shape=group), alpha=.8) +
  geom_smooth(method=lm, color="black", alpha=.4) +
  theme_bw() +
  theme(axis.ticks=element_blank(), 
        panel.grid.minor = element_blank(), 
        panel.grid.major.x = element_blank(),
        legend.position='top')

---

[2] Descriptive table (CV1). We see a large difference in the structure coefficients between the full-sample and present analyis. Before, only one variable contributed past threshold to the Y set. Now, we see several contributions past threshold.

#for table, take just first column
x.weights.std <- x.weights.std[, 1]
x.weights.std <- t(data.frame(lapply(x.weights.std, type.convert), stringsAsFactors = F)) #get into df

y.weights.std <- y.weights.std[, 1]
y.weights.std <- t(data.frame(lapply(y.weights.std, type.convert), stringsAsFactors = F)) #get into df

#bind together standardized weights
weights <- rbind(x.weights.std, y.weights.std)

#bind together structure coefficients
structures <- rbind(x.structures, y.structures)

#make table
tbl <- as.data.frame(cbind(weights, structures))

#add rownames as a variable
tbl <- tibble::rownames_to_column(tbl, "variable")

#redefine row names with full names
tbl$variable <- c(names(X), names(Y))

#set up dynamic names for table
label <- paste("$R_c$=", round(cca_cca$cor[1], 5))
myHeader <- c(" " = 1, label = 2) # create vector with colspan
names(myHeader) <- c(" ", label) #set names

#make into pretty table
tbl_plot <- tbl

#take only first 3 columns -- adjust if want more variates
tbl_plot <- tbl_plot[, 1:3]

#add a column for the communality coefficient
tbl_plot$communalities <- tbl_plot$V2 ^ 2

#round
tbl_plot[,2:4] <- round(tbl_plot[,2:4], 3) #shorten here, because formatting makes characters (?!)

tbl_plot %>%
  mutate(V2 = ifelse((V2 > .45 | V2 < -.45),
          cell_spec(V2, color='black', underline = T, bold = T), #no green
          cell_spec(V2, color='black'))) %>%
  kable(escape=F,
        align = c('l', 'c', 'c', 'c'),
        col.names = c('_Variable_', '$\\beta$', '$r_s$', '$r_s^2$ ($h^2$)')) %>%
kable_styling(bootstrap_options=c('striped', 'hover', 'condensed')) %>%
add_header_above(c(" " = 1, "_Function 1_" = 3)) %>%
pack_rows('X', 1, 23) %>%
pack_rows('Y', 24, 35) %>%
footnote(general = "$r_s$ values > .450 are emphasized, following convention;
         $\\beta$ = standardized canonical function coefficient;
         $r_s$ = structure coefficient;
         $r_s^2$ = squared structure coefficient, here also communality $h^2$")

	Function 1
Variable	\(\beta\)	\(r_s\)	\(r_s^2\) (\(h^2\))
X
AF_left	0.027	-0.272	0.074
CB_left	-0.073	-0.011	0.000
ILF_left	0.196	0.069	0.005
IOFF_left	-0.724	-0.495	0.245
UF_left	0.361	0.008	0.000
AF_right	0.159	-0.025	0.001
CB_right	0.207	-0.172	0.029
ILF_right	-0.085	-0.068	0.005
IOFF_right	0.151	-0.456	0.208
UF_right	-0.253	-0.338	0.115
CR_F_left	0.631	-0.063	0.004
CR_P_left	0.090	0.046	0.002
TF_left	-0.517	-0.477	0.228
CR_F_right	-0.678	-0.362	0.131
CR_P_right	0.035	-0.173	0.030
TF_right	-0.085	-0.378	0.143
CC1_commissural	0.001	-0.156	0.024
CC2_commissural	0.229	-0.1	0.010
CC3_commissural	0.300	0.098	0.010
CC4_commissural	0.167	0.315	0.099
CC5_commissural	0.126	0.311	0.097
CC6_commissural	-0.251	0.038	0.001
CC7_commissural	-0.283	-0.276	0.076
Y
Attention_vigilance	0.028	0.256	0.065
Working_memory	0.316	0.481	0.232
Verbal_learning	0.297	0.488	0.238
Processing_speed	-0.178	0.344	0.118
Visual_learning	-0.511	-0.03	0.001
Problem_solving	-0.095	0.099	0.010
TASIT_1	0.184	0.557	0.310
TASIT_2	0.495	0.645	0.417
TASIT_3	-0.617	0.374	0.140
RMET	0.610	0.72	0.519
RAD	0.173	0.519	0.269
ER_40	0.108	0.279	0.078
Note:
\(r_s\) values > .450 are emphasized, following convention; \(\beta\) = standardized canonical function coefficient; \(r_s\) = structure coefficient; \(r_s^2\) = squared structure coefficient, here also communality \(h^2\)

---

[3] Structure coefficients (CV1). This is the same data as in the table above, but perhaps more clearly presented in this visualization.

#make a new dataframe
tbl_sc <- tbl_plot

#keep only the structure coefficients variable
tbl_sc <- subset(tbl_sc, select = c(variable, V2))

#make a logical variable for colour of bars
tbl_sc$colour <- as.factor(ifelse(abs(tbl_sc$V2) >= .45, 1, 0))

#Lock in factor-level order
tbl_sc$variable <- factor(tbl_sc$variable, levels=tbl_sc$variable)

#make plot
ggplot(tbl_sc, aes(variable, V2)) + 
  geom_bar(stat='identity', aes(fill = colour)) + 
  theme_bw() +
  theme(axis.ticks=element_blank(), 
        panel.grid.minor = element_blank(), 
        panel.grid.major.x = element_blank(),
        axis.text.x = element_text(angle = 90, hjust = 1, size=5),
        legend.position = "none") + 
  xlab('') + 
  ylab('structure coefficient') + 
  geom_hline(yintercept = .45, linetype="dashed", color = "darkgreen", size=1.5) +
  #annotate(geom='text', x=.87, y = 10, color='green', label = "rs = abs(.45)") +
  geom_hline(yintercept = -.45, linetype="dashed", color = "darkgreen", size=1.5) +
  scale_fill_manual(values=c('grey', 'darkgreen'))

---

[4] Model redundancy. Here, we calculate the Stewart-Love redundancy index for each variate (i.e., variance in one set that can be explained by the other set, and vice versa). In the full-sample analysis, we saw that the X set explains ~15% of the variance in Y; here, we have a value of 23%. I believe both values are quite high for similar brain-behaviour relationships.

#proportion of variance in X explained by Y, and vice versa
redundancy(cca_candisc)

## 
## Redundancies for the X variables & total X canonical redundancy
## 
##     Xcan1     Xcan2     Xcan3     Xcan4     Xcan5     Xcan6     Xcan7     Xcan8 
##  0.024941  0.011687  0.009288  0.005992  0.008972  0.006271  0.007492  0.004013 
##     Xcan9    Xcan10    Xcan11    Xcan12 total X|Y 
##  0.005938  0.002187  0.001333  0.001098  0.089210 
## 
## Redundancies for the Y variables & total Y canonical redundancy
## 
##     Ycan1     Ycan2     Ycan3     Ycan4     Ycan5     Ycan6     Ycan7     Ycan8 
##  0.074593  0.036846  0.032499  0.045861  0.008517  0.010943  0.008046  0.003407 
##     Ycan9    Ycan10    Ycan11    Ycan12 total Y|X 
##  0.004034  0.003943  0.002251  0.001633  0.232572

---

[5] Model stability / validation. Lastly, we validated our model via iterative feature removal. We see stability, as was the case in the full-sample analysis.

#run the original correlation 
original <- cc(X, Y)$cor

#write a nested function to leave out a column, and calculate correlation
cc_df_one_out <- function(X, Y){ #takes X and Y dataframe or matrix
  f <- function(x) combn(ncol(x), ncol(x) - 1) #base combn function leave out a column ...
  
  X_inx <- f(X) #from X
  Y_inx <- f(Y) #from Y
  
  ccX <- t(apply(X_inx, 2, function(i) cc(X[, i], Y)$cor)) #run the correlation between dropped X and Y
  ccY <- t(apply(Y_inx, 2, function(i) cc(X, Y[, i])$cor)) #run the correlation between dropped Y and X
  
  list(XVALS = as.data.frame(ccX), YVALS = as.data.frame(ccY)) #add the values to a dataframe
}

#run the function
values <- cc_df_one_out(X, Y)

#extract from list and combine the dataframes 
values <- data.frame(Reduce(rbind.fill, values))

#add the original observations as the last row
values <- rbind(values, original)

#add a column indicating that the last row is original -- used to colour the plot
values$original <- c(rep(FALSE, ncol(X)), rep(FALSE, ncol(Y)), TRUE)

#melt for ggplot
values <- reshape2::melt(values)

#boxplot
ggplot(values, aes(x=variable, y=value)) +
  geom_boxplot(outlier.shape = NA) +
  geom_violin(fill='grey', alpha=.5) +
  geom_jitter(aes(colour = original, size= original)) +
  theme_bw() +
  xlab('') +
  theme(axis.ticks = element_blank()) +
  scale_colour_manual(values = c("FALSE" = "black", "TRUE" = "red")) +
  scale_size_manual(values = c("FALSE" = .5, "TRUE" = 5)) +
  guides(color = FALSE, size = FALSE)

---

SSD subset

We see very similar results when evaluating just the SSD subset.

Correlations

corr_SSD <- matcor(X_SSD, Y_SSD)

plotX <- plotCorr_fn(corr_SSD$Xcor)
plotY <- plotCorr_fn(corr_SSD$Ycor)
plotXY <-  plotCorr_fn(corr_SSD$XYcor)

#arrange side by side
grid.arrange(plotX, plotY, plotXY, ncol=3)

Canonical correlations

Significance testing (F)

n <- nrow(df_SSD)
p_wilks <- p_fn('Wilks')

## Wilks' Lambda, using F-approximation (Rao's F):
##                  stat    approx df1      df2   p.value
## 1 to 12:   0.02512490 0.9909368 276 663.8860 0.5302275
## 2 to 12:   0.04787518 0.9242267 242 620.9453 0.7626370
## 3 to 12:   0.08256254 0.8685558 210 576.2655 0.8856869
## 4 to 12:   0.13779184 0.7992125 180 529.7333 0.9625367
## 5 to 12:   0.20395379 0.7569308 152 481.2149 0.9793041
## 6 to 12:   0.30067516 0.6861351 126 430.5532 0.9939427
## 7 to 12:   0.40726469 0.6317330 102 377.5658 0.9969913
## 8 to 12:   0.53143672 0.5649229  80 322.0447 0.9986770
## 9 to 12:   0.65764659 0.4982351  60 263.7623 0.9991708
## 10 to 12:  0.76887282 0.4466424  42 202.4860 0.9986954
## 11 to 12:  0.87989824 0.3506539  26 138.0000 0.9986038
## 12 to 12:  0.95032442 0.3049213  12  70.0000 0.9865270

p_hotelling <- p_fn('Hotelling')

##  Hotelling-Lawley Trace, using F-approximation:
##                  stat    approx df1 df2   p.value
## 1 to 12:   4.59607462 0.9519647 276 686 0.6815410
## 2 to 12:   3.69058755 0.9023131 242 710 0.8290320
## 3 to 12:   2.96605005 0.8639209 210 734 0.9001368
## 4 to 12:   2.29711096 0.8061158 180 758 0.9617974
## 5 to 12:   1.81695228 0.7789784 152 782 0.9717485
## 6 to 12:   1.34272052 0.7157624 126 806 0.9902222
## 7 to 12:   0.98821989 0.6701164 102 830 0.9941388
## 8 to 12:   0.68332720 0.6078765  80 854 0.9971500
## 9 to 12:   0.44583914 0.5436761  60 878 0.9981581
## 10 to 12:  0.27671148 0.4952257  42 902 0.9971872
## 11 to 12:  0.13231124 0.3926930  26 926 0.9974399
## 12 to 12:  0.05227223 0.3448515  12 950 0.9805562

p_pillai <- p_fn('Pillai')

##  Pillai-Bartlett Trace, using F-approximation:
##                  stat    approx df1  df2   p.value
## 1 to 12:   3.02413824 1.0254056 276  840 0.3925676
## 2 to 12:   2.54893844 0.9628910 242  864 0.6356511
## 3 to 12:   2.12880406 0.9119260 210  888 0.7935194
## 4 to 12:   1.72798710 0.8523290 180  912 0.9088947
## 5 to 12:   1.40359034 0.8156689 152  936 0.9425690
## 6 to 12:   1.08190974 0.7549967 126  960 0.9765632
## 7 to 12:   0.82018920 0.7077413 102  984 0.9861475
## 8 to 12:   0.58653574 0.6475116  80 1008 0.9926666
## 9 to 12:   0.39462435 0.5848616  60 1032 0.9951193
## 10 to 12:  0.24996294 0.5348734  42 1056 0.9936210
## 11 to 12:  0.12378310 0.4329459  26 1080 0.9943903
## 12 to 12:  0.04967558 0.3824292  12 1104 0.9700149

plt.asym(p.asym(rho, n, p, q, tstat = "Roy"), rhostart=1)

##  Roy's Largest Root, using F-approximation:
##               stat   approx df1 df2       p.value
## 1 to 1:  0.3734833 4.023854  12  81 0.00006949685
## 
##  F statistic for Roy's Greatest Root is an upper bound.

#permutation
#CCP::p.perm(X_SSD, Y_SSD, nboot=500, rhostart=1, type='Roy')

Significance testing (perm)

Scores

#calculate correlations
rxx <- cor(X_SSD, use = "p") #pairwise complete observations
ryy <- cor(Y_SSD, use = "p")
rxy <- cor(X_SSD,Y_SSD, use = "p")

#calculate omega
omega = t(solve(chol(rxx))) %*% rxy %*% solve(chol(ryy)) 

#decompose matrix
usv <- svd(omega)

#standardized weights (beta)
x.weights.std <- solve(chol(rxx)) %*% usv$u
y.weights.std <- solve(chol(ryy)) %*% usv$v

#structure coefficients
x.structures <- rxx %*% x.weights.std
y.structures <- ryy %*% y.weights.std

#canonical cross-loadings/redundancy
x.scores <- as.matrix(X_SSD) %*% x.weights.std
y.scores <- as.matrix(Y_SSD) %*% y.weights.std

#redundancy (?)
red.x.yscores <- cor(X_SSD, y.scores) 
red.y.xscores <- cor(Y_SSD, x.scores)

#bind together data
scores.plot <- as.data.frame(cbind(x.scores[, 1], y.scores[, 1]))

#make sure numeric
scores.plot[1:2] <- sapply(scores.plot[1:2],as.numeric)

#change names for clarity
names(scores.plot) <- c('Xset', 'Yset')

#make a scatter plot
ggplot(scores.plot, aes(x=Xset, y=Yset)) +
  geom_point(size = 4, alpha=.8) +
  geom_smooth(method=lm, color="black", alpha=.4) +
  theme_bw() +
  theme(axis.ticks=element_blank(), 
        panel.grid.minor = element_blank(), 
        panel.grid.major.x = element_blank(),
        legend.position='top')

Descriptives

        #for table, take just first column
        x.weights.std <- x.weights.std[, 1]
        x.weights.std <- t(data.frame(lapply(x.weights.std, type.convert), stringsAsFactors = F)) #get into df

        y.weights.std <- y.weights.std[, 1]
        y.weights.std <- t(data.frame(lapply(y.weights.std, type.convert), stringsAsFactors = F)) #get into df

        #bind together standardized weights
        weights <- rbind(x.weights.std, y.weights.std)

        #bind together structure coefficients
        structures <- rbind(x.structures, y.structures)

        #make table
        tbl <- as.data.frame(cbind(weights, structures))

        #add rownames as a variable
        tbl <- tibble::rownames_to_column(tbl, "variable")

        #redefine row names with full names
        tbl$variable <- c(names(X_SSD), names(Y_SSD))

        #set up dynamic names for table
        label <- paste("$R_c$=", round(cca_cca$cor[1], 5))
        myHeader <- c(" " = 1, label = 2) # create vector with colspan
        names(myHeader) <- c(" ", label) #set names

        #make into pretty table
        tbl_plot <- tbl

        #take only first 3 columns -- adjust if want more variates
        tbl_plot <- tbl_plot[, 1:3]

        #add a column for the communality coefficient
        tbl_plot$communalities <- tbl_plot$V2 ^ 2

        #round
        tbl_plot[,2:4] <- round(tbl_plot[,2:4], 3) #shorten here, because formatting makes characters (?!)

        tbl_plot %>%
          mutate(V2 = ifelse((V2 > .45 | V2 < -.45),
                  cell_spec(V2, color='black', underline = T, bold = T), #no green
                  cell_spec(V2, color='black'))) %>%
          kable(escape=F,
                align = c('l', 'c', 'c', 'c'),
                col.names = c('_Variable_', '$\\beta$', '$r_s$', '$r_s^2$ ($h^2$)')) %>%
        kable_styling(bootstrap_options=c('striped', 'hover', 'condensed')) %>%
        add_header_above(c(" " = 1, "_Function 1_" = 3)) %>%
        pack_rows('X', 1, 23) %>%
        pack_rows('Y', 24, 35) %>%
        footnote(general = "$r_s$ values > .450 are emphasized, following convention;
                 $\\beta$ = standardized canonical function coefficient;
                 $r_s$ = structure coefficient;
                 $r_s^2$ = squared structure coefficient, here also communality $h^2$")

	Function 1
Variable	\(\beta\)	\(r_s\)	\(r_s^2\) (\(h^2\))
X
AF_left	-0.563	-0.467	0.218
CB_left	0.054	-0.016	0.000
ILF_left	0.071	0.121	0.015
IOFF_left	-0.353	-0.22	0.048
UF_left	0.184	0.016	0.000
AF_right	0.266	-0.084	0.007
CB_right	0.281	-0.12	0.014
ILF_right	0.120	-0.065	0.004
IOFF_right	0.500	-0.204	0.041
UF_right	0.028	-0.284	0.080
CR_F_left	0.336	-0.12	0.014
CR_P_left	0.241	0.312	0.097
TF_left	-0.067	-0.387	0.149
CR_F_right	-0.241	-0.207	0.043
CR_P_right	0.182	-0.04	0.002
TF_right	-0.558	-0.552	0.305
CC1_commissural	-0.350	-0.476	0.227
CC2_commissural	-0.119	-0.41	0.168
CC3_commissural	-0.004	-0.196	0.038
CC4_commissural	0.700	0.253	0.064
CC5_commissural	-0.367	0.235	0.055
CC6_commissural	-0.739	-0.138	0.019
CC7_commissural	-0.221	-0.008	0.000
Y
Attention_vigilance	-0.005	0.112	0.013
Working_memory	0.390	0.419	0.176
Verbal_learning	0.204	0.537	0.288
Processing_speed	0.036	0.107	0.012
Visual_learning	-0.336	0.002	0.000
Problem_solving	-0.585	-0.43	0.185
TASIT_1	0.346	0.516	0.266
TASIT_2	0.350	0.396	0.157
TASIT_3	-0.117	0.333	0.111
RMET	-0.423	0.028	0.001
RAD	0.447	0.52	0.271
ER_40	-0.153	0.167	0.028
Note:
\(r_s\) values > .450 are emphasized, following convention; \(\beta\) = standardized canonical function coefficient; \(r_s\) = structure coefficient; \(r_s^2\) = squared structure coefficient, here also communality \(h^2\)

Redundancy

cca_candisc <- cancor(X_SSD, Y_SSD)
redundancy(cca_candisc)

## 
## Redundancies for the X variables & total X canonical redundancy
## 
##     Xcan1     Xcan2     Xcan3     Xcan4     Xcan5     Xcan6     Xcan7     Xcan8 
##  0.033298  0.017007  0.011580  0.010947  0.011234  0.008196  0.005255  0.005869 
##     Xcan9    Xcan10    Xcan11    Xcan12 total X|Y 
##  0.008416  0.004652  0.003240  0.004792  0.124485 
## 
## Redundancies for the Y variables & total Y canonical redundancy
## 
##     Ycan1     Ycan2     Ycan3     Ycan4     Ycan5     Ycan6     Ycan7     Ycan8 
##  0.059631  0.049249  0.080219  0.027579  0.011766  0.011239  0.024706  0.015119 
##     Ycan9    Ycan10    Ycan11    Ycan12 total Y|X 
##  0.005434  0.005410  0.004023  0.003645  0.298019