Module 6: Biplots and Networks

#Biplots

library("MPsychoR")
data("BrainIQ")
BrainIQ <- na.omit(BrainIQ[,-1]) ## we omit NAs and gender
rownames(BrainIQ) <- 1:nrow(BrainIQ) ## relabel persons
BrainIQ1 <- as.data.frame(scale(BrainIQ)) ## standardize
regfit <- lm(cbind(VIQ, PIQ, Weight) ~
-1 + Height + MRI_Count, data = BrainIQ1)
colnames(regfit$coef) <- c("VIQ", "PIQ", "Weight")
round(regfit$coef, 3) ## vector coordinates

            VIQ   PIQ Weight
Height    -0.45 -0.48   0.61
MRI_Count  0.57  0.66   0.16

## VIQ PIQ Weight
## Height -0.452 -0.482 0.608
## MRI_Count 0.566 0.662 0.156
R2vec <- sapply(summary(regfit), `[[`, "r.squared")

round(cor(BrainIQ[, 2:4]), 3)

          VIQ   PIQ Weight
VIQ     1.000 0.776 -0.076
PIQ     0.776 1.000  0.003
Weight -0.076 0.003  1.000

#Here are the orthogonal projections.

library("calibrate")

package 㤼㸱calibrate㤼㸲 was built under R version 4.0.5Loading required package: MASS

Attaching package: 㤼㸱MASS㤼㸲

The following object is masked from 㤼㸱package:dplyr㤼㸲:

    select

plot(BrainIQ1$Height, BrainIQ1$MRI_Count, pch = 20, cex = 0.8,
xlab = "Height", ylab = "MRI", col = "darkblue", asp = 1,
main = "Orthogonal Projections")
text(BrainIQ1$Height, BrainIQ1$MRI_Count, labels = 1:nrow(BrainIQ1),
cex = 0.7, pos = 3, col = "darkblue")

abline(h = 0, v = 0, lty = 2, col = "darkgray")
calibrate.Z <- calibrate(regfit$coef[,1], BrainIQ1$VIQ, seq(-2,2, by = 0.5),
cbind(BrainIQ1$Height, BrainIQ1$MRI_Count), dp = TRUE, axiscol = "brown",
axislab = "VIQ", labpos = 3, verb = FALSE)

VIQcal <- calibrate.Z$yt*sd(BrainIQ$VIQ) + mean(BrainIQ$VIQ)

summary(BrainIQ[, 2]) 

   Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
     71      90     113     112     129     150 

BrainIQ[c(20, 11), -c(1, 4)] 

library("MPsychoR")

package 㤼㸱MPsychoR㤼㸲 was built under R version 4.0.3

data("BrainIQ")
BrainIQ1 <- na.omit(BrainIQ[,-1])
head(BrainIQ1, 3)

pca_biq1 <- prcomp(BrainIQ1)
pca_biq2 <- prcomp(BrainIQ1, scale = TRUE)

library(psych)
biplot(pca_biq2)

Here is where I am focusing on multidimension scaling biplots

library("MPsychoR")
library("smacof")

package 㤼㸱smacof㤼㸲 was built under R version 4.0.5Loading required package: plotrix
package 㤼㸱plotrix㤼㸲 was built under R version 4.0.5
Attaching package: 㤼㸱plotrix㤼㸲

The following object is masked from 㤼㸱package:psych㤼㸲:

    rescale

Loading required package: colorspace
package 㤼㸱colorspace㤼㸲 was built under R version 4.0.3Loading required package: e1071
package 㤼㸱e1071㤼㸲 was built under R version 4.0.3Registered S3 methods overwritten by 'car':
  method                          from
  influence.merMod                lme4
  cooks.distance.influence.merMod lme4
  dfbeta.influence.merMod         lme4
  dfbetas.influence.merMod        lme4

Attaching package: 㤼㸱smacof㤼㸲

The following object is masked from 㤼㸱package:psych㤼㸲:

    Procrustes

The following object is masked from 㤼㸱package:base㤼㸲:

    transform

data("NeuralActivity")
delta <- Reduce("+", NeuralActivity)/length(NeuralActivity)
fit_neural <- mds(delta, type = "interval")
fit_neural

Call:
mds(delta = delta, type = "interval")

Model: Symmetric SMACOF 
Number of objects: 60 
Stress-1 value: 0.28 
Number of iterations: 68 

data("NeuralScales")
mdsbi <- biplotmds(fit_neural, NeuralScales)

#Here is a multidimensional scaling biplot

plot(mdsbi)

X <- fit_neural$conf
Y <- scale(NeuralScales, scale = TRUE)
plot(X, pch = 20, cex = 0.6, xlab = "Dimension 1", ylab = "Dimension 2",
col = "darkblue", asp = 1, main = "Biplot MDS Emotion Axis")
text(X, labels = rownames(X), cex = 0.7, pos = 3, col = "darkblue")

abline(h = 0, v = 0, lty = 2, col = "gray")
calEm <- calibrate(mdsbi$coef[,"Emotion"], Y[,"Emotion"],
tm = seq(-2, 1.5, by = 0.5), Fr = X, dp = TRUE,
axiscol = "brown", axislab = "Emotion", labpos = 3, verb = FALSE)

#Below is a corresponding anlaysis biplot.

library("anacor")

package 㤼㸱anacor㤼㸲 was built under R version 4.0.5

superfan <- as.table(matrix(c(9, 12, 8, 1, 13, 1, 6, 20, 15, 4, 23, 18),
ncol = 3))
attr(superfan, "dimnames") <- list(c("Slayer", "Iron Maiden", "Metallica",
"Judas Priest"), c("Horst", "Helga", "Klaus"))
fit_fans <- anacor(superfan, scaling = c("Benzecri", "standard"))
plot(fit_fans, main = "Asymmetric Superfan CA Map")

Srole <- as.table(matrix(c(64, 94, 58, 46,
57, 94, 54, 40,
57, 105, 65, 60,
72, 141, 77, 94,
36, 97, 54, 78,
21, 71, 54, 71), nrow = 4))
attr(Srole, "dimnames") <- list(mhealth = c("well", "mild",
"moderate", "impaired"), ses = LETTERS[1:6])
Srole

          ses
mhealth      A   B   C   D   E   F
  well      64  57  57  72  36  21
  mild      94  94 105 141  97  71
  moderate  58  54  65  77  54  54
  impaired  46  40  60  94  78  71

fit_ses <- anacor(Srole, scaling = c("standard", "Benzecri"))
plot(fit_ses, arrows = c(T, F), main = "Asymmetric CA Map")

#Networks

library("igraph")
album_df <- matrix(c("Helga", "Gertrud", "Horst", "Klaus",
"Horst", "Helga", "Gertrud", "Horst", "Klaus", "Horst",
"Horst", "Helga", "Helga", "Klaus", "Klaus", "Klaus",
"Gertrud", "Gertrud"), ncol = 2)
album_el <- graph.edgelist(album_df, directed = TRUE)
E(album_el)$weight = c(2, 3, 3, 4, 5, 1, 7, 7, 2)

album_ad <- get.adjacency(album_el, sparse = FALSE,
attr = "weight")
album_ad

        Helga Horst Gertrud Klaus
Helga       0     2       0     1
Horst       3     0       7     5
Gertrud     0     3       0     7
Klaus       4     0       2     0

#Below is how you make a network graph.

plot(album_el, vertex.size = 0, edge.arrow.size = 0.5,
vertex.label.dist = 0.8, vertex.color = "black",
vertex.label.color = "black", vertex.label.cex = 1.5)

library("MPsychoR")
data(Rogers)
cormat <- cor(Rogers)

#Below you can set the lower absolute correlation threshold. The thicker an edge, the higher the correlation. Negative correlations are colored differently.

library("qgraph")

package 㤼㸱qgraph㤼㸲 was built under R version 4.0.4

cornet <- qgraph(cormat, layout = "spring", minimum = 0.2,
graph = "cor", groups = list(Depression = 1:16, OCD = 17:26),
color = c("white", "gray"), labels = colnames(Rogers),
title = "Depression/OCD Correlation Network")

centralityPlot(cornet)

Note: z-scores are shown on x-axis rather than raw centrality indices.

library("eigenmodel")

package 㤼㸱eigenmodel㤼㸲 was built under R version 4.0.3

diag(cormat) <- NA ## NA diagonals required
fitEM <- eigenmodel_mcmc(cormat, R = 2, S = 1000,
burn = 200, seed = 123)

5 percent done (iteration 60) Tue Dec 14 15:57:54 2021
10 percent done (iteration 120) Tue Dec 14 15:57:55 2021
15 percent done (iteration 180) Tue Dec 14 15:57:56 2021
20 percent done (iteration 240) Tue Dec 14 15:57:57 2021
25 percent done (iteration 300) Tue Dec 14 15:57:58 2021
30 percent done (iteration 360) Tue Dec 14 15:57:59 2021
35 percent done (iteration 420) Tue Dec 14 15:58:01 2021
40 percent done (iteration 480) Tue Dec 14 15:58:02 2021
45 percent done (iteration 540) Tue Dec 14 15:58:03 2021
50 percent done (iteration 600) Tue Dec 14 15:58:04 2021
55 percent done (iteration 660) Tue Dec 14 15:58:05 2021
60 percent done (iteration 720) Tue Dec 14 15:58:06 2021
65 percent done (iteration 780) Tue Dec 14 15:58:08 2021
70 percent done (iteration 840) Tue Dec 14 15:58:09 2021
75 percent done (iteration 900) Tue Dec 14 15:58:10 2021
80 percent done (iteration 960) Tue Dec 14 15:58:11 2021
85 percent done (iteration 1020) Tue Dec 14 15:58:12 2021
90 percent done (iteration 1080) Tue Dec 14 15:58:14 2021
95 percent done (iteration 1140) Tue Dec 14 15:58:15 2021
100 percent done (iteration 1200) Tue Dec 14 15:58:16 2021

Eigenvalues

evals <- colMeans(fitEM$L_postsamp)
evals

[1] 0.30 0.66

#Below you can see a clear separation between OCD and depression itemsets along dimension 2. This shows that second dimension is more important one. Also, it shows that dimension 1 only discriminates between the “core” depression items and the “periphery” depression items.

evecs <- eigen(fitEM$ULU_postmean)$vec[, 1:2]
cols <- c("coral", "cadetblue")
plot(evecs, type = "n", xlab = "Dimension 1",
ylab = "Dimension 2", xlim = c(-0.30, 0),
main = "Depression/OCD Eigenmodel Network")
corthresh <- 0.2 ## correlation threshold
addlines(evecs, abs(cormat) > corthresh, col = "gray")

ind <- c(rep(1, 16), rep(2, 10))
text(evecs, labels = rownames(cormat), col = cols[ind],
cex = 0.8)
legend("topright", legend = c("Depression", "OCD"),
col = cols, pch = 19)

#Next, lets assume that we do not know which variables belong to which latent model. Here the algorithm will try to determine.

thresh <- 0.2
cormat01 <- ifelse(abs(cormat) > thresh, 1, 0)

library("network")
cornet <- network(cormat01, matrix.type = "adjacency",
directed = FALSE)

library(ergm)
library("latentnet")

set.seed(111)
fitLN1 <- ergmm(cornet ~ euclidean(d = 2, G = 1))
summary(fitLN1)$bic$Z

NOTE: It is not certain whether it is appropriate to use latentnet's BIC to select latent space dimension, whether or not to include actor-specific random effects, and to compare clustered models with the unclustered model.

[1] 248

fitLN2 <- ergmm(cornet ~ euclidean(d = 2, G = 2))
summary(fitLN2)$bic$Z

[1] 236

fitLN3 <- ergmm(cornet ~ euclidean(d = 2, G = 3))
summary(fitLN3)$bic$Z

[1] 242

#Here clusters do not overlap

plot(fitLN2)

clusmemb2 <- fitLN2$mkl$mbc$Z.pZK
dimnames(clusmemb2) <- list(colnames(cormat01),
paste("Cluster", 1:2))
clusmemb2[c("comptime", "suicide", "weightgain"), ]

           Cluster 1 Cluster 2
comptime       0.034    0.9660
suicide        0.963    0.0370
weightgain     1.000    0.0005

#Below you can see that clusters overlap.

plot(fitLN3)

#Here is a cool density map, which gives additional information on the cluster separation.

plot(fitLN3, main = "Latent Class Network",
cluster.col = c("coral", "cadetblue", "darkgoldenrod"),
what = "density")