R/smacofSphere.R In smacof: Multidimensional Scaling

Documented in smacofSphere

```smacofSphere <- function(delta, ndim = 2, type = c("ratio", "interval", "ordinal","mspline"),
algorithm = c("dual", "primal"),
weightmat = NULL, init = "torgerson",  ties = "primary", verbose = FALSE,
penalty = 100, relax = FALSE, modulus = 1, itmax = 1000, eps = 1e-6,
spline.degree = 2, spline.intKnots = 2)
{
# penalty ... penalty term kappa >0, 100 is reasonable

alg <- match.arg(algorithm, c("dual", "primal"))
#type <- match.arg(type, c("ratio", "interval", "ordinal"))
type <- match.arg(type, c("ratio", "interval", "ordinal","mspline"), several.ok = FALSE)

## --------------------------- dual algorithm --------------------------------
if (alg == "dual") {
diss <- delta
if ((is.matrix(diss)) || (is.data.frame(diss))) diss <- strucprep(diss)  #if data are provided as dissimilarity matrix
checkdiss(diss)

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)

if (is.null(weightmat)) {
wgths <- initWeights(diss)
}  else  wgths <- as.dist(weightmat)

dhat <- normDissN(diss,wgths,1)            #normalize dissimilarities
dhat[is.na(dhat)] <- 1     ## in case of missing dissimilarities, pseudo value for dhat

## --- starting values
x <- initConf(init, diss, n, p)
xstart <- x

if (relax) relax <- 2 else relax <- 1

mn <- c(1,rep(0,n))

wgths1 <- as.dist(rbind(0,cbind(0,as.matrix(wgths))))      #distances weights
wgths2 <- as.dist(outer(mn,mn,function(x,y) abs(x-y)))
dhat1 <- as.dist(rbind(0,cbind(0,as.matrix(dhat))))        #0's in the first column
dhat2 <- mean(sqrt(rowSums(x^2)))*wgths2
n1 <- attr(dhat1,"Size")
nn1 <- n1*(n1-1)/2

x <- rbind(0,x)
w <- vmat(wgths1+penalty*wgths2); v<-myGenInv(w); itel<-1;
d <- dist(x)
lb <- sum(wgths1*d*dhat1)/sum(wgths1*d^2)
x <- lb*x
d <- lb*d
sold1 <- sum(wgths1*(dhat1-d)^2)
sold2 <- sum(wgths2*(dhat2-d)^2)
sold <- sold1+penalty*sold2

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)

#---------------- begin majorization ----------------
repeat
{
b <- bmat(dhat1,wgths1,d)+penalty*bmat(dhat2,wgths2,d)
y <- v %*% (b %*% x)
y <- x+relax*(y-x)
e <- dist(y)
ssma1 <- sum(wgths1*(dhat1-e)^2)                       #stress for delta
ssma2 <- sum(wgths2*(dhat2-e)^2)                       #penalty term
ssma <- ssma1+penalty*ssma2                            #joint stress value

#---- nonmetric MDS --------
# if ((type == "ordinal") && ((itel%%modulus) == 0)) {
#   if (ties=="primary") daux<-monregP(diss,e,wgths1)
#   if (ties=="secondary") daux<-monregS(diss,e,wgths1)
#   if (ties=="tertiary") daux<-monregT(diss,e,wgths1)
#   daux<-vecAsDist(daux); dhat1<-normDissN(daux,wgths1,1)
# }
#
# ## --- interval MDS
# if (type == "interval") {
#   Amat <- cbind(1, as.vector(diss), as.vector(diss)^2)
#   daux <- nnlsPred(Amat, as.vector(e), as.vector(wgths1))\$pred
#   daux <- vecAsDist(daux)
#   dhat1 <- normDissN(daux,wgths1,1)
# }

## --- fast OS
dhat3 <- transform(e, disobj, w = wgths1, normq = nn1)  ## dhat update
dhat1 <- vecAsDist(dhat3\$res)
# ## --- end fast OS

dhat2 <- mean(e[1:n])*wgths2
snon1 <- sum(wgths1*(dhat1-e)^2)
snon2 <- sum(wgths2*(dhat2-e)^2)
snon <- snon1+penalty*snon2                            #nonmetric joint stress

if (verbose) cat("Iteration: ",formatC(itel,width=3, format="d")," Stress (not normalized): ",
formatC(c(snon),digits=8,width=12,format="f"),"\n")

if (((sold-snon)<eps) || (itel == itmax)) break()      #convergence

d <- e
sold <- snon
itel <- itel+1

if (itel == itmax) warning("Iteration limit reached! Increase itmax argument!")

}
#-------------- end majorization ---------------

colnames(y) <- paste("D",1:(dim(y)[2]),sep="")
rownames(y) <- labels(diss)
ss <- y[1,]
y <- t(apply(y, 1, function(xx) xx - ss)[,-1] )  ## correct the configurations for plotting
confdiss <- dist(y)

attr(dhat1, "Labels") <- labels(diss)
attr(e, "Labels") <- labels(diss)

stress <- sqrt(snon/nn)                   #stress normalization

e.temp <- as.dist(as.matrix(e)[,-1][-1,])      #remove dummy vector
dummyvec <- as.matrix(e)[,1][-1]
#confdiss <- normDissN(e.temp, wgths, 1)        #final normalization to n(n-1)/2

# point stress
dhat <- as.dist(as.matrix(dhat1)[-1,-1])
spoint <- spp(dhat, confdiss, wgths)
rss <- sum(spoint\$resmat[lower.tri(spoint\$resmat)])  ## residual sum-of-squares

dhat3\$iord.prim <- order(as.vector(as.dist(delta)), as.vector(as.dist(dhat)),  na.last = TRUE)

## -------------------------------------- primal algorithm ------------------------------------------
} else {
diss <- delta
if ((is.matrix(diss)) || (is.data.frame(diss))) diss <- strucprep(diss)  #if data are provided as dissimilarity matrix
checkdiss(diss)

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)

if (is.null(weightmat)) {
wgths <- initWeights(diss)
}  else  wgths <- as.dist(weightmat)

dhat <- normDissN(diss,wgths,1)            #normalize dissimilarities
dhat[is.na(dhat)] <- 1     ## in case of missing dissimilarities, pseudo value for dhat

## --- starting values
x <- initConf(init, diss, n, p)
xstart <- x

w <- vmat(wgths)
v <- myGenInv(w)
itel<-1;

x <- x/sqrt(rowSums(x^2))
d <- dist(x)                           #distance computation (to be extended with geodesic)

lb <- sum(wgths*d*dhat)/sum(wgths*d^2)
x <- lb*x
d <- lb*d
sold <- sum(wgths*(dhat-d)^2)          #initial stress

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)

#------------------------- begin majorization ---------------------------------
repeat {
b <- bmat(dhat,wgths,d)
y <- v%*%b%*%x                       #Guttman transform
y <- sphereProj(y,w)                 #projection on the sphere: slow!

e <- dist(y)                         #extension: distances for Y to be enhanced with geodesics)
ssma <- sum(wgths*(dhat-e)^2)        #metric stress

# if (type == "ordinal") {                       #nonmetric versions
#   if ((itel%%modulus) == 0) {
#     if (ties=="primary") daux <- monregP(diss,e,wgths)
#     if (ties=="secondary") daux <- monregS(diss,e,wgths)
#     if (ties=="tertiary") daux <- monregT(diss,e,wgths)
#     daux <- vecAsDist(daux)
#     dhat <- normDissN(daux,wgths,1)
#   }
# }
#
# ## --- interval MDS
# if (type == "interval") {
#   Amat <- cbind(1, as.vector(diss), as.vector(diss)^2)
#   daux <- nnlsPred(Amat, as.vector(e), as.vector(wgths))\$pred
#   daux <- vecAsDist(daux)
#   dhat <- normDissN(daux,wgths,1)
# }

## --- fast OS
dhat3 <- transform(e, disobj, w = wgths, normq = nn)  ## dhat update
dhat <- dhat3\$res
## --- end fast OS
#
snon <- sum(wgths*(dhat-e)^2)        #nonmetric stress
if (verbose) cat("Iteration: ",formatC(itel,width=3, format="d")," Stress (not normalized): ",
formatC(c(snon),digits=8,width=12,format="f"),"\n")
if (((sold-snon)<eps) || (itel == itmax)) break()

d <- e
sold <- snon
itel <- itel+1
}
#----------------------------- end majorization -------------------------------

colnames(y) <- paste("D",1:(dim(y)[2]),sep="")
rownames(y) <- labels(diss)
attr(dhat, "Labels") <- labels(diss)
attr(e, "Labels") <- labels(diss)
dhat[is.na(diss)] <- NA                 ## in case of NA's

stress <- sqrt(snon/nn)                   #stress normalization

confdiss <- normDissN(e, wgths, 1)        #final normalization to n(n-1)/2

# point stress
spoint <- spp(dhat, dist(y), wgths)
rss <- sum(spoint\$resmat[lower.tri(spoint\$resmat)])  ## residual sum-of-squares

dummyvec <- NA
if (itel == itmax) warning("Iteration limit reached! Increase itmax argument!")

}

result <- list(delta = delta, dhat = dhat, confdist = dist(y), iord = dhat3\$iord.prim, conf = y,
stress = stress, spp = spoint\$spp, ndim = p, weightmat = wgths, resmat = spoint\$resmat, rss = rss,
dummyvec = dummyvec, init = xstart,
model = "Spherical SMACOF", niter = itel, nobj = n, type = type,
algorithm = alg, call = match.call())

class(result) <- c("smacofSP", "smacof")
result
}
```

Try the smacof package in your browser

Any scripts or data that you put into this service are public.

smacof documentation built on July 21, 2021, 3 a.m.