Canonical Correlation Analysis: Liver vs Biochemistry

Looking deeper into the liver

This notebook presents a demonstration canonical correlation analysis - a method for comparing inter-relationships between two sets of variables. In this case, we have three measures of liver health derived from MRI; liver iron, liver inflammation factor (LIF) and liver proton density fat fraction (PDFF). We also have a small set of biochemistry tests measured on the same participants.

Source the data and summarise the variables

source("R/0_setup.R")

library(CCA)

IP<- readRDS("DemoSet2.rds") # 502505 by 71

# Variable sets
biochem<- c("red_blood_cell_erythrocyte_count", "high_light_scatter_reticulocyte_count",                                     "alanine_aminotransferase", "triglycerides", "urate", "c_reactive_protein")

liver<- c("liver_iron", "LIF", "PDFF")

# Make the set

sett<- IP %>% select(all_of(liver), all_of(biochem)) %>% 
              filter(complete.cases(.)) 

sk<- skim(sett) %>% select(variable= skim_variable, mean= numeric.mean, min= numeric.p0, median= numeric.p50,
                           max= numeric.p100, numeric.hist)
sk[, 2:4]<- apply(sk[, 2:4], 2, function(x) round(x, 2))
print(sk)

Pairwise Panels for X and Y

Let’s begin by visualising the variables in X and Y and how they relate to each other.

# Separate and standardise the X and Y matrices
X<- sett[, biochem] %>% scale()
Y<- sett[, liver] %>% scale()

# Pairwise panel for the X matrix (biochemistry)
psych::pairs.panels(X, method = "pearson", hist.col = "#b3d9ff", density = TRUE)

# Pairwise panel for the Y matrix (liver)
psych::pairs.panels(Y, method = "pearson", hist.col = "#b3d9ff", density = TRUE)

The pairwise panels show that for both X and Y matrices, there is some visible correlation between variables. In this situation, canonical correlation analysis is ideal for exploring the overall relationship between the two sets.

Cross-Correlation Panel

correl<-matcor(X, Y)
img.matcor(correl,type = 3)

The cross-correlation panel gives us a preview of the range in which the cross-correlations are likely to sit. Nothing is strongly standing out, so we can expect cross-correlations in the range of 0-0.4 or so.

Applying the CCA method

The main CCA method is given here via the {CCA} package.

can.model<- cc(X,Y)
can.model$cor

[1] 0.42335127 0.06836749 0.04168428

These are the “canonical correlations”, which indicate the overall relationship between X and Y in 3 dimensions. We have three since our smallest input matrix has 3 columns.

Which canonical variates are significant?

rho<- can.model$cor
n<-   nrow(X)
p<-   ncol(X)
q<-   ncol(Y)
CCP::p.asym(rho, n, p, q, tstat="Wilks")

Wilks' Lambda, using F-approximation (Rao's F):
              stat    approx df1      df2     p.value
1 to 3:  0.8155178 47.522623  18 11441.47 0.000000000
2 to 3:  0.9935964  2.603395  10  8092.00 0.003733462
3 to 3:  0.9982624  1.761056   4  4047.00 0.133797814

An assymptotic test on the size of the correlations indicates that only the first two canonical correlations are significantly different from zero. Since the first canoncial correlation is so much larger than the second, the remainder of the analysis will focus on the first component.
In simple terms, the rest of the analysis is geared towards understanding how the matrices relate to each other, via this variate.

Standardised Canonical Weights

# Standardised Canonical Weights

sw1<-diag(sqrt(diag(cov(X)))) %*% can.model$xcoef %>% 
       data.frame() %>% 
       mutate(Var= biochem, Set= "Biochemistry") %>% 
       select(Var, X1, Set) 

sw2<-diag(sqrt(diag(cov(Y)))) %*% can.model$ycoef %>% 
      data.frame() %>% 
      mutate(Var= liver, Set= "Liver") %>% 
      select(Var, X1, Set) 

t0<- bind_rows(sw1, sw2) %>%
       mutate(Direction= ifelse(X1 > 0, "Positive", "Negative"))

ggplot(t0, aes(x= X1, y= fct_reorder(Var, X1), fill= Direction)) +
  geom_col(size= 0.2, color= "grey50") +
  geom_text(aes(label= sprintf("%2.1f", X1*100)), hjust= ifelse(t0$X1 > 0, 1.2, -0.3))+
  facet_wrap(~ Set, scales= "free") +
  theme_bw() +
  theme(legend.position = "None",
         strip.background = element_rect(fill= "grey30"),
         strip.text = element_text(color="white", face= "bold"))  +
  scale_fill_manual(values= c("#a3dca3", "#b3d9ff"))

sw1

sw2

Canonical loadings and cross-loadings

Canonical weights (standardised or not) have fallen out of use in favour of canonical loadings. These are the correlations between the original variables and the scores of the canonical variate.

# Canonical loadings and cross-loadings

can_loadings<- comput(X,Y,cancor)
can_loadings[3:6]

$corr.X.xscores
                                            [,1]        [,2]         [,3]
red_blood_cell_erythrocyte_count      -0.4108776  0.35522284 -0.163735110
high_light_scatter_reticulocyte_count -0.7215862 -0.09267242 -0.009252762
alanine_aminotransferase              -0.6873395  0.20784334 -0.110979461
triglycerides                         -0.6427031  0.35672320 -0.236499349
urate                                 -0.5187563 -0.47980607 -0.588432085
c_reactive_protein                    -0.4918877 -0.38984102  0.675329116

$corr.Y.xscores
                 [,1]         [,2]         [,3]
liver_iron -0.1825606 -0.050259395 -0.021804439
LIF        -0.2537582 -0.034332977  0.025982655
PDFF       -0.4214048 -0.005667556 -0.002000076

$corr.X.yscores
                                            [,1]        [,2]          [,3]
red_blood_cell_erythrocyte_count      -0.1739455  0.02428569 -0.0068251808
high_light_scatter_reticulocyte_count -0.3054844 -0.00633578 -0.0003856948
alanine_aminotransferase              -0.2909861  0.01420973 -0.0046260994
triglycerides                         -0.2720892  0.02438827 -0.0098583060
urate                                 -0.2196162 -0.03280314 -0.0245283701
c_reactive_protein                    -0.2082413 -0.02665245  0.0281506106

$corr.Y.yscores
                 [,1]       [,2]        [,3]
liver_iron -0.4312272 -0.7351359 -0.52308537
LIF        -0.5994034 -0.5021828  0.62332017
PDFF       -0.9954022 -0.0828984 -0.04798153

can_redunds<-  candisc::redundancy(candisc::cancor(X,Y))
can_redunds

Redundancies for the X variables & total X canonical redundancy

    Xcan1     Xcan2     Xcan3 total X|Y 
0.0623131 0.0005355 0.0002599 0.0631086 

Redundancies for the Y variables & total Y canonical redundancy

    Ycan1     Ycan2     Ycan3 total Y|X 
0.0917679 0.0012456 0.0003848 0.0933983 

# X Plot
t0<- tibble(Var= c(biochem, biochem),
            X1= c(can_loadings$corr.X.xscores[,1], can_loadings$corr.X.yscores[,1]),
            Loading= c(rep("X-X Loading", 6), rep("X-Y Cross-loading", 6)))
ggplot(t0, aes(x= X1, y= fct_reorder(Var, -X1), fill= Loading)) +
  geom_col(size= 0.2, color= "grey50") +
  geom_text(aes(label= sprintf("%2.1f", X1*100)), hjust= -0.6)+
  facet_wrap(~ Loading, scales= "free_x") +
  theme_bw() +
  theme(legend.position = "None",
         strip.background = element_rect(fill= "grey30"),
         strip.text = element_text(color="white", face= "bold")) +
  scale_fill_manual(values= c("#a3dca3", "#b3d9ff")) +
  labs(x="", y= "Variable", title= "Loadings and Cross-Loadings for X (Biochemistry)")

# Y Plot
t0<- tibble(Var= c(liver, liver),
            X1= c(can_loadings$corr.Y.yscores[,1], can_loadings$corr.Y.xscores[,1]),
            Loading= c(rep("Y-Y Loading", 3), rep("Y-X Cross-loading", 3))) %>%
     mutate(Loading= ordered(Loading, levels=c("Y-Y Loading", "Y-X Cross-loading")))
ggplot(t0, aes(x= X1, y= fct_reorder(Var, -X1), fill= Loading)) +
  geom_col(size= 0.2, color= "grey50") +
  geom_text(aes(label= sprintf("%2.1f", X1*100)), hjust= -0.6)+
  facet_wrap(~ Loading, scales= "free_x") +
  theme_bw() +
  theme(legend.position = "None",
         strip.background = element_rect(fill= "grey30"),
         strip.text = element_text(color="white", face= "bold")) +
  scale_fill_manual(values= c("#a3dca3", "#b3d9ff")) +
  labs(x="", y= "Variable", title= "Loadings and Cross-Loadings for Y (Liver)")

Examination of the loadings to the first canonical variate show that all variables included in our analysis have a part to play in explaining the relation between biochemistry and liver health. From the liver side, PDFF is the strongest variable, while from the biochemistry side the three most influential variables are light scatter reticulocyte count, alanine aminotransferase and triglycerides.

plt.cc(cc(X,Y),var.label = TRUE)

…

…

…