Advanced usage of gipsDA
Overview
This vignette gives a more detailed overview of the main modeling
options in gipsDA.
It covers:
- data preparation,
- the differences between
gipslda(),
gipsqda(), and gipsmultqda(),
- the
MAP, optimizer, max_iter,
prior, and weighted_avg arguments,
- interpretation of permutation output,
- prediction methods,
- leave-one-out prediction,
- model inspection and diagnostics.
For a shorter first example, see the Getting started vignette.
Example data
We use the built-in iris data set.
set.seed(42)
train_id <- unlist(
lapply(split(seq_len(nrow(iris)), iris$Species), sample, size = 35),
use.names = FALSE
)
train <- iris[train_id, ]
test <- iris[-train_id, ]
table(train$Species)
#>
#> setosa versicolor virginica
#> 35 35 35
table(test$Species)
#>
#> setosa versicolor virginica
#> 15 15 15
Preparing data
gipsDA models assume numeric predictors and a
categorical grouping variable.
Before fitting a model, it is usually useful to:
- remove identifier columns,
- remove or impute missing values,
- remove constant or almost constant columns,
- encode categorical predictors,
- check whether predictors are on comparable scales.
The last point is especially important for gipsDA. The
method searches for permutation symmetries between variables. Such
symmetries are most meaningful when variables are comparable, for
example when they are measured in the same units or represent analogous
sensor readings.
Checking predictor types
str(train)
#> 'data.frame': 105 obs. of 5 variables:
#> $ Sepal.Length: num 5.3 5.5 5.1 4.8 4.9 5 5.1 5.1 4.6 5.1 ...
#> $ Sepal.Width : num 3.7 3.5 3.5 3.4 3.1 3.2 3.5 3.3 3.4 3.8 ...
#> $ Petal.Length: num 1.5 1.3 1.4 1.9 1.5 1.2 1.4 1.7 1.4 1.9 ...
#> $ Petal.Width : num 0.2 0.2 0.2 0.2 0.1 0.2 0.3 0.5 0.3 0.4 ...
#> $ Species : Factor w/ 3 levels "setosa","versicolor",..: 1 1 1 1 1 1 1 1 1 1 ...
For formula-based usage, the response variable should be a factor or
a categorical variable.
is.factor(train$Species)
#> [1] TRUE
The predictors in iris are already numeric.
vapply(train[, 1:4], is.numeric, logical(1))
#> Sepal.Length Sepal.Width Petal.Length Petal.Width
#> TRUE TRUE TRUE TRUE
Scaling predictors
Scaling may be useful when predictors are measured on very different
scales. However, scaling should be done carefully: the center and scale
should be estimated only on the training data and then applied to new
data.
x_train <- train[, 1:4]
x_test <- test[, 1:4]
train_center <- vapply(x_train, mean, numeric(1))
train_scale <- vapply(x_train, sd, numeric(1))
train_scaled <- train
test_scaled <- test
train_scaled[, 1:4] <- scale(
x_train,
center = train_center,
scale = train_scale
)
test_scaled[, 1:4] <- scale(
x_test,
center = train_center,
scale = train_scale
)
Then fit the model on the scaled data.
fit_scaled <- gipslda(Species ~ ., data = train_scaled)
pred_scaled <- predict(fit_scaled, test_scaled)
mean(pred_scaled$class == test_scaled$Species)
#> [1] 0.9777778
Scaling is not always necessary. It depends on whether the original
variables are already comparable and whether scaling is meaningful for
the application.
Model choice
The package provides three main classifiers.
gipslda() |
all classes share one projected covariance matrix |
classes differ mainly in their means |
gipsqda() |
each class has its own projected covariance matrix and its own
permutation structure |
classes may have different covariance patterns |
gipsmultqda() |
each class has its own covariance matrix, but all classes share one
permutation structure |
classes may differ in scale, but share a dependency pattern |
Fit all three models on the same data.
lda_fit <- gipslda(Species ~ ., data = train)
qda_fit <- gipsqda(Species ~ ., data = train)
joint_qda_fit <- gipsmultqda(Species ~ ., data = train)
Compare test-set accuracy.
lda_pred <- predict(lda_fit, test)
qda_pred <- predict(qda_fit, test)
joint_qda_pred <- predict(joint_qda_fit, test)
c(
gipslda = mean(lda_pred$class == test$Species),
gipsqda = mean(qda_pred$class == test$Species),
gipsmultqda = mean(joint_qda_pred$class == test$Species)
)
#> gipslda gipsqda gipsmultqda
#> 0.9777778 0.9777778 1.0000000
Model hierarchy
The inclusion relations between the model classes are shown
below.
The figure shows that gipsqda() is contained in QDA,
gipslda() is contained in LDA, and
gipsmultqda() lies between gipsqda() and
gipslda().
In practical terms, moving inward in the diagram means imposing
stronger assumptions on the covariance structure. Stronger assumptions
can reduce estimation variance, especially when the number of
observations is small, but they may be too restrictive if the assumed
structure is not present in the data.
MAP: Maximum A Posteriori
MAP stands for Maximum A
Posteriori.
The MAP argument controls how the covariance matrix is
projected after the permutation search.
When MAP = TRUE, the model selects the single most
probable permutation structure and projects the covariance matrix onto
the invariant space determined by that structure.
lda_map <- gipslda(
Species ~ .,
data = train,
MAP = TRUE
)
lda_map
#> Call:
#> gipslda(Species ~ ., data = train, MAP = TRUE)
#>
#> Model: gipslda
#> Number of observations: 105
#> Number of groups: 3
#> Number of predictors: 4
#>
#> Fitting options:
#> $MAP
#> [1] TRUE
#>
#> $optimizer
#> [1] "BF"
#>
#> $max_iter
#> NULL
#>
#> $weighted_avg
#> [1] FALSE
#>
#> $store_probabilities
#> [1] TRUE
#>
#>
#> Prior probabilities of groups:
#> setosa versicolor virginica
#> 0.3333333 0.3333333 0.3333333
#>
#> Class counts:
#> setosa versicolor virginica
#> 35 35 35
#>
#> Group means:
#> Sepal.Length Sepal.Width Petal.Length Petal.Width
#> setosa 5.020000 3.420000 1.482857 0.2428571
#> versicolor 5.911429 2.771429 4.302857 1.3371429
#> virginica 6.725714 3.020000 5.654286 2.0685714
#>
#> Selected MAP permutation: (1,2,4,3)
#>
#> Posterior probabilities of retained permutations:
#> (1,2,4,3) (1,3)(2,4) (1,2,3,4) (1,2)(3,4) (1,4)(2,3)
#> 0.549732245 0.423533797 0.018304175 0.004073983 0.003658725
#>
#> Coefficients of linear discriminants:
#> LD1 LD2
#> Sepal.Length -0.2534381 -0.3641984
#> Sepal.Width 2.3946991 -2.3105323
#> Petal.Length -1.3688435 0.8934580
#> Petal.Width -3.5168812 -2.3906963
#>
#> Proportion of trace:
#> LD1 LD2
#> 0.9888 0.0112
When MAP = FALSE, the model uses posterior probabilities
over retained permutation structures and computes a posterior-weighted
projection.
lda_avg <- gipslda(
Species ~ .,
data = train,
MAP = FALSE
)
lda_avg
#> Call:
#> gipslda(Species ~ ., data = train, MAP = FALSE)
#>
#> Model: gipslda
#> Number of observations: 105
#> Number of groups: 3
#> Number of predictors: 4
#>
#> Fitting options:
#> $MAP
#> [1] FALSE
#>
#> $optimizer
#> [1] "BF"
#>
#> $max_iter
#> NULL
#>
#> $weighted_avg
#> [1] FALSE
#>
#> $store_probabilities
#> [1] TRUE
#>
#>
#> Prior probabilities of groups:
#> setosa versicolor virginica
#> 0.3333333 0.3333333 0.3333333
#>
#> Class counts:
#> setosa versicolor virginica
#> 35 35 35
#>
#> Group means:
#> Sepal.Length Sepal.Width Petal.Length Petal.Width
#> setosa 5.020000 3.420000 1.482857 0.2428571
#> versicolor 5.911429 2.771429 4.302857 1.3371429
#> virginica 6.725714 3.020000 5.654286 2.0685714
#>
#> Posterior probabilities of retained permutations:
#> (1,2,4,3) (1,3)(2,4) (1,2,3,4) (1,2)(3,4) (1,4)(2,3)
#> 0.549732245 0.423533797 0.018304175 0.004073983 0.003658725
#>
#> Coefficients of linear discriminants:
#> LD1 LD2
#> Sepal.Length 0.06710214 -0.5349658
#> Sepal.Width 2.15454818 -2.2243646
#> Petal.Length -1.59385584 0.9770055
#> Petal.Width -3.36967333 -2.4060085
#>
#> Proportion of trace:
#> LD1 LD2
#> 0.9886 0.0114
Conceptually, if \(S\) is an
empirical covariance matrix and \(S_c\)
is its projection under permutation structure \(c\), then the two approaches can be
summarized as follows.
For MAP = TRUE:
\[
\hat{S} = S_{c^*},
\]
where
\[
c^* = \arg\max_c P(c \mid X).
\]
For MAP = FALSE:
\[
\hat{S} =
\sum_c P(c \mid X) S_c.
\]
The second option averages over several possible symmetry structures
instead of using only one selected structure.
By default, gipsDA stores posterior probabilities of
retained permutations. For faster MAP-only fitting, set
store_probabilities = FALSE. In that case, the selected MAP
permutation is still stored and shown, but posterior probabilities are
not stored in the fitted model object.
With one predictor, gipslda() and gipsqda()
require MAP = TRUE. The identity permutation
() is the only possible permutation and has probability
1, which is stored when
store_probabilities = TRUE. gipsmultqda()
requires at least two predictors.
Both QDA fitters reject unused levels in the grouping factor with an
error listing those levels. Use droplevels(grouping) before
fitting to remove them.
Interpreting permutation output
Printed model output may contain permutations such as:
(1,2)
(1,2)(3,4)
(1,2,4,3)
()
This notation is called cycle notation.
For example, the cycle
(1,2,4,3)
means that the permutation maps:
1 -> 2
2 -> 4
4 -> 3
3 -> 1
For cycles of length greater than two, this should be understood as
invariance under the cyclic permutation and its repeated applications.
It does not necessarily mean full exchangeability under every possible
pairwise swap of features in the cycle.
A product of cycles such as:
(1,2)(3,4)
means that feature 1 is swapped with feature 2, and feature 3 is
swapped with feature 4.
The empty permutation:
()
is the identity permutation. It means that no non-trivial permutation
symmetry was selected.
In the context of gipsDA, a selected permutation
structure describes invariance constraints imposed on the covariance
estimator.
For gipslda(), the permutation search is performed after
centering observations by their class means and scaling the resulting
within-class residuals to unit marginal variance. Therefore, the
selected permutation describes symmetry in the standardized within-class
covariance structure, not necessarily symmetry of the raw covariance
matrix in the original units.
For gipsqda() and gipsmultqda(), the class
covariance matrices are projected on the original predictor scale.
Because of this difference, selected permutations from LDA and QDA
should not always be interpreted in exactly the same way when predictors
are measured on different scales.
Optimizer
The optimizer argument controls how permutation
structures are searched.
"BF" |
brute-force search |
small number of dimensions, default for p <= 10 |
"MH" |
Metropolis-Hastings search |
larger number of dimensions, default for p > 10 |
The brute-force optimizer searches the relevant permutation space
exhaustively. It is deterministic, but its cost grows quickly with the
number of features.
fit_bf <- gipsqda(
Species ~ .,
data = train,
optimizer = "BF"
)
fit_bf
#> Call:
#> gipsqda(Species ~ ., data = train, optimizer = "BF")
#>
#> Model: gipsqda
#> Number of observations: 105
#> Number of groups: 3
#> Number of predictors: 4
#>
#> Fitting options:
#> $MAP
#> [1] TRUE
#>
#> $optimizer
#> [1] "BF"
#>
#> $max_iter
#> NULL
#>
#> $store_probabilities
#> [1] TRUE
#>
#>
#> Prior probabilities of groups:
#> setosa versicolor virginica
#> 0.3333333 0.3333333 0.3333333
#>
#> Class counts:
#> setosa versicolor virginica
#> 35 35 35
#>
#> Group means:
#> Sepal.Length Sepal.Width Petal.Length Petal.Width
#> setosa 5.020000 3.420000 1.482857 0.2428571
#> versicolor 5.911429 2.771429 4.302857 1.3371429
#> virginica 6.725714 3.020000 5.654286 2.0685714
#>
#> Group: setosa
#> Selected MAP permutation: (1,2)
#>
#> Group: versicolor
#> Selected MAP permutation: (2,3)
#>
#> Group: virginica
#> Selected MAP permutation: (1,3)(2,4)
#>
#> Posterior probabilities of retained permutations:
#>
#> Group: setosa
#> (1,2) (1,2)(3,4) (3,4) ()
#> 0.562114561 0.404501310 0.024221480 0.008570546
#>
#> Group: versicolor
#> (2,3) () (1,3) (1,2,3) (1,2) (1,3)(2,4)
#> 0.889698934 0.061618087 0.029058419 0.015073718 0.002647551 0.001632264
#>
#> Group: virginica
#> (1,3)(2,4) (1,3) (2,4) ()
#> 0.849678403 0.089183371 0.058826178 0.001555686
#>
#> Log determinants of projected covariance matrices:
#> [1] -10.197338 -8.531728 -8.239145
For larger problems, use the Metropolis-Hastings optimizer.
fit_mh <- gipsqda(
Species ~ .,
data = train,
optimizer = "MH",
max_iter = 1000
)
max_iter
max_iter controls the number of Metropolis-Hastings
iterations when optimizer = "MH".
Increasing max_iter gives the stochastic search more
time to explore the permutation space, but also increases runtime.
fit_mh_100 <- gipsqda(
Species ~ .,
data = train,
optimizer = "MH",
max_iter = 100
)
fit_mh_1000 <- gipsqda(
Species ~ .,
data = train,
optimizer = "MH",
max_iter = 1000
)
For optimizer = "BF", max_iter is
ignored.
Class priors
The prior argument specifies prior probabilities of
classes.
By default, priors are estimated from the training data.
lda_fit$prior
#> setosa versicolor virginica
#> 0.3333333 0.3333333 0.3333333
You can set them manually.
equal_prior <- rep(1 / length(levels(train$Species)), length(levels(train$Species)))
names(equal_prior) <- levels(train$Species)
lda_equal_prior <- gipslda(
Species ~ .,
data = train,
prior = equal_prior
)
lda_equal_prior$prior
#> setosa versicolor virginica
#> 0.3333333 0.3333333 0.3333333
The prior vector should contain one value per class and should sum to
one.
sum(equal_prior)
#> [1] 1
Class priors affect posterior probabilities and may affect predicted
classes, especially when classes overlap.
weighted_avg in gipslda()
The weighted_avg argument is specific to
gipslda().
It controls how the pooled covariance matrix is constructed before
the gips projection is applied.
Let:
- \(K\) be the number of
classes,
- \(n_k\) be the number of
observations in class \(k\),
- \(n = \sum_{k=1}^{K} n_k\) be the
total number of observations,
- \(S_k\) be the sample covariance
matrix in class \(k\).
With weighted_avg = FALSE, gipslda() uses
the classic pooled covariance estimator:
\[
S_{\mathrm{classic}}
=
\frac{1}{n - K}
\sum_{k = 1}^{K}
(n_k - 1) S_k.
\]
With weighted_avg = TRUE, gipslda()
uses:
\[
S_{\mathrm{weighted}}
=
\frac{1}{n}
\sum_{k = 1}^{K}
n_k S_k.
\]
Fit both variants.
lda_classic <- gipslda(
Species ~ .,
data = train,
weighted_avg = FALSE
)
lda_weighted <- gipslda(
Species ~ .,
data = train,
weighted_avg = TRUE
)
Compare predictions.
pred_classic <- predict(lda_classic, test)
pred_weighted <- predict(lda_weighted, test)
c(
classic = mean(pred_classic$class == test$Species),
weighted = mean(pred_weighted$class == test$Species)
)
#> classic weighted
#> 0.9777778 0.9777778
The two variants can behave differently when class sizes are
imbalanced or when class-specific covariance estimates differ
substantially.
Matrix interface
The matrix interface separates predictors from class labels.
x <- as.matrix(iris[, 1:4])
grouping <- iris$Species
fit_matrix_lda <- gipslda(x, grouping)
fit_matrix_qda <- gipsqda(x, grouping)
fit_matrix_joint <- gipsmultqda(x, grouping)
Prediction can then be performed on a matrix with the same
columns.
predict(fit_matrix_lda, x[1:5, ])$class
#> [1] setosa setosa setosa setosa setosa
#> Levels: setosa versicolor virginica
predict(fit_matrix_qda, x[1:5, ])$class
#> [1] setosa setosa setosa setosa setosa
#> Levels: setosa versicolor virginica
predict(fit_matrix_joint, x[1:5, ])$class
#> [1] setosa setosa setosa setosa setosa
#> Levels: setosa versicolor virginica
Data-frame interface
Predictors can also be passed as a data frame, with labels supplied
separately.
x_df <- iris[, 1:4]
y <- iris$Species
fit_df_lda <- gipslda(x_df, y)
fit_df_qda <- gipsqda(x_df, y)
fit_df_joint <- gipsmultqda(x_df, y)
predict(fit_df_lda, x_df[1:5, ])$class
#> [1] setosa setosa setosa setosa setosa
#> Levels: setosa versicolor virginica
predict(fit_df_qda, x_df[1:5, ])$class
#> [1] setosa setosa setosa setosa setosa
#> Levels: setosa versicolor virginica
predict(fit_df_joint, x_df[1:5, ])$class
#> [1] setosa setosa setosa setosa setosa
#> Levels: setosa versicolor virginica
Prediction output
Prediction returns a list.
pred <- predict(lda_fit, test)
names(pred)
#> [1] "class" "posterior" "x"
The most important components are:
class |
predicted class labels |
posterior |
posterior class probabilities |
x |
discriminant coordinates, when available |
head(pred$class)
#> [1] setosa setosa setosa setosa setosa setosa
#> Levels: setosa versicolor virginica
head(pred$posterior)
#> setosa versicolor virginica
#> 6 1 4.955555e-22 1.956129e-41
#> 9 1 1.894514e-16 8.017336e-36
#> 14 1 8.720623e-21 8.034755e-42
#> 16 1 1.868380e-28 2.037539e-49
#> 17 1 1.869751e-24 1.343241e-44
#> 19 1 5.185944e-22 1.062783e-41
For gipsqda() and gipsmultqda(), the output
has the same main structure.
names(qda_pred)
#> [1] "class" "posterior"
names(joint_qda_pred)
#> [1] "class" "posterior"
Prediction methods for gipslda()
For gipslda() objects, predict() supports
the same main prediction-method names as
MASS::predict.lda():
"plug-in","predictive","debiased".
The default method is "plug-in".
pred_plugin <- predict(lda_fit, test, method = "plug-in")
pred_predictive <- predict(lda_fit, test, method = "predictive")
pred_debiased <- predict(lda_fit, test, method = "debiased")
The "plug-in" method uses estimated parameters directly
in the discriminant rule.
The "predictive" and "debiased" methods are
alternative LDA prediction rules following the
MASS::predict.lda() interface. See
?MASS::predict.lda for the original description of these
prediction rules.
For easy data sets, the predicted classes may be identical across
methods.
c(
plugin = mean(pred_plugin$class == test$Species),
predictive = mean(pred_predictive$class == test$Species),
debiased = mean(pred_debiased$class == test$Species)
)
#> plugin predictive debiased
#> 0.9777778 0.9555556 0.9777778
Differences may be more visible in posterior probabilities.
head(pred_plugin$posterior)
#> setosa versicolor virginica
#> 6 1 4.955555e-22 1.956129e-41
#> 9 1 1.894514e-16 8.017336e-36
#> 14 1 8.720623e-21 8.034755e-42
#> 16 1 1.868380e-28 2.037539e-49
#> 17 1 1.869751e-24 1.343241e-44
#> 19 1 5.185944e-22 1.062783e-41
head(pred_predictive$posterior)
#> setosa versicolor virginica
#> 6 1 3.321301e-15 3.513263e-23
#> 9 1 3.298368e-12 4.408852e-21
#> 14 1 1.246057e-14 2.184951e-23
#> 16 1 3.198616e-17 1.370273e-24
#> 17 1 3.437129e-16 5.177073e-24
#> 19 1 2.773375e-15 2.170032e-23
head(pred_debiased$posterior)
#> setosa versicolor virginica
#> 6 1 2.097526e-21 3.081126e-40
#> 9 1 5.494451e-16 8.635043e-35
#> 14 1 3.392590e-20 1.299122e-40
#> 16 1 1.221829e-27 5.510501e-48
#> 17 1 9.325508e-24 2.621213e-43
#> 19 1 2.192110e-21 1.704312e-40
Leave-one-out prediction for QDA models
For gipsqda() and gipsmultqda() objects,
leave-one-out cross-validation can be requested by omitting
newdata and using method = "looCV".
In leave-one-out prediction, each training observation is classified
as if it had not been used to fit the model. This provides an internal
estimate of classification performance without creating a separate test
set.
qda_loo <- predict(qda_fit, method = "looCV")
joint_qda_loo <- predict(joint_qda_fit, method = "looCV")
c(
gipsqda_loo_accuracy = mean(qda_loo$class == train$Species),
gipsmultqda_loo_accuracy = mean(joint_qda_loo$class == train$Species)
)
#> gipsqda_loo_accuracy gipsmultqda_loo_accuracy
#> 0.9714286 0.9714286
Use leave-one-out results as an internal diagnostic, not as a
replacement for a proper independent test set when one is available.
Inspecting fitted models
The most readable way to inspect a fitted model is to print it.
print(lda_fit)
#> Call:
#> gipslda(Species ~ ., data = train)
#>
#> Model: gipslda
#> Number of observations: 105
#> Number of groups: 3
#> Number of predictors: 4
#>
#> Fitting options:
#> $MAP
#> [1] TRUE
#>
#> $optimizer
#> [1] "BF"
#>
#> $max_iter
#> NULL
#>
#> $weighted_avg
#> [1] FALSE
#>
#> $store_probabilities
#> [1] TRUE
#>
#>
#> Prior probabilities of groups:
#> setosa versicolor virginica
#> 0.3333333 0.3333333 0.3333333
#>
#> Class counts:
#> setosa versicolor virginica
#> 35 35 35
#>
#> Group means:
#> Sepal.Length Sepal.Width Petal.Length Petal.Width
#> setosa 5.020000 3.420000 1.482857 0.2428571
#> versicolor 5.911429 2.771429 4.302857 1.3371429
#> virginica 6.725714 3.020000 5.654286 2.0685714
#>
#> Selected MAP permutation: (1,2,4,3)
#>
#> Posterior probabilities of retained permutations:
#> (1,2,4,3) (1,3)(2,4) (1,2,3,4) (1,2)(3,4) (1,4)(2,3)
#> 0.549732245 0.423533797 0.018304175 0.004073983 0.003658725
#>
#> Coefficients of linear discriminants:
#> LD1 LD2
#> Sepal.Length -0.2534381 -0.3641984
#> Sepal.Width 2.3946991 -2.3105323
#> Petal.Length -1.3688435 0.8934580
#> Petal.Width -3.5168812 -2.3906963
#>
#> Proportion of trace:
#> LD1 LD2
#> 0.9888 0.0112
print(qda_fit)
#> Call:
#> gipsqda(Species ~ ., data = train)
#>
#> Model: gipsqda
#> Number of observations: 105
#> Number of groups: 3
#> Number of predictors: 4
#>
#> Fitting options:
#> $MAP
#> [1] TRUE
#>
#> $optimizer
#> [1] "BF"
#>
#> $max_iter
#> NULL
#>
#> $store_probabilities
#> [1] TRUE
#>
#>
#> Prior probabilities of groups:
#> setosa versicolor virginica
#> 0.3333333 0.3333333 0.3333333
#>
#> Class counts:
#> setosa versicolor virginica
#> 35 35 35
#>
#> Group means:
#> Sepal.Length Sepal.Width Petal.Length Petal.Width
#> setosa 5.020000 3.420000 1.482857 0.2428571
#> versicolor 5.911429 2.771429 4.302857 1.3371429
#> virginica 6.725714 3.020000 5.654286 2.0685714
#>
#> Group: setosa
#> Selected MAP permutation: (1,2)
#>
#> Group: versicolor
#> Selected MAP permutation: (2,3)
#>
#> Group: virginica
#> Selected MAP permutation: (1,3)(2,4)
#>
#> Posterior probabilities of retained permutations:
#>
#> Group: setosa
#> (1,2) (1,2)(3,4) (3,4) ()
#> 0.562114561 0.404501310 0.024221480 0.008570546
#>
#> Group: versicolor
#> (2,3) () (1,3) (1,2,3) (1,2) (1,3)(2,4)
#> 0.889698934 0.061618087 0.029058419 0.015073718 0.002647551 0.001632264
#>
#> Group: virginica
#> (1,3)(2,4) (1,3) (2,4) ()
#> 0.849678403 0.089183371 0.058826178 0.001555686
#>
#> Log determinants of projected covariance matrices:
#> [1] -10.197338 -8.531728 -8.239145
print(joint_qda_fit)
#> Call:
#> gipsmultqda(Species ~ ., data = train)
#>
#> Model: gipsmultqda
#> Number of observations: 105
#> Number of groups: 3
#> Number of predictors: 4
#>
#> Fitting options:
#> $MAP
#> [1] TRUE
#>
#> $optimizer
#> [1] "BF"
#>
#> $max_iter
#> NULL
#>
#> $store_probabilities
#> [1] TRUE
#>
#>
#> Prior probabilities of groups:
#> setosa versicolor virginica
#> 0.3333333 0.3333333 0.3333333
#>
#> Class counts:
#> setosa versicolor virginica
#> 35 35 35
#>
#> Group means:
#> Sepal.Length Sepal.Width Petal.Length Petal.Width
#> setosa 5.020000 3.420000 1.482857 0.2428571
#> versicolor 5.911429 2.771429 4.302857 1.3371429
#> virginica 6.725714 3.020000 5.654286 2.0685714
#>
#> Selected MAP permutation: (1,3)
#>
#> Posterior probabilities of retained permutations:
#> (1,3) ()
#> 0.6386674 0.3609011
#>
#> Log determinants of projected covariance matrices:
#> [1] -10.093275 -8.701035 -8.306090
The printed output shows the model call, prior probabilities, group
means, and information about selected or averaged permutation
structures.
For a structural view of the object, use summary().
summary(lda_fit)
#> Call:
#> gipslda(Species ~ ., data = train)
#>
#> Model: gipslda
#> Number of observations: 105
#> Number of groups: 3
#> Number of predictors: 4
#>
#> Fitting options:
#> $MAP
#> [1] TRUE
#>
#> $optimizer
#> [1] "BF"
#>
#> $max_iter
#> NULL
#>
#> $weighted_avg
#> [1] FALSE
#>
#> $store_probabilities
#> [1] TRUE
#>
#>
#> Class counts:
#> setosa versicolor virginica
#> 35 35 35
#>
#> Prior probabilities of groups:
#> setosa versicolor virginica
#> 0.3333333 0.3333333 0.3333333
#>
#> Group means:
#> Sepal.Length Sepal.Width Petal.Length Petal.Width
#> setosa 5.020000 3.420000 1.482857 0.2428571
#> versicolor 5.911429 2.771429 4.302857 1.3371429
#> virginica 6.725714 3.020000 5.654286 2.0685714
#>
#> Proportion of trace:
#> LD1 LD2
#> 0.9888 0.0112
#>
#> Selected MAP permutation: (1,2,4,3)
#>
#> Posterior probabilities of retained permutations:
#> (1,2,4,3) (1,3)(2,4) (1,2,3,4) (1,2)(3,4) (1,4)(2,3)
#> 0.549732245 0.423533797 0.018304175 0.004073983 0.003658725
summary(qda_fit)
#> Call:
#> gipsqda(Species ~ ., data = train)
#>
#> Model: gipsqda
#> Number of observations: 105
#> Number of groups: 3
#> Number of predictors: 4
#>
#> Fitting options:
#> $MAP
#> [1] TRUE
#>
#> $optimizer
#> [1] "BF"
#>
#> $max_iter
#> NULL
#>
#> $store_probabilities
#> [1] TRUE
#>
#>
#> Class counts:
#> setosa versicolor virginica
#> 35 35 35
#>
#> Prior probabilities of groups:
#> setosa versicolor virginica
#> 0.3333333 0.3333333 0.3333333
#>
#> Group means:
#> Sepal.Length Sepal.Width Petal.Length Petal.Width
#> setosa 5.020000 3.420000 1.482857 0.2428571
#> versicolor 5.911429 2.771429 4.302857 1.3371429
#> virginica 6.725714 3.020000 5.654286 2.0685714
#>
#> Log determinants of projected covariance matrices:
#> [1] -10.197338 -8.531728 -8.239145
#>
#> Scaling array dimensions:
#> [1] 4 4 3
#>
#> Group: setosa
#> Selected MAP permutation: (1,2)
#>
#> Group: versicolor
#> Selected MAP permutation: (2,3)
#>
#> Group: virginica
#> Selected MAP permutation: (1,3)(2,4)
#>
#> Posterior probabilities of retained permutations:
#>
#> Group: setosa
#> (1,2) (1,2)(3,4) (3,4) ()
#> 0.562114561 0.404501310 0.024221480 0.008570546
#>
#> Group: versicolor
#> (2,3) () (1,3) (1,2,3) (1,2) (1,3)(2,4)
#> 0.889698934 0.061618087 0.029058419 0.015073718 0.002647551 0.001632264
#>
#> Group: virginica
#> (1,3)(2,4) (1,3) (2,4) ()
#> 0.849678403 0.089183371 0.058826178 0.001555686
summary(joint_qda_fit)
#> Call:
#> gipsmultqda(Species ~ ., data = train)
#>
#> Model: gipsmultqda
#> Number of observations: 105
#> Number of groups: 3
#> Number of predictors: 4
#>
#> Fitting options:
#> $MAP
#> [1] TRUE
#>
#> $optimizer
#> [1] "BF"
#>
#> $max_iter
#> NULL
#>
#> $store_probabilities
#> [1] TRUE
#>
#>
#> Class counts:
#> setosa versicolor virginica
#> 35 35 35
#>
#> Prior probabilities of groups:
#> setosa versicolor virginica
#> 0.3333333 0.3333333 0.3333333
#>
#> Group means:
#> Sepal.Length Sepal.Width Petal.Length Petal.Width
#> setosa 5.020000 3.420000 1.482857 0.2428571
#> versicolor 5.911429 2.771429 4.302857 1.3371429
#> virginica 6.725714 3.020000 5.654286 2.0685714
#>
#> Log determinants of projected covariance matrices:
#> [1] -10.093275 -8.701035 -8.306090
#>
#> Scaling array dimensions:
#> [1] 4 4 3
#>
#> Selected MAP permutation: (1,3)
#>
#> Posterior probabilities of retained permutations:
#> (1,3) ()
#> 0.6386674 0.3609011
Fitted models are list-like S3 objects, so you can inspect their
component names.
names(lda_fit)
#> [1] "prior" "counts"
#> [3] "means" "scaling"
#> [5] "lev" "svd"
#> [7] "N" "optimization_info"
#> [9] "selected_map_permutation" "fit_info"
#> [11] "call" "terms"
#> [13] "xlevels"
names(qda_fit)
#> [1] "prior" "counts"
#> [3] "means" "scaling"
#> [5] "ldet" "lev"
#> [7] "N" "call"
#> [9] "optimization_info" "selected_map_permutation"
#> [11] "fit_info" "terms"
#> [13] "xlevels"
names(joint_qda_fit)
#> [1] "prior" "counts"
#> [3] "means" "scaling"
#> [5] "ldet" "lev"
#> [7] "N" "call"
#> [9] "optimization_info" "selected_map_permutation"
#> [11] "fit_info" "terms"
#> [13] "xlevels"
The exact set of components depends on the model family, but the most
useful components are usually:
prior |
prior probabilities of classes |
counts |
number of observations in each class |
means |
class-wise feature means |
scaling |
scaling or decomposition information used for prediction |
ldet |
log-determinant information for QDA-type models |
lev |
class labels |
call |
original function call |
optimization_info |
information returned by the permutation optimization step |
Inspect the optimization information directly.
lda_fit$optimization_info
#> (1,2,4,3) (1,3)(2,4) (1,2,3,4) (1,2)(3,4) (1,4)(2,3)
#> 0.549732245 0.423533797 0.018304175 0.004073983 0.003658725
qda_fit$optimization_info
#> $setosa
#> (1,2) (1,2)(3,4) (3,4) ()
#> 0.562114561 0.404501310 0.024221480 0.008570546
#>
#> $versicolor
#> (2,3) () (1,3) (1,2,3) (1,2) (1,3)(2,4)
#> 0.889698934 0.061618087 0.029058419 0.015073718 0.002647551 0.001632264
#>
#> $virginica
#> (1,3)(2,4) (1,3) (2,4) ()
#> 0.849678403 0.089183371 0.058826178 0.001555686
joint_qda_fit$optimization_info
#> (1,3) ()
#> 0.6386674 0.3609011
A compact helper can be useful when debugging fitted objects.
inspect_model <- function(object) {
data.frame(
component = names(object),
class = vapply(
object,
function(x) paste(class(x), collapse = ", "),
character(1)
),
length = vapply(object, length, integer(1)),
dim = vapply(
object,
function(x) {
d <- dim(x)
if (is.null(d)) "" else paste(d, collapse = " x ")
},
character(1)
),
row.names = NULL
)
}
inspect_model(lda_fit)
#> component class length dim
#> 1 prior numeric 3
#> 2 counts integer 3
#> 3 means matrix, array 12 3 x 4
#> 4 scaling matrix, array 8 4 x 2
#> 5 lev character 3
#> 6 svd numeric 2
#> 7 N integer 1
#> 8 optimization_info numeric 5
#> 9 selected_map_permutation gips_perm 1
#> 10 fit_info list 5
#> 11 call call 3
#> 12 terms terms, formula 3
#> 13 xlevels list 0
inspect_model(qda_fit)
#> component class length dim
#> 1 prior numeric 3
#> 2 counts integer 3
#> 3 means matrix, array 12 3 x 4
#> 4 scaling array 48 4 x 4 x 3
#> 5 ldet numeric 3
#> 6 lev character 3
#> 7 N integer 1
#> 8 call call 3
#> 9 optimization_info list 3
#> 10 selected_map_permutation list 3
#> 11 fit_info list 4
#> 12 terms terms, formula 3
#> 13 xlevels list 0
inspect_model(joint_qda_fit)
#> component class length dim
#> 1 prior numeric 3
#> 2 counts integer 3
#> 3 means matrix, array 12 3 x 4
#> 4 scaling array 48 4 x 4 x 3
#> 5 ldet numeric 3
#> 6 lev character 3
#> 7 N integer 1
#> 8 call call 3
#> 9 optimization_info numeric 2
#> 10 selected_map_permutation gips_perm 3
#> 11 fit_info list 4
#> 12 terms terms, formula 3
#> 13 xlevels list 0
LDA diagnostics
For gipslda() objects, standard LDA-style diagnostics
are available.
coef(lda_fit)
#> LD1 LD2
#> Sepal.Length -0.2534381 -0.3641984
#> Sepal.Width 2.3946991 -2.3105323
#> Petal.Length -1.3688435 0.8934580
#> Petal.Width -3.5168812 -2.3906963
pairs(lda_fit, type = "std")
These methods are useful for inspecting the fitted discriminant
directions and class separation.
Practical workflow
A typical workflow is:
- prepare numeric predictors and a categorical response,
- split data into training and test sets,
- start with
gipslda(),
- compare with
gipsqda() or gipsmultqda() if
class-specific covariance structure may matter,
- inspect
optimization_info,
- evaluate performance on held-out data.
fit_lda <- gipslda(Species ~ ., data = train)
fit_qda <- gipsqda(Species ~ ., data = train)
fit_joint <- gipsmultqda(Species ~ ., data = train)
pred_lda <- predict(fit_lda, test)
pred_qda <- predict(fit_qda, test)
pred_joint <- predict(fit_joint, test)
c(
gipslda = mean(pred_lda$class == test$Species),
gipsqda = mean(pred_qda$class == test$Species),
gipsmultqda = mean(pred_joint$class == test$Species)
)
#> gipslda gipsqda gipsmultqda
#> 0.9777778 0.9777778 1.0000000
Troubleshooting
The model is slow
Use optimizer = "MH" for larger numbers of features.
fit <- gipsqda(
Species ~ .,
data = train,
optimizer = "MH",
max_iter = 1000
)
Decrease max_iter for faster exploratory runs, then
increase it for final runs.
The output selects ()
The permutation () is the identity permutation. It means
that the selected structure did not impose a non-trivial permutation
symmetry.
This can happen when:
- the data do not contain strong exchangeability patterns,
- the number of observations is too small,
- predictors are not meaningfully comparable,
- the optimizer did not find a more structured permutation.
Predictions are identical across models
This may happen on simple data sets. Compare posterior probabilities
and inspect optimization_info for more detail.
head(pred_lda$posterior)
#> setosa versicolor virginica
#> 6 1 4.955555e-22 1.956129e-41
#> 9 1 1.894514e-16 8.017336e-36
#> 14 1 8.720623e-21 8.034755e-42
#> 16 1 1.868380e-28 2.037539e-49
#> 17 1 1.869751e-24 1.343241e-44
#> 19 1 5.185944e-22 1.062783e-41
head(pred_qda$posterior)
#> setosa versicolor virginica
#> 6 1 1.983970e-19 1.033310e-27
#> 9 1 2.720258e-16 1.166576e-19
#> 14 1 1.613600e-19 2.068832e-22
#> 16 1 6.597486e-25 3.982454e-36
#> 17 1 1.485936e-23 4.934051e-32
#> 19 1 4.669037e-19 2.042588e-29
head(pred_joint$posterior)
#> setosa versicolor virginica
#> 6 1 1.809368e-18 3.630155e-32
#> 9 1 4.185357e-13 9.716563e-23
#> 14 1 8.534255e-16 4.000877e-26
#> 16 1 2.045692e-24 3.136984e-41
#> 17 1 1.110504e-21 8.419286e-37
#> 19 1 2.676609e-18 3.735764e-34
Training and test data have different preprocessing
When preprocessing uses estimated quantities, such as means and
standard deviations for scaling, estimate them on the training data only
and apply the same transformation to test data.
Summary
The main advanced controls are:
MAP = TRUE |
use the single Maximum A Posteriori permutation |
MAP = FALSE |
average covariance projections using posterior probabilities |
optimizer = "BF" |
exhaustive search, typically for p <= 10 |
optimizer = "MH" |
stochastic search, typically for p > 10 |
max_iter |
number of Metropolis-Hastings iterations |
prior |
manually set class prior probabilities |
weighted_avg |
choose the pooled covariance estimator in
gipslda() |
store_probabilities |
whether to store posterior probabilities of retained
permutations |
The three main models are:
gipslda() |
classes can share one covariance structure |
gipsqda() |
each class may have its own covariance structure |
gipsmultqda() |
classes may have different covariance matrices but one shared
permutation structure |
