Nothing
smacofRect <- function(delta, ndim = 2, type = c("ratio", "interval", "ordinal", "mspline"),
conditionality = "unconditional", lambda = 0.5, omega = 1,
circle = c("none", "row", "column"), weightmat = NULL, init = NULL,
fixed = c("none", "row", "column"), fixed.coord = NULL,
ties = c("primary", "secondary"), verbose = FALSE, relax = TRUE, itmax = 10000,
eps = 1e-6, spline.degree = 2, spline.intKnots = 2, parallelize = FALSE)
# init ... either a list of 2 matrices of dimension n \times p and m \times p with
# starting values. if NULL, svd is used.
{
circle <- match.arg(circle, c("none", "row", "column"), several.ok = FALSE)
type <- match.arg(type, c("ratio", "interval", "ordinal", "mspline"), several.ok = FALSE)
ties <- match.arg(ties, c("primary","secondary"), several.ok = FALSE)
conditionality <- match.arg(conditionality, c("unconditional", "matrix", "row"), several.ok = FALSE)
if (conditionality == "unconditional") conditionality <- "matrix" ## for backward compatibility
fixed <- match.arg(fixed, c("none", "row", "column"), several.ok = FALSE)
diss <- delta
rnames <- rownames(delta)
if (is.data.frame(diss)) diss <- as.matrix(diss)
checkdiss(diss)
n <- dim(diss)[1] #number of individuals
m <- dim(diss)[2] #number of objects
p <- ndim
## --- sanity check external unfolding
if (is.null(fixed.coord) & fixed != "none") stop("fixed.coord is not specified.")
if (fixed == "row"){
if (nrow(fixed.coord) != n | ncol(fixed.coord) < p)
stop("Matrix of fixed row coordinates specified by fixed.coord has an incorrect size")
} else if (fixed == "col"){
if (nrow(fixed.coord) != m | ncol(fixed.coord) < p)
stop("Matrix of fixed column coordinates specified by fixed.coord has an incorrect size")
}
## --- weight init
if (is.null(weightmat)) {
w <- matrix(1,n,m) #initialize weights (as 1)
if (any(is.na(diss))) w[is.na(diss)] <- 0 # blank out missings
} else w <- weightmat
wpp <- sum(w)
## --- Prepare for parallization
if (parallelize) {
noCores <- detectCores()/2
cl <- parallel::makeCluster(noCores)
#cl <- parallel::makeCluster(4)
doParallel::registerDoParallel(cl)
}
delta <- ifelse(is.na(diss),0,diss) #replace NA's by 0
## --- 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=="spline"){
trans <- "mspline"
}
## --- end optimal scaling prep
## --- conditionality prep
disobj <- list()
if (conditionality == "matrix"){
disobj[[1]] <- transPrep(as.vector(delta), trans = trans, spline.intKnots = spline.intKnots, spline.degree = spline.degree)
if (trans == "mspline") disobj[[1]]$base <- cbind(rep(1, nrow(disobj[[1]]$base)), disobj[[1]]$base)
tt <- transform(as.vector(delta), disobj[[1]], w = as.vector(w)) #, normq = wpp
dhat <- matrix(tt$res, n, m) ## dhat update
tt <- vector(mode = "list", length = 1)
} else { ## conditionality == "row"
dhat <- matrix(0, n, m)
for (i in 1:n) {
disobj[[i]] <- transPrep(delta[i, ], trans = trans, spline.intKnots = spline.intKnots, spline.degree = spline.degree)
if (trans == "mspline") disobj[[i]]$base <- cbind(rep(1, nrow(disobj[[i]]$base)), disobj[[i]]$base)
dhat[i, ] <- transform(delta[i, ], disobj[[i]], w = w[i, ], normq = sum(w))$res ## dhat update
}
tt <- vector(mode = "list", length = n)
}
## --- end conditionality prep
itel <- 1
e <- dhat^2
e <- -0.5*(e - outer(rowSums(e)/m, colSums(e)/n, "+") + (sum(e)/(n*m)))
## --------- external unfolding ------------
if (is.list(init) & fixed == "none") { ## Initial coordinates given
x <- init[[1]][, 1:p]
y <- init[[2]][, 1:p]
} else if (is.list(init) & fixed == "row") { ## Initial coordinates given + fixed rows
fixed.coord <- fixed.coord[, 1:p]
x <- fixed.coord
y <- init[[2]][, 1:p]
} else if (is.list(init) & fixed == "column") { ## Initial coordinates given + fixed columns
fixed.coord <- fixed.coord[, 1:p]
x <- init[[1]][, 1:p]
y <- fixed.coord
} else if (!is.list(init) & fixed == "none") { ## No initial coordinates given
z <- svd(e, nu = p, nv = 0) #SVD for e (pos. distances)
x <- z$u #starting value for x
y <- crossprod(e, x) #starting value for y
} else if (!is.list(init) & fixed == "row") { ## No initial coordinates given + fixed rows
fixed.coord <- fixed.coord[, 1:p]
x <- fixed.coord #starting value for x
y <- crossprod(e, x) %*% solve(crossprod(x)) #starting value for y
} else if (!is.list(init) & fixed == "column") { ## No initial coordinates given + fixed cols
fixed.coord <- fixed.coord[, 1:p]
y <- fixed.coord #starting value for y
x <- solve(crossprod(y), crossprod(y, t(e))) #starting value for x
x <- t(x)
}
## Rescale x and y by a constant to fit e
dil <- sum(diag(t(x) %*% (e %*% y))) / sum(diag(crossprod(x) %*% crossprod(y)))
x <- (dil^.5/ sum(x^2)^.5 ) * x
y <- (dil^.5/ sum(y^2)^.5 ) * y
## --------- circular restriction
if (circle != "none"){
r <- projCircle(x,y,x,y,circle=circle)
x <- r$x
y <- r$y
wr <- rowSums(w)
wc <- colSums(w)
lambda.circ <- 2*max(c(wr,wc))
}
d <- distRect(x,y,0) #n times m of reproduced diss
ps <- pstress(dhat, d, w, omega, lambda, wpp, conditionality) #pstress value
lold <- ps$pstress
if (circle == "none") {
ww <- w
wr <- rowSums(ww)
wc <- colSums(ww)
v <- solve(diag(wc) + (1/m) - crossprod(ww, ww/wr)) - (1/m)
}
#------------------- begin majorization -----------------------------
repeat {
d.is.0 <- d < eps
if (circle == "none") {
b <- w * dhat * (!d.is.0) / (d + d.is.0) # B matrix
br <- rowSums(b) #rows B
bc <- colSums(b) #columns W
xraw <- (br * x) - ( b %*% y)
yraw <- (bc * y) - crossprod(b, x)
xold <- x
yold <- y
yunc <- v %*% (yraw + crossprod(ww, xraw/wr))
xunc <- (xraw + (ww %*% yunc))/wr
## x update
if (fixed != "column") {
y <- yunc
} else {
y <- (sum(diag(crossprod(fixed.coord, yunc))/sum(fixed.coord^2))) * fixed.coord
}
## y update
if (fixed != "row") {
x <- xunc
} else {
x <- (sum(diag(crossprod(fixed.coord, xunc))/sum(fixed.coord^2))) * fixed.coord
}
if (relax & itel > 100){
x <- 2 * x - xold
y <- 2 * y - yold
}
} else {
b <- w * (1 - (!d.is.0) * dhat / (d + d.is.0)) #B matrix
br <- rowSums(b) #rows B
bc <- colSums(b) #columns W
xunc <- x - outer(br, rep(1/lambda.circ, p), "*") * x + b %*% (y/lambda.circ)
yunc <- y - outer(bc, rep(1/lambda.circ, p), "*") * y + t(b) %*% (x/lambda.circ)
r <- projCircle(xunc, yunc, x, y, circle = circle)
if (fixed != "row") {
x <- r$x
} else {
x <- (sum(diag(crossprod(fixed.coord, xunc))/sum(fixed.coord^2))) * fixed.coord
}
if (fixed != "column") {
y <- r$y
} else {
y <- (sum(diag(crossprod(fixed.coord, yunc))/sum(fixed.coord^2))) * fixed.coord
}
}
d <- distRect(x, y, 0) #compute distances (update)
dhat.old <- dhat
if ( conditionality == "matrix" ) {
nrmw <- sqrt( ps$nrmw2 )
nrmm <- sqrt( ps$nrmm2 )
nrmv <- sqrt( ps$nrmv2 )
g1 <- ps$nrmd2 / wpp
g2 <- ps$nrmv2 + omega * ps$nrmm2
g3 <- sqrt( ps$nrmv2 )
g <- g1^( lambda / 2 ) * sqrt( g2 ) / g3
alpha2 <- 0.5 * lambda * sqrt( g2 ) * g1^( lambda / 2 - 1 )
alpha3 <- 0.5 * g1^( lambda / 2 ) / sqrt( g2 )
fmin <- min( dhat )
tau1 <- ps$wsum / wpp
tau2 <- 0.5 / eps
if ( fmin >= eps ) tau2 <- 0.5 / fmin
tau4 <- tau1 / nrmv
beta1 <- 1 / wpp
beta3 <- 1 + omega + 2 * omega * tau1 * tau2
beta5 <- tau2 * tau4
lower <- alpha2 * beta1 + alpha3 * beta3 + g * beta5
t1 <- tau2 * as.vector( dhat ) - 0.5
t2 <- as.vector( dhat ) / ( 2 * nrmv )
b1 <- as.vector( d ) / wpp
#b2 <- ifelse( as.vector( dhat ) <= eps, tau1, tau1 + omega * as.vector( dhat ) + ( 2 * omega * tau1 ) * t1 )
b2 <- tau1 + (as.vector(dhat) > eps) * omega * (dhat + ( 2 * tau1 ) * t1)
#b3 <- ifelse( as.vector( dhat ) <= eps, 0, t2 + tau4 * t1 )
b3 <- (as.vector(dhat) > eps) * (t2 + tau4 * t1)
upper <- alpha2 * b1 + alpha3 * b2 + g * b3
ksi <- upper / lower
##tt[[1]] <- transform( ksi, disobj[[1]], w = as.vector( w ) )
tt[[1]] <- smacof::transform( ksi, disobj[[1]], w = as.vector( w ) )
dhat <- matrix( tt[[1]]$res, n, m ) ## dhat update
} else { ## conditionality == "row"
if (parallelize){
ii <- c(0, floor(n * (1:noCores)/noCores))
s <- 2:(noCores + 1)
tt <- foreach(s = s, .combine = c, .verbose = FALSE,
.maxcombine = noCores, .multicombine = TRUE, export = ps) %dopar% {
dhatRowUpdateBlock(ps, s, ii, n, omega, lambda, wpp, eps, dhat, d, w, disobj)
}
for (i in 1:n) {
dhat[i, ] <- tt[[i]]$res
}
} else {
for ( i in 1:n ) {
mnc1 <- sum( ps$c1 )
mnc2 <- sum( ps$c2 )
nrmw <- sqrt( ps$nrmw2[i] )
nrmm <- sqrt( ps$nrmm2[i] )
nrmv <- sqrt( ps$nrmv2[i] )
sqrtn <- sqrt( n )
sumc1 <- mnc1 - ps$nrmd2[i]
sumc2 <- mnc2 - ( ps$nrmv2[i] + omega * ps$nrmm2[i] ) / ps$nrmv2[i]
g1 = ( sumc1 + ps$nrmd2[i] ) / wpp
g2 = ( 1 + sumc2 ) * ps$nrmv2[i] + omega * ps$nrmm2[i]
g3 = n * ps$nrmv2[i]
g = sqrt ( g1^lambda * g2 / g3 )
alpha2 <- 0.5 * lambda * sqrt( g2 ) * g1^( lambda/2 - 1 )
alpha3 <- 0.5 * g1^( lambda / 2 ) / sqrt( g2 )
fmin <- min( dhat[i, ] )
tau1 <- ps$wsum[i] / ps$sumw[i]
tau2 <- 0.5 / eps
if ( fmin >= eps ) tau2 <- 0.5 / fmin
tau4 <- tau1 / nrmv
beta1 <- 1 / wpp
beta3 <- 1 + sumc2 + omega + 2 * omega * tau1 * tau2
beta5 <- sqrtn * tau2 * tau4
lower <- alpha2 * beta1 + alpha3 * beta3 + g * beta5
t1 <- tau2 * dhat[i, ] - 0.5
t2 <- dhat[i, ] / ( 2 * nrmv )
b1 <- d[i, ] / wpp
#b2 <- ifelse( dhat[i, ] <= eps, tau1 + tau1 * sumc2, tau1 + tau1 * sumc2 + omega * dhat[i, ] + ( 2 * omega * tau1 ) * t1 )
b2 <- (tau1 + tau1 * sumc2) + (dhat[i, ] > eps) * omega * (dhat[i, ] + ( 2 * tau1 ) * t1)
#b2 <- ifelse( dhat[i, ] <= eps, tau1 + tau1 * sumc2, tau1 + tau1 * sumc2 + omega * dhat[i, ] + ( 2 * omega * tau1 ) * t1 )
#b3 <- ifelse( dhat[i, ] <= eps, 0, sqrtn * (t2 + tau4 * t1) )
b3 <- (dhat[i, ] > eps) * sqrtn * (t2 + tau4 * t1)
upper <- alpha2 * b1 + alpha3 * b2 + g * b3
ksi <- upper / lower
tt[[i]] <- transform(ksi, disobj[[i]], w = w[i, ])
dhat[i, ] <- tt[[i]]$res ## dhat update
# ps$c1[i] <- ps$nrmd2[i] <- sum( w[i, ] * ( dhat[i, ] - d[i, ] )^2 )
# ps$wsum[i] <- sum( w[i, ] * dhat[i, ] )
# ps$nrmw2[i] <- sum( w[i, ] * dhat[i, ]^2 )
# ps$nrmm2[i] <- ps$wsum[i]^2 / ps$sumw[i]
# ps$nrmv2[i] <- ps$nrmw2[i] - ps$nrmm2[i]
# ps$c2[i] <- ( ps$nrmv2[i] + omega * ps$nrmm2[i] ) / ps$nrmv2[i]
}
}
}
ps <- pstress( dhat, d, w, omega, lambda, wpp, conditionality ) #pstress value
lnew <- ps$pstress
if (verbose) cat("Iteration: ", formatC(itel, digits=6, width=6),
" P-Stress:", formatC(lnew, digits=6, width=12, format="f"),
" Dif:", formatC(lold - lnew, digits=6,width=12, format="f"),
"\n")
#if ( ( (lold-lnew) < eps & itel > 1) || (itel==itmax)) break()
if ( itel == itmax ) break()
if ( 2 * ( lold - lnew ) <= eps * ( lold + lnew + 1e-15 ) ) break()
lold <- lnew #update stress
itel <- itel+1
}
#-------------------- end majorization --------------------------
## --- Finish for parallization
if (parallelize) parallel::stopCluster(cl)
# point stress
resmat <- as.matrix(d - dhat)^2 #point stress
spp.col <- colSums(w * resmat, na.rm = TRUE) / colSums(w, na.rm = TRUE)
spp.col <- spp.col/sum(spp.col)*100
spp.row <- rowSums(w * resmat, na.rm = TRUE) / rowSums(w, na.rm = TRUE)
spp.row <- spp.row/sum(spp.row)*100
# Add names
colnames(y) <- colnames(x) <- paste("D",1:(dim(y)[2]),sep="")
names(spp.row) <- rownames(x) <- rownames(diss) <- rownames(d) <- rnames
names(spp.col) <- rownames(y) <- colnames(delta)
# final dilation
# first re-scale distances (unconditional)
ssqd <- sum( w * d^2 )
scale <- sqrt( sum( w ) / ssqd )
d <- scale * d
x <- scale * x
y <- scale * y
dhat <- scale * dhat
if (itel == itmax) warning("Iteration limit reached! Increase itmax argument!")
## stress normalization
lnew <- sqrt(sum(w*(dhat - d)^2, na.rm = TRUE)/wpp)
lnew <- sqrt(1 - sum(w * dhat * d, na.rm = TRUE)^2 /
(sum(w * dhat^2, na.rm = TRUE) * sum(w * d^2, na.rm = TRUE)) )
## congruence coefficients
diss0 <- diss
diss0[is.na(diss0)] <- 0
congnum <- diag(diss0 %*% t(d))
congdenom <- sqrt(diag(diss0 %*% t(diss0)) * diag(d %*% t(d)))
congvec <- congnum/congdenom
##
iord <- vector(mode = "list", length = length(tt))
if (conditionality == "matrix"){
iord[[1]] <- tt[[1]]$iord.prim
} else {
iord <- vector(mode = "list", n)
for (i in 1:n) iord[[i]] <- tt[[i]]$iord.prim
}
#return configuration distances, row and column configurations, stress
result <- list(obsdiss = diss, confdist = d, dhat = dhat, iord = iord, conf.row = x, conf.col = y, stress = lnew,
pstress = ps$pstress, spp.row = spp.row, spp.col = spp.col, congvec = congvec, weightmat = w,
ndim = p, model = "Rectangular smacof", niter = itel, nind = n, nobj = m, trans = trans, conditionality = conditionality, call = match.call())
class(result) <- c("smacofR")
return(result)
}
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.