Skript mit Erklärung der Verfahren und Vorgehensweisen anhand von Beispielen.

MDS.pdf (10 MB)

Für eine und mehrere Unähnlichkeitsmatrizen

Datenstruktur: Die hier dargestellten Verfahren setzen voraus, dass die Daten in Form einer symmetrischen Distanzmatrix vorliegen. Die Einschränkung auf symmetrische Distanzmatrizen hat zur Folge, dass mit der Ankerpunktmethode erhobene asymmetrische zeilenkonditionale Unähnlichkeitsmatrizen nicht analysisert werden können.
Dateneingabe: Um Arbeit bei der Dateneingabe zu sparen, kann man z.B, in Excel eine “Dreiecksmatrix” eingeben, welche in der unteren Hälfte die Unähnlichkeiten und sonst Nullen enthält.

Erforderliche Pakete laden

library(smacof)     # MDS
## Loading required package: plotrix
## Loading required package: colorspace
## Loading required package: e1071
## 
## Attaching package: 'smacof'
## The following object is masked from 'package:base':
## 
##     transform
library(ggplot2)    # Plot
library(gridExtra)  # Plot
library(plotrix)    # Plot

Eine Unähnlichkeitsmatrix

Datensatz einlesen und Variablen spezifizieren

# Dateneingabe
data(ekman)
data <- sim2diss(ekman)

# Variablen spezifizieren
Skalenniveau <- "ordinal" # "ratio", "interval", "ordinal" oder "mspline"
k <- 2                    # Anzahl Dimensionen eingeben
ties <- "primary"         # verbundene Ränge: "primary" break ties, "secondary" keep ties

MDS

library(smacof)
d <- as.dist(as.matrix(data))
res <- smacofSym(d, ndim=k, type=Skalenniveau, ties=ties)
# Stress
# Der vom Programm gelieferte Stress entspricht dem Stress nach Kruskal 1 und wird mit
# sqrt(sum((res$confdiss-res$dhat)^2)/sum(res$dhat^2)) berechnet.
# quadrierte Korrelation zwischen Disparit??ten und Distanzen
RSQ <- cor(res$delta,res$dhat)^2

# Permutationstest
# p-value = Wahrscheinlichkeit für Stresswerte kleiner als Stress der obigen MDS
res.perm <- permtest(res, verbose=FALSE)

# Stabilitätsmass
res.jk <- jackmds(res)

# Mittelwert des Stress von 100 Zufallskonfigurationen
Stress.Zufallskonf <- mean(randomstress(n = res$nobj, ndim = k, nrep = 100, type = Skalenniveau))

list("MDS-Output"=res, "quadrierte Korrelation zwischen Disparitäten und Distanzen"=RSQ, Stabilität=res.jk, Permutationstest=res.perm, "Mittelwert des Stress von 100 Zufallskonfigurationen"=Stress.Zufallskonf)
## $`MDS-Output`
## 
## Call:
## smacofSym(delta = d, ndim = k, type = Skalenniveau, ties = ties)
## 
## Model: Symmetric SMACOF 
## Number of objects: 14 
## Stress-1 value: 0.023 
## Number of iterations: 27 
## 
## 
## $`quadrierte Korrelation zwischen Disparitäten und Distanzen`
## [1] 0.8981498
## 
## $Stabilität
## 
## Call: jackmds.smacofB(object = res)
## 
## SMACOF Jackknife
## Number of objects: 14 
## Value loss function: 0.0487 
## Number of iterations: 5 
## 
## Stability measure: 0.9996 
## Cross validity: 1 
## Dispersion: 4e-04 
## 
## 
## $Permutationstest
## 
## Call: permtest.smacof(object = res, verbose = FALSE)
## 
## SMACOF Permutation Test
## Number of objects: 14 
## Number of replications (permutations): 100 
## 
## Observed stress value: 0.023 
## p-value: <0.001 
## 
## 
## $`Mittelwert des Stress von 100 Zufallskonfigurationen`
## [1] 0.2488249

Weitere Output-Bestandteile:
Die MDS-Koordinaten können mit res$conf abgerufen werden.
Die MDS-Distanzen können mit res[[3]] abgerufen werden.
Die Proximitäten können mit res$dhat abgerufen werden.
Weitere Outputmöglichkeiten finden Sie, wenn Sie ?smacofSym eingeben unter “Value”.

Grafiken

# Konfigurationsplot
plot(res, plot.type = "confplot")
abline(h=0, lty=2); abline(v=0, lty=2)

# Diagnostikplots
par(mfrow=c(2, 2))
# Jacknife Plot
plot(res.jk, hclpar = list(c = 0, l = 0))
# Shepard-Diagramm
plot(res, plot.type = "Shepard")
# Rsiduenplot
plot(res, plot.type = "resplot", main="Residual Plot")
# Stressplot: Aufteilung des Stress auf die Objekte
plot(res, plot.type = "stressplot")


Hinweis 1: Wenn Sie eine sphärische Punktekonfiguration erwarten, können Sie einen Kreis in Ihren 2D-Konfigurationsplot fitten.

# plot(res, plot.type = "confplot")
# abline(h=0, lty=2); abline(v=0, lty=2)
# circle <- fitCircle(res$conf[,1], res$conf[,2])
# draw.circle(circle$cx, circle$cy, radius = circle$radius, border = "gray")


Hinweis 2: Wenn Sie eine sphärische Punktekonfiguration erwarten, können Sie eine konfirmatorische MDS rechnen, d.h. die Punkte der MDS-Konfiguration liegen exakt auf einem Kreis (2D) oder einer Kugel (3D).

# res.circle <- smacofSphere(d, ndim=k, type=Skalenniveau, ties=ties)
# res.circle
# plot(res.circle, plot.type = "confplot")
# abline(h=0, lty=2); abline(v=0, lty=2)


Hinweis 3: Wenn Sie externe Skalen (Eigenschaften oder Präferenzen) als Vektoren in Ihre MDS-L??sung fitten wollen, müssen Sie diese Skalen als n x m - Matrix laden: n = Anzahl Objekte, m = Anzahl Skalen.

# res.biplot <- biplotmds(res, externe.variablen)
# plot(res.biplot)
# abline(h=0, lty=2); abline(v=0, lty=2)
# Vektor mit den R2-Werten. Niedrige R2-Werte deuten auf eine schlechte Vorhersagbarkeit hin.
# res.biplot$R2vec 


Hinweis 4: Vergleich zweier MDS-Lösungen. Zuerst wird die zweite Konfiguration (Y) durch Translation, Drehung und Streckung der ersten (X) möglichst ähnlich gemacht. Anschliessend wird die Ähnlichkeit beider Lösungen mit Hilfe des Kongruenzkoeffizients ausgedrückt.

# res.proc <- Procrustes(X=res1$conf, Y=res2$conf)
# res.proc
# plot(res.proc, plot.type="transplot", col.X = "red", col.Y = "blue")
# abline(h=0, lty=2); abline(v=0, lty=2)

Mehrere Unähnlichkeitsmatrizen

Datensatz einlesen und Variablen spezifizieren

Datenstruktur: Wir gehen davon aus, dass die Unähnlichkeitsmatrizen direkt untereinader eingegeben werden. Dieser Datensatz wird anschliessend in die benötigte Form gebracht.

Bedeutung der Constraints:
Mit constraint=“identity” müssen die matrixspezifischen Lösungen mit der Gruppenlösung übereinstimmen.
Mit constraint=“indscal” wird jede matrixspezifische Lösung als eine pro Dimension individuell gestreckte oder gestauchte Variante einer Gruppenlösung betrachtet.
Mit constraint=“idioscal” wird jede matrixspezifische Lösung als rotierte und gestreckte oder gestauchte Variante der Gesamtlösung betrachtet.

# Dateneingabe
data<-read.table(file="https://mmi.psycho.unibas.ch/r-toolbox/data/INDSCAL.txt", header=TRUE)

# Variablen spezifizieren
Skalenniveau <- "ordinal" # "ratio", "interval" oder "ordinal"
constraint <- "indscal"   # "indscal", "idioscal" oder "identity"
n.mat <- 16               # Anzahl untereinanderstehender Matrizen eingeben
k <- 2                    # Anzahl Dimensionen eingeben
ties <- "primary"         # verbundene Ränge: "primary" break ties, "secondary" keep ties

MDS

# Umwandeln des Datensatzes mit untereinderstehenden Matrizen in eine Liste von Matrizen
d <- vector("list", length=n.mat)
for (i in 1:n.mat) {
matr <- as.matrix(data[((i-1)*length(data)+1):(i*length(data)),])
rownames(matr) <- names(data)
d[[i]] <- as.dist(matr)}

# MDS
library(smacof)
res <- smacofIndDiff(d, ndim=k, constraint=constraint, type=Skalenniveau, ties=ties)

# Berechnung des Stress (Kruskal 1) und RSQ pro Matrix
Stress.pro.Matrix <- vector(length=n.mat)
RSQ.pro.Matrix <- vector(length=n.mat)
for (i in 1:n.mat) {
Stress.pro.Matrix[i] <- sqrt(sum((res[[3]][i][[1]]-res$dhat[[i]])^2)/sum(res$dhat[[i]]^2))
RSQ.pro.Matrix[i] <- cor(res$delta[[i]],res$dhat[[i]])^2
}
stress.rsq <- data.frame(Stress.pro.Matrix, RSQ.pro.Matrix)
rownames(stress.rsq) <- paste("M",1:n.mat,sep="")
colnames(stress.rsq) <- c("Stress 1", "RSQ")

# Falls INDSCAL gerechnet wurde, Tabelle mit den Gewichtsfaktoren für die Dimensionen erstellen
if (constraint == "indscal") {
   Faktoren <- matrix(ncol=k, nrow=n.mat)
   for (i in 1:n.mat) Faktoren[i,] <- diag(array(unlist(res$cweights), dim = c(nrow(res$cweights[[1]]), ncol(res$cweights[[1]]), length(res$cweights)))[,,i])
colnames(Faktoren) <- paste("D",1:k,sep=""); rownames(Faktoren) <- paste("M",1:n.mat,sep="")
}

if (constraint == "idioscal"){
list("MDS-Output"=res, "Stress (Kruskal 1) und RSQ pro Matrix"= stress.rsq)
} else {
   list("MDS-Output"=res, "Stress (Kruskal 1) und RSQ pro Matrix"= stress.rsq, Streckungsfaktoren=Faktoren)
}
## $`MDS-Output`
## 
## Call: smacofIndDiff(delta = d, ndim = k, type = Skalenniveau, constraint = constraint, 
##     ties = ties)
## 
## Model: Three-way SMACOF 
## Number of objects: 10 
## Stress-1 value: 0.086 
## Number of iterations: 229 
## 
## 
## $`Stress (Kruskal 1) und RSQ pro Matrix`
##       Stress 1       RSQ
## M1  0.03879643 0.9687276
## M2  0.05582684 0.9460614
## M3  0.06861372 0.9783981
## M4  0.05661686 0.9844709
## M5  0.12902000 0.9463517
## M6  0.06485966 0.9621078
## M7  0.06842392 0.9589968
## M8  0.10636993 0.9646128
## M9  0.05905280 0.9712914
## M10 0.08006403 0.9694182
## M11 0.08847646 0.9760974
## M12 0.09993405 0.9794522
## M13 0.04938955 0.9863941
## M14 0.09661412 0.9632224
## M15 0.09804244 0.9718232
## M16 0.14048467 0.9094536
## 
## $Streckungsfaktoren
##            D1        D2
## M1  0.9904678 1.0453607
## M2  0.8976326 1.1649039
## M3  0.9352340 1.1168441
## M4  0.9470295 1.1032860
## M5  1.0889909 0.8601671
## M6  0.9524200 1.0949721
## M7  0.8638650 1.2016757
## M8  1.0487493 0.9409445
## M9  1.2614702 0.3913095
## M10 0.9624682 1.0787905
## M11 0.9678666 1.0696543
## M12 0.9145211 1.1368119
## M13 0.8598956 1.2084367
## M14 1.0847193 0.8791954
## M15 0.9968371 1.0258582
## M16 1.2258462 0.4968163

Weitere Output-Bestandteile:
Die MDS-Koordinaten der Gruppenlösung können mit res$gspace abgerufen werden.
Die MDS-Koordinaten der einzelnen Matrizen können mit res$conf abgerufen werden.
Die MDS-Distanzen der einzelnen Matrizen können mit res[[3]] abgerufen werden.
Die Proximitäten der einzelnen Matrizen können mit res$dhat abgerufen werden.
Die Konfigurationsgewichte (configuration weights) können mit res$cweights abfgerufen werden.
Weitere Outputmöglichkeiten finden Sie, wenn Sie ?smacofIndDiff eingeben unter “Value”.

Grafiken

Das Shepard-Diagramm (plot.type=“Shepard”) und der Residuenplot (plot.type=“resplot”) werden nicht ausgegeben. Deshalb erstellen wir diese Plots aus dem Ouput der MDS.
Da die Matrizen der normalisierten Dissimilaritäten (res$obsdiss) leer sind, verwenden wir die Rohwerte der Dissimilaritäten (res$delta).

# Konfigurationsplot
plot(res, plot.type = "confplot")
abline(h=0, lty=2); abline(v=0, lty=2)

# Falls INDSCAL gewählt wurde: Plot mit den Gewichtsfaktoren für die ersten beiden Diemensionen
library(ggplot2)
if (constraint == "indscal") {
  ggplot(data=data.frame(Faktoren)) + geom_segment(mapping=aes(x=0,y=0,xend=D1,yend=D2)) + geom_text(aes(x=D1,y=D2, label=1:n.mat), size=4, hjust=0, vjust=0) + xlab("D1") + ylab("D2") + xlim(0,max(Faktoren[,1:2])) + ylim(0,max(Faktoren[,1:2])) + ggtitle("Gewichtsfaktoren") + theme_bw()
}

# Diagnostikplots
# Matrixspezifische Shepard-Diagramme und Residuenplots
library(ggplot2); library(gridExtra); library(plotrix)
pltList.Shepard <- list()
pltList.Residuen <- list()
for (i in 1:n.mat) {
    dat <- data.frame(Dissimilarities=as.matrix(res$delta[[i]])[lower.tri(as.matrix(res$delta[[i]]))], Distances=as.matrix(res[[3]][i][[1]])[lower.tri(as.matrix(res[[3]][i][[1]]))], Disparities=as.matrix(res$dhat[[i]])[lower.tri(as.matrix(res$dhat[[i]]))])
    pltList.Shepard[[i]] <- ggplot(dat, aes(x=Dissimilarities, y=Distances)) + geom_point(size=1) + geom_line(aes(x=Dissimilarities, y=Disparities), colour="grey") + theme(axis.title.x=element_blank(), axis.title.y=element_blank()) + theme_bw()
    pltList.Residuen[[i]] <- ggplot(dat, aes(x=Disparities, y=Distances)) + geom_point(size=1) + theme(axis.title.x=element_blank(), axis.title.y=element_blank()) + geom_smooth(method=lm, se=FALSE, col=1, size=1) + theme_bw()
}
# Matrixspezifische Shepard-Diagramme: Nummerierung zeilenweise
grid.arrange(grobs=pltList.Shepard, top="Separd Diagram: x=Dissimilarities, y=Distances and Disparities (grey line)")

# Matrixspezifische Residuenplots: Nummerierung zeilenweise
grid.arrange(grobs=pltList.Residuen, top="Residual Plot: x=Disparities, y=Distances")
## `geom_smooth()` using formula 'y ~ x'
## `geom_smooth()` using formula 'y ~ x'
## `geom_smooth()` using formula 'y ~ x'
## `geom_smooth()` using formula 'y ~ x'
## `geom_smooth()` using formula 'y ~ x'
## `geom_smooth()` using formula 'y ~ x'
## `geom_smooth()` using formula 'y ~ x'
## `geom_smooth()` using formula 'y ~ x'
## `geom_smooth()` using formula 'y ~ x'
## `geom_smooth()` using formula 'y ~ x'
## `geom_smooth()` using formula 'y ~ x'
## `geom_smooth()` using formula 'y ~ x'
## `geom_smooth()` using formula 'y ~ x'
## `geom_smooth()` using formula 'y ~ x'
## `geom_smooth()` using formula 'y ~ x'
## `geom_smooth()` using formula 'y ~ x'

# Stressplot: Aufteilung des Stress auf die Objekte
plot(res, plot.type = "stressplot")


Hinweis 1: Wenn Sie eine sphärische Punktekonfiguration erwarten, können Sie einen Kreis in Ihren 2D-Konfigurationsplot fitten.

# plot(res, plot.type = "confplot")
# abline(h=0, lty=2); abline(v=0, lty=2)
# circle <- fitCircle(res$gspace[,1], res$gspace[,2])
# draw.circle(circle$cx, circle$cy, radius = circle$radius, border = "gray")


Hinweis 2: Wenn Sie externe Skalen (Eigenschaften oder Präferenzen) als Vektoren in Ihre MDS-Lösung fitten wollen, müssen Sie diese Skalen als n x m - Matrix laden: n = Anzahl Objekte, m = Anzahl Skalen.

# res.biplot <- biplotmds(res, externe.variablen)
# plot(res.biplot)
# abline(h=0, lty=2); abline(v=0, lty=2)
# Vektor mit den R2-Werten. Niedrige R2-Werte deuten auf eine schlechte Vorhersagbarkeit hin.
# res.biplot$R2vec 


Hinweis 3: Vergleich zweier MDS-Lösungen. Zuerst wird die zweite Konfiguration (Y) durch Translation, Drehung und Streckung der ersten (X) möglichst ähnlich gemacht. Anschliessend wird die Änlichkeit beider Lösungen mit Hilfe des Kongruenzkoeffizients ausgedrückt.

# res.proc <- Procrustes(X=res1$conf, Y=res2$conf)
# res.proc
# plot(fit.proc, plot.type="transplot", col.X = "red", col.Y = "blue")
# abline(h=0, lty=2); abline(v=0, lty=2)