# below is a function that converts a distance structure to a matrix

dist2full <- function(dis)

{

        n <- attr(dis, "Size")

        full <- matrix(0, n, n)

        full[lower.tri(full)] <- dis

        full + t(full)

}


# below is a function that calculates two measures of stress

# it is fairly slow for datasets of more than 50 or so.

cmdscale.gof <- function(dis, k = 4)

{

        amat <- -0.5 * (dist2full(dis))^2       # see dist help file

        bmat <- sweep(amat, 1, apply(amat, 1, mean))

        bmat <- sweep(bmat, 2, apply(bmat, 2, mean))

        eigs <- svd(bmat, 0, 0)

        gof1 <- 1 - (cumsum(abs(eigs$d[1:k]))/sum(abs(eigs$d)))

        gof2 <- 1 - (cumsum(eigs$d[1:k]^2)/sum(eigs$d^2))

        list(gof1 = gof1, gof2 = gof2, eig = eigs$d)

}