par(mfrow=c(1,1)) #install.packages("CCA") library(CCA) data(nutrimouse) # http://www6.toulouse.inra.fr/toxalim/pharmaco-moleculaire/acceuil.html # R package, CCA # The nutrimouse dataset comes from a nutrition study in the mouse. It was provided by Pascal Martin from the Toxicology and Pharmacology Laboratory (French National Institute for Agronomic Research). # n = 40, r (# of genes) = 120, s (# varianble in 2nd set) = 21 -- concentrations of 21 hepatic fatty acids. #boxplot(nutrimouse$lipid) #boxplot(nutrimouse$gene) # take the first 10 genes X=as.matrix(nutrimouse$gene[,1:10]) Y=as.matrix(nutrimouse$lipid) correl=matcor(X,Y) img.matcor(correl) img.matcor(correl,type=2) # canonical correlation res.cc=cc(X,Y) # correlations plot(res.cc$cor,type="b") # another view cancorrel = matcor(res.cc$scores$xscores ,res.cc$scores$yscores) img.matcor(cancorrel,type="b") par(mfrow=c(2,2)) plot(res.cc$scores$xscores[,1], res.cc$scores$yscores[,1],xlab="xi_1",ylab='omega_1', main="first canonical covariates") plot(res.cc$scores$xscores[,2], res.cc$scores$yscores[,2],xlab="xi_1",ylab='omega_1', main="second canonical covariates") plot(res.cc$scores$xscores[,3], res.cc$scores$yscores[,3],xlab="xi_1",ylab='omega_1', main="third canonical covariates") plot(res.cc$scores$xscores[,4], res.cc$scores$yscores[,4],xlab="xi_1",ylab='omega_1', main="fourth canonical covariates") par(mfrow=c(1,1)) # take all into account X=as.matrix(nutrimouse$gene) Y=as.matrix(nutrimouse$lipid) correl=matcor(X,Y) img.matcor(correl) img.matcor(correl,type=2) # regularized CCA # cross-validation by # http://www.math.univ-toulouse.fr/~igonzal/smpgd07-slides.pdf par(mfrow=c(3,3)) damp = 0.0001 res.cc=rcc(X,Y,damp,damp) plot(res.cc$cor,main = "Regularized by lambda = 0.0001") barplot(res.cc$xcoef[,1]) barplot(res.cc$ycoef[,1]) damp = 0.01 res.cc=rcc(X,Y,damp,damp) plot(res.cc$cor,main = "Regularized by lambda = 0.01") barplot(res.cc$xcoef[,1]) barplot(res.cc$ycoef[,1]) damp = 1 res.cc=rcc(X,Y,damp,damp) plot(res.cc$cor,main = "Regularized by lambda = 1") barplot(res.cc$xcoef[,1]) barplot(res.cc$ycoef[,1])