knitr::opts_chunk$set(
  collapse = TRUE,
  comment = "#>"
)

Importing the package

library(clusterwiselmm)
library(ggplot2)

Defining the random effets correlation matrices

lvl2_corr_mat = matrix(c(1, 0,
                         0, 1), 2, 2)
lvl3_corr_mat = matrix(c(1000, 0,
                         0, 50), 2, 2)

Generate 3 levels data with 3 groups that contains each 8 subgroups, and each subgroup contains 25 individuals

set.seed(50)
data = simulate_data(
  n1 = 25,
  n2 = 8,
  n3 = 3,
  fixed_slope = -2,
  fixed_intercept = 2,
  lvl2_corr_mat = lvl2_corr_mat,
  lvl3_corr_mat = lvl3_corr_mat,
  sigma2 = 1
)
head(data, 5)

Random effets coefficients

print(toString(c("random effet lvl3 intercepts :", unique(data$lvl3_intercepts))))
print(toString(c("random effet lvl3 slopes :", unique(data$lvl3_slopes))))
print(toString(c("random effet lvl2 intercepts :", unique(data$lvl2_intercepts))))
print(toString(c("random effet lvl2 slopes :", unique(data$lvl2_slopes))))

Visualise level 3 groups and level 3 groups for each level 2 group

print(ggplot2::ggplot(data = data) + geom_point(aes(x = x, y= y, col=lvl3)) + ggtitle("level 3 groups"))
print(ggplot2::ggplot(data = data) + geom_point(aes(x = x, y= y, col=lvl3)) + facet_wrap(data$lvl2)+ ggtitle("level 3 groups for each level 2 group")) 

Fit a clusterwise lmm with 3 classes that tries to find level 3 clustering.

formula = formula(y ~ x + (0 + x || lvl2))
final_model = clusterwiselmm(
  data = data,
  target = data$y,
  K = 3,
  formula = formula,
  nb_trials = 5
)
print(toString(c("Final mse :", final_model$cost)))

Visualize final clustering

print(ggplot2::ggplot(data = data) + geom_point(aes(x = x, y= y, col=as.factor(final_model$clusters))) + ggtitle("level 3 groups"))
print(ggplot2::ggplot(data = data) + geom_point(aes(x = x, y= y, col=as.factor(final_model$clusters))) + facet_wrap(data$lvl2)+ ggtitle("level 3 groups for each level 2 group")) 


Redha-ALLICHE/clusterwiselmm documentation built on Aug. 3, 2020, 1:10 a.m.