ds <- readr::read_csv("2 Longterm/assigndata.csv", na = c("", "NA")) |> na.omit()

ds <- ds |> dplyr::mutate(mrs_1 = mrs_1 > 2)


library(tidymodels)
library(tidyverse)

ds |>
  na.omit() |>
  kmeans(centers = 3)


ds_num <- ds |> mutate(
  across(where(is.double), ~ scale(.x)),
  across(is.character, ~ as.numeric(factor(.x)))
)

ds_num |> kmeans(centers = 3)

set.seed(321)
kclusts <-
  tibble(k = 1:9) %>%
  mutate(
    kclust = map(k, ~ kmeans(ds_num, .x)),
    tidied = map(kclust, tidy),
    glanced = map(kclust, glance),
    augmented = map(kclust, augment, ds_num)
  )

kclusts

clusters <-
  kclusts %>%
  unnest(cols = c(tidied))

assignments <-
  kclusts %>%
  unnest(cols = c(augmented))

clusterings <-
  kclusts %>%
  unnest(cols = c(glanced))

ggplot(assignments, aes(x = pase_0, y = age)) +
  geom_point(aes(color = .cluster), alpha = 0.8) +
  facet_wrap(~k)

ggplot(clusterings, aes(k, tot.withinss)) +
  geom_line() +
  geom_point()

PCA

pairs(assignments[-c(2:4)],
  gap = 0,
  bg = c("red", "yellow", "blue")[assignments$k],
  pch = 21
)


pc.out <- prcomp(ds_num, center = TRUE, scale = TRUE)

pc.sum <- summary(pc.out)

pscr <- tibble(
  x = 1:dim(pc.sum$importance)[2],
  Proportion = pc.sum$importance[2, ],
  Cumulative = pc.sum$importance[3, ]
) %>%
  pivot_longer(cols = -x) %>%
  ggplot(aes(x = x, y = value, color = name)) +
  geom_line() +
  geom_point() +
  ylim(0, 1) +
  labs(
    title = "Scree plot",
    color = "Variance"
  ) +
  ylab("Variance") +
  xlab("Principal components")


###
###

install.packages("VarSelLCM")
library(VarSelLCM)

# Please indicate the number of cores you wan to use for parallelization
nb.CPU <- 4
# clustering without variable selection (about than 10/20 sec on 4 CPU)
res_without <- ds |> select(!mrs_1) |> 
  mutate(across(is.character, ~ factor(.x))) |> as.data.frame() |> 
  VarSelCluster(
    gvals = 1:5,
    crit.varsel = "BIC",
    vbleSelec = FALSE,
    nbcores = nb.CPU
  )

summary(res_without)

plot(res_without)

plot(x=res_without, y="mdi_1")

plot(x=res_without, y="sex")

print(res_without)

coef(res_without)

VarSelShiny(res_without)