Nothing
smacofSym <- function(delta, ndim = 2, type = c("ratio", "interval", "ordinal","mspline"),
weightmat = NULL, init = "torgerson", ties = "primary", principal = FALSE, verbose = FALSE,
relax = FALSE, modulus = 1, itmax = 1000, eps = 1e-6,
spline.degree = 2, spline.intKnots = 2) {
# delta ... dissimilarity matrix
# wghts ... weight structure. if not specified, weights is 1-structure
# p ... number of dimensions
# init ... matrix with starting values of dimension n \times p
# type ... either "ratio", "interval", "ordinal", "spline" (replaces metric)
# ties ... ties for pava (primary, secondary, tertiary)
# relax ... relaxation factor
# modulus ... modulus for nonmetric update
# itmax ... maximum number of iterations
# eps ... change in loss function
# spline.degree ... degree of the spline in case a spline transformation is chosen
# spline.intKnots ... number of interior knots for the spline in case a spline transformation is chosen
## --- sanity checks
type <- match.arg(type, c("ratio", "interval", "ordinal","mspline"), several.ok = FALSE)
diss <- delta
if ((is.matrix(diss)) || (is.data.frame(diss))) {
diss <- strucprep(diss) #if data are provided as dissimilarity matrix
attr(diss, "Labels") <- rownames(delta)
}
checkdiss(diss) ## check whether dissimilarities are all positive
p <- ndim
n <- attr(diss,"Size")
if (p > (n - 1)) stop("Maximum number of dimensions is n-1!")
nn <- n*(n-1)/2
m <- length(diss)
if (is.null(attr(diss, "Labels"))) attr(diss, "Labels") <- paste(1:n)
## --- weights
if (is.null(weightmat)) {
wgths <- initWeights(diss)
} else {
wgths <- as.dist(weightmat)
}
## --- starting values
x <- initConf(init, diss, n, p)
xstart <- x
## --- Prepare for optimal scaling
trans <- type
if (trans=="ratio"){
trans <- "none"
} else if (trans=="ordinal" & ties=="primary"){
trans <- "ordinalp"
} else if(trans=="ordinal" & ties=="secondary"){
trans <- "ordinals"
} else if(trans=="ordinal" & ties=="tertiary"){
trans <- "ordinalt"
} else if(trans=="spline"){
trans <- "mspline"
}
disobj <- transPrep(diss,trans = trans, spline.intKnots = spline.intKnots, spline.degree = spline.degree)
## Add an intercept to the spline base transformation
if (trans == "mspline") disobj$base <- cbind(rep(1, nrow(disobj$base)), disobj$base)
## dhats and missings
dhat <- normDissN(diss, wgths, 1) ## normalize to n(n-1)/2
dhat[is.na(dhat)] <- 1 ## in case of missing dissimilarities, pseudo value for dhat
if (relax) relax <- 2 else relax <- 1
w <- vmat(wgths) #matrix V of weights and unit vectors
v <- myGenInv(w) #Moore-Penrose inverse
itel <- 1 #iteration number
d <- dist(x) #Euclidean distances d(X)
lb <- sum(wgths*d*dhat, na.rm = TRUE)/sum(wgths*d^2) #denominator: normalization tr(X'VX);
x <- lb*x #modify x with lb-factor
d <- lb*d #modify d with lb-factor
sold <- sum(wgths*(dhat-d)^2, na.rm = TRUE)/nn #stress (to be minimized in repeat loop)
#--------------- begin majorization --------------------
repeat { #majorization loop
b <- bmat(dhat,wgths,d)
y <- v%*%(b%*%x) #apply Guttman transform denoted as \bar(Y) in the paper
y <- x+relax*(y-x) #n \times p matrix of Guttman transformed distances x's
e <- dist(y) #new distance matrix for Y
ssma <- sum(wgths*(dhat-e)^2) ## stress
dhat2 <- transform(e, disobj, w = wgths, normq = nn) ## dhat update
dhat <- dhat2$res
snon <- sum(wgths*(dhat-e)^2)/nn
#print out intermediate stress
if (verbose) cat("Iteration: ",formatC(itel,width=3, format="d"),
" Stress (raw):", formatC(c(snon),digits=8,width=12,format="f"),
" Difference: ", formatC(sold-snon,digits=8,width=12,format="f"),"\n")
if (((sold-snon)<eps) || (itel == itmax)) break()
x <- y #update configurations
d <- e #update configuration distances
sold <- snon #update stress
itel <- itel+1 #increase iterations
}
#------------------ end majorization ---------------
stress <- sqrt(snon) #stress normalization
colnames(y) <- paste("D",1:(dim(y)[2]),sep="")
rownames(y) <- labels(diss)
dhat <- structure(dhat, Size = n, call = quote(as.dist.default(m=b)), class = "dist", Diag = FALSE, Upper = FALSE)
attr(dhat, "Labels") <- labels(diss)
attr(e, "Labels") <- labels(diss)
dhat[is.na(diss)] <- NA ## in case of NA's
confdiss <- normDissN(e, wgths, 1) #final normalization to n(n-1)/2
## stress-per-point
spoint <- spp(dhat, dist(y), wgths)
rss <- sum(spoint$resmat[lower.tri(spoint$resmat)]) ## residual sum-of-squares
if (itel == itmax) warning("Iteration limit reached! You may want to increase the itmax argument!")
## principal axis transformation
if (principal) {
y_svd <- svd(y)
y <- y %*% y_svd$v ## rotation
}
#return configurations, configuration distances, normalized observed distances
result <- list(delta = diss, dhat = dhat, confdist = dist(y), iord = dhat2$iord.prim, conf = y, stress = stress,
spp = spoint$spp, ndim = p, weightmat = wgths, resmat = spoint$resmat, rss = rss, init = xstart, model = "Symmetric SMACOF", niter = itel, nobj = n,
type = type, call = match.call())
class(result) <- c("smacofB","smacof")
result
}
mds <- smacofSym
Any scripts or data that you put into this service are public.
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.