set.seed(1)
library('tidyverse')
library('tidyLPA')
library('patchwork')
library('ggforce')Latent profile analysis
R
Latent profile analysis
- Underlying model: finite Gaussian mixture model
- We observe \(n\) cases on \(m\) numeric variables and want to identify latent subgroups
- Assume \(k\) a priori unknown profiles, each characterized by an \(m\)-dimensional mean vector and covariance structure
- Conditional on profile membership, observations follow a multivariate normal distribution
- LPA estimates the profile parameters, mixing proportions, and posterior membership probabilities for each case
- \(k\) is typically chosen by comparing candidate models using criteria such as BIC or AIC
- Different constraints can be imposed on within-profile variances and covariances
Iris example
- Iris dataset
- 4 numeric variables, one categorical variable (
Species) defining class membership - For our example let us assume we don’t have the
Speciescolumn and want to infer class membership by latent profile analysis
- 4 numeric variables, one categorical variable (
fitted_models <- iris |>
select(-Species) |>
estimate_profiles(2:5)
fitted_models |> get_fit()# A tibble: 4 × 20
Model Classes LogLik parameters n AIC AWE BIC CAIC CLC KIC
<dbl> <int> <dbl> <dbl> <int> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl>
1 1 2 -489. 13 150 1004. 1145. 1043. 1056. 980. 1020.
2 1 3 -361. 18 150 759. 955. 813. 831. 725. 780.
3 1 4 -356. 23 150 758. 1010. 827. 850. 714. 784.
4 1 5 -301. 28 150 658. 964. 742. 770. 603. 689.
# ℹ 9 more variables: SABIC <dbl>, ICL <dbl>, Entropy <dbl>, prob_min <dbl>,
# prob_max <dbl>, n_min <dbl>, n_max <dbl>, BLRT_val <dbl>, BLRT_p <dbl>est <- get_estimates(fitted_models) |>
filter(Category == 'Means') |>
pivot_wider(id_cols = c('Class', 'Model', 'Classes'),
names_from = 'Parameter',
values_from = 'Estimate')
Extract probabilities
get_data(fitted_models) |>
filter(classes_number == 3) |>
mutate(Class_prob = as_factor(as.integer(Class_prob))) |>
select(Class, Class_prob, Probability) |>
ggplot(aes(Class_prob, Probability, color = factor(round(Class)))) +
geom_jitter(height = 0, width = 0.1) +
scale_y_log10(breaks = 10^(seq(0, -100, by = -10)))
Decision boundary
ii <- iris |>
select(Sepal.Width, Sepal.Length, true_species = Species)
prof_3 <- ii |>
select(-true_species) |>
estimate_profiles(3)
iii <- bind_cols(
ii,
get_data(prof_3) |> select(-c(Sepal.Width, Sepal.Length))
) |>
mutate(
p = pmax(CPROB1, CPROB2, CPROB3)
)
fit <- prof_3[[1]]$model
grid <-
expand_grid(
Sepal.Width = seq(min(ii$Sepal.Width), max(ii$Sepal.Width), length.out = 100),
Sepal.Length = seq(min(ii$Sepal.Length), max(ii$Sepal.Length), length.out = 100)
)
pred <- predict(fit, newdata = grid)
gg <- grid |>
mutate(predicted_species = factor(pred$classification)) |>
mutate(
predicted_species = fct_recode(
predicted_species,
setosa = "1",
versicolor = "2",
virginica = "3"
)
)
###
ellipse_dat <- get_estimates(prof_3) |>
filter(Category == "Means") |>
select(Class, Parameter, Estimate, se) |>
pivot_wider(
names_from = Parameter,
values_from = c(Estimate, se)
)
###
ellipse_dat |> glimpse()Rows: 3
Columns: 5
$ Class <int> 1, 2, 3
$ Estimate_Sepal.Width <dbl> 3.441080, 2.687226, 3.060952
$ Estimate_Sepal.Length <dbl> 5.071534, 5.743185, 6.741152
$ se_Sepal.Width <dbl> 0.08457775, 0.09339336, 0.07713682
$ se_Sepal.Length <dbl> 0.0751666, 0.2122229, 0.2444219ggplot() +
geom_raster(
data = gg,
aes(Sepal.Length, Sepal.Width, fill = predicted_species),
alpha = 0.25
) +
geom_ellipse(
data = ellipse_dat,
aes(
x0 = Estimate_Sepal.Length,
y0 = Estimate_Sepal.Width,
a = se_Sepal.Length * 2,
b = se_Sepal.Width * 2,
angle = 0,
group = Class
),
inherit.aes = FALSE,
fill = NA
) +
geom_point(
data = iii,
aes(
Sepal.Length,
Sepal.Width,
shape = true_species,
colour = Class |> factor(),
alpha = p),
size = 3
) +
guides(color = 'none') +
theme_minimal()
Graphical methods to check model adequacy
Conceptual illustration of the 6 tidyLPA covariance models
- I am not fully convinced the following plot is correct…
Try out different variance / covariance models
Model 1: Variances equal, covariance zero
fm <- iris |>
select(-Species) |>
estimate_profiles(3, variances = 'equal', covariances = 'zero')
get_estimates(fm) |> print(n = Inf) # A tibble: 24 × 8
Category Parameter Estimate se p Class Model Classes
<chr> <chr> <dbl> <dbl> <dbl> <int> <dbl> <dbl>
1 Means Sepal.Length 5.01 0.0457 0 1 1 3
2 Means Sepal.Width 3.43 0.0581 0 1 1 3
3 Means Petal.Length 1.46 0.0280 0 1 1 3
4 Means Petal.Width 0.246 0.0155 1.03e- 56 1 1 3
5 Variances Sepal.Length 0.235 0.0329 8.88e- 13 1 1 3
6 Variances Sepal.Width 0.107 0.0143 6.04e- 14 1 1 3
7 Variances Petal.Length 0.187 0.0252 1.22e- 13 1 1 3
8 Variances Petal.Width 0.0379 0.00731 2.15e- 7 1 1 3
9 Means Sepal.Length 5.92 0.0777 0 2 1 3
10 Means Sepal.Width 2.75 0.0501 0 2 1 3
11 Means Petal.Length 4.33 0.110 0 2 1 3
12 Means Petal.Width 1.35 0.0615 1.53e-107 2 1 3
13 Variances Sepal.Length 0.235 0.0329 8.88e- 13 2 1 3
14 Variances Sepal.Width 0.107 0.0143 6.04e- 14 2 1 3
15 Variances Petal.Length 0.187 0.0252 1.22e- 13 2 1 3
16 Variances Petal.Width 0.0379 0.00731 2.15e- 7 2 1 3
17 Means Sepal.Length 6.68 0.142 0 3 1 3
18 Means Sepal.Width 3.02 0.0557 0 3 1 3
19 Means Petal.Length 5.61 0.133 0 3 1 3
20 Means Petal.Width 2.07 0.0543 3.36e-318 3 1 3
21 Variances Sepal.Length 0.235 0.0329 8.88e- 13 3 1 3
22 Variances Sepal.Width 0.107 0.0143 6.04e- 14 3 1 3
23 Variances Petal.Length 0.187 0.0252 1.22e- 13 3 1 3
24 Variances Petal.Width 0.0379 0.00731 2.15e- 7 3 1 3Model 2: Variances equal, covariance zero
fm <- iris |>
select(-Species) |>
estimate_profiles(3, variances = 'varying', covariances = 'zero')
get_estimates(fm) |>
print(n = Inf) # A tibble: 24 × 8
Category Parameter Estimate se p Class Model Classes
<chr> <chr> <dbl> <dbl> <dbl> <int> <dbl> <dbl>
1 Means Sepal.Length 5.01 0.0492 0 1 2 3
2 Means Sepal.Width 3.43 0.0552 0 1 2 3
3 Means Petal.Length 1.46 0.0260 0 1 2 3
4 Means Petal.Width 0.246 0.0159 4.62e- 54 1 2 3
5 Variances Sepal.Length 0.122 0.0211 8.43e- 9 1 2 3
6 Variances Sepal.Width 0.141 0.0325 1.49e- 5 1 2 3
7 Variances Petal.Length 0.0296 0.00591 5.73e- 7 1 2 3
8 Variances Petal.Width 0.0109 0.00258 2.52e- 5 1 2 3
9 Means Sepal.Length 5.93 0.131 0 2 2 3
10 Means Sepal.Width 2.75 0.0709 0 2 2 3
11 Means Petal.Length 4.40 0.180 1.96e-132 2 2 3
12 Means Petal.Width 1.41 0.0952 8.07e- 50 2 2 3
13 Variances Sepal.Length 0.232 0.0497 3.06e- 6 2 2 3
14 Variances Sepal.Width 0.0874 0.0161 5.22e- 8 2 2 3
15 Variances Petal.Length 0.276 0.0593 3.41e- 6 2 2 3
16 Variances Petal.Width 0.0687 0.0237 3.68e- 3 2 2 3
17 Means Sepal.Length 6.81 0.148 0 3 2 3
18 Means Sepal.Width 3.07 0.0535 0 3 2 3
19 Means Petal.Length 5.72 0.186 1.18e-208 3 2 3
20 Means Petal.Width 2.10 0.0878 7.00e-127 3 2 3
21 Variances Sepal.Length 0.286 0.0603 2.10e- 6 3 2 3
22 Variances Sepal.Width 0.0822 0.0217 1.54e- 4 3 2 3
23 Variances Petal.Length 0.250 0.0691 2.95e- 4 3 2 3
24 Variances Petal.Width 0.0605 0.0215 4.95e- 3 3 2 3Model 3: Variances equal, covariances equal
fm <- iris |>
select(-Species) |>
estimate_profiles(3, variances = 'equal', covariances = 'equal')
est <- get_estimates(fm) |>
arrange(Category)
est |> print(n = Inf)# A tibble: 30 × 8
Category Parameter Estimate se p Class Model Classes
<chr> <chr> <dbl> <dbl> <dbl> <int> <dbl> <dbl>
1 Covariances Sepal.Length.WITH… 0.0898 0.0159 1.68e- 8 1 3 3
2 Covariances Sepal.Length.WITH… 0.170 0.0329 2.50e- 7 1 3 3
3 Covariances Sepal.Length.WITH… 0.0393 0.0129 2.27e- 3 1 3 3
4 Covariances Sepal.Width.WITH.… 0.0511 0.0140 2.58e- 4 1 3 3
5 Covariances Sepal.Width.WITH.… 0.0299 0.00729 4.04e- 5 1 3 3
6 Covariances Petal.Length.WITH… 0.0420 0.0129 1.12e- 3 1 3 3
7 Means Sepal.Length 5.01 0.0433 0 1 3 3
8 Means Sepal.Width 3.43 0.0533 0 1 3 3
9 Means Petal.Length 1.46 0.0254 0 1 3 3
10 Means Petal.Width 0.246 0.0134 2.70e- 75 1 3 3
11 Means Sepal.Length 5.94 0.0809 0 2 3 3
12 Means Sepal.Width 2.76 0.0470 0 2 3 3
13 Means Petal.Length 4.26 0.105 0 2 3 3
14 Means Petal.Width 1.32 0.0424 1.09e-212 2 3 3
15 Means Sepal.Length 6.57 0.107 0 3 3 3
16 Means Sepal.Width 2.98 0.0433 0 3 3 3
17 Means Petal.Length 5.54 0.0910 0 3 3 3
18 Means Petal.Width 2.03 0.0565 1.94e-281 3 3 3
19 Variances Sepal.Length 0.264 0.0342 1.26e- 14 1 3 3
20 Variances Sepal.Width 0.112 0.0147 2.69e- 14 1 3 3
21 Variances Petal.Length 0.187 0.0390 1.67e- 6 1 3 3
22 Variances Petal.Width 0.0397 0.00651 1.13e- 9 1 3 3
23 Variances Sepal.Length 0.264 0.0342 1.26e- 14 2 3 3
24 Variances Sepal.Width 0.112 0.0147 2.69e- 14 2 3 3
25 Variances Petal.Length 0.187 0.0390 1.67e- 6 2 3 3
26 Variances Petal.Width 0.0397 0.00651 1.13e- 9 2 3 3
27 Variances Sepal.Length 0.264 0.0342 1.26e- 14 3 3 3
28 Variances Sepal.Width 0.112 0.0147 2.69e- 14 3 3 3
29 Variances Petal.Length 0.187 0.0390 1.67e- 6 3 3 3
30 Variances Petal.Width 0.0397 0.00651 1.13e- 9 3 3 3xtabs(~ Category, data = est)Category
Covariances Means Variances
6 12 12 Model 5: Variances varying, covariances varying
fm <- iris |>
select(-Species) |>
estimate_profiles(3, variances = 'varying', covariances = 'varying')
est <- get_estimates(fm) |>
arrange(Category)
est |> print(n = Inf)# A tibble: 42 × 8
Category Parameter Estimate se p Class Model Classes
<chr> <chr> <dbl> <dbl> <dbl> <int> <dbl> <dbl>
1 Covariances Sepal.Length.WITH… 0.0972 0.0215 6.08e- 6 1 6 3
2 Covariances Sepal.Length.WITH… 0.0160 0.0101 1.13e- 1 1 6 3
3 Covariances Sepal.Length.WITH… 0.0101 0.00398 1.11e- 2 1 6 3
4 Covariances Sepal.Width.WITH.… 0.0115 0.00794 1.49e- 1 1 6 3
5 Covariances Sepal.Width.WITH.… 0.00911 0.00513 7.58e- 2 1 6 3
6 Covariances Petal.Length.WITH… 0.00595 0.00272 2.87e- 2 1 6 3
7 Covariances Sepal.Length.WITH… 0.0969 0.0333 3.59e- 3 2 6 3
8 Covariances Sepal.Length.WITH… 0.185 0.0565 1.08e- 3 2 6 3
9 Covariances Sepal.Length.WITH… 0.0544 0.0182 2.85e- 3 2 6 3
10 Covariances Sepal.Width.WITH.… 0.0911 0.0280 1.13e- 3 2 6 3
11 Covariances Sepal.Width.WITH.… 0.0430 0.0103 3.09e- 5 2 6 3
12 Covariances Petal.Length.WITH… 0.0610 0.0179 6.30e- 4 2 6 3
13 Covariances Sepal.Length.WITH… 0.0922 0.0395 1.97e- 2 3 6 3
14 Covariances Sepal.Length.WITH… 0.303 0.0670 6.29e- 6 3 6 3
15 Covariances Sepal.Length.WITH… 0.0615 0.0296 3.77e- 2 3 6 3
16 Covariances Sepal.Width.WITH.… 0.0842 0.0409 3.93e- 2 3 6 3
17 Covariances Sepal.Width.WITH.… 0.0560 0.0173 1.21e- 3 3 6 3
18 Covariances Petal.Length.WITH… 0.0743 0.0439 9.07e- 2 3 6 3
19 Means Sepal.Length 5.01 0.0487 0 1 6 3
20 Means Sepal.Width 3.43 0.0526 0 1 6 3
21 Means Petal.Length 1.46 0.0224 0 1 6 3
22 Means Petal.Width 0.246 0.0159 7.53e- 54 1 6 3
23 Means Sepal.Length 5.92 0.0940 0 2 6 3
24 Means Sepal.Width 2.78 0.0591 0 2 6 3
25 Means Petal.Length 4.20 0.0816 0 2 6 3
26 Means Petal.Width 1.30 0.0298 0 2 6 3
27 Means Sepal.Length 6.54 0.0879 0 3 6 3
28 Means Sepal.Width 2.95 0.0530 0 3 6 3
29 Means Petal.Length 5.48 0.101 0 3 6 3
30 Means Petal.Width 1.98 0.0572 1.96e-263 3 6 3
31 Variances Sepal.Length 0.122 0.0210 6.57e- 9 1 6 3
32 Variances Sepal.Width 0.141 0.0326 1.54e- 5 1 6 3
33 Variances Petal.Length 0.0296 0.00676 1.22e- 5 1 6 3
34 Variances Petal.Width 0.0109 0.00291 1.84e- 4 1 6 3
35 Variances Sepal.Length 0.275 0.0577 1.86e- 6 2 6 3
36 Variances Sepal.Width 0.0926 0.0180 2.61e- 7 2 6 3
37 Variances Petal.Length 0.201 0.0594 7.29e- 4 2 6 3
38 Variances Petal.Width 0.0320 0.00648 7.77e- 7 2 6 3
39 Variances Sepal.Length 0.387 0.0737 1.53e- 7 3 6 3
40 Variances Sepal.Width 0.110 0.0276 6.48e- 5 3 6 3
41 Variances Petal.Length 0.328 0.0847 1.09e- 4 3 6 3
42 Variances Petal.Width 0.0857 0.0218 8.37e- 5 3 6 3xtabs(~ Category, data = est)Category
Covariances Means Variances
18 12 12 Different models
| Model | Variances | Covariances | Number of free parameters |
|---|---|---|---|
| 1 | Equal across classes | Zero | \((kn + n + (k-1))\) |
| 2 | Varying across classes | Zero | \((kn + kn + (k-1) = 2kn + k - 1)\) |
| 3 | Equal across classes | Equal across classes | \((kn + \frac{n(n+1)}{2} + (k-1))\) |
| 4 | Varying across classes | Equal across classes | \((kn + kn + \frac{n(n-1)}{2} + (k-1))\) |
| 5 | Equal across classes | Varying across classes | \((kn + n + k\frac{n(n-1)}{2} + (k-1))\) |
| 6 | Varying across classes | Varying across classes | \((kn + k\frac{n(n+1)}{2} + (k-1))\) |
