Profile Analysis via Multidimensional Scaling with pams

Se-Kang Kim and Donghoh Kim

2026-10-05

What pams adds

Profile Analysis via Multidimensional Scaling (PAMS) separates each person’s scores into an overall level and an ipsatized pattern. Without pams, an R user would need to program the standardization, level–pattern decomposition, inter-variable dissimilarities, MDS fitting, person-level bootstrap, sign alignment, BCa intervals, and person-weight regressions as separate steps. BootSmacof() integrates those operations and returns a fitted object with summary() and plot() methods.

Let \(z_i\) be person \(i\)’s \(J\)-variable score vector. The level and pattern are

\[ \ell_i = J^{-1}\sum_j z_{ij}, \qquad p_i = z_i - \ell_i\mathbf{1}. \]

If the retained MDS coordinates form the \(J \times K\) matrix \(C\), PAMS represents the pattern as

\[ p_i = Cw_i + e_i. \]

The values w1, …, wK are the unstandardized no-intercept ordinary least-squares coefficients in \(w_i\). The corDim values are partial correlations between the person pattern and each coordinate vector after controlling for the remaining vectors.

Workflow

The recommended workflow is:

  1. prepare the persons-by-variables matrix and standardize when measurement scales differ;
  2. inspect preliminary MDS solutions and dimensionality evidence;
  3. choose interpretable axis directions;
  4. run BootSmacof() with at least 1,000 bootstrap samples for an analysis;
  5. use summary(), plot(), and the person-level output; and
  6. treat downstream group comparisons as optional analyses of returned quantities rather than part of the PAMS estimator.

The small bootstrap count below keeps vignette-building fast. It is not a recommended value for substantive work.

Reproducible example

The built-in USArrests data are used only to demonstrate the software workflow. States are treated as persons and the four variables as related measurements.

library(pams)

example_data <- as.data.frame(USArrests[1:12, ])
example_names <- colnames(example_data)

Because the variables use different units, the preliminary dissimilarities are computed after column standardization. Stress can be inspected across candidate dimensions; with four variables, two dimensions are used here.

standardized <- scale(example_data)
proximity <- dist(t(standardized))
preliminary <- smacof::smacofSym(proximity, ndim = 2, type = "ordinal")
preliminary$stress
#> [1] 2.631087e-05
preliminary$conf
#>                  D1          D2
#> Murder   -0.4703018 -0.44644956
#> Assault  -0.2884457  0.05335699
#> UrbanPop  0.8249120 -0.14355135
#> Rape     -0.0661645  0.53664391

Axis signs are arbitrary. Inspect the preliminary coordinates and select 1 or -1 for each axis so that prespecified anchor variables appear on the desired side. The example retains both signs.

set.seed(2026)
fit <- suppressWarnings(BootSmacof(
  testdata = example_data,
  participant = 1:3,
  mds = "smacof",
  type = "ordinal",
  distance = "euclid",
  scale = TRUE,
  nprofile = 2,
  direction = c(1, 1),
  cl = 0.95,
  nBoot = 10,
  testname = example_names
))

Summaries and plots

summary(fit)
#> Profile Analysis via Multidimensional Scaling
#> 
#> Call:
#> BootSmacof(testdata = example_data, participant = 1:3, mds = "smacof", 
#>     type = "ordinal", distance = "euclid", scale = TRUE, nprofile = 2, 
#>     direction = c(1, 1), cl = 0.95, nBoot = 10, testname = example_names)
#> 
#> Model: smacof( ordinal )with euclid distances
#> Persons: 12  | Variables: 4  | Core profiles: 2 
#> Bootstrap samples: 10  | Confidence level: 0.95 
#> Original-sample stress: 0 
#> Mean person-level R-squared: 0.793 
#> Core-profile collinearity R-squared:
#>    D1    D2 
#> 0.003 0.003 
#> Coordinates with pointwise BCa intervals excluding zero:
#> Profile1 Profile2 
#>        2        1
round(fit$MDSsummary[[1]], 3)
#>             Ori   Mean    SE  Lower  Upper BCaLower BCaUpper
#> Murder   -0.470 -0.470 0.211 -0.782 -0.170   -0.787   -0.175
#> Assault  -0.288 -0.217 0.195 -0.463  0.130   -0.480    0.100
#> UrbanPop  0.825  0.771 0.068  0.668  0.877    0.769    0.892
#> Rape     -0.066 -0.084 0.341 -0.577  0.325   -0.588    0.324
round(fit$Weight[1:5, ], 3)
#>        w1     w2  level   R^2 corDim1 corDim2
#> #1 -1.116 -0.960  0.005 0.999  -0.999  -0.998
#> #2 -1.721  2.148  0.265 0.982  -0.985   0.981
#> #3  0.363  0.478  0.501 0.337   0.449   0.430
#> #4 -0.979 -0.191 -0.551 0.998  -0.999  -0.950
#> #5  0.888  0.942  0.903 0.991   0.992   0.987
plot(fit, profiles = 1:2, interval = "BCa")

The original-sample MDS fit stored in fit$MDS is the result of applying the selected directions to the same SMACOF analysis. The stress therefore agrees with the preliminary fit. Coordinate signs should be resolved before direct comparison.

c(preliminary = preliminary$stress, BootSmacof = fit$MDS$stress)
#>  preliminary   BootSmacof 
#> 2.631087e-05 2.631087e-05

signs <- sign(diag(cor(preliminary$conf, fit$MDS$conf)))
signs[signs == 0] <- 1
aligned_preliminary <- sweep(preliminary$conf, 2, signs, `*`)
max(abs(aligned_preliminary - fit$MDS$conf))
#> [1] 0

Tidyverse-style handling

The fitted object remains a named list for backward compatibility, so it can be passed through a base-R pipe and its tabular components can be converted directly to tibbles. The tibble package is optional rather than a required dependency.

fit |> summary()
#> Profile Analysis via Multidimensional Scaling
#> 
#> Call:
#> BootSmacof(testdata = example_data, participant = 1:3, mds = "smacof", 
#>     type = "ordinal", distance = "euclid", scale = TRUE, nprofile = 2, 
#>     direction = c(1, 1), cl = 0.95, nBoot = 10, testname = example_names)
#> 
#> Model: smacof( ordinal )with euclid distances
#> Persons: 12  | Variables: 4  | Core profiles: 2 
#> Bootstrap samples: 10  | Confidence level: 0.95 
#> Original-sample stress: 0 
#> Mean person-level R-squared: 0.793 
#> Core-profile collinearity R-squared:
#>    D1    D2 
#> 0.003 0.003 
#> Coordinates with pointwise BCa intervals excluding zero:
#> Profile1 Profile2 
#>        2        1

if (requireNamespace("tibble", quietly = TRUE)) {
  weights_tbl <- tibble::as_tibble(fit$Weight, rownames = "participant")
  weights_tbl
}
#> # A tibble: 12 × 7
#>    participant     w1     w2    level `R^2` corDim1 corDim2
#>    <chr>        <dbl>  <dbl>    <dbl> <dbl>   <dbl>   <dbl>
#>  1 #1          -1.12  -0.960  0.00454 0.999  -0.999  -0.998
#>  2 #2          -1.72   2.15   0.265   0.982  -0.985   0.981
#>  3 #3           0.363  0.478  0.501   0.337   0.449   0.430
#>  4 #4          -0.979 -0.191 -0.551   0.998  -0.999  -0.950
#>  5 #5           0.888  0.942  0.903   0.991   0.992   0.987
#>  6 #6           0.508  1.16   0.360   0.842   0.759   0.887
#>  7 #7           1.50  -0.759 -0.816   0.995   0.997  -0.980
#>  8 #8           0.457 -0.469 -0.270   0.256   0.436  -0.337
#>  9 #9          -0.498 -0.715  1.04    0.748  -0.759  -0.769
#> 10 #10         -1.36  -1.38   0.299   0.840  -0.874  -0.795
#> 11 #11          1.80  -0.496 -0.588   0.755   0.867  -0.326
#> 12 #12          0.157  0.237 -1.15    0.775   0.774   0.798

Sign and alignment limitation

For every bootstrap and jackknife sample, BootSmacof() aligns each dimension to its original-sample counterpart using the sign of their correlation. It does not perform general rotational alignment or dimension-permutation matching. Coordinate-wise intervals should therefore be interpreted cautiously when dimensions are weak or nearly interchangeable. Reversing an axis changes only its reported sign; it does not change interpoint distances, stress, overall fit, or whether an interval excludes zero.