"dudi.acm" <- function (df, row.w = rep(1, nrow(df)), scannf = TRUE, nf = 2) {
if (!all(unlist(lapply(df, is.factor))))
stop("All variables must be factors")
df <- as.data.frame(df)
X <- acm.disjonctif(df)
lig <- nrow(X)
col <- ncol(X)
var <- ncol(df)
if (length(row.w) != lig)
stop("Non convenient row weights")
if (any(row.w < 0))
stop("row weight < 0")
row.w <- row.w/sum(row.w)
col.w <- apply(X, 2, function(x) sum(x*row.w))
if (any(col.w == 0))
stop("One category with null weight")
X <- t(t(X)/col.w) - 1
col.w <- col.w/var
X <- as.dudi(data.frame(X), col.w, row.w, scannf = scannf,
nf = nf, call = match.call(), type = "acm")
rcor <- matrix(0, ncol(df), X$nf)
rcor <- row(rcor) + 0 + (0+1i) * col(rcor)
floc <- function(x) {
i <- Re(x)
j <- Im(x)
x <- X$l1[, j] * X$lw
qual <- df[, i]
poicla <- unlist(tapply(X$lw, qual, sum))
z <- unlist(tapply(x, qual, sum))/poicla
return(sum(poicla * z * z))
}
rcor <- apply(rcor, c(1, 2), floc)
rcor <- data.frame(rcor)
row.names(rcor) <- names(df)
names(rcor) <- names(X$l1)
X$cr <- rcor
return(X)
}
"boxplot.acm" <- function (x, xax = 1, ...) {
# correction d'un bug par P. Cornillon 29/10/2004
if (!inherits(x, "acm"))
stop("Object of class 'acm' expected")
if ((xax < 1) || (xax > x$nf))
stop("non convenient axe number")
def.par <- par(no.readonly = TRUE)
on.exit(par(def.par))
oritab <- eval.parent(as.list(x$call)[[2]])
nvar <- ncol(oritab)
if (nvar <= 7)
sco.boxplot(x$l1[, xax], oritab[, 1:nvar], clabel = 1)
else if (nvar <= 14) {
par(mfrow = c(1, 2))
sco.boxplot(x$l1[, xax], oritab[, 1:(nvar%/%2)], clabel = 1.3)
sco.boxplot(x$l1[, xax], oritab[, (nvar%/%2 + 1):nvar],
clabel = 1.3)
}
else {
par(mfrow = c(1, 3))
if ((a0 <- nvar%/%3) < nvar/3)
a0 <- a0 + 1
sco.boxplot(x$l1[, xax], oritab[, 1:a0], clabel = 1.6)
sco.boxplot(x$l1[, xax], oritab[, (a0 + 1):(2 * a0)],
clabel = 1.6)
sco.boxplot(x$l1[, xax], oritab[, (2 * a0 + 1):nvar],
clabel = 1.6)
}
}
"acm.burt" <- function (df1, df2, counts = rep(1, nrow(df1))) {
if (!all(unlist(lapply(df1, is.factor))))
stop("All variables must be factors")
if (!all(unlist(lapply(df2, is.factor))))
stop("All variables must be factors")
if (nrow(df1) != nrow(df2))
stop("non convenient row numbers")
if (length(counts) != nrow(df2))
stop("non convenient row numbers")
g1 <- acm.disjonctif(df1)
g1 <- g1 * counts
g2 <- acm.disjonctif(df2)
burt <- as.matrix(t(g1)) %*% as.matrix(g2)
burt <- data.frame(burt)
names(burt) <- names(g2)
row.names(burt) <- names(g1)
return(burt)
}
"acm.disjonctif" <- function (df) {
acm.util.df <- function(i) {
cl <- df[,i]
cha <- names(df)[i]
n <- length(cl)
cl <- as.factor(cl)
x <- matrix(0, n, length(levels(cl)))
x[(1:n) + n * (unclass(cl) - 1)] <- 1
dimnames(x) <- list(row.names(df), paste(cha,levels(cl),sep="."))
return(x)
}
G <- lapply(1:ncol(df), acm.util.df)
G <- data.frame (G, check.names = FALSE)
return(G)
}
fac2disj<- function(fac, drop = FALSE) {
## Returns the disjunctive table corrseponding to a factor
n <- length(fac)
fac <- as.factor(fac)
if(drop)
fac <- factor(fac)
x <- matrix(0, n, nlevels(fac))
x[(1:n) + n * (unclass(fac) - 1)] <- 1
dimnames(x) <- list(names(fac), as.character(levels(fac)))
return(data.frame(x, check.names = FALSE))
}
scatterutil.eti <- function (x, y, label, clabel, boxes = TRUE, coul = rep(1, length(x)),
horizontal = TRUE, bg = "white")
{
if (length(label) == 0)
return(invisible())
if (is.null(label))
return(invisible())
if (any(label == ""))
return(invisible())
cex0 <- par("cex") * clabel
for (i in 1:(length(x))) {
cha <- as.character(label[i])
cha <- paste(" ", cha, " ", sep = "")
x1 <- x[i]
y1 <- y[i]
xh <- strwidth(cha, cex = cex0)
yh <- strheight(cha, cex = cex0) * 5/3
if (!horizontal) {
tmp <- scatterutil.convrot90(xh, yh)
xh <- tmp[1]
yh <- tmp[2]
}
if (boxes) {
rect(x1 - xh/2, y1 - yh/2, x1 + xh/2, y1 + yh/2,
col = bg, border = coul[i])
}
if (horizontal) {
text(x1, y1, cha, cex = cex0, col = coul[i])
}
else {
text(x1, y1, cha, cex = cex0, col = coul[i], srt = 90)
}
}
}
scatterutil.sub <- function (cha, csub, possub = "bottomleft")
{
cha <- as.character(cha)
if (length(cha) == 0)
return(invisible())
if (is.null(cha))
return(invisible())
if (is.na(cha))
return(invisible())
if (any(cha == ""))
return(invisible())
if (csub == 0)
return(invisible())
cex0 <- par("cex") * csub
cha <- paste(" ", cha, " ", sep = "")
xh <- strwidth(cha, cex = cex0)
yh <- strheight(cha, cex = cex0) * 5/3
if (possub == "bottomleft") {
x1 <- par("usr")[1]
y1 <- par("usr")[3]
rect(x1, y1, x1 + xh, y1 + yh, col = "white", border = 0)
text(x1 + xh/2, y1 + yh/2, cha, cex = cex0)
}
else if (possub == "topleft") {
x1 <- par("usr")[1]
y1 <- par("usr")[4]
rect(x1, y1, x1 + xh, y1 - yh, col = "white", border = 0)
text(x1 + xh/2, y1 - yh/2, cha, cex = cex0)
}
else if (possub == "bottomright") {
x1 <- par("usr")[2]
y1 <- par("usr")[3]
rect(x1, y1, x1 - xh, y1 + yh, col = "white", border = 0)
text(x1 - xh/2, y1 + yh/2, cha, cex = cex0)
}
else if (possub == "topright") {
x1 <- par("usr")[2]
y1 <- par("usr")[4]
rect(x1, y1, x1 - xh, y1 - yh, col = "white", border = 0)
text(x1 - xh/2, y1 - yh/2, cha, cex = cex0)
}
}
scatterutil.grid <- function (cgrid)
{
col <- "lightgray"
lty <- 1
xaxp <- par("xaxp")
ax <- (xaxp[2] - xaxp[1])/xaxp[3]
yaxp <- par("yaxp")
ay <- (yaxp[2] - yaxp[1])/yaxp[3]
a <- min(ax, ay)
v0 <- seq(xaxp[1], xaxp[2], by = a)
h0 <- seq(yaxp[1], yaxp[2], by = a)
abline(v = v0, col = col, lty = lty)
abline(h = h0, col = col, lty = lty)
if (cgrid <= 0)
return(invisible())
cha <- paste(" d = ", a, " ", sep = "")
cex0 <- par("cex") * cgrid
xh <- strwidth(cha, cex = cex0)
yh <- strheight(cha, cex = cex0) * 5/3
x1 <- par("usr")[2]
y1 <- par("usr")[4]
rect(x1 - xh, y1 - yh, x1 + xh, y1 + yh, col = "white", border = 0)
text(x1 - xh/2, y1 - yh/2, cha, cex = cex0)
}
scatterutil.base <- function (dfxy, xax, yax, xlim, ylim, grid, addaxes, cgrid, include.origin,
origin, sub, csub, possub, pixmap, contour, area, add.plot)
{
df <- data.frame(dfxy)
if (!is.data.frame(df))
stop("Non convenient selection for df")
if ((xax < 1) || (xax > ncol(df)))
stop("Non convenient selection for xax")
if ((yax < 1) || (yax > ncol(df)))
stop("Non convenient selection for yax")
x <- df[, xax]
y <- df[, yax]
if (is.null(xlim)) {
x1 <- x
if (include.origin)
x1 <- c(x1, origin[1])
x1 <- c(x1 - diff(range(x1)/10), x1 + diff(range(x1))/10)
xlim <- range(x1)
}
if (is.null(ylim)) {
y1 <- y
if (include.origin)
y1 <- c(y1, origin[2])
y1 <- c(y1 - diff(range(y1)/10), y1 + diff(range(y1))/10)
ylim <- range(y1)
}
if (!is.null(pixmap)) {
if (is.null(class(pixmap)))
pixmap <- NULL
if (is.na(charmatch("pixmap", class(pixmap))))
pixmap <- NULL
}
if (!is.null(contour)) {
if (!is.data.frame(contour))
contour <- NULL
if (ncol(contour) != 4)
contour <- NULL
}
if (!is.null(area)) {
if (!is.data.frame(area))
area <- NULL
if (!is.factor(area[, 1]))
area <- NULL
if (ncol(area) < 3)
area <- NULL
}
if (!add.plot)
plot.default(0, 0, type = "n", asp = 1, xlab = "", ylab = "",
xaxt = "n", yaxt = "n", xlim = xlim, ylim = ylim,
xaxs = "i", yaxs = "i", frame.plot = FALSE)
if (!is.null(pixmap)) {
plot(pixmap, add = TRUE)
}
if (!is.null(contour)) {
apply(contour, 1, function(x) segments(x[1], x[2], x[3],
x[4], lwd = 1))
}
if (grid & !add.plot)
scatterutil.grid(cgrid)
if (addaxes & !add.plot)
abline(h = 0, v = 0, lty = 1, lwd=0.4)
if (!is.null(area)) {
nlev <- nlevels(area[, 1])
x1 <- area[, 2]
x2 <- area[, 3]
for (i in 1:nlev) {
lev <- levels(area[, 1])[i]
a1 <- x1[area[, 1] == lev]
a2 <- x2[area[, 1] == lev]
polygon(a1, a2)
}
}
if (csub > 0)
scatterutil.sub(sub, csub, possub)
return(list(x = x, y = y))
}
scatterutil.star <- function (x, y, z, cstar, coul = rep(1, length(x)), star.lwd=0.5)
{
z <- z/sum(z)
x1 <- sum(x * z)
y1 <- sum(y * z)
for (i in which(z > 0)) {
hx <- cstar * (x[i] - x1)
hy <- cstar * (y[i] - y1)
segments(x1, y1, x1 + hx, y1 + hy, col = coul, lwd=star.lwd)
}
}
scatterutil.ellipse <- function (x, y, z, cellipse, axesell, coul = rep(1, length(x)), ellipse.lwd=2, axesell.lwd=1, star.lwd=0.5)
{
if (any(is.na(z)))
return(invisible())
if (sum(z * z) == 0)
return(invisible())
util.ellipse <- function(mx, my, vx, cxy, vy, coeff) {
lig <- 100
epsi <- 1e-10
x <- 0
y <- 0
if (vx < 0)
vx <- 0
if (vy < 0)
vy <- 0
if (vx == 0 && vy == 0)
return(NULL)
delta <- (vx - vy) * (vx - vy) + 4 * cxy * cxy
delta <- sqrt(delta)
l1 <- (vx + vy + delta)/2
l2 <- vx + vy - l1
if (l1 < 0)
l1 <- 0
if (l2 < 0)
l2 <- 0
l1 <- sqrt(l1)
l2 <- sqrt(l2)
test <- 0
if (vx == 0) {
a0 <- 0
b0 <- 1
test <- 1
}
if ((vy == 0) && (test == 0)) {
a0 <- 1
b0 <- 0
test <- 1
}
if (((abs(cxy)) < epsi) && (test == 0)) {
a0 <- 1
b0 <- 0
test <- 1
}
if (test == 0) {
a0 <- 1
b0 <- (l1 * l1 - vx)/cxy
norm <- sqrt(a0 * a0 + b0 * b0)
a0 <- a0/norm
b0 <- b0/norm
}
a1 <- 2 * pi/lig
c11 <- coeff * a0 * l1
c12 <- (-coeff) * b0 * l2
c21 <- coeff * b0 * l1
c22 <- coeff * a0 * l2
angle <- 0
for (i in 1:lig) {
cosinus <- cos(angle)
sinus <- sin(angle)
x[i] <- mx + c11 * cosinus + c12 * sinus
y[i] <- my + c21 * cosinus + c22 * sinus
angle <- angle + a1
}
return(list(x = x, y = y, seg1 = c(mx + c11, my + c21,
mx - c11, my - c21), seg2 = c(mx + c12, my + c22,
mx - c12, my - c22)))
}
z <- z/sum(z)
m1 <- sum(x * z)
m2 <- sum(y * z)
v1 <- sum((x - m1) * (x - m1) * z)
v2 <- sum((y - m2) * (y - m2) * z)
cxy <- sum((x - m1) * (y - m2) * z)
ell <- util.ellipse(m1, m2, v1, cxy, v2, cellipse)
if (is.null(ell))
return(invisible())
polygon(ell$x, ell$y, border = coul, lwd=ellipse.lwd)
if (axesell)
segments(ell$seg1[1], ell$seg1[2], ell$seg1[3], ell$seg1[4],
lty = 2, col = coul, lwd=axesell.lwd)
if (axesell)
segments(ell$seg2[1], ell$seg2[2], ell$seg2[3], ell$seg2[4],
lty = 2, col = coul, lwd=axesell.lwd)
}
s.class2 <- function (dfxy, fac, wt = rep(1, length(fac)), xax = 1, yax = 2,
cstar = 1, cellipse = 1.5, axesell = TRUE, label = levels(fac),
clabel = 1, cpoint = 1, pch = 20, col = rep(1, length(levels(fac))),
xlim = NULL, ylim = NULL, grid = TRUE, addaxes = TRUE, origin = c(0,
0), include.origin = TRUE, sub = "", csub = 1, possub = "bottomleft",
cgrid = 1, pixmap = NULL, contour = NULL, area = NULL, add.plot = FALSE,
ellipse.lwd=2, axesell.lwd=1, star.lwd=0.5, axes=TRUE, xlab='NMDS1', ylab='NMDS2')
{
opar <- par(mar = par("mar"))
if (axes) {
par(mar = c(0.1, 0.1, 0.1, 0.1) * 40)
} else {
par(mar = c(0.1, 0.1, 0.1, 0.1))
}
on.exit(par(opar))
dfxy <- data.frame(dfxy)
if (!is.data.frame(dfxy))
stop("Non convenient selection for dfxy")
if (any(is.na(dfxy)))
stop("NA non implemented")
if (!is.factor(fac))
stop("factor expected for fac")
dfdistri <- fac2disj(fac) * wt
coul <- col
w1 <- unlist(lapply(dfdistri, sum))
dfdistri <- t(t(dfdistri)/w1)
coox <- as.matrix(t(dfdistri)) %*% dfxy[, xax]
cooy <- as.matrix(t(dfdistri)) %*% dfxy[, yax]
if (nrow(dfxy) != nrow(dfdistri))
stop(paste("Non equal row numbers", nrow(dfxy), nrow(dfdistri)))
coo <- scatterutil.base(dfxy = dfxy, xax = xax, yax = yax,
xlim = xlim, ylim = ylim, grid = grid, addaxes = addaxes,
cgrid = cgrid, include.origin = include.origin, origin = origin,
sub = sub, csub = csub, possub = possub, pixmap = pixmap,
contour = contour, area = area, add.plot = add.plot)
if (cpoint > 0)
for (i in 1:ncol(dfdistri)) {
pch <- rep(pch, length = nrow(dfxy))
points(coo$x[dfdistri[, i] > 0], coo$y[dfdistri[,
i] > 0], pch = pch[dfdistri[, i] > 0], cex = par("cex") *
cpoint, col = coul[i])
}
if (cstar > 0)
for (i in 1:ncol(dfdistri)) {
scatterutil.star(coo$x, coo$y, dfdistri[, i], cstar = cstar,
coul[i], star.lwd=star.lwd)
}
if (cellipse > 0)
for (i in 1:ncol(dfdistri)) {
scatterutil.ellipse(coo$x, coo$y, dfdistri[, i],
cellipse = cellipse, axesell = axesell, coul[i], ellipse.lwd=ellipse.lwd, axesell.lwd=axesell.lwd)
}
if (clabel > 0)
scatterutil.eti(coox, cooy, label, clabel, coul = col)
box(lwd=0.4)
if (axes) {
title(xlab=xlab, ylab=ylab)
axis(side=1, lwd=0.4)
axis(side=2, las=1, lwd=0.4)
}
invisible(match.call())
}
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.