Nothing
R/lpostmds.R
In cops: Cluster Optimized Proximity Scaling
Defines functions
localmds
#' Local Power Stress MDS either as self-organizing or not
#'
#' A local version of power stress MDS (POST-)MDS by quasi-majorization with ratio, interval and spline optimal scaling for dissimilarities and optional power transformations. This marries POST-MDS with the local multidimensional scaling idea of Venna and Karski (2006). It works like this: Find an initial weight matrix where w_ij(X^0)=0 if d_ij(X^0)> tau and delta_ij > tau, w_ij(X^0)= rho if d_ij(X^0)> tau and delta_ij <= tau, w_ij(X^0)= 1-rho if d_ij(X^0)<= tau and delta_ij > tau, w_ij(X^0)=1 if d_ij(X^0)<= tau and delta_ij <= tau.
#'
#' Then fit a POST-MDS with these weights. Then in each successive iteration step, the weightmat is recalculated so that w_ij(X^(n+1))=0 if d_ij(X^(n+1))> tau and delta_ij > tau, w_ij(X^(n+1))= rho if d_ij(X^(n+1))> tau and delta_ij <= tau, w_ij(X^(n+1))= 1-rho if d_ij(X^(n+1))<= tau and delta_ij > tau, w_ij(X^(n+1))=1 if d_ij(X^(n+1))<= tau and delta_ij <= tau.
#'
#' If rho = 0, this is eCLPCA/SPMDS see \code{\link[smacofx]{spmds}}.
#'
#' There is a wrappers 'lmdsvk' where the exponents are 1 and 'type=ratio'. The neighborhood parameter tau is kept fixed in 'lpostmds' and 'lmdsvk'. The functions 'so_lpostmds', 'so_lmdsvk' implement a self-organising principle, where the model is repeatedly fitted for a decreasing sequence of taus.
#'
#' @param delta dist object or a symmetric, numeric data.frame or matrix of distances
#' @param lambda exponent of the power transformation of the dissimilarities; defaults to 1., which is also the setup of 'lmdsvk'
#' @param kappa exponent of the power transformation of the fitted distances; defaults to 1, which is also the setup of 'lmdsvk'.
#' @param nu exponent of the power of the weighting matrix; defaults to 1 which is also the setup for 'lmdsvk'.
#' @param rho is the continuity parameter and must be a value in [0,1]. This trades-off trustworthiness with continuity of the mapping. rho = 0 will give a weighting as in curvilinear component analysis. If rho=1 it will give an MDS with all inputs > tau treated as missing.
#' @param tau the boundary/neighbourhood parameter(s) (called sigma in the original paper). For 'lpostmds' and 'lmdsvk' it is supposed to be a numeric scalar (if a sequence is supplied the maximum is taken as tau) and all the transformed fitted distances exceeding tau are set to 0 via the weightmat (assignment can change between iterations). It defaults to the 25\% quantile of the explictly normalized (power transformed) delta. We recommend not to use the defaults!
#' @param type what type of MDS to fit. Currently one of "ratio", "interval" or "mspline". Default is "ratio".
#' @param spline.degree Degree of the spline for ‘mspline’ MDS type
#' @param spline.intKnots Number of interior knots of the spline for ‘mspline’ MDS type
#' @param weightmat a matrix of finite weights.
#' @param init starting configuration. If NULL (default) we fit a full rstress model.
#' @param ndim dimension of the configuration; defaults to 2
#' @param acc numeric accuracy of the iteration. Default is 1e-8.
#' @param itmax maximum number of iterations. Default is 10000.
#' @param verbose should fitting information be printed; if > 0 then yes
#' @param principal If 'TRUE', principal axis transformation is applied to the final configuration
#' @param traceIt save the iteration progress in a vector (stress values)
#'
#' @return a 'smacofP' object (inheriting from 'smacofB', see \code{\link[smacof]{smacofSym}}). It is a list with the components
#' \itemize{
#' \item delta: Observed, untransformed dissimilarities
#' \item tdelta: Observed explicitly transformed dissimilarities, normalized
#' \item dhat: Explicitly transformed dissimilarities (dhats), optimally scaled and normalized
#' \item confdist: Transformed configuration distances
#' \item conf: Matrix of fitted configuration
#' \item stress: Default stress (stress1; sqrt of explicitly normalized stress)
#' \item spp: Stress per point
#' \item ndim: Number of dimensions
#' \item model: Name of smacof model
#' \item niter: Number of iterations
#' \item nobj: Number of objects
#' \item type: Type of MDS model
#' \item weightmat: weighting matrix as supplied
#' \item stress.m: Default stress (square of stress1, explicitly normalized stress)
#' \item tweightmat: transformed weighting matrix; it is weightmat^nu and containing all the 0s for the distances/dissimilartities set to 0.
#' \item trace: if 'traceIt=TRUE' a vector with the iteration progress
#'}
#'
#'
#' @details
#' The solution is found by "quasi-majorization", which means that the majorization is only real majorization once the weightmat no longer changes. This typically happens after a few iterations. Due to that it can be that in the beginning the stress may not decrease monotonically and that there's a chance it might never.
#'
#' If tau is too small it may happen that all distances for one i to all j are zero and then there will be an error, so make sure to set a larger tau.
#'
#'
#'
#' @importFrom stats dist as.dist quantile
#' @importFrom smacof transform transPrep
#'
#' @examples
#' dis<-smacof::morse
#' res<-lpostmds(dis, type = "interval", rho = 0.5, kappa = 2, lambda = 2,tau = 0.3, itmax = 100) #use higher itmax
#' res
#' summary(res)
#' plot(res)
#'
#' ##which d_ij(X)^kappa exceeded tau at convergence (i.e., have been set to 0)?
#' res$tweightmat
#'
#' # We use Quasi-Majorization
#' res$trace
#'
#'
#'
#' @export
lpostmds <- function (delta, lambda=1, kappa=1, nu=1, tau = NULL, rho = 0, type="ratio", weightmat=1-diag(nrow(delta)), init=NULL, ndim = 2, acc= 1e-8, itmax = 10000, verbose = FALSE, principal=FALSE, spline.degree = 2, spline.intKnots = 2, traceIt=FALSE) {
## TODO: Should we allow for different tau for input and output?
if(inherits(delta,"dist") || is.data.frame(delta)) delta <- as.matrix(delta)
if(!isSymmetric(delta)) stop("delta is not symmetric.\n")
if(inherits(weightmat,"dist") || is.data.frame(weightmat)) weightmat <- as.matrix(weightmat)
if(!isSymmetric(weightmat)) stop("weightmat is not symmetric.\n")
r <- kappa/2
if(rho < 0 && rho > 1) stop("rho must be in the interval [0, 1].")
## -- Setup for MDS type
if(missing(type)) type <- "ratio"
type <- match.arg(type, c("ratio", "interval", "mspline","spline"),several.ok = FALSE) #no ordinal allowed
trans <- type
typo <- type
if (trans=="ratio"){
trans <- "none"
}
# else if (trans=="ordinal" & ties=="primary"){
# trans <- "ordinalp"
# typo <- "ordinal (primary)"
# } else if(trans=="ordinal" & ties=="secondary"){
# trans <- "ordinals"
# typo <- "ordinal (secondary)"
# } else if(trans=="ordinal" & ties=="tertiary"){
# trans <- "ordinalt"
# typo <- "ordinal (tertiary)"
# } else
if(trans=="spline"){
trans <- "mspline"
type <- "mspline"
typo <- "mspline"
}
if(type %in% c("mspline", "spline")) lambda <- 1 #We dont allow powers for dissimilarities in nonmetric and spline MDS
n <- nrow (delta)
normi <- 0.5
##normi <- n #if normi=n we can use the iord structure in plot.smacofP
## but the problem is we don't get the correct stress then anymore.
p <- ndim
if (p > (n - 1)) stop("Maximum number of dimensions is n-1!")
if(is.null(rownames(delta))) rownames(delta) <- 1:n
labos <- rownames(delta) #labels
deltaorig <- delta
delta <- delta^lambda ## Transfrom the deltas
weightmato <- weightmat #3 so we can retunr this as supplied
weightmat <- weightmat^nu ##transform the weightmat
weightmat[!is.finite(weightmat)] <- 0 ## all infinite or NA weights get a weight of 0
delta <- delta / enorm (delta, weightmat) ## explicit normalization of input
if(is.null(tau)) tau <- stats::quantile(delta,0.25) ## tau default
if(length(tau)>1)
{
warning("Supplied tau is of length > 1. The min(tau) was used as tau.")
tau <- min(tau)
}
if(tau<=0) stop("tau must be positive.")
if(verbose>0) cat(paste("Fitting",type,"Local POST-MDS with lambda=",lambda, "kappa=",kappa,"nu=",nu, "rho=",rho,"and tau=",tau,"\n"))
disobj <- smacof::transPrep(as.dist(delta), 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)
deltaold <- delta
itel <- 1
##Starting Configs
xold <- init
# if(is.null(init)) xold <- smacof::torgerson (delta,r=kappa/2,type=type,ties=ties,weightmat=weightmat,ndim=ndim,init=init,itmax=itmax,principal=principal)$conf
if(is.null(init)) xold <- smacof::torgerson(delta, p = p)
xstart <- xold
xold <- xold / enorm (xold) #explicit normalization of output
nn <- diag (n)
dold <- sqdist (xold) #squared distances
doldpow <- mkPower(dold, kappa / 2)# distances^kappa
## --- Make local weightmat
trustwn <- (1-rho)*(doldpow <= tau) ##CCA penalty
contin <- rho*(deltaold <= tau)
tcmat <- trustwn+contin ##
lweightmat <- weightmat*tcmat
diag(lweightmat) <- 0 ## sanitizing
lweightmat[!is.finite(lweightmat)] <- 0 ## all infinite or NA weights get a weight of 0
##first optimal scaling
eold <- as.dist(mkPower(dold,r))
dhat <- smacof::transform(eold, disobj, w = as.dist(lweightmat), normq = normi)
dhatt <- dhat$res
dhatd <- structure(dhatt, Size = n, call = quote(as.dist.default(m=b)), class = "dist", Diag = FALSE, Upper = FALSE)
delta <- as.matrix(dhatd)
## first stress
rold <- sum (lweightmat * delta * mkPower (dold, r))
nold <- sum (lweightmat * mkPower (dold, 2 * r))
aold <- rold / nold
sold <- 1 - 2 * aold * rold + (aold ^ 2) * nold
tracev <- NULL
if(isTRUE(traceIt)) tracev <- rep(NA,itmax)
## Optimizing
## Loop Logic for iteration i: X^(i)|W^(i-1) -> W^(i)|X^(i) -> stress^(i) -> X^(i+1)|W^(i) -> etc.
repeat {
if(tau <= min(doldpow[lower.tri(doldpow)]) && min(delta[lower.tri(delta)])) stop("Current tau is lower than the smallest fitted distance and the smallest delta (therefore all distances are set to 0). Increase tau.")
p1 <- mkPower (dold, r - 1)
p2 <- mkPower (dold, (2 * r) - 1)
by <- mkBmat (lweightmat * delta * p1)
cy <- mkBmat (lweightmat * p2)
ga <- 2 * sum (lweightmat * p2)
be <- (2 * r - 1) * (2 ^ r) * sum (lweightmat * delta)
de <- (4 * r - 1) * (4 ^ r) * sum (lweightmat)
if (r >= 0.5) {
my <- by - aold * (cy - de * nn)
}
if (r < 0.5) {
my <- (by - be * nn) - aold * (cy - ga * nn)
}
xnew <- my %*% xold
xnew <- xnew / enorm (xnew)
dnew <- sqdist (xnew)
dnewpow <- mkPower(dnew,kappa/2)
### We now reset the weight matrix, so it is possible that a w_{ij} can change from 0 to >0 again
weightmat <- weightmato^nu #the supplied weightmat and carrying the nu with us
weightmat[!is.finite(weightmat)] <- 0
## --- Make local weightmat
trustwn <- (1-rho)*(dnewpow <= tau) ##trustworthiness
contin <- rho*(delta <= tau) #continuity; delta are now the dhats
tcmat <- trustwn+contin ##convex combo
lweightmat <- weightmat*tcmat #multiply with given weights
diag(lweightmat) <- 0 ## sanitizing
lweightmat[!is.finite(lweightmat)] <- 0 ## all infinite or NA weights get a weight of 0
##optimal scaling
e <- as.dist(mkPower(dnew,r)) #I need the dist(x) here for interval
dhat2 <- smacof::transform(e, disobj, w = as.dist(lweightmat), normq = normi) ## dhat update
dhatt <- dhat2$res
dhatd <- structure(dhatt, Size = n, call = quote(as.dist.default(m=b)), class = "dist", Diag = FALSE, Upper = FALSE)
delta <- as.matrix(dhatd)
rnew <- sum (lweightmat * delta * mkPower (dnew, r))
nnew <- sum (lweightmat * mkPower (dnew, 2 * r))
anew <- rnew / nnew
snew <- 1 - 2 * anew * rnew + (anew ^ 2) * nnew
if(is.na(snew)) #if there are issues with the values
{
snew <- sold
dnew <- dold
anew <- aold
xnew <- xold
}
if (verbose>2) {
cat (
formatC (itel, width = 4, format = "d"),
formatC (
sold,
digits = 10,
width = 13,
format = "f"
),
formatC (
snew,
digits = 10,
width = 13,
format = "f"
),
"\n"
)
}
if(isTRUE(traceIt)) tracev[itel] <- snew
if ((itel == itmax) || (abs(sold - snew) < acc)) #new
{
if(sold < snew) itel <- itel-1
snew <- min(sold,snew) #to make sure we use the lowest stress as it is quasi majorization and sold might be < snew at this point
break ()
}
itel <- itel + 1
xold <- xnew
dold <- dnew
sold <- snew
aold <- anew
}
xnew <- xnew/enorm(xnew)
## relabeling as they were removed in the optimal scaling
rownames(delta) <- labos
attr(xnew,"dimnames")[[1]] <- rownames(delta)
attr(xnew,"dimnames")[[2]] <- paste("D",1:p,sep="")
doutm <- mkPower(sqdist(xnew),r)
deltam <- delta
#delta <- structure(delta, Size = n, call = quote(as.dist.default(m=b)),
# class = "dist", Diag = FALSE, Upper = FALSE)
delta <- stats::as.dist(delta)
deltaorig <- stats::as.dist(deltaorig)
deltaold <- stats::as.dist(deltaold)
#doute <- doutm/enorm(doutm) #this is an issue here!
#doute <- stats::as.dist(doute)
dout <- stats::as.dist(doutm)
weightmatm <-lweightmat ##the local weightmat used
weightmat <- stats::as.dist(weightmatm)
spoint <- spp(delta, dout, weightmat)
resmat<-spoint$resmat
rss <- sum(spoint$resmat[lower.tri(spoint$resmat)])
spp <- spoint$spp
if (verbose > 0 && itel == itmax) warning("Iteration limit reached! You may want to increase the itmax argument!")
if (principal) {
xnew_svd <- svd(xnew)
xnew <- xnew %*% xnew_svd$v
}
if(isTRUE(traceIt)) tracev <- tracev[!is.na(tracev)]
if(verbose>1) cat("*** Stress:",snew, "; Stress-1 (default reported):",sqrt(snew),"\n")
#delta is input delta, tdelta is input delta with explicit transformation and normalized, dhat is dhats
out <- list(delta=deltaorig, dhat=delta, confdist=dout, iord=dhat2$iord.prim, conf = xnew, stress=sqrt(snew), spp=spp, ndim=p, weightmat=weightmato, resmat=resmat, rss=rss, init=xstart, model="local POST-MDS", niter = itel, nobj = dim(xnew)[1], type = type, call=match.call(), stress.m=snew, alpha = anew, sigma = snew, tdelta=deltaold, parameters=c(kappa=kappa,lambda=lambda,nu=nu,rho=rho,tau=tau), pars=c(kappa=kappa,lambda=lambda,nu=nu,rho=rho,tau=tau), theta=c(kappa=kappa,lambda=lambda,nu=nu,rho=rho,tau=tau), tweightmat=weightmat, trace = tracev)
class(out) <- c("smacofP","smacofB","smacof")
out
}
#' Local MDS (Venna and Kaski Style) either as self-organizing or not
#'
#' Fits Local MDS (Venna and Kaski, 2006) by quasi-majorization. It works like this: Find an initial weight matrix where w_ij(X^0)=0 if d_ij(X^0)> tau and delta_ij > tau, w_ij(X^0)= rho if d_ij(X^0)> tau and delta_ij <= tau, w_ij(X^0)= 1-rho if d_ij(X^0)<= tau and delta_ij > tau, w_ij(X^0)=1 if d_ij(X^0)<= tau and delta_ij <= tau.
#'
#' If rho = 0, this is CLCA/SMDS see \code{\link[smacofx]{spmds}}.
#'
#' This should not be confused with Local MDS of Chen and Buja (2009) as implemented in \code{\link[smacofx]{lmds}}.
#'
#' The neighborhood parameter tau is kept fixed in 'localmd'. The function 'so_localmds' implement a self-organising principle, where the model is repeatedly fitted for a decreasing sequence of taus.
#'
#' @param delta dist object or a symmetric, numeric data.frame or matrix of distances
#' @param rho is the continuity parameter and must be a value in [0,1]. This trades-off trustworthiness with continuity of the mapping. rho = 0 will give a weighting as in curvilinear component analysis. If rho=1 it will give an MDS with all inputs > tau treated as missing.
#' @param tau the boundary/neighbourhood parameter(s) (called sigma in the original paper). For 'lpostmds' and 'lmdsvk' it is supposed to be a numeric scalar (if a sequence is supplied the maximum is taken as tau) and all the transformed fitted distances exceeding tau are set to 0 via the weightmat (assignment can change between iterations). It defaults to the 25\% quantile of the explictly normalized (power transformed) delta. We recommend not to use the defaults!
#' @param weightmat a matrix of finite weights.
#' @param init starting configuration. If NULL (default) we fit a full rstress model.
#' @param ndim dimension of the configuration; defaults to 2
#' @param acc numeric accuracy of the iteration. Default is 1e-8.
#' @param itmax maximum number of iterations. Default is 10000.
#' @param verbose should fitting information be printed; if > 0 then yes
#' @param principal If 'TRUE', principal axis transformation is applied to the final configuration
#' @param traceIt save the iteration progress in a vector (stress values)
#'
#' @return a 'smacofP' object (inheriting from 'smacofB', see \code{\link[smacof]{smacofSym}}). It is a list with the components
#' \itemize{
#' \item delta: Observed, untransformed dissimilarities
#' \item tdelta: Observed explicitly transformed dissimilarities, normalized
#' \item dhat: Explicitly transformed dissimilarities (dhats), optimally scaled and normalized
#' \item confdist: Transformed configuration distances
#' \item conf: Matrix of fitted configuration
#' \item stress: Default stress (stress1; sqrt of explicitly normalized stress)
#' \item spp: Stress per point
#' \item ndim: Number of dimensions
#' \item model: Name of smacof model
#' \item niter: Number of iterations
#' \item nobj: Number of objects
#' \item type: Type of MDS model
#' \item weightmat: weighting matrix as supplied
#' \item stress.m: Default stress (square of stress1, explicitly normalized stress)
#' \item tweightmat: transformed weighting matrix; it is weightmat^nu and containing all the 0s for the distances/dissimilartities set to 0.
#' \item trace: if 'traceIt=TRUE' a vector with the iteration progress
#'}
#'
#'
#' @details
#' The solution is found by "quasi-majorization", which means that the majorization is only real majorization once the weightmat no longer changes. This typically happens after a few iterations. Due to that it can be that in the beginning the stress may not decrease monotonically and that there's a chance it might never.
#'
#' If tau is too small it may happen that all distances for one i to all j are zero and then there will be an error, so make sure to set a larger tau.
#'
#' In the standard functions 'localmds' we keep tau fixed throughout. This means that if tau is large enough, then the result is the same as the corresponding MDS.
#'
#' @importFrom stats dist as.dist quantile
#' @importFrom smacof transform transPrep
#'
#' @examples
#' dis<-smacof::morse
#' res2<-localmds(dis, rho = 0.5, tau = 0.3, itmax = 500, traceIt = TRUE) #use higher itmax#'
#' res2
#' summary(res2)
#' plot(res2)
#'
#' ##which d_{ij}(X)^kappa exceeded tau at convergence (i.e., have been set to 0)?
#' res2$tweightmat
#'
#' # We use Quasi-Majorization
#' res2$trace
#'
#'
#'
#' @export
localmds <- function(delta, rho = 0, tau = NULL, weightmat=1-diag(nrow(delta)), init=NULL, ndim = 2, acc= 1e-8, itmax = 10000, verbose = FALSE, principal=FALSE, traceIt=FALSE) {
cc <- match.call()
type <- "ratio"
#type <- match.arg(type, c("ratio", "interval","mspline","spline"), several.ok = FALSE)
if(inherits(delta,"dist") || is.data.frame(delta)) delta <- as.matrix(delta)
if(!isSymmetric(delta)) stop("delta is not symmetric.\n")
weightmato <- weightmat
weightmat[!is.finite(weightmat)] <- 0
if(is.null(tau)) tau <- stats::quantile((delta/enorm(delta,weightmat)),0.25)
out <- lpostmds(delta=delta, lambda=1, kappa=1, nu=1, rho=rho, tau=tau, type=type, weightmat=weightmato, init=init, ndim=ndim, acc=acc, itmax=itmax, verbose=verbose, principal=principal,traceIt=traceIt)
out$model <- "Local MDS (VK)"
out$call <- cc
out$parameters <- out$theta <- out$pars <- c(rho=rho,tau=tau)
out
}
#' @rdname localmds
#' @export
lmdsvk <- localmds
#' @rdname lpostmds
#' @export
localpostmds <- lpostmds
### TODO THIS
## #' @rdname spmds
## #' @param epochs for 'so_spmds' and tau being scalar, it gives the number of passes through the data. The sequence of taus created is 'seq(tau,tau/epochs,length.out=epochs)'. If tau is of length >1, this argument is ignored. For 'so_lpostmds' tau is supposed to be either a user supplied decreasing sequence of taus or if a scalar the maximum tau from which a decreasing sequence of taus is generated automatically as 'seq(from=tau,to=tau/epochs,length.out=epochs)' and then used in sequence.
## #' In the standard functions 'lpostmds' and 'lmdsvk' we keep tau fixed throughout. This means that if tau is large enough, then the result is the same as the corresponding MDS. In the orginal publication the idea was that of a self-organizing map which decreased tau over epochs (i.e., passes through the data). This can be achieved with our function 'so_lpostmds' 'so_lmdsvk' which creates a vector of decreasing tau values, calls the function 'lpostmds' with the first tau, then supplies the optimal configuration obtained as the init for the next call with the next tau and so on.
## #' \donttest{
## #' ## Self-organizing map style (as in the clca publication)
## #' #run the som-style (p)smds
## #' sommod1<-so_lpostmds(dis,tau=1,kappa=0.5,lambda=2,epochs=10,verbose=1)
## #' sommod1
## #' }
## #' @export
## so_spmds <- function(delta, kappa=1, lambda=1, nu=1, tau=max(delta), epochs=10, type="ratio", ties="primary", weightmat=1-diag(nrow(delta)), init=NULL, ndim = 2, acc= 1e-8, itmax = 10000, verbose = FALSE, principal=FALSE, spline.degree = 2, spline.intKnots = 2) {
## cc <- match.call()
## if(inherits(delta,"dist") || is.data.frame(delta)) delta <- as.matrix(delta)
## if(!isSymmetric(delta)) stop("delta is not symmetric.\n")
## if(length(tau)<2)
## {
## taumax <- tau
## taumin <- tau/epochs
## taus <- seq(taumax,taumin,length.out=epochs)
## } else taus <- tau
## if(any(diff(taus)>0)) taus <- sort(taus,decreasing=TRUE)
## finconf <- init
## for(i in 1:length(taus))
## {
## if(verbose>0) cat(paste0("Epoch ",i,": tau=",taus[i],"\n"))
## tmp<-spmds(delta=delta, lambda=lambda, kappa=kappa, nu=nu, tau=taus[i], type=type, ties=ties, weightmat=weightmat, init=finconf, ndim=ndim, verbose=verbose-1, acc=acc, itmax=itmax, principal=principal, spline.degree=spline.degree, spline.intKnots=spline.intKnots)
## finconf<-tmp$conf
## finmod<-tmp
## }
## finmod$call <- cc
## finmod$model <- "SO-SPMDS"
## return(finmod)
## }
## #' @rdname spmds
## #' \donttest{
## #' ## Self-organizing map style (as in the clca publication)
## #' #run the som-style (p)smds
## #' sommod2<-so_localmds(dis,tau=1,epochs=10,verbose=1)
## #' sommod2
## #' }
## #' In the orginal publication the idea was that of a self-organizing map which decreased tau over epochs (i.e., passes through the data). This can be achieved with our function 'so_localmds' which creates a vector of decreasing tau values, calls the function 'localmds' with the first tau, then supplies the optimal configuration obtained as the init for the next call with the next tau and so on.
## #'For 'so_lpostmds' tau is supposed to be either a user supplied decreasing sequence of taus or if a scalar the maximum tau from which a decreasing sequence of taus is generated automatically as 'seq(from=tau,to=tau/epochs,length.out=epochs)' and then used in sequence
## #' @export
## #' @param epochs for 'so_spmds' and tau being scalar, it gives the number of passes through the data. The sequence of taus created is 'seq(tau,tau/epochs,length.out=epochs)'. If tau is of length >1, this argument is ignored.
## so_smds <- function(delta, tau=max(delta), epochs=10, type="ratio", ties="primary", weightmat=1-diag(nrow(delta)), init=NULL, ndim = 2, acc= 1e-8, itmax = 10000, verbose = FALSE, principal=FALSE, spline.degree = 2, spline.intKnots = 2) {
## cc <- match.call()
## if(inherits(delta,"dist") || is.data.frame(delta)) delta <- as.matrix(delta)
## if(!isSymmetric(delta)) stop("delta is not symmetric.\n")
## if(length(tau)<2)
## {
## taumax <- tau
## taumin <- tau/epochs
## taus <- seq(taumax,taumin,length.out=epochs)
## } else taus <- tau
## if(any(diff(taus)>0)) taus <- sort(taus,decreasing=TRUE)
## finconf <- init
## for(i in 1:length(taus))
## {
## if(verbose>0) cat(paste0("Epoch ",i,": tau=",taus[i],"\n"))
## tmp<-smds(delta=delta, tau=taus[i], type=type, ties=ties, weightmat=weightmat, init=finconf, ndim=ndim, verbose=verbose, acc=acc, itmax=itmax, principal=principal, spline.degree=spline.degree, spline.intKnots=spline.intKnots)
## finconf<-tmp$conf
## finmod<-tmp
## }
## finmod$call <- cc
## finmod$model <- "SO-SMDS"
## return(finmod)
## }
Try the cops package in your browser
Any scripts or data that you put into this service are public.
cops documentation built on Oct. 22, 2026, 3 a.m.
R Package Documentation
Browse R Packages
We want your feedback!
Note that we can't provide technical support on individual packages. You should contact the package authors for that.
Defines functions localmds
#' Local Power Stress MDS either as self-organizing or not
#'
#' A local version of power stress MDS (POST-)MDS by quasi-majorization with ratio, interval and spline optimal scaling for dissimilarities and optional power transformations. This marries POST-MDS with the local multidimensional scaling idea of Venna and Karski (2006). It works like this: Find an initial weight matrix where w_ij(X^0)=0 if d_ij(X^0)> tau and delta_ij > tau, w_ij(X^0)= rho if d_ij(X^0)> tau and delta_ij <= tau, w_ij(X^0)= 1-rho if d_ij(X^0)<= tau and delta_ij > tau, w_ij(X^0)=1 if d_ij(X^0)<= tau and delta_ij <= tau.
#'
#' Then fit a POST-MDS with these weights. Then in each successive iteration step, the weightmat is recalculated so that w_ij(X^(n+1))=0 if d_ij(X^(n+1))> tau and delta_ij > tau, w_ij(X^(n+1))= rho if d_ij(X^(n+1))> tau and delta_ij <= tau, w_ij(X^(n+1))= 1-rho if d_ij(X^(n+1))<= tau and delta_ij > tau, w_ij(X^(n+1))=1 if d_ij(X^(n+1))<= tau and delta_ij <= tau.
#'
#' If rho = 0, this is eCLPCA/SPMDS see \code{\link[smacofx]{spmds}}.
#'
#' There is a wrappers 'lmdsvk' where the exponents are 1 and 'type=ratio'. The neighborhood parameter tau is kept fixed in 'lpostmds' and 'lmdsvk'. The functions 'so_lpostmds', 'so_lmdsvk' implement a self-organising principle, where the model is repeatedly fitted for a decreasing sequence of taus.
#'
#' @param delta dist object or a symmetric, numeric data.frame or matrix of distances
#' @param lambda exponent of the power transformation of the dissimilarities; defaults to 1., which is also the setup of 'lmdsvk'
#' @param kappa exponent of the power transformation of the fitted distances; defaults to 1, which is also the setup of 'lmdsvk'.
#' @param nu exponent of the power of the weighting matrix; defaults to 1 which is also the setup for 'lmdsvk'.
#' @param rho is the continuity parameter and must be a value in [0,1]. This trades-off trustworthiness with continuity of the mapping. rho = 0 will give a weighting as in curvilinear component analysis. If rho=1 it will give an MDS with all inputs > tau treated as missing.
#' @param tau the boundary/neighbourhood parameter(s) (called sigma in the original paper). For 'lpostmds' and 'lmdsvk' it is supposed to be a numeric scalar (if a sequence is supplied the maximum is taken as tau) and all the transformed fitted distances exceeding tau are set to 0 via the weightmat (assignment can change between iterations). It defaults to the 25\% quantile of the explictly normalized (power transformed) delta. We recommend not to use the defaults!
#' @param type what type of MDS to fit. Currently one of "ratio", "interval" or "mspline". Default is "ratio".
#' @param spline.degree Degree of the spline for ‘mspline’ MDS type
#' @param spline.intKnots Number of interior knots of the spline for ‘mspline’ MDS type
#' @param weightmat a matrix of finite weights.
#' @param init starting configuration. If NULL (default) we fit a full rstress model.
#' @param ndim dimension of the configuration; defaults to 2
#' @param acc numeric accuracy of the iteration. Default is 1e-8.
#' @param itmax maximum number of iterations. Default is 10000.
#' @param verbose should fitting information be printed; if > 0 then yes
#' @param principal If 'TRUE', principal axis transformation is applied to the final configuration
#' @param traceIt save the iteration progress in a vector (stress values)
#'
#' @return a 'smacofP' object (inheriting from 'smacofB', see \code{\link[smacof]{smacofSym}}). It is a list with the components
#' \itemize{
#' \item delta: Observed, untransformed dissimilarities
#' \item tdelta: Observed explicitly transformed dissimilarities, normalized
#' \item dhat: Explicitly transformed dissimilarities (dhats), optimally scaled and normalized
#' \item confdist: Transformed configuration distances
#' \item conf: Matrix of fitted configuration
#' \item stress: Default stress (stress1; sqrt of explicitly normalized stress)
#' \item spp: Stress per point
#' \item ndim: Number of dimensions
#' \item model: Name of smacof model
#' \item niter: Number of iterations
#' \item nobj: Number of objects
#' \item type: Type of MDS model
#' \item weightmat: weighting matrix as supplied
#' \item stress.m: Default stress (square of stress1, explicitly normalized stress)
#' \item tweightmat: transformed weighting matrix; it is weightmat^nu and containing all the 0s for the distances/dissimilartities set to 0.
#' \item trace: if 'traceIt=TRUE' a vector with the iteration progress
#'}
#'
#'
#' @details
#' The solution is found by "quasi-majorization", which means that the majorization is only real majorization once the weightmat no longer changes. This typically happens after a few iterations. Due to that it can be that in the beginning the stress may not decrease monotonically and that there's a chance it might never.
#'
#' If tau is too small it may happen that all distances for one i to all j are zero and then there will be an error, so make sure to set a larger tau.
#'
#'
#'
#' @importFrom stats dist as.dist quantile
#' @importFrom smacof transform transPrep
#'
#' @examples
#' dis<-smacof::morse
#' res<-lpostmds(dis, type = "interval", rho = 0.5, kappa = 2, lambda = 2,tau = 0.3, itmax = 100) #use higher itmax
#' res
#' summary(res)
#' plot(res)
#'
#' ##which d_ij(X)^kappa exceeded tau at convergence (i.e., have been set to 0)?
#' res$tweightmat
#'
#' # We use Quasi-Majorization
#' res$trace
#'
#'
#'
#' @export
lpostmds <- function (delta, lambda=1, kappa=1, nu=1, tau = NULL, rho = 0, type="ratio", weightmat=1-diag(nrow(delta)), init=NULL, ndim = 2, acc= 1e-8, itmax = 10000, verbose = FALSE, principal=FALSE, spline.degree = 2, spline.intKnots = 2, traceIt=FALSE) {
## TODO: Should we allow for different tau for input and output?
if(inherits(delta,"dist") || is.data.frame(delta)) delta <- as.matrix(delta)
if(!isSymmetric(delta)) stop("delta is not symmetric.\n")
if(inherits(weightmat,"dist") || is.data.frame(weightmat)) weightmat <- as.matrix(weightmat)
if(!isSymmetric(weightmat)) stop("weightmat is not symmetric.\n")
r <- kappa/2
if(rho < 0 && rho > 1) stop("rho must be in the interval [0, 1].")
## -- Setup for MDS type
if(missing(type)) type <- "ratio"
type <- match.arg(type, c("ratio", "interval", "mspline","spline"),several.ok = FALSE) #no ordinal allowed
trans <- type
typo <- type
if (trans=="ratio"){
trans <- "none"
}
# else if (trans=="ordinal" & ties=="primary"){
# trans <- "ordinalp"
# typo <- "ordinal (primary)"
# } else if(trans=="ordinal" & ties=="secondary"){
# trans <- "ordinals"
# typo <- "ordinal (secondary)"
# } else if(trans=="ordinal" & ties=="tertiary"){
# trans <- "ordinalt"
# typo <- "ordinal (tertiary)"
# } else
if(trans=="spline"){
trans <- "mspline"
type <- "mspline"
typo <- "mspline"
}
if(type %in% c("mspline", "spline")) lambda <- 1 #We dont allow powers for dissimilarities in nonmetric and spline MDS
n <- nrow (delta)
normi <- 0.5
##normi <- n #if normi=n we can use the iord structure in plot.smacofP
## but the problem is we don't get the correct stress then anymore.
p <- ndim
if (p > (n - 1)) stop("Maximum number of dimensions is n-1!")
if(is.null(rownames(delta))) rownames(delta) <- 1:n
labos <- rownames(delta) #labels
deltaorig <- delta
delta <- delta^lambda ## Transfrom the deltas
weightmato <- weightmat #3 so we can retunr this as supplied
weightmat <- weightmat^nu ##transform the weightmat
weightmat[!is.finite(weightmat)] <- 0 ## all infinite or NA weights get a weight of 0
delta <- delta / enorm (delta, weightmat) ## explicit normalization of input
if(is.null(tau)) tau <- stats::quantile(delta,0.25) ## tau default
if(length(tau)>1)
{
warning("Supplied tau is of length > 1. The min(tau) was used as tau.")
tau <- min(tau)
}
if(tau<=0) stop("tau must be positive.")
if(verbose>0) cat(paste("Fitting",type,"Local POST-MDS with lambda=",lambda, "kappa=",kappa,"nu=",nu, "rho=",rho,"and tau=",tau,"\n"))
disobj <- smacof::transPrep(as.dist(delta), 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)
deltaold <- delta
itel <- 1
##Starting Configs
xold <- init
# if(is.null(init)) xold <- smacof::torgerson (delta,r=kappa/2,type=type,ties=ties,weightmat=weightmat,ndim=ndim,init=init,itmax=itmax,principal=principal)$conf
if(is.null(init)) xold <- smacof::torgerson(delta, p = p)
xstart <- xold
xold <- xold / enorm (xold) #explicit normalization of output
nn <- diag (n)
dold <- sqdist (xold) #squared distances
doldpow <- mkPower(dold, kappa / 2)# distances^kappa
## --- Make local weightmat
trustwn <- (1-rho)*(doldpow <= tau) ##CCA penalty
contin <- rho*(deltaold <= tau)
tcmat <- trustwn+contin ##
lweightmat <- weightmat*tcmat
diag(lweightmat) <- 0 ## sanitizing
lweightmat[!is.finite(lweightmat)] <- 0 ## all infinite or NA weights get a weight of 0
##first optimal scaling
eold <- as.dist(mkPower(dold,r))
dhat <- smacof::transform(eold, disobj, w = as.dist(lweightmat), normq = normi)
dhatt <- dhat$res
dhatd <- structure(dhatt, Size = n, call = quote(as.dist.default(m=b)), class = "dist", Diag = FALSE, Upper = FALSE)
delta <- as.matrix(dhatd)
## first stress
rold <- sum (lweightmat * delta * mkPower (dold, r))
nold <- sum (lweightmat * mkPower (dold, 2 * r))
aold <- rold / nold
sold <- 1 - 2 * aold * rold + (aold ^ 2) * nold
tracev <- NULL
if(isTRUE(traceIt)) tracev <- rep(NA,itmax)
## Optimizing
## Loop Logic for iteration i: X^(i)|W^(i-1) -> W^(i)|X^(i) -> stress^(i) -> X^(i+1)|W^(i) -> etc.
repeat {
if(tau <= min(doldpow[lower.tri(doldpow)]) && min(delta[lower.tri(delta)])) stop("Current tau is lower than the smallest fitted distance and the smallest delta (therefore all distances are set to 0). Increase tau.")
p1 <- mkPower (dold, r - 1)
p2 <- mkPower (dold, (2 * r) - 1)
by <- mkBmat (lweightmat * delta * p1)
cy <- mkBmat (lweightmat * p2)
ga <- 2 * sum (lweightmat * p2)
be <- (2 * r - 1) * (2 ^ r) * sum (lweightmat * delta)
de <- (4 * r - 1) * (4 ^ r) * sum (lweightmat)
if (r >= 0.5) {
my <- by - aold * (cy - de * nn)
}
if (r < 0.5) {
my <- (by - be * nn) - aold * (cy - ga * nn)
}
xnew <- my %*% xold
xnew <- xnew / enorm (xnew)
dnew <- sqdist (xnew)
dnewpow <- mkPower(dnew,kappa/2)
### We now reset the weight matrix, so it is possible that a w_{ij} can change from 0 to >0 again
weightmat <- weightmato^nu #the supplied weightmat and carrying the nu with us
weightmat[!is.finite(weightmat)] <- 0
## --- Make local weightmat
trustwn <- (1-rho)*(dnewpow <= tau) ##trustworthiness
contin <- rho*(delta <= tau) #continuity; delta are now the dhats
tcmat <- trustwn+contin ##convex combo
lweightmat <- weightmat*tcmat #multiply with given weights
diag(lweightmat) <- 0 ## sanitizing
lweightmat[!is.finite(lweightmat)] <- 0 ## all infinite or NA weights get a weight of 0
##optimal scaling
e <- as.dist(mkPower(dnew,r)) #I need the dist(x) here for interval
dhat2 <- smacof::transform(e, disobj, w = as.dist(lweightmat), normq = normi) ## dhat update
dhatt <- dhat2$res
dhatd <- structure(dhatt, Size = n, call = quote(as.dist.default(m=b)), class = "dist", Diag = FALSE, Upper = FALSE)
delta <- as.matrix(dhatd)
rnew <- sum (lweightmat * delta * mkPower (dnew, r))
nnew <- sum (lweightmat * mkPower (dnew, 2 * r))
anew <- rnew / nnew
snew <- 1 - 2 * anew * rnew + (anew ^ 2) * nnew
if(is.na(snew)) #if there are issues with the values
{
snew <- sold
dnew <- dold
anew <- aold
xnew <- xold
}
if (verbose>2) {
cat (
formatC (itel, width = 4, format = "d"),
formatC (
sold,
digits = 10,
width = 13,
format = "f"
),
formatC (
snew,
digits = 10,
width = 13,
format = "f"
),
"\n"
)
}
if(isTRUE(traceIt)) tracev[itel] <- snew
if ((itel == itmax) || (abs(sold - snew) < acc)) #new
{
if(sold < snew) itel <- itel-1
snew <- min(sold,snew) #to make sure we use the lowest stress as it is quasi majorization and sold might be < snew at this point
break ()
}
itel <- itel + 1
xold <- xnew
dold <- dnew
sold <- snew
aold <- anew
}
xnew <- xnew/enorm(xnew)
## relabeling as they were removed in the optimal scaling
rownames(delta) <- labos
attr(xnew,"dimnames")[[1]] <- rownames(delta)
attr(xnew,"dimnames")[[2]] <- paste("D",1:p,sep="")
doutm <- mkPower(sqdist(xnew),r)
deltam <- delta
#delta <- structure(delta, Size = n, call = quote(as.dist.default(m=b)),
# class = "dist", Diag = FALSE, Upper = FALSE)
delta <- stats::as.dist(delta)
deltaorig <- stats::as.dist(deltaorig)
deltaold <- stats::as.dist(deltaold)
#doute <- doutm/enorm(doutm) #this is an issue here!
#doute <- stats::as.dist(doute)
dout <- stats::as.dist(doutm)
weightmatm <-lweightmat ##the local weightmat used
weightmat <- stats::as.dist(weightmatm)
spoint <- spp(delta, dout, weightmat)
resmat<-spoint$resmat
rss <- sum(spoint$resmat[lower.tri(spoint$resmat)])
spp <- spoint$spp
if (verbose > 0 && itel == itmax) warning("Iteration limit reached! You may want to increase the itmax argument!")
if (principal) {
xnew_svd <- svd(xnew)
xnew <- xnew %*% xnew_svd$v
}
if(isTRUE(traceIt)) tracev <- tracev[!is.na(tracev)]
if(verbose>1) cat("*** Stress:",snew, "; Stress-1 (default reported):",sqrt(snew),"\n")
#delta is input delta, tdelta is input delta with explicit transformation and normalized, dhat is dhats
out <- list(delta=deltaorig, dhat=delta, confdist=dout, iord=dhat2$iord.prim, conf = xnew, stress=sqrt(snew), spp=spp, ndim=p, weightmat=weightmato, resmat=resmat, rss=rss, init=xstart, model="local POST-MDS", niter = itel, nobj = dim(xnew)[1], type = type, call=match.call(), stress.m=snew, alpha = anew, sigma = snew, tdelta=deltaold, parameters=c(kappa=kappa,lambda=lambda,nu=nu,rho=rho,tau=tau), pars=c(kappa=kappa,lambda=lambda,nu=nu,rho=rho,tau=tau), theta=c(kappa=kappa,lambda=lambda,nu=nu,rho=rho,tau=tau), tweightmat=weightmat, trace = tracev)
class(out) <- c("smacofP","smacofB","smacof")
out
}
#' Local MDS (Venna and Kaski Style) either as self-organizing or not
#'
#' Fits Local MDS (Venna and Kaski, 2006) by quasi-majorization. It works like this: Find an initial weight matrix where w_ij(X^0)=0 if d_ij(X^0)> tau and delta_ij > tau, w_ij(X^0)= rho if d_ij(X^0)> tau and delta_ij <= tau, w_ij(X^0)= 1-rho if d_ij(X^0)<= tau and delta_ij > tau, w_ij(X^0)=1 if d_ij(X^0)<= tau and delta_ij <= tau.
#'
#' If rho = 0, this is CLCA/SMDS see \code{\link[smacofx]{spmds}}.
#'
#' This should not be confused with Local MDS of Chen and Buja (2009) as implemented in \code{\link[smacofx]{lmds}}.
#'
#' The neighborhood parameter tau is kept fixed in 'localmd'. The function 'so_localmds' implement a self-organising principle, where the model is repeatedly fitted for a decreasing sequence of taus.
#'
#' @param delta dist object or a symmetric, numeric data.frame or matrix of distances
#' @param rho is the continuity parameter and must be a value in [0,1]. This trades-off trustworthiness with continuity of the mapping. rho = 0 will give a weighting as in curvilinear component analysis. If rho=1 it will give an MDS with all inputs > tau treated as missing.
#' @param tau the boundary/neighbourhood parameter(s) (called sigma in the original paper). For 'lpostmds' and 'lmdsvk' it is supposed to be a numeric scalar (if a sequence is supplied the maximum is taken as tau) and all the transformed fitted distances exceeding tau are set to 0 via the weightmat (assignment can change between iterations). It defaults to the 25\% quantile of the explictly normalized (power transformed) delta. We recommend not to use the defaults!
#' @param weightmat a matrix of finite weights.
#' @param init starting configuration. If NULL (default) we fit a full rstress model.
#' @param ndim dimension of the configuration; defaults to 2
#' @param acc numeric accuracy of the iteration. Default is 1e-8.
#' @param itmax maximum number of iterations. Default is 10000.
#' @param verbose should fitting information be printed; if > 0 then yes
#' @param principal If 'TRUE', principal axis transformation is applied to the final configuration
#' @param traceIt save the iteration progress in a vector (stress values)
#'
#' @return a 'smacofP' object (inheriting from 'smacofB', see \code{\link[smacof]{smacofSym}}). It is a list with the components
#' \itemize{
#' \item delta: Observed, untransformed dissimilarities
#' \item tdelta: Observed explicitly transformed dissimilarities, normalized
#' \item dhat: Explicitly transformed dissimilarities (dhats), optimally scaled and normalized
#' \item confdist: Transformed configuration distances
#' \item conf: Matrix of fitted configuration
#' \item stress: Default stress (stress1; sqrt of explicitly normalized stress)
#' \item spp: Stress per point
#' \item ndim: Number of dimensions
#' \item model: Name of smacof model
#' \item niter: Number of iterations
#' \item nobj: Number of objects
#' \item type: Type of MDS model
#' \item weightmat: weighting matrix as supplied
#' \item stress.m: Default stress (square of stress1, explicitly normalized stress)
#' \item tweightmat: transformed weighting matrix; it is weightmat^nu and containing all the 0s for the distances/dissimilartities set to 0.
#' \item trace: if 'traceIt=TRUE' a vector with the iteration progress
#'}
#'
#'
#' @details
#' The solution is found by "quasi-majorization", which means that the majorization is only real majorization once the weightmat no longer changes. This typically happens after a few iterations. Due to that it can be that in the beginning the stress may not decrease monotonically and that there's a chance it might never.
#'
#' If tau is too small it may happen that all distances for one i to all j are zero and then there will be an error, so make sure to set a larger tau.
#'
#' In the standard functions 'localmds' we keep tau fixed throughout. This means that if tau is large enough, then the result is the same as the corresponding MDS.
#'
#' @importFrom stats dist as.dist quantile
#' @importFrom smacof transform transPrep
#'
#' @examples
#' dis<-smacof::morse
#' res2<-localmds(dis, rho = 0.5, tau = 0.3, itmax = 500, traceIt = TRUE) #use higher itmax#'
#' res2
#' summary(res2)
#' plot(res2)
#'
#' ##which d_{ij}(X)^kappa exceeded tau at convergence (i.e., have been set to 0)?
#' res2$tweightmat
#'
#' # We use Quasi-Majorization
#' res2$trace
#'
#'
#'
#' @export
localmds <- function(delta, rho = 0, tau = NULL, weightmat=1-diag(nrow(delta)), init=NULL, ndim = 2, acc= 1e-8, itmax = 10000, verbose = FALSE, principal=FALSE, traceIt=FALSE) {
cc <- match.call()
type <- "ratio"
#type <- match.arg(type, c("ratio", "interval","mspline","spline"), several.ok = FALSE)
if(inherits(delta,"dist") || is.data.frame(delta)) delta <- as.matrix(delta)
if(!isSymmetric(delta)) stop("delta is not symmetric.\n")
weightmato <- weightmat
weightmat[!is.finite(weightmat)] <- 0
if(is.null(tau)) tau <- stats::quantile((delta/enorm(delta,weightmat)),0.25)
out <- lpostmds(delta=delta, lambda=1, kappa=1, nu=1, rho=rho, tau=tau, type=type, weightmat=weightmato, init=init, ndim=ndim, acc=acc, itmax=itmax, verbose=verbose, principal=principal,traceIt=traceIt)
out$model <- "Local MDS (VK)"
out$call <- cc
out$parameters <- out$theta <- out$pars <- c(rho=rho,tau=tau)
out
}
#' @rdname localmds
#' @export
lmdsvk <- localmds
#' @rdname lpostmds
#' @export
localpostmds <- lpostmds
### TODO THIS
## #' @rdname spmds
## #' @param epochs for 'so_spmds' and tau being scalar, it gives the number of passes through the data. The sequence of taus created is 'seq(tau,tau/epochs,length.out=epochs)'. If tau is of length >1, this argument is ignored. For 'so_lpostmds' tau is supposed to be either a user supplied decreasing sequence of taus or if a scalar the maximum tau from which a decreasing sequence of taus is generated automatically as 'seq(from=tau,to=tau/epochs,length.out=epochs)' and then used in sequence.
## #' In the standard functions 'lpostmds' and 'lmdsvk' we keep tau fixed throughout. This means that if tau is large enough, then the result is the same as the corresponding MDS. In the orginal publication the idea was that of a self-organizing map which decreased tau over epochs (i.e., passes through the data). This can be achieved with our function 'so_lpostmds' 'so_lmdsvk' which creates a vector of decreasing tau values, calls the function 'lpostmds' with the first tau, then supplies the optimal configuration obtained as the init for the next call with the next tau and so on.
## #' \donttest{
## #' ## Self-organizing map style (as in the clca publication)
## #' #run the som-style (p)smds
## #' sommod1<-so_lpostmds(dis,tau=1,kappa=0.5,lambda=2,epochs=10,verbose=1)
## #' sommod1
## #' }
## #' @export
## so_spmds <- function(delta, kappa=1, lambda=1, nu=1, tau=max(delta), epochs=10, type="ratio", ties="primary", weightmat=1-diag(nrow(delta)), init=NULL, ndim = 2, acc= 1e-8, itmax = 10000, verbose = FALSE, principal=FALSE, spline.degree = 2, spline.intKnots = 2) {
## cc <- match.call()
## if(inherits(delta,"dist") || is.data.frame(delta)) delta <- as.matrix(delta)
## if(!isSymmetric(delta)) stop("delta is not symmetric.\n")
## if(length(tau)<2)
## {
## taumax <- tau
## taumin <- tau/epochs
## taus <- seq(taumax,taumin,length.out=epochs)
## } else taus <- tau
## if(any(diff(taus)>0)) taus <- sort(taus,decreasing=TRUE)
## finconf <- init
## for(i in 1:length(taus))
## {
## if(verbose>0) cat(paste0("Epoch ",i,": tau=",taus[i],"\n"))
## tmp<-spmds(delta=delta, lambda=lambda, kappa=kappa, nu=nu, tau=taus[i], type=type, ties=ties, weightmat=weightmat, init=finconf, ndim=ndim, verbose=verbose-1, acc=acc, itmax=itmax, principal=principal, spline.degree=spline.degree, spline.intKnots=spline.intKnots)
## finconf<-tmp$conf
## finmod<-tmp
## }
## finmod$call <- cc
## finmod$model <- "SO-SPMDS"
## return(finmod)
## }
## #' @rdname spmds
## #' \donttest{
## #' ## Self-organizing map style (as in the clca publication)
## #' #run the som-style (p)smds
## #' sommod2<-so_localmds(dis,tau=1,epochs=10,verbose=1)
## #' sommod2
## #' }
## #' In the orginal publication the idea was that of a self-organizing map which decreased tau over epochs (i.e., passes through the data). This can be achieved with our function 'so_localmds' which creates a vector of decreasing tau values, calls the function 'localmds' with the first tau, then supplies the optimal configuration obtained as the init for the next call with the next tau and so on.
## #'For 'so_lpostmds' tau is supposed to be either a user supplied decreasing sequence of taus or if a scalar the maximum tau from which a decreasing sequence of taus is generated automatically as 'seq(from=tau,to=tau/epochs,length.out=epochs)' and then used in sequence
## #' @export
## #' @param epochs for 'so_spmds' and tau being scalar, it gives the number of passes through the data. The sequence of taus created is 'seq(tau,tau/epochs,length.out=epochs)'. If tau is of length >1, this argument is ignored.
## so_smds <- function(delta, tau=max(delta), epochs=10, type="ratio", ties="primary", weightmat=1-diag(nrow(delta)), init=NULL, ndim = 2, acc= 1e-8, itmax = 10000, verbose = FALSE, principal=FALSE, spline.degree = 2, spline.intKnots = 2) {
## cc <- match.call()
## if(inherits(delta,"dist") || is.data.frame(delta)) delta <- as.matrix(delta)
## if(!isSymmetric(delta)) stop("delta is not symmetric.\n")
## if(length(tau)<2)
## {
## taumax <- tau
## taumin <- tau/epochs
## taus <- seq(taumax,taumin,length.out=epochs)
## } else taus <- tau
## if(any(diff(taus)>0)) taus <- sort(taus,decreasing=TRUE)
## finconf <- init
## for(i in 1:length(taus))
## {
## if(verbose>0) cat(paste0("Epoch ",i,": tau=",taus[i],"\n"))
## tmp<-smds(delta=delta, tau=taus[i], type=type, ties=ties, weightmat=weightmat, init=finconf, ndim=ndim, verbose=verbose, acc=acc, itmax=itmax, principal=principal, spline.degree=spline.degree, spline.intKnots=spline.intKnots)
## finconf<-tmp$conf
## finmod<-tmp
## }
## finmod$call <- cc
## finmod$model <- "SO-SMDS"
## return(finmod)
## }
Try the cops package in your browser
Any scripts or data that you put into this service are public.
R Package Documentation
Browse R Packages
We want your feedback!
Note that we can't provide technical support on individual packages. You should contact the package authors for that.
Embedding an R snippet on your website
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.