Code
# Paths
library(here)
# Multivariate analysis
library(ade4)
library(adegraphics)
# Matrix algebra
library(expm)
# Plots
library(ggplot2)
library(CAnetwork)
library(patchwork)Canonical Correlation Analysis (CCorA)
Code
# Define bound for comparison to zero
zero <- 10e-10Introduction
Canonical correlation analysis is a symmetrical analysis to examine the link between 2 matrices \(Y_R\) and \(Y_C\). It is intended to analyze quantitative continuous data: but here, we use it to analyze count data from inflated matrices.
Here, we will analyze matrix \(Y\), representing the occurrences of species in sites:
Inflate matrices
We compute the inflated tables \(R\) (\(n_\bar{0} \times r\)) and \(C\) (\(n_\bar{0} \times c\)) from \(Y\) (\(r \times c\)).
# Transform matrix to count table
Yfreq <- as.data.frame(as.table(Y))
colnames(Yfreq) <- c("row", "col", "Freq")
# Remove the cells with no observation
Yfreq0 <- Yfreq[-which(Yfreq$Freq == 0),]
Yfreq0$colind <- match(Yfreq0$col, colnames(Y)) # match index and species names
# Create indicator tables
YR <- acm.disjonctif(as.data.frame(Yfreq0$row))
colnames(YR) <- rownames(Y)
YR <- as.matrix(YR)
YC <- acm.disjonctif(as.data.frame(Yfreq0$col))
colnames(YC) <- colnames(Y)
YC <- as.matrix(YC)P <- Y/sum(Y)
Pfreq <- as.data.frame(as.table(P))
colnames(Pfreq) <- c("row", "col", "Freq")
# Remove the cells with no observation
Pfreq0 <- Pfreq[-which(Pfreq$Freq == 0),]
Pfreq0$colind <- match(Pfreq0$col, colnames(P)) # match index and species namesClassical canonical correlation analysis
We perform CCorA with the R function:
cancorRC <- cancor(YR,
YC,
xcenter = FALSE, ycenter = FALSE)
neig <- min(r-1, c-1)Now, we do it by hand. First, we compute matrix \[K = {S_{RR}^\top}^{-1/2} S_{RC} {S_{CC}^\top}^{-1/2}\] Where \(S_{AB}\) is the variance covariance matrix of matrices \(A\) and \(B\).
SRR <- t(YR) %*% YR
SRC <- t(YR) %*% YC
SCC <- t(YC) %*% YC
K <- solve(chol(t(SRR))) %*% SRC %*% solve(chol(SCC))
svdK <- svd(K)Compare singular values to canonical correlations:
all((svdK$d[1:neig]/cancorRC$cor[1:neig] - 1) < zero)[1] TRUEThen, we compute canonical coefficients:
# Canonical coefficients
xcoefK <- solve(chol(SRR)) %*% svdK$u
ycoefK <- solve(chol(t(SCC))) %*% svdK$v
all(abs(abs(xcoefK[, 1:neig]/cancorRC$xcoef[, 1:neig]) - 1) < zero)[1] TRUEall(abs(abs(ycoefK[, 1:neig]/cancorRC$ycoef[, 1:neig]) - 1) < zero)[1] TRUEAnd the canonical variates:
# Canonical variates
scoreRK <- YR %*% xcoefK
scoreCK <- YC %*% ycoefK
scoreR <- YR %*% cancorRC$xcoef
scoreC <- YC %*% cancorRC$ycoef
all(abs(abs(scoreRK[, 1:neig]/scoreR[, 1:neig]) - 1) < zero)[1] TRUEall(abs(abs(scoreCK[, 1:neig]/scoreC[, 1:neig]) - 1) < zero)[1] TRUEWeighted canonical correlation analysis
# Get weights
wt <- Pfreq0$Freq
Dw <- diag(wt)First, we scale the matrices with custom weights:
# Center tables
YR_scaled <- scalewt(YR, wt, scale = FALSE)
YC_scaled <- scalewt(YC, wt, scale = FALSE)Now, we perform a weighted version with ADE4.
cancorRC <- cancor(diag(sqrt(wt)) %*% YR_scaled,
diag(sqrt(wt)) %*% YC_scaled,
xcenter = FALSE, ycenter = FALSE)Now, we do it by hand. First, we compute matrix \(K\):
SRR <- t(YR) %*% Dw %*% YR
SRC <- t(YR_scaled) %*% Dw %*% YC_scaled
SCC <- t(YC) %*% Dw %*% YC
K <- solve(chol(t(SRR))) %*% SRC %*% solve(chol(SCC))
svdK <- svd(K, nu = 25, nv = 20)Compare singular values to canonical correlations:
all((svdK$d[1:neig]/cancorRC$cor[1:neig] - 1) < zero)[1] TRUEThen, we compute canonical coefficients:
# Canonical coefficients
xcoefK <- solve(chol(SRR)) %*% svdK$u
ycoefK <- solve(chol(t(SCC))) %*% svdK$v
xcoefK[1:(r-1), ]/cancorRC$xcoef [,1] [,2] [,3] [,4] [,5] [,6]
1 0.02856273 -1.4279321 0.28175678 -0.30578706 2.981125219 -1.0733512
2 -0.68164197 -0.8302571 0.19316330 -2.43709227 -8.277800302 -0.7146887
3 0.02647019 -1.6225095 -0.16059216 0.13168986 -0.107103007 -2.2693713
4 -0.26009508 -1.5104118 -0.35755537 -0.42652154 -0.263473731 -0.9307892
5 -0.21151429 -1.7756576 -0.35889045 0.45752661 0.279372125 -1.7411333
6 -0.36482860 -3.6255560 -0.24419283 -0.33392424 -0.340050834 -0.7462876
7 -0.29679540 -0.6792450 -0.34363000 -0.58953651 -0.147895938 -0.1436754
8 1.84806898 -0.8227405 2.72718615 1.97318399 -0.606594445 -1.1932838
9 -0.09089339 -1.5641257 -0.44942699 -0.40471671 0.373864141 -1.2105851
10 -0.38350610 -0.6677276 0.22254515 -0.56201016 0.158216977 -1.0413845
11 -0.48151529 19.9796383 -0.41520339 0.10266166 -0.657341396 -0.8681612
12 -0.14342694 -0.5650092 -0.16266620 -0.57762970 4.168794211 -1.1457626
13 -5.34243891 -0.9065452 -0.07898875 -0.08688740 0.008964368 -0.8394007
14 1.34542164 -1.3262371 -0.15576563 2.13837992 -4.207172775 -2.5321872
15 5.79418273 -0.8681254 0.43525069 -0.53249301 -1.387522789 -0.9682030
16 0.14368795 -1.3505090 -0.38759196 0.10296499 1.422970140 -0.8840121
17 -0.09577067 -1.7053198 0.21611608 -0.44223632 0.221509198 -0.8325864
18 0.20299484 -1.3071248 -0.11796715 -0.03264285 1.733884367 -0.9133820
19 1.58127716 0.6487321 0.48252719 0.54231535 -0.519331012 -0.8993666
20 0.06093118 -1.8048533 0.01665392 -0.36884287 0.722577358 -0.9354730
21 1.44649963 -1.2390520 -0.60709717 3.18028540 -4.000778869 -1.9779093
22 0.26522088 -1.5047232 0.07271057 0.26676779 -0.307931653 -0.6860382
23 1.09703385 -0.4309448 0.54154815 0.35274069 -0.618132642 -0.8515017
24 2.55136234 -4.8417623 0.52748847 36.18812571 -0.012843839 -1.0543019
25 62.12195151 0.1388428 -2.69958370 -8.88729369 1.165523354 -1.1246051
[,7] [,8] [,9] [,10] [,11] [,12]
1 -0.4878032 0.97649296 1.11063181 0.28396721 2.86376219 0.47914792
2 14.4606672 0.88191535 0.28039479 0.08966871 -0.07206365 -0.18444939
3 -0.3558283 1.13010958 0.24038957 2.04235361 0.52160813 0.33895844
4 -13.1162605 0.94975040 0.73948188 5.10052793 0.21040563 0.77947799
5 -6.5973674 1.18241492 0.45877257 1.27788823 0.71260630 0.60134395
6 -0.4747574 1.14544086 0.45893584 0.76634908 1.43774605 0.87642842
7 1.4127476 1.13856225 0.53368568 0.60538299 0.74258618 3.65088420
8 -0.5233280 0.94462675 0.09005177 -1.85943316 0.66666151 0.53116340
9 -1.7095934 1.04475742 0.51448798 1.20308305 4.21088640 0.78265676
10 0.2474649 0.90278231 -0.10435726 5.21426849 0.10066761 -0.02506802
11 -0.7372381 1.04772596 1.72249328 0.14164086 -2.58983837 1.73109122
12 -2.6833721 1.05205193 0.28978551 0.63496093 2.09555129 1.43561700
13 3.1014514 0.81001866 0.07517764 0.60923928 7.00127811 0.14663020
14 0.2883166 1.07295245 5.32021987 0.59272974 0.04358356 4.26937024
15 -0.8605384 1.03223123 1.49413300 1.26711241 0.70496989 0.73377342
16 -2.1076658 3.52741398 1.47586147 0.85129592 0.68684280 0.54671978
17 -0.6665257 0.97803485 1.56728278 2.09913974 0.26109694 -4.31671862
18 -0.3408633 0.93813776 0.75187532 -0.01500123 1.37269779 1.85986371
19 -1.5250799 1.11425991 3.88815493 1.32585972 2.98565616 2.20427367
20 -0.5089793 0.96090805 0.81072968 0.57326607 0.30308278 1.72857814
21 -2.1734407 1.07278201 3.08468981 0.69626683 0.62893432 -0.20240914
22 -0.1044161 0.90562351 0.70599250 -1.50599388 0.61548603 3.20418918
23 -1.7954181 0.70732071 -11.26391387 0.33062537 0.39835748 0.46426227
24 -0.7757194 1.13647649 0.76003873 0.68179912 0.41972064 0.41300486
25 -0.3408509 0.01677946 -0.96560152 1.88194131 1.67741512 -14.89054034
[,13] [,14] [,15] [,16] [,17] [,18]
1 0.47483627 -2.97570665 0.61738444 -74.34121804 -0.60121512 0.59109788
2 -3.48184482 -0.20154831 -1.54717251 -0.02803428 5.64198559 0.12307293
3 -4.30610650 0.03763125 0.86161131 -6.18680505 -1.20093504 0.12892607
4 0.52187287 2.97195260 0.65796501 0.73869671 0.42680438 0.49230721
5 0.74196976 0.06297276 0.83356404 -2.28940918 -0.88396676 6.21874368
6 2.33511326 0.50588778 0.24721017 1.74247436 0.03786166 -0.32475278
7 1.40710171 -2.00072615 0.69631429 1.91427122 -2.99228477 0.10917849
8 1.29813046 0.38753345 1.84152073 0.52233678 -0.54384503 0.43081655
9 0.66289928 1.90352951 1.13284050 0.29508101 -2.00164895 0.14057387
10 0.67528307 0.61564082 -0.49605777 0.54898942 -0.32032988 -0.04825296
11 0.72710236 -0.28208113 -5.32793198 0.20735034 -0.50584029 -0.09426923
12 1.89069940 -0.85282039 -1.65464296 0.55565673 -0.18484450 -0.11720708
13 0.82885819 -2.15029399 0.62141131 5.19699878 -0.43276860 0.03139150
14 -1.21096270 0.41390252 0.56795836 1.62210349 -1.77036168 -4.35625655
15 2.85252450 0.57787774 2.72684255 0.72076097 2.23955212 0.29382495
16 -2.97395391 -0.52664285 0.52255725 0.79995263 -0.55973017 -0.68135577
17 2.59170676 3.90477833 2.33280940 1.88222129 -4.32719549 8.09593418
18 0.44526448 0.29048338 0.81066557 0.61307231 -0.56099690 0.37786031
19 1.19236721 0.33466126 0.46868979 4.36729385 -0.77826842 0.55932609
20 132.14078704 0.43433604 1.20619918 10.30991184 -0.75981426 -1.72107851
21 0.76443664 0.43795971 1.29198322 1.65917378 -0.59716231 0.54179762
22 0.82937268 0.42432897 -1.46188754 -0.07400184 -1.25071072 0.43037406
23 -0.04502113 0.36304249 2.94268721 0.52252737 -1.90729352 -5.40948873
24 3.98561919 4.38989880 -0.10569180 0.25771647 -0.65815909 0.21912943
25 1.99337708 0.25446321 -0.04873034 0.60135152 -1.59277480 14.41526153
[,19] [,20] [,21] [,22] [,23] [,24]
1 -0.349344610 -0.41020999 -0.27257603 -1.99027470 0.2898110 -0.04153318
2 -0.031892622 0.02609518 0.88509473 0.29638159 -17.1260930 -1.70344828
3 -0.182155560 -0.29962135 0.47261060 1.27439145 -4.8907105 0.74059445
4 0.875084265 -0.54193570 -1.83598171 0.03658023 -1.8575407 -0.49288098
5 -0.190033062 0.28031799 -4.78828857 -11.04932299 0.2192397 10.05245941
6 0.053219846 0.48334576 2.92252652 -1.12936411 -0.2616219 0.73195351
7 0.041244449 -0.25524133 -17.28287951 0.96624634 5.3898282 -0.29357864
8 -0.217599494 1.00648471 1.39135997 12.87462358 -1.0391273 1.97313247
9 -0.241790859 0.34709774 -0.31158600 -0.20649307 0.4241271 -1.25448799
10 0.524259640 -0.07964137 0.21085236 1.41617988 -0.3242136 -2.78391261
11 0.100540504 -0.15005335 0.71246489 0.70114695 -1.3811112 -2.17393044
12 -0.221705252 0.73961865 -8.57362208 -2.81071533 0.4221634 -0.73194173
13 0.028398676 -0.07598651 0.65556954 0.33006757 4.5211127 -4.16643808
14 -0.034522726 0.04511354 -1.28512702 -0.87478938 -0.7236337 100.62567577
15 -0.003325141 0.24934659 3.93942649 0.37152564 14.5940992 0.24452783
16 0.192722611 0.62160191 0.24310145 -0.03251290 -3.4062737 -1.57339846
17 0.216858429 6.28956953 0.55610121 -2.12499626 3.1471710 -0.05492724
18 -0.196115921 1.03687188 -0.86476378 -0.83927038 -1.1152451 -0.77825172
19 0.624864261 0.19466050 1.21240596 -0.19489976 -2.2359949 -3.06313664
20 -0.130736358 -0.48170187 1.10153048 -0.35169042 -2.8814839 39.85728892
21 0.067975188 -0.31039241 0.08784137 -1.13924424 1.9441293 -1.29676754
22 0.053874374 11.27145430 -0.92009320 -1.74335359 -0.3675612 1.03686211
23 -0.194346850 -0.11035707 0.03444202 -0.40076599 -3.0108416 0.33537646
24 0.037245096 0.07366010 0.88509473 0.29638159 -17.1260930 -1.70344828
25 0.114008395 -0.51317679 -3.46139298 1.09607572 -0.2779958 0.28015362
[,25]
1 -0.011568783
2 -0.080015957
3 0.250183880
4 0.018641378
5 -0.150872794
6 -0.153408989
7 -0.309069918
8 0.015139288
9 -0.062147042
10 -0.008399421
11 -0.061165449
12 -0.021392366
13 -0.060263505
14 0.147244547
15 -0.102850032
16 -0.110241089
17 -0.169822351
18 -0.141179287
19 -0.100916197
20 0.147909566
21 -0.180921155
22 -0.242599584
23 -0.051590025
24 -0.080015957
25 -0.354832746ycoefK[1:(c-1), ]/cancorRC$ycoef [,1] [,2] [,3] [,4] [,5] [,6]
sp1 -0.3206920 -0.3918323 -0.2053836 -0.07896689 -0.6843283 -0.2765287
sp2 -0.6223393 -1.5253529 0.7797516 -0.44323344 -1.1951936 -0.5397811
sp3 0.9829681 -0.5416906 0.6348328 -1.58275444 -0.7535566 -1.2142380
sp4 -0.1144668 -0.3724289 0.6258095 -0.42808183 -3.2286382 -1.5231839
sp5 -0.6886700 -1.1490643 0.6558947 -2.00202077 -0.9391458 -0.8970107
sp6 -0.2353731 -0.5758130 -6.1838709 -0.26184873 -0.5893347 -0.4165000
sp7 -0.6040801 -0.5440551 0.8591842 -0.75907885 0.1259773 -1.1179251
sp8 1.8331627 -0.3046293 0.2942855 5.08906572 141.5749623 -0.3446304
sp9 0.1088741 -0.7290569 -0.2393746 -5.16790319 -0.5106421 -0.8179855
sp10 -0.2285032 4.8950199 0.6341077 0.45810516 -0.5737832 -0.5248353
sp11 -0.3012959 -0.6462378 -0.8355436 -0.33394363 -0.4970987 1.9089822
sp12 0.1585970 -1.3534852 -0.3836982 -1.88275932 -0.7543297 -3.7388062
sp13 -1.2819071 -1.2754031 0.5668180 -0.84593634 -0.9072413 -0.5501736
sp14 -0.5466955 0.2080719 0.8524276 -0.49783796 -1.0942941 -0.7817410
sp15 -1.5423338 -9.1712585 -0.8023357 -0.66742969 -1.2418855 -0.6919100
sp16 -1.4747134 -1.5794849 0.4473012 -0.72728601 -0.8662902 -2.0494508
sp17 22.1509449 0.7527273 -1.0882740 -1.46139822 0.9320075 -1.3583809
sp18 2.3485268 -0.1210567 13.8514375 0.47127422 -1.2592610 -0.4551580
sp19 -0.2387618 -0.5551409 -6.9633955 -0.47809033 -0.5600797 0.6642413
sp21 -0.6694132 -1.2661550 1.2881805 -0.62553089 -1.4757146 -7.9283583
[,7] [,8] [,9] [,10] [,11] [,12]
sp1 0.2959162 -0.2769511 -0.5138021 1.0368710 0.2541210 0.0887245548
sp2 -1.2655302 0.1415636 -0.3197066 0.9588997 -6.1598085 -0.0790847669
sp3 -0.8257572 0.9199730 0.8937208 1.0184481 2.1678595 -0.0210570528
sp4 -0.7435170 0.4048965 -2.2754349 2.1496216 0.4457759 -0.0002848192
sp5 -0.9332524 1.2487255 0.7421351 1.0105397 0.7606758 12.9131924281
sp6 -2.2626495 0.3112946 0.4280161 0.9868145 0.7306161 -0.1721994504
sp7 -0.9469877 1.4970668 1.1924591 0.9923087 0.7899175 -1.0873887945
sp8 -2.1463116 0.4862592 0.2710334 0.7743712 -0.7240926 -0.0822760584
sp9 -0.7650638 0.8673584 1.2950009 0.9426169 1.2555752 0.0616131382
sp10 -0.7602276 -0.9041519 0.3521552 0.9250803 -0.6029831 -0.0646770646
sp11 -1.7721359 3.1683072 0.1653756 1.0388643 -0.6981898 -0.0561522312
sp12 -2.6389967 -1.0190020 -2.9186849 0.9706587 0.4299928 0.3601474823
sp13 -1.2610451 0.7298176 -3.2057448 0.8557230 -0.5371215 0.1033786059
sp14 -1.0811060 4.0627160 0.6889667 1.0093784 1.8423115 0.2629667464
sp15 -0.8575001 1.4323253 0.6477168 0.9695854 1.3184301 0.2084798032
sp16 -1.3011255 0.6458378 -3.0322722 1.0152158 0.7130269 -0.1772533147
sp17 -2.1724796 1.3618953 -0.4635685 0.3082314 1.3423466 -0.0293750531
sp18 -0.6645233 -0.2318054 -0.1198323 1.0528794 0.7127659 0.1264080987
sp19 -0.4825079 4.6021851 0.3926907 0.9657913 1.4668998 0.1539695012
sp21 -0.4546468 0.7318890 0.1801397 1.0885759 2.7166412 -0.1674831945
[,13] [,14] [,15] [,16] [,17] [,18]
sp1 0.002470018 0.77095493 -0.05017785 -0.34503625 -4.00561471 0.45435166
sp2 0.821287319 0.81269278 -0.40760478 -1.10952707 2.32307382 0.37281958
sp3 0.693060486 0.88310740 0.29601025 2.16037117 -0.74438628 1.35786141
sp4 0.445316629 1.14801603 0.09795254 0.21424458 -1.51949815 -0.16370021
sp5 0.891936101 1.06737441 -0.09159134 0.66338350 -0.31566513 1.11273630
sp6 0.019659567 1.49185030 0.40024659 -0.11962761 -0.32930306 0.07928034
sp7 1.176333322 0.98354642 0.20476546 0.50471400 -0.96123400 1.08809204
sp8 0.075050996 0.79015477 -0.81630743 0.46175312 0.92511305 0.14904314
sp9 0.574904643 1.15889179 0.52083075 -0.39243356 -0.78423843 2.69662474
sp10 1.433598052 -0.04519393 0.24431031 -0.11878771 -0.34965583 0.74068866
sp11 0.567270842 -0.10417639 0.47881819 0.05066364 -0.70620285 3.41625822
sp12 0.285850612 2.07477837 -0.50128554 -0.41081565 -0.60022028 2.35136591
sp13 1.346947164 0.88562168 0.62066515 -0.08700197 -1.10060668 1.19314680
sp14 1.068515550 1.04208064 0.46005554 -0.42316623 -0.90064739 1.11962754
sp15 0.804264900 0.67470173 0.62078386 -1.02896785 -0.70849161 0.48137277
sp16 0.723990619 1.10243893 2.85921934 0.20505978 -0.93452348 0.87788113
sp17 -0.702845737 2.23892531 0.74459729 0.30124642 -2.25281730 0.65629233
sp18 -0.395698696 0.10067506 -0.37086792 -0.15876033 -0.04479285 -2.02237483
sp19 -0.723035446 1.57991743 2.03987835 0.03398848 0.09780847 -8.92740948
sp21 2.579096215 0.42610603 0.42706830 0.37002862 -0.39676232 0.46555329
[,19] [,20]
sp1 -1.8949223 -0.1777219
sp2 -0.6677453 -0.5724865
sp3 -0.4572508 -0.6069470
sp4 -13.9076394 1.7416979
sp5 -1.9266635 -0.7709239
sp6 -4.0054196 -0.3849816
sp7 -1.1941697 -0.7973706
sp8 3.0564047 1.6828834
sp9 -0.4744434 -0.3741881
sp10 -0.7439925 3.1885706
sp11 -0.7764606 -3.0636462
sp12 0.3761583 -2.7331721
sp13 -0.8494519 -0.5932805
sp14 -1.0886923 -0.5272949
sp15 -1.1727279 -1.8475689
sp16 0.5613588 -0.7053780
sp17 -0.2916728 -0.8416820
sp18 -0.6780628 -0.4516535
sp19 -0.5632682 -0.5347594
sp21 2.0013854 0.3667632all(abs(abs(xcoefK[1:(r-1), ]/cancorRC$xcoef) - 1) < zero)[1] FALSEall(abs(abs(ycoefK[1:(c-1), ]/cancorRC$ycoef) - 1) < zero)[1] FALSEAnd the canonical variates:
# Canonical variates with ade4
scoreR <- YR_scaled[, 1:(r-1)] %*% cancorRC$xcoef
scoreC <- YC_scaled[, 1:(c-1)] %*% cancorRC$ycoef
# Check R C are normed at wt
res <- t(scoreR) %*% diag(wt) %*% scoreR # R^t Dw R = 1
all(abs(res - diag(1, ncol(scoreR))) < zero) [1] TRUEres <- t(scoreC) %*% diag(wt) %*% scoreC # R^t Dw R = 1
all(abs(res - diag(1, ncol(scoreC))) < zero)[1] TRUE# Canonical variates with function
scoreRK <- YR %*% xcoefK
scoreCK <- YC %*% ycoefK
# Check R C are normed at wt
res <- t(scoreRK) %*% diag(wt) %*% scoreRK # R^t Dw R = 1
all(abs(res - diag(1, ncol(scoreRK))) < zero) [1] TRUEres <- t(scoreCK) %*% diag(wt) %*% scoreCK # R^t Dw R = 1
all(abs(res - diag(1, ncol(scoreCK))) < zero)[1] TRUEall(abs(abs(scoreR[,1:neig]/scoreRK[,1:neig]) - 1) < zero)[1] TRUEall(abs(abs(scoreC[,1:neig]/scoreCK[,1:neig]) - 1) < zero)[1] TRUE