Vijayan Padmanabhan
2011-Jan-11 06:45 UTC
[R] modified FAST Script from package SensoMineR for the R community - Reg
###Dear R users
###I have been using SensoMineR package from CRAN for most of my work in
sensory data analysis and from my usage experience, I encountered some
areas for improvement and considered ###modifying the function in
SensoMineR package for my personal use. I felt that it could be useful to
share this to the community for enabling adoption by other users where
they might require a similar ###functionality.
###Thanks to Dr.Husson and Dr.Le for making available the above package in
the first place.
###Basic utility delivered through the modified fast function: (Refer the
function "fastmod" below:)
##Plotting the most relevant 6 plots from "fast" function of
SensoMiner
Package in a single layout (Very useful when working with Sweave or
odfWeave, wherein the figure returned as result from the ##function
execution can't be more than one per code chunk). The nearest workaround
would mean saving each graphical output as a image and reinserting the
stored image of plots through multiple ##code chunks or a code chunk
executing image insert in loops.. a bit tedious i suppose!).
## Decluttering overlapping labels produced by plot.MCA and plot.MFA
function in SensoMineR (Refer the functions plot.MCA2 and plot.MFA2 for
this feature)
##The fastmod function includes a call to plot.MCA2 and plot.MFA2 as well
so the final output is also decluttered of overlapping labels, however
before calling the fastmod function a source or call to ##plot.MFA2 and
plot.MCA2 is required.
## library(maptools) need to be called before calling the functions.
###plot.MFA2 Function
plot.MFA2<-
function (x, axes = c(1, 2), choix = "ind", ellipse = NULL,
ellipse.par =
NULL,
lab.grpe = TRUE, lab.var = TRUE, lab.ind.moy = TRUE, lab.par = FALSE,
habillage = "ind", col.hab = NULL, invisible = NULL, partial =
NULL,
lim.cos2.var = 0, chrono = FALSE, xlim = NULL, ylim = NULL,
cex = 1, title = NULL, palette = NULL, new.plot = TRUE, ...)
{
res.mfa <- x
if (!inherits(res.mfa, "MFA"))
stop("non convenient data")
if (is.null(palette))
palette(c("black", "red", "green3",
"blue", "cyan", "magenta",
"darkgray", "darkgoldenrod",
"darkgreen", "violet",
"turquoise", "orange", "lightpink",
"lavender", "yellow",
"lightgreen", "lightgrey",
"lightblue", "darkkhaki",
"darkmagenta", "darkolivegreen",
"lightcyan", "darkorange",
"darkorchid", "darkred", "darksalmon",
"darkseagreen",
"darkslateblue", "darkslategray",
"darkslategrey",
"darkturquoise", "darkviolet",
"lightgray", "lightsalmon",
"lightyellow", "maroon"))
lab.x <- paste("Dim ", axes[1], " (",
signif(res.mfa$eig[axes[1],
2], 4), " %)", sep = "")
lab.y <- paste("Dim ", axes[2], " (",
signif(res.mfa$eig[axes[2],
2], 4), " %)", sep = "")
group <- res.mfa$call$group
nbre.grpe <- length(group)
type <- res.mfa$call$type
num.group.sup = NULL
if (!is.null(res.mfa$group$coord.sup)) {
num.group.sup <- res.mfa$call$num.group.sup
nbre.grpe.sup <- length(num.group.sup)
type.sup <- res.mfa$call$type[num.group.sup]
type.act <- type[-num.group.sup]
nbre.grpe <- nbre.grpe - length(num.group.sup)
}
if (choix == "axes") {
if (new.plot)
dev.new()
if (is.null(title))
title <- "Partial axes"
coord.axes <- res.mfa$partial.axes$coord[, axes, drop = FALSE]
plot(0, 0, xlab = lab.x, ylab = lab.y, xlim = c(-1.1,
1.1), ylim = c(-1.1, 1.1), col = "white", asp = 1,
cex = cex, main = title)
x.cercle <- seq(-1, 1, by = 0.01)
y.cercle <- sqrt(1 - x.cercle^2)
lines(x.cercle, y = y.cercle)
lines(x.cercle, y = -y.cercle)
abline(v = 0, lty = 2, cex = cex)
abline(h = 0, lty = 2, cex = cex)
if (habillage == "group") {
if (is.null(col.hab) | length(col.hab) < length(group)) {
col.hab <- 2:(length(group) + 1)
}
i = 1
couleur.axes <- col.hab[i]
auxil = strsplit(rownames(res.mfa$partial.axes$coord)[1],
".", fixed = TRUE)[[1]]
auxil2 = auxil[length(auxil)]
for (j in 2:nrow(res.mfa$partial.axes$coord)) {
auxil = strsplit(rownames(res.mfa$partial.axes$coord)[j],
".", fixed = TRUE)[[1]]
if (auxil2 != auxil[length(auxil)]) {
i = i + 1
auxil2 = auxil[length(auxil)]
}
couleur.axes <- c(couleur.axes, col.hab[i])
}
}
else {
couleur.axes <- NULL
for (i in 1:length(group)) couleur.axes <- c(couleur.axes,
rep("black", ncol(res.mfa$partial.axes$coord)))
}
for (v in 1:nrow(coord.axes)) {
arrows(0, 0, coord.axes[v, 1], coord.axes[v, 2],
length = 0.1, angle = 15, code = 2, col = couleur.axes[v])
if (abs(coord.axes[v, 1]) > abs(coord.axes[v, 2])) {
if (coord.axes[v, 1] >= 0)
pos <- 4
else pos <- 2
}
else {
if (coord.axes[v, 2] >= 0)
pos <- 3
else pos <- 1
}
text(coord.axes[v, 1], y = coord.axes[v, 2], labels =
rownames(coord.axes)[v],
pos = pos, col = couleur.axes[v])
}
if (habillage == "group") {
legend("topleft", legend =
rownames(res.mfa$group$Lg)[-length(rownames(res.mfa$group$Lg))],
text.col = unique(couleur.axes), cex = 0.8)
}
}
if (choix == "group") {
if (new.plot)
dev.new()
if (is.null(title))
title <- "Groups representation"
coord.actif <- res.mfa$group$coord[, axes, drop = FALSE]
if (!is.null(res.mfa$group$coord.sup))
coord.illu <- res.mfa$group$coord.sup[, axes, drop = FALSE]
if (is.null(col.hab)) {
col.hab = rep("darkred", nrow(coord.actif))
if (!is.null(res.mfa$group$coord.sup))
col.hab = c(col.hab, rep("darkolivegreen",
nrow(coord.illu)))
}
if (habillage == "group")
col.hab <- (2:(length(group) + 1))
plot(coord.actif, xlab = lab.x, ylab = lab.y, xlim = c(0,
1), ylim = c(0, 1), pch = 17, col =
col.hab[1:nrow(coord.actif)],
cex = cex, main = title, cex.main = cex * 1.2, asp = 1)
if (lab.grpe)
pointLabel(coord.actif[, 1], y = coord.actif[, 2], labels =
rownames(coord.actif),
allowSmallOverlap=FALSE,trace=FALSE,doPlot=TRUE, col =
col.hab[1:nrow(coord.actif)])
if (!is.null(res.mfa$group$coord.sup)) {
points(coord.illu, pch = 17, col = col.hab[(nrow(coord.actif)
+
1):(nrow(coord.actif) + nrow(coord.illu))])
if (lab.grpe)
pointLabel(coord.illu[, 1], y = coord.illu[, 2], labels =
rownames(coord.illu),
allowSmallOverlap=FALSE,trace=FALSE,doPlot=TRUE,col =
col.hab[(nrow(coord.actif) +
1):(nrow(coord.actif) + nrow(coord.illu))])
}
}
if (choix == "var") {
if (new.plot)
dev.new()
test.invisible <- vector(length = 2)
if (!is.null(invisible)) {
test.invisible[1] <- match("actif", invisible)
test.invisible[2] <- match("sup", invisible)
}
else test.invisible <- rep(NA, 2)
col <- NULL
if (habillage == "group") {
if (is.null(col.hab) | length(col.hab) < length(group[type !=
"n"]))
col.hab <- 2:(length(group[type != "n"]) + 1)
for (i in 1:length(group[type != "n"])) col <- c(col,
rep(col.hab[i], group[type != "n"][i]))
}
else {
if (is.null(col.hab) | length(col.hab) < sum(group[type !=
"n"]))
col <- rep(1, sum(group[type != "n"]))
else col <- col.hab
}
if (is.null(title))
title <- "Correlation circle"
plot(0, 0, main = title, xlab = lab.x, ylab = lab.y,
xlim = c(-1.1, 1.1), ylim = c(-1.1, 1.1), col = "white",
asp = 1, cex = cex)
x.cercle <- seq(-1, 1, by = 0.01)
y.cercle <- sqrt(1 - x.cercle^2)
lines(x.cercle, y = y.cercle)
lines(x.cercle, y = -y.cercle)
abline(v = 0, lty = 2, cex = cex)
abline(h = 0, lty = 2, cex = cex)
if (habillage == "group" & is.na(test.invisible[1]) &
is.na(test.invisible[2]))
legend("topleft", legend =
rownames(res.mfa$group$Lg[-nrow(res.mfa$group$Lg),
])[type != "n"], text.col = col.hab, cex = 0.8)
if (habillage == "group" & is.na(test.invisible[1]) &
!is.na(test.invisible[2]))
legend("topleft", legend =
rownames(res.mfa$group$Lg[-c(num.group.sup,
nrow(res.mfa$group$Lg)), ])[type.act != "n"],
text.col = col.hab, cex = 0.8)
if (habillage == "group" & !is.na(test.invisible[1]) &
is.na(test.invisible[2]))
legend("topleft", legend =
rownames(res.mfa$group$Lg[num.group.sup,
])[type.sup != "n"], text.col = col.hab, cex = 0.8)
nrow.coord.var <- 0
if (!is.null(res.mfa["quanti.var"]$quanti.var$coord) &
is.na(test.invisible[1])) {
coord.var <- res.mfa$quanti.var$cor[, axes, drop = FALSE]
nrow.coord.var <- nrow(coord.var)
for (v in 1:nrow(coord.var)) {
if (sum(res.mfa$quanti.var$cos2[v, axes], na.rm = TRUE) >=
lim.cos2.var && !is.na(sum(res.mfa$quanti.var$cos2[v,
axes], na.rm = TRUE))) {
arrows(0, 0, coord.var[v, 1], coord.var[v,
2], length = 0.1, angle = 15, code = 2, col = col[v])
if (lab.var) {
if (abs(coord.var[v, 1]) > abs(coord.var[v,
2])) {
if (coord.var[v, 1] >= 0)
pos <- 4
else pos <- 2
}
else {
if (coord.var[v, 2] >= 0)
pos <- 3
else pos <- 1
}
pointLabel(coord.var[v, 1], y = coord.var[v, 2],
labels = rownames(coord.var)[v],
allowSmallOverlap=FALSE,trace=FALSE,doPlot=TRUE,
col = col[v])
}
}
}
}
if (!is.null(res.mfa$quanti.var.sup$coord) &
is.na(test.invisible[2])) {
coord.var.sup <- res.mfa$quanti.var.sup$cor[, axes,
drop = FALSE]
for (q in 1:nrow(coord.var.sup)) {
print(sum(res.mfa$quanti.var.sup$cos2[q, axes],
na.rm = TRUE))
if (sum(res.mfa$quanti.var.sup$cos2[q, axes],
na.rm = TRUE) >= lim.cos2.var &&
!is.na(sum(res.mfa$quanti.var.sup$cos2[q,
axes], na.rm = TRUE))) {
arrows(0, 0, coord.var.sup[q, 1], coord.var.sup[q,
2], length = 0.1, angle = 15, code = 2, lty = 2,
col = col[nrow.coord.var + q])
if (lab.var) {
if (abs(coord.var.sup[q, 1]) > abs(coord.var.sup[q,
2])) {
if (coord.var.sup[q, 1] >= 0)
pos <- 4
else pos <- 2
}
else {
if (coord.var.sup[q, 2] >= 0)
pos <- 3
else pos <- 1
}
pointLabel(coord.var.sup[q, 1], y = coord.var.sup[q,
2], labels = rownames(coord.var.sup)[q],
allowSmallOverlap=FALSE,trace=FALSE,doPlot=TRUE, col
= col[nrow.coord.var + q])
}
}
}
}
par(mar = c(5, 4, 4, 2) + 0.1)
}
if (choix == "ind") {
test.invisible <- vector(length = 3)
if (!is.null(invisible)) {
test.invisible[1] <- match("ind", invisible)
test.invisible[2] <- match("ind.sup", invisible)
test.invisible[3] <- match("quali", invisible)
}
else test.invisible <- rep(NA, 3)
nb.ind.actif <- nrow(res.mfa$ind$coord)
nb.ind.illu <- 0
if (!is.null(res.mfa$ind.sup))
nb.ind.illu <- nrow(res.mfa$ind.sup$coord)
nb.ind <- nb.ind.actif + nb.ind.illu
coord.ind <- res.mfa$ind$coord[, axes, drop = FALSE]
coord.ind.partiel <- res.mfa$ind$coord.partiel[, axes,
drop = FALSE]
coord.ind.sup <- NULL
if (!is.null(res.mfa$ind.sup)) {
coord.ind.sup <- res.mfa$ind.sup$coord[, axes, drop = FALSE]
coord.ind.partiel.sup <- res.mfa$ind.sup$coord.partiel[,
axes, drop = FALSE]
}
coord.quali <- coord.quali.sup <- coord.quali.partiel <-
coord.quali.sup.partiel <- NULL
nrow.coord.quali <- 0
if (!is.null(res.mfa["quali.var"]$quali.var)) {
coord.quali <- res.mfa$quali.var$coord[, axes, drop = FALSE]
coord.quali.partiel <- res.mfa$quali.var$coord.partiel[,
axes, drop = FALSE]
nrow.coord.quali <- nrow(coord.quali)
}
if (!is.null(res.mfa["quali.var.sup"])) {
coord.quali.sup <- res.mfa$quali.var.sup$coord[,
axes, drop = FALSE]
coord.quali.partiel.sup <-
res.mfa$quali.var.sup$coord.partiel[,
axes, drop = FALSE]
}
group.ind.actif <- group.ind.sup <- group.quali <-
group.quali.sup
<- NULL
if (!is.null(partial)) {
if (length(partial) == 1) {
if (partial == "all") {
group.ind.actif <- 1:nrow(coord.ind)
if (!is.null(res.mfa$ind.sup))
group.ind.sup <- 1:nrow(coord.ind.sup)
if (!is.null(res.mfa["quali.var"]$quali.var))
group.quali <- 1:nrow(coord.quali)
if
(!is.null(res.mfa["quali.var.sup"]$quali.var.sup))
group.quali.sup <- 1:nrow(coord.quali.sup)
}
else {
for (i in 1:length(partial)) {
if (partial[i] %in% rownames(coord.ind))
group.ind.actif <- c(group.ind.actif,
grep(partial[i],
rownames(coord.ind)))
if (partial[i] %in% rownames(coord.ind.sup))
group.ind.sup <- c(group.ind.sup, grep(partial[i],
rownames(coord.ind.sup)))
if (partial[i] %in% rownames(coord.quali))
group.quali <- c(group.quali, grep(partial[i],
rownames(coord.quali)))
if (partial[i] %in% rownames(coord.quali.sup))
group.quali.sup <- c(group.quali.sup,
grep(partial[i],
rownames(coord.quali.sup)))
}
}
}
else {
for (i in 1:length(partial)) {
if (partial[i] %in% rownames(coord.ind))
group.ind.actif <- c(group.ind.actif, grep(partial[i],
rownames(coord.ind)))
if (partial[i] %in% rownames(coord.ind.sup))
group.ind.sup <- c(group.ind.sup, grep(partial[i],
rownames(coord.ind.sup)))
if (partial[i] %in% rownames(coord.quali))
group.quali <- c(group.quali, grep(partial[i],
rownames(coord.quali)))
if (partial[i] %in% rownames(coord.quali.sup))
group.quali.sup <- c(group.quali.sup, grep(partial[i],
rownames(coord.quali.sup)))
}
}
}
if (!is.null(ellipse)) {
coord.ellipse <- ellipse$res
npoint.ellipse <- ellipse$call
}
else coord.ellipse <- NULL
if (!is.null(ellipse.par)) {
coord.ellipse.par <- ellipse.par$res
npoint.ellipse.par <- ellipse.par$call
}
else coord.ellipse.par <- NULL
if (is.null(xlim)) {
xmin <- xmax <- 0
if (is.na(test.invisible[1]))
xmin <- min(xmin, coord.ind[, 1])
if (is.na(test.invisible[1]))
xmax <- max(xmax, coord.ind[, 1])
if (!is.null(res.mfa$ind.sup) & is.na(test.invisible[2]))
xmin <- min(xmin, coord.ind.sup[, 1])
if (!is.null(res.mfa$ind.sup) & is.na(test.invisible[2]))
xmax <- max(xmax, coord.ind.sup[, 1])
if (is.na(test.invisible[1]))
xmin <- min(xmin,
coord.ind.partiel[unlist(lapply(group.ind.actif,
function(k) seq(nbre.grpe * (k - 1) + 1, length =
nbre.grpe))),
1])
if (is.na(test.invisible[1]))
xmax <- max(xmax,
coord.ind.partiel[unlist(lapply(group.ind.actif,
function(k) seq(nbre.grpe * (k - 1) + 1, length =
nbre.grpe))),
1])
if (!is.null(res.mfa$ind.sup) & is.na(test.invisible[2]))
xmin <- min(xmin,
coord.ind.partiel.sup[unlist(lapply(group.ind.sup,
function(k) seq(nbre.grpe * (k - 1) + 1, length =
nbre.grpe))),
1])
if (!is.null(res.mfa$ind.sup) & is.na(test.invisible[2]))
xmax <- max(xmax,
coord.ind.partiel.sup[unlist(lapply(group.ind.sup,
function(k) seq(nbre.grpe * (k - 1) + 1, length =
nbre.grpe))),
1])
if (!is.null(res.mfa["quali.var"]$quali.var) &
is.na(test.invisible[3]))
xmin <- min(xmin, coord.quali[, 1])
if (!is.null(res.mfa["quali.var"]$quali.var) &
is.na(test.invisible[3]))
xmax <- max(xmax, coord.quali[, 1])
if (!is.null(res.mfa["quali.var"]$quali.var) &
is.na(test.invisible[3]))
xmin <- min(xmin,
coord.quali.partiel[unlist(lapply(group.quali,
function(k) seq(nbre.grpe * (k - 1) + 1, length =
nbre.grpe))),
1])
if (!is.null(res.mfa["quali.var"]$quali.var) &
is.na(test.invisible[3]))
xmax <- max(xmax,
coord.quali.partiel[unlist(lapply(group.quali,
function(k) seq(nbre.grpe * (k - 1) + 1, length =
nbre.grpe))),
1])
if (!is.null(res.mfa$quali.var.sup) &
is.na(test.invisible[3]))
xmin <- min(xmin, coord.quali[, 1], coord.quali.sup[,
1])
if (!is.null(res.mfa$quali.var.sup) &
is.na(test.invisible[3]))
xmax <- max(xmax, coord.quali[, 1], coord.quali.sup[,
1])
if (!is.null(res.mfa$quali.var.sup) &
is.na(test.invisible[3]))
xmin <- min(xmin,
coord.quali.partiel.sup[unlist(lapply(group.quali.sup,
function(k) seq(nbre.grpe * (k - 1) + 1, length =
nbre.grpe))),
1])
if (!is.null(res.mfa$quali.var.sup) &
is.na(test.invisible[3]))
xmax <- max(xmax,
coord.quali.partiel.sup[unlist(lapply(group.quali.sup,
function(k) seq(nbre.grpe * (k - 1) + 1, length =
nbre.grpe))),
1])
xlim <- c(xmin, xmax) * 1.1
}
else {
xmin = xlim[1]
xmax = xlim[2]
}
if (is.null(ylim)) {
ymin <- ymax <- 0
if (is.na(test.invisible[1]))
ymin <- min(ymin, coord.ind[, 2])
if (is.na(test.invisible[1]))
ymax <- max(ymax, coord.ind[, 2])
if (!is.null(res.mfa$ind.sup) & is.na(test.invisible[2]))
ymin <- min(ymin, coord.ind.sup[, 2])
if (!is.null(res.mfa$ind.sup) & is.na(test.invisible[2]))
ymax <- max(ymax, coord.ind.sup[, 2])
if (is.na(test.invisible[1]))
ymin <- min(ymin,
coord.ind.partiel[unlist(lapply(group.ind.actif,
function(k) seq(nbre.grpe * (k - 1) + 1, length =
nbre.grpe))),
2])
if (is.na(test.invisible[1]))
ymax <- max(ymax,
coord.ind.partiel[unlist(lapply(group.ind.actif,
function(k) seq(nbre.grpe * (k - 1) + 1, length =
nbre.grpe))),
2])
if (!is.null(res.mfa$ind.sup) & is.na(test.invisible[2]))
ymin <- min(ymin,
coord.ind.partiel.sup[unlist(lapply(group.ind.sup,
function(k) seq(nbre.grpe * (k - 1) + 1, length =
nbre.grpe))),
2])
if (!is.null(res.mfa$ind.sup) & is.na(test.invisible[2]))
ymax <- max(ymax,
coord.ind.partiel.sup[unlist(lapply(group.ind.sup,
function(k) seq(nbre.grpe * (k - 1) + 1, length =
nbre.grpe))),
2])
if (!is.null(res.mfa["quali.var"]$quali.var) &
is.na(test.invisible[3]))
ymin <- min(ymin, coord.quali[, 2])
if (!is.null(res.mfa["quali.var"]$quali.var) &
is.na(test.invisible[3]))
ymax <- max(ymax, coord.quali[, 2])
if (!is.null(res.mfa["quali.var"]$quali.var) &
is.na(test.invisible[3]))
ymin <- min(ymin,
coord.quali.partiel[unlist(lapply(group.quali,
function(k) seq(nbre.grpe * (k - 1) + 1, length =
nbre.grpe))),
2])
if (!is.null(res.mfa["quali.var"]$quali.var) &
is.na(test.invisible[3]))
ymax <- max(ymax,
coord.quali.partiel[unlist(lapply(group.quali,
function(k) seq(nbre.grpe * (k - 1) + 1, length =
nbre.grpe))),
2])
if (!is.null(res.mfa$quali.var.sup) &
is.na(test.invisible[3]))
ymin <- min(ymin, coord.quali[, 1], coord.quali.sup[,
2])
if (!is.null(res.mfa$quali.var.sup) &
is.na(test.invisible[3]))
ymax <- max(ymax, coord.quali[, 1], coord.quali.sup[,
2])
if (!is.null(res.mfa$quali.var.sup) &
is.na(test.invisible[3]))
ymin <- min(ymin,
coord.quali.partiel.sup[unlist(lapply(group.quali.sup,
function(k) seq(nbre.grpe * (k - 1) + 1, length =
nbre.grpe))),
2])
if (!is.null(res.mfa$quali.var.sup) &
is.na(test.invisible[3]))
ymax <- max(ymax,
coord.quali.partiel.sup[unlist(lapply(group.quali.sup,
function(k) seq(nbre.grpe * (k - 1) + 1, length =
nbre.grpe))),
2])
ylim <- c(ymin, ymax) * 1.1
}
else {
ymin = ylim[1]
ymax = ylim[2]
}
if (habillage == "group") {
if (is.null(col.hab) | length(col.hab) != (nbre.grpe))
col.hab <- 2:(nbre.grpe + 1)
col.ind <- c(rep(1, nb.ind.actif), rep(col.hab, nb.ind.actif))
if (!is.null(res.mfa$ind.sup))
col.ind.sup <- c(rep(1, nb.ind - nb.ind.actif),
rep(col.hab, nb.ind - nb.ind.actif))
if (length(group[type == "n"]) != 0)
col.quali <- c(rep(1, sum(res.mfa$call$group.mod[type ==
"n"])), rep(col.hab, sum(res.mfa$call$group.mod[type
==
"n"])))
if (!is.null(res.mfa$quali.var.sup))
col.quali.sup <- c(rep(1,
sum(res.mfa$call$group.mod[num.group.sup][type.sup ==
"n"])), rep(col.hab,
sum(res.mfa$call$group.mod[num.group.sup][type.sup ==
"n"])))
if (!is.null(ellipse))
col.ellipse <- rep(1, nb.ind.actif)
if (!is.null(ellipse.par))
col.ellipse.par <- rep(col.hab, nb.ind.actif)
}
if (habillage == "ind") {
if (is.null(col.hab) | length(col.hab) != nb.ind) {
col.hab <- 1:nb.ind
}
col.ind <- c(col.hab[1:nb.ind.actif],
rep(col.hab[1:nb.ind.actif],
each = nbre.grpe))
if (!is.null(res.mfa$ind.sup))
col.ind.sup <- c(col.hab[(nb.ind.actif + 1):nb.ind],
rep(col.hab[(nb.ind.actif + 1):nb.ind], each =
nbre.grpe))
if (length(group[type == "n"]) != 0)
col.quali <- col.quali.sup <- rep("black", (1 +
nbre.grpe) * sum(res.mfa$call$group.mod[type ==
"n"]))
if (!is.null(ellipse))
col.ellipse <- col.hab[1:nb.ind.actif]
if (!is.null(ellipse.par))
col.ellipse.par <- rep(col.hab[1:nb.ind.actif],
each = nbre.grpe)
}
if ((habillage != "none") & (habillage != "ind")
& (habillage !=
"group")) {
group.act <- (1:length(group))
if (!is.null(num.group.sup))
group.act <- group.act[-num.group.sup]
nbre.modalite <- NULL
liste.quali <- NULL
for (i in group.act) {
if (type[i] == "n") {
for (k in 1:ncol(res.mfa$separate.analyses[[i]]$call$X))
nbre.modalite <- c(nbre.modalite,
nlevels(res.mfa$separate.analyses[[i]]$call$X[,
k]))
if (i == 1)
liste.quali <- c(liste.quali,
colnames(res.mfa$call$X[1:group[1]]))
else liste.quali <- c(liste.quali,
colnames(res.mfa$call$X[(sum(group[1:(i -
1)]) + 1):sum(group[1:i])]))
}
}
if (!is.null(num.group.sup)) {
for (i in num.group.sup) {
if (type[i] == "n") {
if (i == 1)
liste.quali <- c(liste.quali,
colnames(res.mfa$call$X[1:group[1]]))
else liste.quali <- c(liste.quali,
colnames(res.mfa$call$X[(sum(group[1:(i -
1)]) + 1):sum(group[1:i])]))
for (k in
1:ncol(res.mfa$separate.analyses[[i]]$call$X)) nbre.modalite <-
c(nbre.modalite,
nlevels(res.mfa$separate.analyses[[i]]$call$X[,
k]))
}
}
}
if (is.double(habillage))
nom.quali <- colnames(res.mfa$call$X)[habillage]
else nom.quali = habillage
if (!(nom.quali %in% liste.quali))
stop("The variable ", habillage, " is not
qualitative")
modalite <- levels(as.factor(res.mfa$call$X[, nom.quali]))
col.ind <- as.numeric(as.factor(res.mfa$call$X[,
nom.quali]))
if (is.null(col.hab) | length(col.hab) != length(modalite))
col.hab <- 2:(1 + length(modalite))
col.ind <- col.hab[col.ind]
if (!is.null(res.mfa$call$ind.sup)) {
col.ind.sup <- col.ind[res.mfa$call$ind.sup]
col.ind <- col.ind[-res.mfa$call$ind.sup]
col.ind.sup <- c(col.ind.sup, rep(col.ind.sup,
each = nbre.grpe))
}
col.ind <- c(col.ind, rep(col.ind, each = nbre.grpe))
col.ellipse <- col.ind[1:nb.ind.actif]
col.ellipse.par <- col.ind[-c(1:nb.ind.actif)]
indice.inf <- sum(nbre.modalite[0:(match(nom.quali,
liste.quali) - 1)]) + 1
indice.sup <- indice.inf + length(modalite) - 1
col.quali <- rep("black",
sum(res.mfa$call$group.mod[type ==
"n"]))
if (length(group[type == "n"]) != 0) {
for (i in 1:length(liste.quali)) {
if (liste.quali[i] == nom.quali)
col.quali[indice.inf:indice.sup] <- col.hab
}
}
col.quali <- c(col.quali, rep(col.quali, each = nbre.grpe))
col.quali.sup <- col.quali
}
if (habillage == "none")
col.ind <- col.ind.sup <- col.quali.sup <- col.quali <-
col.ellipse <- col.ellipse.par <- rep("black",
nb.ind * (nbre.grpe + 1))
if (new.plot)
dev.new(width = min(14, max(8, 8 * (xmax - xmin)/(ymax -
ymin))), height = 8)
if (is.null(title))
title <- "Individual factor map"
plot(0, 0, main = title, xlab = lab.x, ylab = lab.y,
xlim = xlim, ylim = ylim, col = "white", asp = 1,
cex = cex)
abline(v = 0, lty = 2, cex = cex)
abline(h = 0, lty = 2, cex = cex)
if (is.na(test.invisible[1])) {
points(coord.ind, pch = 20, col = col.ind[1:nb.ind.actif],
cex = cex)
if (lab.ind.moy)
pointLabel(coord.ind[, 1], y = coord.ind[, 2], labels =
rownames(coord.ind),
allowSmallOverlap=FALSE,trace=FALSE,doPlot=TRUE, col =
col.ind[1:nb.ind.actif])
for (i in group.ind.actif) {
for (j in 1:nbre.grpe) {
points(coord.ind.partiel[(i - 1) * nbre.grpe +
j, ], cex = 0.8 * cex, col = col.ind[nb.ind.actif +
(i - 1) * nbre.grpe + j], pch = 20)
if (lab.par)
pointLabel(coord.ind.partiel[(i - 1) * nbre.grpe +
j, 1], y = coord.ind.partiel[(i - 1) *
nbre.grpe + j, 2], labels =
rownames(coord.ind.partiel)[(i -
1) * nbre.grpe + j],
allowSmallOverlap=FALSE,trace=FALSE,doPlot=TRUE, col =
col.ind[nb.ind.actif +
(i - 1) * nbre.grpe + j])
if (chrono) {
if (j > 1)
lines(c(coord.ind.partiel[(i - 1) * nbre.grpe +
(j - 1), 1], coord.ind.partiel[(i - 1) *
nbre.grpe + j, 1]), c(coord.ind.partiel[(i -
1) * nbre.grpe + (j - 1), 2], coord.ind.partiel[(i
-
1) * nbre.grpe + j, 2]), col = col.ind[i])
}
else lines(c(coord.ind[i, 1], coord.ind.partiel[(i -
1) * nbre.grpe + j, 1]), c(coord.ind[i, 2],
coord.ind.partiel[(i - 1) * nbre.grpe + j,
2]), col = col.ind[nb.ind.actif + (i -
1) * nbre.grpe + j], lty = j)
}
}
}
if (!is.null(res.mfa$ind.sup) & is.na(test.invisible[2])) {
points(coord.ind.sup, pch = 21, col = col.ind.sup[1:(nb.ind -
nb.ind.actif)], cex = cex)
if (lab.ind.moy)
pointLabel(coord.ind.sup[, 1], y = coord.ind.sup[,
2], labels = rownames(coord.ind.sup),
allowSmallOverlap=FALSE,trace=FALSE,doPlot=TRUE,
col = col.ind.sup[1:(nb.ind - nb.ind.actif)])
for (i in group.ind.sup) {
for (j in 1:nbre.grpe) {
points(coord.ind.partiel.sup[(i - 1) * nbre.grpe +
j, ], cex = 0.8 * cex, col = col.ind.sup[nb.ind -
nb.ind.actif + (i - 1) * nbre.grpe + j],
pch = 21)
if (lab.par)
pointLabel(coord.ind.partiel.sup[(i - 1) * nbre.grpe +
j, 1], y = coord.ind.partiel.sup[nb.ind +
(i - 1) * nbre.grpe + j, 2], labels =
rownames(coord.ind.partiel.sup)[(i -
1) * nbre.grpe + j],
allowSmallOverlap=FALSE,trace=FALSE,doPlot=TRUE, col = col.ind.sup[nb.ind
-
nb.ind.actif + (i - 1) * nbre.grpe + j])
if (chrono) {
if (j > 1)
lines(c(coord.ind.partiel.sup[(i - 1) *
nbre.grpe + (j - 1), 1], coord.ind.partiel.sup[(i
-
1) * nbre.grpe + j, 1]),
c(coord.ind.partiel.sup[(i -
1) * nbre.grpe + (j - 1), 2],
coord.ind.partiel.sup[(i -
1) * nbre.grpe + j, 2]), col =
col.ind[nb.ind.actif +
i])
}
else lines(c(coord.ind.sup[i, 1],
coord.ind.partiel.sup[(i -
1) * nbre.grpe + j, 1]), c(coord.ind.sup[i,
2], coord.ind.partiel.sup[(i - 1) * nbre.grpe +
j, 2]), col = col.ind.sup[nb.ind - nb.ind.actif +
(i - 1) * nbre.grpe + j], lty = j)
}
}
}
if (!is.null(coord.quali) & is.na(test.invisible[3])) {
points(coord.quali, pch = 15, col =
col.quali[1:nrow.coord.quali],
cex = cex)
if (lab.var)
pointLabel(coord.quali[, 1], y = coord.quali[, 2],
labels = rownames(coord.quali),
allowSmallOverlap=FALSE,trace=FALSE,doPlot=TRUE, col =
col.quali[1:nrow.coord.quali])
for (i in group.quali) {
for (j in 1:nbre.grpe) {
points(coord.quali.partiel[(i - 1) * nbre.grpe +
j, ], pch = 15, col = col.quali[nrow.coord.quali +
(i - 1) * nbre.grpe + j], cex = cex * 0.8)
if (lab.var & lab.par)
pointLabel(coord.quali.partiel[(i - 1) * nbre.grpe +
j, 1], y = coord.quali.partiel[(i - 1) *
nbre.grpe + j, 2], labels =
rownames(coord.quali.partiel)[(i -
1) * nbre.grpe +
j],allowSmallOverlap=FALSE,trace=FALSE,doPlot=TRUE, col =
col.quali[nrow.coord.quali +
(i - 1) * nbre.grpe + j])
if (chrono) {
if (j > 1)
lines(c(coord.quali.partiel[(i - 1) * nbre.grpe +
(j - 1), 1], coord.quali.partiel[(i -
1) * nbre.grpe + j, 1]), c(coord.quali.partiel[(i
-
1) * nbre.grpe + (j - 1), 2],
coord.quali.partiel[(i -
1) * nbre.grpe + j, 2]), col = col.quali[i])
}
else lines(c(coord.quali[i, 1], coord.quali.partiel[(i -
1) * nbre.grpe + j, 1]), c(coord.quali[i,
2], coord.quali.partiel[(i - 1) * nbre.grpe +
j, 2]), col = col.quali[nrow.coord.quali +
(i - 1) * nbre.grpe + j], lty = j)
}
}
}
if (!is.null(coord.quali.sup) & is.na(test.invisible[3])) {
points(coord.quali.sup, pch = 22, col =
col.quali.sup[1:nrow(coord.quali.sup)],
cex = cex)
if (lab.var)
pointLabel(coord.quali.sup[, 1], y = coord.quali.sup[,
2], labels = rownames(coord.quali.sup),
allowSmallOverlap=FALSE,trace=FALSE,doPlot=TRUE,
col = col.quali.sup[1:nrow(coord.quali.sup)])
for (i in group.quali.sup) {
for (j in 1:nbre.grpe) {
points(coord.quali.partiel.sup[(i - 1) * nbre.grpe +
j, ], pch = 22, col =
col.quali.sup[nrow(coord.quali.sup) +
(i - 1) * nbre.grpe + j], cex = cex * 0.8)
if (lab.var & lab.par)
pointLabel(coord.quali.partiel.sup[(i - 1) * nbre.grpe
+
j, 1], y = coord.quali.partiel.sup[(i -
1) * nbre.grpe + j, 2], labels =
rownames(coord.quali.partiel.sup)[(i -
1) * nbre.grpe + j],
allowSmallOverlap=FALSE,trace=FALSE,doPlot=TRUE, col =
col.quali.sup[nrow(coord.quali.sup) +
(i - 1) * nbre.grpe + j])
if (chrono) {
if (j > 1)
lines(c(coord.quali.partiel.sup[(i - 1) *
nbre.grpe + (j - 1), 1],
coord.quali.partiel.sup[(i -
1) * nbre.grpe + j, 1]),
c(coord.quali.partiel.sup[(i -
1) * nbre.grpe + (j - 1), 2],
coord.quali.partiel.sup[(i -
1) * nbre.grpe + j, 2]), col =
col.quali[nrow.coord.quali +
i])
}
else lines(c(coord.quali.sup[i, 1],
coord.quali.partiel.sup[(i -
1) * nbre.grpe + j, 1]), c(coord.quali.sup[i,
2], coord.quali.partiel.sup[(i - 1) * nbre.grpe +
j, 2]), col = col.quali.sup[nrow(coord.quali.sup) +
(i - 1) * nbre.grpe + j], lty = j)
}
}
}
if ((!is.null(partial)) & (habillage == "group"))
legend("topleft", legend =
rownames(res.mfa$group$Lg)[-c(num.group.sup,
length(rownames(res.mfa$group$Lg)))], lty =
1:length(rownames(res.mfa$group$Lg)[-c(num.group.sup,
length(rownames(res.mfa$group$Lg)))]), text.col = col.hab,
col = col.hab, cex = 0.8)
if ((!is.null(partial)) & (habillage != "group"))
legend("topleft", legend =
rownames(res.mfa$group$Lg)[-c(num.group.sup,
length(rownames(res.mfa$group$Lg)))], lty =
1:length(rownames(res.mfa$group$Lg)[-c(num.group.sup,
length(rownames(res.mfa$group$Lg)))]), cex = 0.8)
if ((habillage != "none") & (habillage != "ind")
& (habillage !=
"group"))
legend("topleft", legend = levels(res.mfa$call$X[,
habillage]), text.col = col.hab, cex = 0.8)
if (!is.null(coord.ellipse) & is.na(test.invisible[2])) {
for (e in 1:nb.ind.actif) {
debut <- ((nb.ind.actif - 1) * npoint.ellipse) +
1
fin <- debut + npoint.ellipse - 1
data.elli <- coord.ellipse[debut:fin, -1]
lines(data.elli[, 1], y = data.elli[, 2], col =
col.ellipse[e])
}
}
if (!is.null(coord.ellipse)) {
for (e in 1:nlevels(coord.ellipse[, 1])) {
data.elli <- coord.ellipse[(npoint.ellipse *
(e - 1) + 1):(npoint.ellipse * e), -1]
lines(data.elli[, 1], y = data.elli[, 2], col =
col.ellipse[e])
}
}
if (!is.null(coord.ellipse.par)) {
for (i in group.ind.actif) {
for (j in 1:nbre.grpe) {
ind.e <- (i - 1) * nbre.grpe + j
data.elli <- coord.ellipse.par[(npoint.ellipse.par *
(ind.e - 1) + 1):(npoint.ellipse.par * ind.e),
-1]
lines(data.elli[, 1], y = data.elli[, 2], col =
col.ellipse.par[ind.e],
lty = 2)
}
}
}
}
###plot.MCA2 Function
plot.MCA2<-
function (x, axes = c(1, 2), choix = "ind", xlim = NULL, ylim = NULL,
invisible = NULL, col.ind = "blue", col.var = "red",
col.quali.sup =
"darkgreen",
col.ind.sup = "darkblue", col.quanti.sup = "black",
label = "all",
cex = 1, title = NULL, habillage = "none", palette = NULL,
new.plot = TRUE, ...)
{
res.mca <- x
if (!inherits(res.mca, "MCA"))
stop("non convenient data")
if (is.null(palette))
palette(c("black", "red", "green3",
"blue", "cyan", "magenta",
"darkgray", "darkgoldenrod",
"darkgreen", "violet",
"turquoise", "orange", "lightpink",
"lavender", "yellow",
"lightgreen", "lightgrey",
"lightblue", "darkkhaki",
"darkmagenta", "darkolivegreen",
"lightcyan", "darkorange",
"darkorchid", "darkred", "darksalmon",
"darkseagreen",
"darkslateblue", "darkslategray",
"darkslategrey",
"darkturquoise", "darkviolet",
"lightgray", "lightsalmon",
"lightyellow", "maroon"))
lab.x = paste("Dim ", axes[1], " (",
signif(res.mca$eig[axes[1],
2], 4), "%)", sep = "")
lab.y = paste("Dim ", axes[2], " (",
signif(res.mca$eig[axes[2],
2], 4), "%)", sep = "")
if (choix == "ind") {
lab.ind <- lab.var <- lab.quali.sup <- lab.ind.sup <- FALSE
if (length(label) == 1 && label == "all")
lab.ind <- lab.var <- lab.quali.sup <- lab.ind.sup <-
TRUE
if ("ind" %in% label)
lab.ind <- TRUE
if ("var" %in% label)
lab.var <- TRUE
if ("quali.sup" %in% label)
lab.quali.sup <- TRUE
if ("ind.sup" %in% label)
lab.ind.sup <- TRUE
test.invisible <- vector(length = 5)
if (!is.null(invisible)) {
test.invisible[1] <- match("ind", invisible)
test.invisible[2] <- match("var", invisible)
test.invisible[3] <- match("quanti.sup", invisible)
test.invisible[4] <- match("ind.sup", invisible)
test.invisible[5] <- match("quali.sup", invisible)
}
else test.invisible <- rep(NA, 5)
coord.var <- res.mca$var$coord[, axes]
coord.ind <- res.mca$ind$coord[, axes]
coord.ind.sup <- coord.quali.sup <- NULL
if (!is.null(res.mca$ind.sup))
coord.ind.sup <- res.mca$ind.sup$coord[, axes]
if (!is.null(res.mca$quali.sup))
coord.quali.sup <- res.mca$quali.sup$coord[, axes]
if (is.null(xlim)) {
xmin <- xmax <- 0
if (is.na(test.invisible[1]))
xmin <- min(xmin, coord.ind[, 1])
if (is.na(test.invisible[1]))
xmax <- max(xmax, coord.ind[, 1])
if (is.na(test.invisible[4]))
xmin <- min(xmin, coord.ind.sup[, 1])
if (is.na(test.invisible[4]))
xmax <- max(xmax, coord.ind.sup[, 1])
if (is.na(test.invisible[2]))
xmin <- min(xmin, coord.var[, 1])
if (is.na(test.invisible[2]))
xmax <- max(xmax, coord.var[, 1])
if (is.na(test.invisible[5]))
xmin <- min(xmin, coord.quali.sup[, 1])
if (is.na(test.invisible[5]))
xmax <- max(xmax, coord.quali.sup[, 1])
xlim <- c(xmin, xmax) * 1.2
}
else {
xmin = xlim[1]
xmax = xlim[2]
}
if (is.null(ylim)) {
ymin <- ymax <- 0
if (is.na(test.invisible[1]))
ymin <- min(ymin, coord.ind[, 2])
if (is.na(test.invisible[1]))
ymax <- max(ymax, coord.ind[, 2])
if (is.na(test.invisible[4]))
ymin <- min(ymin, coord.ind.sup[, 2])
if (is.na(test.invisible[4]))
ymax <- max(ymax, coord.ind.sup[, 2])
if (is.na(test.invisible[2]))
ymin <- min(ymin, coord.var[, 2])
if (is.na(test.invisible[2]))
ymax <- max(ymax, coord.var[, 2])
if (is.na(test.invisible[5]))
ymin <- min(ymin, coord.quali.sup[, 2])
if (is.na(test.invisible[5]))
ymax <- max(ymax, coord.quali.sup[, 2])
ylim <- c(ymin, ymax) * 1.2
}
else {
ymin = ylim[1]
ymax = ylim[2]
}
if (habillage == "quali") {
aux = 1
col.var = NULL
for (j in res.mca$call$quali) {
col.var <- c(col.var, rep(aux, nlevels(res.mca$call$X[,
j])))
aux = aux + 1
}
if (!is.null(res.mca$call$quali.sup)) {
col.quali.sup = NULL
for (j in res.mca$call$quali.sup) {
col.quali.sup <- c(col.quali.sup, rep(aux,
nlevels(res.mca$call$X[, j])))
aux = aux + 1
}
}
}
if ((habillage != "none") & (habillage !=
"quali")) {
if (!is.factor(res.mca$call$X[, habillage]))
stop("The variable ", habillage, " is not
qualitative")
col.ind <- as.numeric(as.factor(res.mca$call$X[,
habillage]))
n.mod <- nlevels(as.factor(res.mca$call$X[, habillage]))
col.ind.sup <- col.ind[res.mca$call$ind.sup]
if (!is.null(res.mca$call$ind.sup))
col.ind <- col.ind[-res.mca$call$ind.sup]
}
titre = title
if (is.null(title))
titre <- "MCA factor map"
if (is.na(test.invisible[1]) | is.na(test.invisible[2]) |
is.na(test.invisible[4]) | is.na(test.invisible[5])) {
if (new.plot)
dev.new(width = min(14, max(8, 8 * (xmax - xmin)/(ymax -
ymin))), height = 8)
plot(0, 0, main = titre, xlab = lab.x, ylab = lab.y,
xlim = xlim, ylim = ylim, col = "white", asp = 1,
cex = cex)
abline(v = 0, lty = 2, cex = cex)
abline(h = 0, lty = 2, cex = cex)
if (is.na(test.invisible[1])) {
points(coord.ind, pch = 16, col = col.ind, cex = cex)
if (lab.ind)
pointLabel(coord.ind[, 1], y = coord.ind[, 2], labels =
rownames(coord.ind),
allowSmallOverlap=FALSE,trace=FALSE,doPlot=TRUE, col =
col.ind, cex = cex)
}
if (is.na(test.invisible[2])) {
points(coord.var[, 1], y = coord.var[, 2], pch = 17,
col = col.var, cex = cex)
if (lab.var)
pointLabel(coord.var[, 1], y = coord.var[, 2], labels =
rownames(coord.var),
allowSmallOverlap=FALSE,trace=FALSE,doPlot=TRUE,col =
col.var, cex = cex)
}
if (!is.null(res.mca$quali.sup) & is.na(test.invisible[5])) {
points(coord.quali.sup[, 1], y = coord.quali.sup[,
2], pch = 17, col = col.quali.sup, cex = cex)
if (lab.quali.sup)
pointLabel(coord.quali.sup[, 1], y = coord.quali.sup[,
2], labels = rownames(coord.quali.sup),
allowSmallOverlap=FALSE,trace=FALSE,doPlot=TRUE,
col = col.quali.sup, cex = cex)
}
if (!is.null(res.mca$ind.sup) & is.na(test.invisible[4])) {
points(coord.ind.sup[, 1], y = coord.ind.sup[,
2], pch = 16, col = col.ind.sup, cex = cex)
if (lab.ind.sup)
pointLabel(coord.ind.sup[, 1], y = coord.ind.sup[,
2], labels = rownames(coord.ind.sup),
allowSmallOverlap=FALSE,trace=FALSE,doPlot=TRUE,
col = col.ind.sup, cex = cex)
}
if ((habillage != "none") & (habillage !=
"quali") &
(is.na(test.invisible[1]) | is.na(test.invisible[2])))
legend("topleft", legend = levels(res.mca$call$X[,
habillage]), text.col = 1:n.mod, cex = 0.8)
}
}
if (choix == "quanti.sup") {
if (!is.null(res.mca$quanti.sup)) {
if (new.plot)
dev.new()
if (is.null(title))
title <- "Supplementary variables on the MCA factor
map"
plot(0, 0, main = title, xlab = paste("Dim ", axes[1],
" (", signif(res.mca$eig[axes[1], 2], 4),
"%)",
sep = ""), ylab = paste("Dim ", axes[2],
" (",
signif(res.mca$eig[axes[2], 2], 4), "%)", sep =
""),
xlim = c(-1.1, 1.1), ylim = c(-1.1, 1.1), col =
"white",
asp = 1, cex = cex)
abline(v = 0, lty = 2, cex = cex)
abline(h = 0, lty = 2, cex = cex)
x.cercle <- seq(-1, 1, by = 0.01)
y.cercle <- sqrt(1 - x.cercle^2)
lines(x.cercle, y = y.cercle)
lines(x.cercle, y = -y.cercle)
for (v in 1:nrow(res.mca$quanti.sup$coord)) {
arrows(0, 0, res.mca$quanti.sup$coord[v, axes[1]],
res.mca$quanti.sup$coord[v, axes[2]], length = 0.1,
angle = 15, code = 2, col = col.quanti.sup)
if (abs(res.mca$quanti.sup$coord[v, axes[1]]) >
abs(res.mca$quanti.sup$coord[v, axes[2]])) {
if (res.mca$quanti.sup$coord[v, axes[1]] >=
0)
pos <- 4
else pos <- 2
}
else {
if (res.mca$quanti.sup$coord[v, axes[2]] >=
0)
pos <- 3
else pos <- 1
}
if ((!is.null(label)) && (label == "all" |
"quanti.sup"
%in%
label)) {
pointLabel(res.mca$quanti.sup$coord[v, axes[1]],
y = res.mca$quanti.sup$coord[v, axes[2]],
labels = rownames(res.mca$quanti.sup$coord)[v],
allowSmallOverlap=FALSE,trace=FALSE,doPlot=TRUE,cex =
cex, col = col.quanti.sup)
}
}
}
}
if (choix == "var") {
lab.var <- lab.quali.sup <- FALSE
if (length(label) == 1 && label == "all")
lab.var <- lab.quali.sup <- TRUE
if ("var" %in% label)
lab.var <- TRUE
if ("quali.sup" %in% label)
lab.quali.sup <- TRUE
test.invisible <- vector(length = 2)
if (!is.null(invisible)) {
test.invisible[1] <- match("var", invisible)
test.invisible[2] <- match("quali.sup", invisible)
}
else test.invisible <- rep(NA, 2)
if (new.plot)
dev.new()
if (is.null(title))
title <- "Variables representation"
coord.actif <- res.mca$var$eta2[, axes]
if (!is.null(res.mca$quali.sup$eta2))
coord.illu <- res.mca$quali.sup$eta2[, axes, drop = FALSE]
if (is.na(test.invisible[1])) {
plot(coord.actif, xlab = lab.x, ylab = lab.y, xlim = c(0,
1), ylim = c(0, 1), pch = 20, col = col.var,
cex = cex, main = title, cex.main = cex * 1.2,
asp = 1)
if (lab.var)
pointLabel(coord.actif[, 1], y = coord.actif[, 2],
labels = rownames(coord.actif),
allowSmallOverlap=FALSE,trace=FALSE,doPlot=TRUE, col = col.var)
}
if ((!is.null(res.mca$quali.sup$eta2)) &&
(is.na(test.invisible[2]))) {
if (!is.na(test.invisible[1]))
plot(coord.illu, xlab = lab.x, ylab = lab.y,
xlim = c(0, 1), ylim = c(0, 1), pch = 20, col =
col.quali.sup,
cex = cex, main = title, cex.main = cex * 1.2,
asp = 1)
else points(coord.illu, pch = 20, col = col.quali.sup)
if (lab.quali.sup)
pointLabel(coord.illu[, 1], y = coord.illu[, 2], labels =
rownames(coord.illu),
allowSmallOverlap=FALSE,trace=FALSE,doPlot=TRUE, col =
col.quali.sup)
}
}
## fastmod function
fastmod<-function (don, alpha = 0.05, mot_min = 2, graph = TRUE, ncp = 5,
B = 200)
{
don = as.data.frame(don)
I = nrow(don)
J = ncol(don)
for (i in 1:J) {
don[, i] = as.factor(don[, i])
}
acm = MCA(don, graph = F, ncp = ncp)
if (graph) {
layout(rbind(c(1,1,2,2,3,3),c(4,4,5,5,6,6)))
##plot.MCA2(acm, choix = "ind", invisible =
"var",new.plot=FALSE)
plot.MCA2(acm, choix = "ind", invisible =
"ind",new.plot=FALSE)
plot.MCA2(acm, choix = "ind",new.plot=FALSE)
}
X = 0.05
tab = as.matrix(acm$var$v.test)
vtestmean1 = mean(tab[, 1])
vtestmean2 = mean(tab[, 2])
vtestsd1 = sd(tab[, 1])
vtestsd2 = sd(tab[, 2])
seuil1 = vtestmean1 + vtestsd1
seuil2 = vtestmean2 + vtestsd2
modext = which(abs(tab[, 1]) >= seuil1 | abs(tab[, 2]) >=
seuil2)
tab2 = which(abs(tab[, 1]) < seuil1 & abs(tab[, 2]) < seuil2)
modmoy = tab2[sample(1:length(tab2), X * length(tab2))]
mod_kept = cbind(t(modext), t(modmoy))
res = acm
res$var$coord = acm$var$coord[mod_kept, ]
res$var$cos2 = acm$var$cos2[mod_kept, ]
res$var$contrib = acm$var$contrib[mod_kept, ]
res$var$v.test = acm$var$v.test[mod_kept, ]
if (graph) {
##plot.MCA2(res, choix = "ind", invisible =
"ind",new.plot=FALSE)
}
compte = matrix(0, I, I)
for (i in 1:J) {
for (j in 1:I) {
for (k in j:I) {
if (don[j, i] == don[k, i]) {
compte[j, k] = compte[j, k] + 1
}
}
}
}
for (i in 1:I) {
for (j in i:I) {
compte[j, i] = compte[i, j]
}
}
colnames(compte) = rownames(compte) = rownames(don)
chi2_t = function(x, y) {
obs = table(x, y)
chi = matrix(NA, length(levels(x)), length(levels(y)))
for (i in 1:length(levels(x))) {
for (j in 1:length(levels(y))) {
chi[i, j] = (obs[i, j] - (sum(obs[i, ]) * sum(obs[,
j])/length(x)))^2/(sum(obs[i, ]) * sum(obs[,
j])/length(x))
}
}
chi2 = sum(rowSums(chi))
return(chi2)
}
cramer = function(x, y) {
chi2 = chi2_t(x, y)
n = length(x)
p = length(levels(x))
q = length(levels(y))
m = min(p - 1, q - 1)
V = sqrt(chi2/(n * m))
return(V)
}
res = matrix(NA, J, J)
for (i in 1:J) {
for (j in i:J) {
res[i, j] = res[j, i] = cramer(don[, i], don[, j])
}
}
colnames(res) = rownames(res) = colnames(don)
afc = CA(res, graph = F)
ord = order(afc$row$coord[, 1])
res2 = res[ord, ]
res2 = res2[, ord]
afm = MFA(don, group = rep(1, J), type = rep("n", J), graph = F,
name.group = colnames(don))
if (graph) {
plot.MFA2(afm, choix = "group",new.plot=FALSE)
}
ordre_prod = order(acm$ind$coord[, 1])
out = matrix(NA, I, J)
lev = rep(NA, J)
for (i in 1:J) {
lev[i] = length(levels(don[, i]))
}
tdc = tab.disjonctif(don)
gp = 0
for (i in 1:J) {
conc = cbind(c(1:lev[i]), acm$var$coord[(1 + gp):(gp +
lev[i]), 1])
o = order(conc[, 2])
conc2 = cbind(conc[o, ], c(1:lev[i]))
o2 = order(conc2[, 1])
conc3 = conc2[o2, ]
out[, i] = tdc[, (1 + gp):(gp + lev[i])] %*% conc3[,
3]
gp = gp + lev[i]
}
out = data.frame(out)
for (i in 1:J) {
out[, i] = as.factor(out[, i])
}
catego_num2 = out[ordre_prod, ]
catego_num2 = catego_num2[, ord]
rownames(catego_num2) = rownames(don)[ordre_prod]
colnames(catego_num2) = colnames(don)[ord]
compte2 = compte[ordre_prod, ]
compte2 = compte2[, ordre_prod]
lev2 = as.factor(lev)
if (graph) {
##x11()
plot(lev2, main = "Number of groups formed from sorting
tasks")
}
out1 = matrix(NA, (I * J), 4)
lev = rep(NA, J)
for (i in 1:J) {
lev[i] = length(levels(don[, i]))
}
tdc = tab.disjonctif(don)
gp = 0
for (i in 1:J) {
out1[((i - 1) * I + 1):(i * I), 1:2] = tdc[, (1 + gp):(gp +
lev[i])] %*% acm$var$coord[(1 + gp):(gp + lev[i]),
1:2]
gp = gp + lev[i]
}
prod = rep(rownames(don), J)
jug = rep(colnames(don), each = I)
out2 = data.frame(out1)
out2[, 3] = prod
out2[, 4] = jug
out2[, 3] = as.factor(out2[, 3])
out2[, 4] = as.factor(out2[, 4])
d12 = acm$ind$coord[, 1:2]
pr = rownames(don)
ju = rep("0", I)
out22 = data.frame(d12, pr, ju)
colnames(out22) = colnames(out2)
out3 = rbind(out22, out2)
out3[, 3] = as.factor(out3[, 3])
out3[, 4] = as.factor(out3[, 4])
out4 = out3
out4[-c(1:I), 1] = out3[-c(1:I), 1]/sqrt(acm$eig[1, 1])
out4[-c(1:I), 2] = out3[-c(1:I), 2]/sqrt(acm$eig[2, 1])
out5 = out4
out5$moyen = out5
sim = simulation(out5, nbsimul = B)
if (graph) {
##x11()
plotellipse(sim, alpha = alpha, eig = signif(acm$eig,
4))
}
texte = matrix(NA, (I * J), 3)
texte = data.frame(texte)
texte[, 1] = rep(rownames(don), J)
texte[, 2] = rep(colnames(don), each = I)
for (i in 1:J) {
texte[((I * (i - 1)) + 1):(I * i), 3] = paste(don[, i])
}
nbp = strsplit(summary(don, maxsum = max(lev)), ":")
agg = rep(0, J * max(lev))
for (i in 1:(J * max(lev))) {
agg[i] = nbp[[i]][2]
}
agg2 = na.omit(agg)
agg2 = as.factor(agg2)
if (graph) {
##x11()
plot(agg2, main = "Number of products per group")
}
acm_call = list(X = acm$call$X, marge.col = acm$call$marge.col,
marge.row = acm$call$marge.row, ncp = acm$call$ncp, quali =
acm$call$quali,
name.group = afm$call$name.group, sim = sim)
group_afm = list(coord = afm$group$coord, cos2 = afm$group$cos2,
contrib = afm$group$contrib)
res = list(eig = acm$eig, var = acm$var, ind = acm$ind, group =
group_afm,
call = acm_call, cooccur = compte2, reord = catego_num2,
cramer = res2)
class(res) <- c("catego", "list ")
return(res)
}
Regards
Vijayan Padmanabhan
"What is expressed without proof can be denied without proof" -
Euclide.
Can you avoid printing this?
Think of the environment before printing the email.
-------------------------------------------------------------------------------
Please visit us at www.itcportal.com
******************************************************************************
This Communication is for the exclusive use of the intended recipient (s) and
shall
not attach any liability on the originator or ITC Ltd./its Subsidiaries/its
Group
Companies. If you are the addressee, the contents of this email are intended for
your
use only and it shall not be forwarded to any third party, without first
obtaining
written authorisation from the originator or ITC Ltd./its Subsidiaries/its Group
Companies. It may contain information which is confidential and legally
privileged
and the same shall not be used or dealt with by any third party in any manner
whatsoever without the specific consent of ITC Ltd./its Subsidiaries/its Group
Companies.
[[alternative HTML version deleted]]
Maybe Matching Threads
Wisdom of the Ancients
