superfan <- as.table(matrix(c(9, 12, 8, 1, 13, 1, 6, 20, 15, 4,
23, 18), ncol = 3))
attr(superfan, "dimnames") <- list(Band = c("Slayer",
"Iron Maiden", "Metallica", "Judas Priest"), Fan = c("Horst",
"Helga", "Klaus"))
superfan
              Fan
Band           Horst Helga Klaus
  Slayer           9    13    15
  Iron Maiden     12     1     4
  Metallica        8     6    23
  Judas Priest     1    20    18
fit_chisq <- chisq.test(superfan)
Chi-squared approximation may be incorrect
fit_chisq

    Pearson's Chi-squared test

data:  superfan
X-squared = 39.523, df = 6, p-value = 5.653e-07
S <- fit_chisq$residuals
round(S, 3)
              Fan
Band            Horst  Helga  Klaus
  Slayer        0.158  0.479 -0.503
  Iron Maiden   4.078 -1.850 -1.373
  Metallica    -0.184 -1.596  1.433
  Judas Priest -2.667  2.309  0.000
library("vcd")
package 㤼㸱vcd㤼㸲 was built under R version 4.0.5Loading required package: grid
mosaic(superfan, shade = TRUE)

plot(mosaic)
Error in dimnames(x) <- lapply(d, seq) : 'dimnames' applied to non-array

This plot shows that Horst and Iron Maiden shows the largest deviation from independence.

P <- prop.table(superfan)
round(P, 3)
              Fan
Band           Horst Helga Klaus
  Slayer       0.069 0.100 0.115
  Iron Maiden  0.092 0.008 0.031
  Metallica    0.062 0.046 0.177
  Judas Priest 0.008 0.154 0.138
r_mass <- margin.table(P, 1)
round(r_mass, 3)
Band
      Slayer  Iron Maiden    Metallica Judas Priest 
       0.285        0.131        0.285        0.300 
c_mass <- margin.table(P, 2)
round(c_mass, 3)
Fan
Horst Helga Klaus 
0.231 0.308 0.462 
r_profile <- prop.table(P, 1)
round(r_profile, 3)
              Fan
Band           Horst Helga Klaus
  Slayer       0.243 0.351 0.405
  Iron Maiden  0.706 0.059 0.235
  Metallica    0.216 0.162 0.622
  Judas Priest 0.026 0.513 0.462
c_profile <- prop.table(P, 2)
round(c_profile, 3)
              Fan
Band           Horst Helga Klaus
  Slayer       0.300 0.325 0.250
  Iron Maiden  0.400 0.025 0.067
  Metallica    0.267 0.150 0.383
  Judas Priest 0.033 0.500 0.300
ar_profile <- t(r_profile) %*% r_mass 
round(as.vector(ar_profile), 3)
[1] 0.231 0.308 0.462
round(as.vector(c_mass), 3)
[1] 0.231 0.308 0.462
ac_profile <- c_profile %*% c_mass 
round(as.vector(ac_profile), 3)
[1] 0.285 0.131 0.285 0.300
round(as.vector(r_mass), 3) 
[1] 0.285 0.131 0.285 0.300
library("plot3D")
package 㤼㸱plot3D㤼㸲 was built under R version 4.0.5
tc <- r_profile
scatter3D(x = tc[,1], y = tc[,2], z = tc[,3], xlab = "Horst", ylab = "Helga", zlab = "Klaus", colkey = FALSE, 
          col = 1, pch = 20, xlim = c(0,1), ylim = c(0,1), zlim = c(0,1), ticktype = "simple", type = "h", 
          phi = 40, theta = 50, main = "Row Profiles", bty = "g")
points3D(x = c(0,0,1), y = c(0,1,0), z = c(1,0,0), col = "red", add = TRUE, pch = 20)
lines3D(x = c(0,0,1,0), y = c(0,1,0,0), z = c(1,0,0,1), col = "red", add = TRUE)
text3D(x = tc[,1], y = tc[,2], z = tc[,3], labels = rownames(tc), pos = 3, add = TRUE, cex = 0.8, adj = -0.1) 


library("ggtern")
package 㤼㸱ggtern㤼㸲 was built under R version 4.0.5Loading required package: ggplot2
package 㤼㸱ggplot2㤼㸲 was built under R version 4.0.5--
Remember to cite, run citation(package = 'ggtern') for further info.
--

Attaching package: 㤼㸱ggtern㤼㸲

The following objects are masked from 㤼㸱package:ggplot2㤼㸲:

    aes, annotate, ggplot,
    ggplot_build, ggplot_gtable,
    ggplotGrob, ggsave, layer_data,
    theme_bw, theme_classic,
    theme_dark, theme_gray,
    theme_light, theme_linedraw,
    theme_minimal, theme_void
tf <- as.data.frame.matrix(superfan)
c_mass <- as.vector(c_mass)
lines <- data.frame(x = c(c_mass[1], 1-c_mass[3], 0),
                    y = c(1-c_mass[1], 0, c_mass[2]),
                    z = c(0, c_mass[3], 1-c_mass[2]),
                    xend = c(c_mass[1], c_mass[1], c_mass[1]),
                    yend = c(c_mass[2], c_mass[2], c_mass[2]),
                    zend = c(c_mass[3], c_mass[3], c_mass[3]), row.names = NULL)
gt <- ggtern(data = tf, aes(Horst, Helga, Klaus))
gt + geom_point() + theme_rgbw() + geom_text(label = rownames(tf), vjust = -0.5) +
geom_point(aes(x = c_mass[1], y = c_mass[2], z = c_mass[3]), colour = "red", size = 4) +
geom_segment(data = lines, aes(x = x, y = y, z = z, xend = xend, yend = yend, zend = zend), 
             color = 'red', size = 0.5) +
labs(title = "Ternary Plot")
Ignoring unknown aesthetics: zIgnoring unknown aesthetics: z, zend

The red dot represents the average row profile with the corresponding projections on the three fan axes. Here you can see that the Iron Maiden scores high on Horst’s axis, but low on Helga’s and Kalus’ axes, as it is there least favorite band.

#Interpretations & use of confidence ellipsoids!!

library("anacor")
package 㤼㸱anacor㤼㸲 was built under R version 4.0.5Registered S3 method overwritten by 'data.table':
  method           from
  print.data.table     
Registered S3 methods overwritten by 'car':
  method                          from
  influence.merMod                lme4
  cooks.distance.influence.merMod lme4
  dfbeta.influence.merMod         lme4
  dfbetas.influence.merMod        lme4
ca_fans <- anacor(superfan, ellipse = TRUE)
ca_fans

CA fit: 

Total chi-square value: 39.523 
Sum of eigenvalues (total inertia): 0.304 
Eigenvalues (principal inertias):
0.256 0.049 

Benzecri RMSE rows:  1.098991e-17 
Benzecri RMSE columns:  1.603251e-18 

Chi-square decomposition: 
             Chisq Proportion Cumulative Proportion
Dimension 1 33.216       0.84                  0.84
Dimension 2  6.307       0.16                  1.00

z-test on singular values: 
NA

#the first dimension explains 84% of the dispersion in the row/column profiles. Note - if it is significant, then the dimension is needed.

plot(ca_fans, main = "Symmetric CA Map")

Here it plots the principal row and column coordinates with confidence ellipsoids. Looking at the row coordinates, Slayer is close to the origin since their row progile is close to the average. Looking at the extremities, Iron Maiden is the most extreme band - as seen in the ternary plot. However, it is important to be careful not to compare row and column distances, as it is mapping two different plots onto the same one.

LS0tDQp0aXRsZTogIk1vZHVsZSAzOkNvcnJlc3BvbmRlbmNlIEFuYWx5c2lzIg0KYXV0aG9yOiBKYWtlIFJleW5vbGRzLCBGYWxsIDIwMjEgLSBJbmRlcGVuZGVudCBTdHVkeQ0Kb3V0cHV0OiBodG1sX25vdGVib29rDQotLS0NCg0KDQpgYGB7cn0NCnN1cGVyZmFuIDwtIGFzLnRhYmxlKG1hdHJpeChjKDksIDEyLCA4LCAxLCAxMywgMSwgNiwgMjAsIDE1LCA0LA0KMjMsIDE4KSwgbmNvbCA9IDMpKQ0KYXR0cihzdXBlcmZhbiwgImRpbW5hbWVzIikgPC0gbGlzdChCYW5kID0gYygiU2xheWVyIiwNCiJJcm9uIE1haWRlbiIsICJNZXRhbGxpY2EiLCAiSnVkYXMgUHJpZXN0IiksIEZhbiA9IGMoIkhvcnN0IiwNCiJIZWxnYSIsICJLbGF1cyIpKQ0Kc3VwZXJmYW4NCmBgYA0KYGBge3J9DQpmaXRfY2hpc3EgPC0gY2hpc3EudGVzdChzdXBlcmZhbikNCmZpdF9jaGlzcQ0KYGBgDQpgYGB7cn0NClMgPC0gZml0X2NoaXNxJHJlc2lkdWFscw0Kcm91bmQoUywgMykNCmBgYA0KYGBge3J9DQpsaWJyYXJ5KCJ2Y2QiKQ0KbW9zYWljKHN1cGVyZmFuLCBzaGFkZSA9IFRSVUUpDQpgYGANCiMgVGhpcyBwbG90IHNob3dzIHRoYXQgSG9yc3QgYW5kIElyb24gTWFpZGVuIHNob3dzIHRoZSBsYXJnZXN0IGRldmlhdGlvbiBmcm9tIGluZGVwZW5kZW5jZS4gDQpgYGB7cn0NClAgPC0gcHJvcC50YWJsZShzdXBlcmZhbikNCmBgYA0KDQoNCmBgYHtyfQ0Kcm91bmQoUCwgMykNCmBgYA0KYGBge3J9DQpyX21hc3MgPC0gbWFyZ2luLnRhYmxlKFAsIDEpDQpyb3VuZChyX21hc3MsIDMpDQpgYGANCmBgYHtyfQ0KY19tYXNzIDwtIG1hcmdpbi50YWJsZShQLCAyKQ0Kcm91bmQoY19tYXNzLCAzKQ0KYGBgDQoNCmBgYHtyfQ0Kcl9wcm9maWxlIDwtIHByb3AudGFibGUoUCwgMSkNCnJvdW5kKHJfcHJvZmlsZSwgMykNCmBgYA0KDQpgYGB7cn0NCmNfcHJvZmlsZSA8LSBwcm9wLnRhYmxlKFAsIDIpDQpyb3VuZChjX3Byb2ZpbGUsIDMpDQpgYGANCg0KYGBge3J9DQphcl9wcm9maWxlIDwtIHQocl9wcm9maWxlKSAlKiUgcl9tYXNzIA0Kcm91bmQoYXMudmVjdG9yKGFyX3Byb2ZpbGUpLCAzKQ0Kcm91bmQoYXMudmVjdG9yKGNfbWFzcyksIDMpDQphY19wcm9maWxlIDwtIGNfcHJvZmlsZSAlKiUgY19tYXNzIA0Kcm91bmQoYXMudmVjdG9yKGFjX3Byb2ZpbGUpLCAzKQ0Kcm91bmQoYXMudmVjdG9yKHJfbWFzcyksIDMpIA0KYGBgDQpgYGB7cn0NCmxpYnJhcnkoInBsb3QzRCIpDQp0YyA8LSByX3Byb2ZpbGUNCnNjYXR0ZXIzRCh4ID0gdGNbLDFdLCB5ID0gdGNbLDJdLCB6ID0gdGNbLDNdLCB4bGFiID0gIkhvcnN0IiwgeWxhYiA9ICJIZWxnYSIsIHpsYWIgPSAiS2xhdXMiLCBjb2xrZXkgPSBGQUxTRSwgDQogICAgICAgICAgY29sID0gMSwgcGNoID0gMjAsIHhsaW0gPSBjKDAsMSksIHlsaW0gPSBjKDAsMSksIHpsaW0gPSBjKDAsMSksIHRpY2t0eXBlID0gInNpbXBsZSIsIHR5cGUgPSAiaCIsIA0KICAgICAgICAgIHBoaSA9IDQwLCB0aGV0YSA9IDUwLCBtYWluID0gIlJvdyBQcm9maWxlcyIsIGJ0eSA9ICJnIikNCnBvaW50czNEKHggPSBjKDAsMCwxKSwgeSA9IGMoMCwxLDApLCB6ID0gYygxLDAsMCksIGNvbCA9ICJyZWQiLCBhZGQgPSBUUlVFLCBwY2ggPSAyMCkNCmxpbmVzM0QoeCA9IGMoMCwwLDEsMCksIHkgPSBjKDAsMSwwLDApLCB6ID0gYygxLDAsMCwxKSwgY29sID0gInJlZCIsIGFkZCA9IFRSVUUpDQp0ZXh0M0QoeCA9IHRjWywxXSwgeSA9IHRjWywyXSwgeiA9IHRjWywzXSwgbGFiZWxzID0gcm93bmFtZXModGMpLCBwb3MgPSAzLCBhZGQgPSBUUlVFLCBjZXggPSAwLjgsIGFkaiA9IC0wLjEpIA0KDQpsaWJyYXJ5KCJnZ3Rlcm4iKQ0KdGYgPC0gYXMuZGF0YS5mcmFtZS5tYXRyaXgoc3VwZXJmYW4pDQpjX21hc3MgPC0gYXMudmVjdG9yKGNfbWFzcykNCmxpbmVzIDwtIGRhdGEuZnJhbWUoeCA9IGMoY19tYXNzWzFdLCAxLWNfbWFzc1szXSwgMCksDQogICAgICAgICAgICAgICAgICAgIHkgPSBjKDEtY19tYXNzWzFdLCAwLCBjX21hc3NbMl0pLA0KICAgICAgICAgICAgICAgICAgICB6ID0gYygwLCBjX21hc3NbM10sIDEtY19tYXNzWzJdKSwNCiAgICAgICAgICAgICAgICAgICAgeGVuZCA9IGMoY19tYXNzWzFdLCBjX21hc3NbMV0sIGNfbWFzc1sxXSksDQogICAgICAgICAgICAgICAgICAgIHllbmQgPSBjKGNfbWFzc1syXSwgY19tYXNzWzJdLCBjX21hc3NbMl0pLA0KICAgICAgICAgICAgICAgICAgICB6ZW5kID0gYyhjX21hc3NbM10sIGNfbWFzc1szXSwgY19tYXNzWzNdKSwgcm93Lm5hbWVzID0gTlVMTCkNCmd0IDwtIGdndGVybihkYXRhID0gdGYsIGFlcyhIb3JzdCwgSGVsZ2EsIEtsYXVzKSkNCmd0ICsgZ2VvbV9wb2ludCgpICsgdGhlbWVfcmdidygpICsgZ2VvbV90ZXh0KGxhYmVsID0gcm93bmFtZXModGYpLCB2anVzdCA9IC0wLjUpICsNCmdlb21fcG9pbnQoYWVzKHggPSBjX21hc3NbMV0sIHkgPSBjX21hc3NbMl0sIHogPSBjX21hc3NbM10pLCBjb2xvdXIgPSAicmVkIiwgc2l6ZSA9IDQpICsNCmdlb21fc2VnbWVudChkYXRhID0gbGluZXMsIGFlcyh4ID0geCwgeSA9IHksIHogPSB6LCB4ZW5kID0geGVuZCwgeWVuZCA9IHllbmQsIHplbmQgPSB6ZW5kKSwgDQogICAgICAgICAgICAgY29sb3IgPSAncmVkJywgc2l6ZSA9IDAuNSkgKw0KbGFicyh0aXRsZSA9ICJUZXJuYXJ5IFBsb3QiKQ0KYGBgDQoNCiMgVGhlIHJlZCBkb3QgcmVwcmVzZW50cyB0aGUgYXZlcmFnZSByb3cgcHJvZmlsZSB3aXRoIHRoZSBjb3JyZXNwb25kaW5nIHByb2plY3Rpb25zIG9uIHRoZSB0aHJlZSBmYW4gYXhlcy4gSGVyZSB5b3UgY2FuIHNlZSB0aGF0IHRoZSBJcm9uIE1haWRlbiBzY29yZXMgaGlnaCBvbiBIb3JzdCdzIGF4aXMsIGJ1dCBsb3cgb24gSGVsZ2EncyBhbmQgS2FsdXMnIGF4ZXMsIGFzIGl0IGlzIHRoZXJlIGxlYXN0IGZhdm9yaXRlIGJhbmQuIA0KDQoNCg0KI0ludGVycHJldGF0aW9ucyAmIHVzZSBvZiBjb25maWRlbmNlIGVsbGlwc29pZHMhISANCg0KYGBge3J9DQpsaWJyYXJ5KCJhbmFjb3IiKQ0KY2FfZmFucyA8LSBhbmFjb3Ioc3VwZXJmYW4sIGVsbGlwc2UgPSBUUlVFKQ0KY2FfZmFucw0KYGBgDQojdGhlIGZpcnN0IGRpbWVuc2lvbiBleHBsYWlucyA4NCUgb2YgdGhlIGRpc3BlcnNpb24gaW4gdGhlIHJvdy9jb2x1bW4gcHJvZmlsZXMuIE5vdGUgLSBpZiBpdCBpcyBzaWduaWZpY2FudCwgdGhlbiB0aGUgZGltZW5zaW9uIGlzIG5lZWRlZC4NCg0KYGBge3J9DQpwbG90KGNhX2ZhbnMsIG1haW4gPSAiU3ltbWV0cmljIENBIE1hcCIpDQpgYGANCiMgSGVyZSBpdCBwbG90cyB0aGUgcHJpbmNpcGFsIHJvdyBhbmQgY29sdW1uIGNvb3JkaW5hdGVzIHdpdGggY29uZmlkZW5jZSBlbGxpcHNvaWRzLiBMb29raW5nIGF0IHRoZSByb3cgY29vcmRpbmF0ZXMsIFNsYXllciBpcyBjbG9zZSB0byB0aGUgb3JpZ2luIHNpbmNlIHRoZWlyIHJvdyBwcm9naWxlIGlzIGNsb3NlIHRvIHRoZSBhdmVyYWdlLiBMb29raW5nIGF0IHRoZSBleHRyZW1pdGllcywgSXJvbiBNYWlkZW4gaXMgdGhlIG1vc3QgZXh0cmVtZSBiYW5kIC0gYXMgc2VlbiBpbiB0aGUgdGVybmFyeSBwbG90LiBIb3dldmVyLCBpdCBpcyBpbXBvcnRhbnQgdG8gYmUgY2FyZWZ1bCBub3QgdG8gY29tcGFyZSByb3cgYW5kIGNvbHVtbiBkaXN0YW5jZXMsIGFzIGl0IGlzIG1hcHBpbmcgdHdvIGRpZmZlcmVudCBwbG90cyBvbnRvIHRoZSBzYW1lIG9uZS4gDQo=