Ever been in a scenario where you needed to come up with pairwise correlation, covariance, or cosine matrices for data on the fly without the help of a function? Probably not. Even so, it’s worth taking a look at how these matrices could be calculated as there are some neat commonalities in their respective calculations. Playing a role in it all is the \(X^tX\) matrix, an old friend we met a few posts back when illustrating multiple regression from an algebraic viewpoint. Here, you can think of \(X\) as some data you’ve collected where every column is a different measure and every row a different subject.

Below, we will illustrate what was just discussed by creating a function that outputs one of the aforementiond matrices based on user input.

Our function

my_fun<-function(x,type =c("ss","cov","cor","cosine")){
  type<-match.arg(type)
  out<-matrix(rep(0,dim(x)[2]*dim(x)[2]),nrow=dim(x)[2])
  if (type == "cov"){
    x<-scale(x,center=TRUE,scale=FALSE)
    out<-t(x)%*%(x)/(nrow(x)-1)
  } else if (type =="cor"){
    x<-scale(x,center=TRUE,scale=TRUE)
    out<-t(x)%*%(x)/(nrow(x)-1)
  } else if (type == "cosine"){
    x<-apply(x, 2,function(x) x / sqrt(sum(x^2)) )
    out<-t(x)%*%(x)
  } else {
     out<-t(x)%*%(x)
  }
  return(out)   
}

Create some data

Now, that we have our function let’s make some data to apply it to.

dat<-cbind(replicate(4,sample(1:20,100,rep=TRUE)))

SSCP

We can now compare the output of our function to the complement functions in base R. Starting with the \(SSCP\) matrix

crossprod(dat)
      [,1]  [,2]  [,3]  [,4]
[1,] 16498 11696 12111 12395
[2,] 11696 12972  9451 10356
[3,] 12111  9451 13220  9592
[4,] 12395 10356  9592 14279
my_fun(dat, type = "ss")
      [,1]  [,2]  [,3]  [,4]
[1,] 16498 11696 12111 12395
[2,] 11696 12972  9451 10356
[3,] 12111  9451 13220  9592
[4,] 12395 10356  9592 14279

Covariance

cov(dat)
          [,1]      [,2]      [,3]      [,4]
[1,] 33.988283  3.772929  8.890909  1.688687
[2,]  3.772929 32.429899 -2.337374 -1.878384
[3,]  8.890909 -2.337374 36.525253 -8.733333
[4,]  1.688687 -1.878384 -8.733333 29.233434
my_fun(dat, type = "cov")
          [,1]      [,2]      [,3]      [,4]
[1,] 33.988283  3.772929  8.890909  1.688687
[2,]  3.772929 32.429899 -2.337374 -1.878384
[3,]  8.890909 -2.337374 36.525253 -8.733333
[4,]  1.688687 -1.878384 -8.733333 29.233434

Correlation

cor(dat)
           [,1]        [,2]        [,3]        [,4]
[1,] 1.00000000  0.11364262  0.25233933  0.05357283
[2,] 0.11364262  1.00000000 -0.06791391 -0.06100587
[3,] 0.25233933 -0.06791391  1.00000000 -0.26726587
[4,] 0.05357283 -0.06100587 -0.26726587  1.00000000
my_fun(dat, type = "cor")
           [,1]        [,2]        [,3]        [,4]
[1,] 1.00000000  0.11364262  0.25233933  0.05357283
[2,] 0.11364262  1.00000000 -0.06791391 -0.06100587
[3,] 0.25233933 -0.06791391  1.00000000 -0.26726587
[4,] 0.05357283 -0.06100587 -0.26726587  1.00000000

Cosine similarity

We’ll load the library “philentropy” to check our work here as it contains many useful distance functions. Note that we are transposing our data as the default behavior of this function is to make pairwise comparisons of all rows.

library(philentropy)
distance(t(dat), method = "cosine")
          v1        v2        v3        v4
v1 1.0000000 0.7994996 0.8200657 0.8075734
v2 0.7994996 1.0000000 0.7217031 0.7609212
v3 0.8200657 0.7217031 1.0000000 0.6981433
v4 0.8075734 0.7609212 0.6981433 1.0000000
my_fun(dat, type = "cosine")
          [,1]      [,2]      [,3]      [,4]
[1,] 1.0000000 0.7994996 0.8200657 0.8075734
[2,] 0.7994996 1.0000000 0.7217031 0.7609212
[3,] 0.8200657 0.7217031 1.0000000 0.6981433
[4,] 0.8075734 0.7609212 0.6981433 1.0000000
LS0tDQp0aXRsZTogIlhgWCwgY292YXJpYW5jZSwgY29ycmVsYXRpb24sIGFuZCBjb3NpbmUgbWF0cmljZXMiDQpvdXRwdXQ6IGh0bWxfbm90ZWJvb2sNCi0tLQ0KDQo8c3R5bGUgdHlwZT0idGV4dC9jc3MiPg0KDQpoMS50aXRsZSB7DQogIGZvbnQtc2l6ZTogMzhweDsNCiAgY29sb3I6IEJsYWNrOw0KICBmb250LWZhY2U6IEJvbGQ7DQp9DQoNCmJvZHksIHRkIHsNCiAgIGZvbnQtc2l6ZTogMThweDsNCn0NCmNvZGUucnsNCiAgZm9udC1zaXplOiAxNnB4Ow0KfQ0KcHJlIHsNCiAgZm9udC1zaXplOiAxNHB4DQp9DQo8L3N0eWxlPg0KPGJyPg0KDQoNCkV2ZXIgYmVlbiBpbiBhIHNjZW5hcmlvIHdoZXJlIHlvdSBuZWVkZWQgdG8gY29tZSB1cCB3aXRoIHBhaXJ3aXNlIGNvcnJlbGF0aW9uLCBjb3ZhcmlhbmNlLCBvciBjb3NpbmUgbWF0cmljZXMgZm9yIGRhdGEgb24gdGhlIGZseSB3aXRob3V0IHRoZSBoZWxwIG9mIGEgZnVuY3Rpb24/IFByb2JhYmx5IG5vdC4gIEV2ZW4gc28sIGl0J3Mgd29ydGggdGFraW5nIGEgbG9vayBhdCBob3cgdGhlc2UgbWF0cmljZXMgY291bGQgYmUgY2FsY3VsYXRlZCBhcyB0aGVyZSBhcmUgc29tZSBuZWF0IGNvbW1vbmFsaXRpZXMgaW4gdGhlaXIgcmVzcGVjdGl2ZSBjYWxjdWxhdGlvbnMuICBQbGF5aW5nIGEgcm9sZSBpbiBpdCBhbGwgaXMgdGhlICRYXnRYJCBtYXRyaXgsIGFuIG9sZCBmcmllbmQgd2UgbWV0IDxhIGhyZWY9Imh0dHBzOi8vZGFzaWx2YWExMC5naXRodWIuaW8vL2IxLyI+PGI+IGEgZmV3IHBvc3RzIGJhY2s8L2I+PC9hPiB3aGVuIGlsbHVzdHJhdGluZyBtdWx0aXBsZSByZWdyZXNzaW9uIGZyb20gYW4gYWxnZWJyYWljIHZpZXdwb2ludC4gSGVyZSwgeW91IGNhbiB0aGluayBvZiAkWCQgYXMgc29tZSBkYXRhIHlvdSd2ZSBjb2xsZWN0ZWQgd2hlcmUgZXZlcnkgY29sdW1uIGlzIGEgZGlmZmVyZW50IG1lYXN1cmUgYW5kIGV2ZXJ5IHJvdyBhIGRpZmZlcmVudCBzdWJqZWN0LiAgDQoNCiogIFRha2luZyB0aGUgZGF0YSAkWCQgYXMgaXMgYW5kIG11bHRpcGx5aW5nIGl0cyB0cmFuc3Bvc2UgYnkgaXRzZWxmICRYXnRYJCByZXN1bHRzIGluIHRoZSBzdW0gb2Ygc3F1YXJlcyBjcm9zcyBwcm9kdWN0cyBtYXRyaXggKFNTQ1ApIHdoZXJlIFNTIGZhbGwgb24gdGhlIGRpYWdvbmFsIGFuZCBjcm9zcyBwcm91ZGN0cyBvbiB0aGUgb2ZmIGRpYWdvbmFsLg0KDQoqICBDZW50ZXJpbmcgJFgkLCBtdWx0aXBseWluZyBpdHMgdHJhbnNwb3NlIGJ5IGl0c2VsZiAkWF50WCQgYW5kIGRpdmlkaW5nIGJ5IG4gLSAxICh3aGVyZSBuID0gIyBvZiByb3dzIGluIFgpIHJlc3VsdHMgaW4gdGhlIHZhcmlhbmNlIGNvdmFyaWFuY2UgbWF0cml4IHdpdGggdmFyaWFuY2VzIG9uIHRoZSBkaWFnb25hbCBhbmQgY292YXJpYW5jZXMgb24gdGhlIG9mZiBkaWFnb25hbC4NCg0KKiAgU3RhbmRhcmRpemluZyAkWCQsIG11bHRpcGx5aW5nIGl0cyB0cmFuc3Bvc2UgYnkgaXRzZWxmICRYXnRYJCBhbmQgZGl2aWRpbmcgYnkgbiAtIDEgKHdoZXJlIG4gPSAjIG9mIHJvd3MgaW4gWCkgcmVzdWx0cyBpbiB0aGUgcGVhcnNvbiBjb3JyZWxhdGlvbiBiZXR3ZWVuIHZhcmlhYmxlIHBhaXJzLg0KDQoqICBVbml0LXNjYWxpbmcgJFgkIGFuZCBtdWx0aXBseWluZyBpdHMgdHJhbnNwb3NlIGJ5IGl0c2VsZiAkWF50WCQgcmVzdWx0cyBpbiB0aGUgY29zaW5lIHNpbWlsYXJpdHkgYmV0d2VlbiB2YXJpYWJsZSBwYWlycy4NCg0KQmVsb3csIHdlIHdpbGwgaWxsdXN0cmF0ZSB3aGF0IHdhcyBqdXN0IGRpc2N1c3NlZCBieSBjcmVhdGluZyBhIGZ1bmN0aW9uIHRoYXQgb3V0cHV0cyBvbmUgb2YgdGhlIGFmb3JlbWVudGlvbmQgbWF0cmljZXMgYmFzZWQgb24gdXNlciBpbnB1dC4NCg0KDQojIyMqKk91ciBmdW5jdGlvbioqDQoNCmBgYHtyfQ0KbXlfZnVuPC1mdW5jdGlvbih4LHR5cGUgPWMoInNzIiwiY292IiwiY29yIiwiY29zaW5lIikpew0KICB0eXBlPC1tYXRjaC5hcmcodHlwZSkNCiAgb3V0PC1tYXRyaXgocmVwKDAsZGltKHgpWzJdKmRpbSh4KVsyXSksbnJvdz1kaW0oeClbMl0pDQogIGlmICh0eXBlID09ICJjb3YiKXsNCiAgICB4PC1zY2FsZSh4LGNlbnRlcj1UUlVFLHNjYWxlPUZBTFNFKQ0KICAgIG91dDwtdCh4KSUqJSh4KS8obnJvdyh4KS0xKQ0KICB9IGVsc2UgaWYgKHR5cGUgPT0iY29yIil7DQogICAgeDwtc2NhbGUoeCxjZW50ZXI9VFJVRSxzY2FsZT1UUlVFKQ0KICAgIG91dDwtdCh4KSUqJSh4KS8obnJvdyh4KS0xKQ0KICB9IGVsc2UgaWYgKHR5cGUgPT0gImNvc2luZSIpew0KICAgIHg8LWFwcGx5KHgsIDIsZnVuY3Rpb24oeCkgeCAvIHNxcnQoc3VtKHheMikpICkNCiAgICBvdXQ8LXQoeCklKiUoeCkNCiAgfSBlbHNlIHsNCiAgICAgb3V0PC10KHgpJSolKHgpDQogIH0NCiAgcmV0dXJuKG91dCkgICANCn0NCg0KYGBgDQoNCiMjIyoqQ3JlYXRlIHNvbWUgZGF0YSoqDQpOb3csIHRoYXQgd2UgaGF2ZSBvdXIgZnVuY3Rpb24gbGV0J3MgbWFrZSBzb21lIGRhdGEgdG8gYXBwbHkgaXQgdG8uDQoNCmBgYHtyfQ0KDQpkYXQ8LWNiaW5kKHJlcGxpY2F0ZSg0LHNhbXBsZSgxOjIwLDEwMCxyZXA9VFJVRSkpKQ0KDQpgYGANCg0KIyMjKipTU0NQKioNCldlIGNhbiBub3cgY29tcGFyZSB0aGUgb3V0cHV0IG9mIG91ciBmdW5jdGlvbiB0byB0aGUgY29tcGxlbWVudCBmdW5jdGlvbnMgaW4gYmFzZSBSLiAgU3RhcnRpbmcgd2l0aCB0aGUgJFNTQ1AkIG1hdHJpeA0KDQpgYGB7cn0NCmNyb3NzcHJvZChkYXQpDQoNCm15X2Z1bihkYXQsIHR5cGUgPSAic3MiKQ0KDQpgYGANCg0KIyMjKipDb3ZhcmlhbmNlKioNCg0KYGBge3J9DQpjb3YoZGF0KQ0KDQpteV9mdW4oZGF0LCB0eXBlID0gImNvdiIpDQoNCmBgYA0KDQojIyMqKkNvcnJlbGF0aW9uKioNCg0KYGBge3J9DQpjb3IoZGF0KQ0KDQpteV9mdW4oZGF0LCB0eXBlID0gImNvciIpDQoNCmBgYA0KIyMjKipDb3NpbmUgc2ltaWxhcml0eSoqDQoNCldlJ2xsIGxvYWQgdGhlIGxpYnJhcnkgInBoaWxlbnRyb3B5IiB0byBjaGVjayBvdXIgd29yayBoZXJlIGFzIGl0IGNvbnRhaW5zIG1hbnkgdXNlZnVsIGRpc3RhbmNlIGZ1bmN0aW9ucy4gIE5vdGUgdGhhdCB3ZSBhcmUgdHJhbnNwb3Npbmcgb3VyIGRhdGEgYXMgdGhlIGRlZmF1bHQgYmVoYXZpb3Igb2YgdGhpcyBmdW5jdGlvbiBpcyB0byBtYWtlIHBhaXJ3aXNlIGNvbXBhcmlzb25zIG9mIGFsbCAqcm93cyouDQoNCmBgYHtyfQ0KbGlicmFyeShwaGlsZW50cm9weSkNCg0KZGlzdGFuY2UodChkYXQpLCBtZXRob2QgPSAiY29zaW5lIikNCg0KbXlfZnVuKGRhdCwgdHlwZSA9ICJjb3NpbmUiKQ0KYGBgDQoNCg0KDQoNCg0KDQo=