2. Get communities of patients with unique symptom
severity/frequency patterns
To run second-order unsupervised clustering, we use the
get_communities function, which has 6 parameters:
data: Dataframe of symptom severities (formatted as required by
PRONA)
nrows: Number of rows in the dataframe (Defaults to the number of
rows in data, you should not change this parameter)
ncols: Number of columns in the dataframe (Defaults to the number
of columns in data, you should not change this parameter)
detectAlgo: The community detection algorithm to use. Options
include FG (fast greedy), IM (infomap), LP (label propagation), LE
(leading eigen), LV (louvain), or WT (walktrap). (Default: WT)
simil_measure: The similarity measure to use, either ‘Euclidean’
or ‘ARI’ (Adjusted Rand Index) (Default: ARI)
simplify_graphs: Boolean that indicates whether to simplify the
graph by removing multi-edges and loops (Default: TRUE)
It returns a new dataframe that is a copy of the original symptom
severities dataframe, but with a new column representing which community
each patient belongs to.
community_df <- get_communities(reduced_df)
Lets see how many patients are in each community
table(community_df$community)
##
## 1 2 3
## 136 218 5
Next, lets plot a heatmap of symptom severities for each community
using the plot_community_heatmap function, which has four
parameters:
df: The output of the get_communities function
cluster_rows: Boolean representing whether or not to cluster the
rows (symptoms) via heirarchical clustering (Default: TRUE)
network: A network object output by construct_ggm. If a network
is passed into this function, the heatmap will automatically order rows
(symptoms) by the symptom clusters identified in the network. This
overrides the cluster_rows parameter. (Default: NULL)
row_label_size: The font size of the row labels (Default:
10)
Here, “PC” stands for “Patient Community”:
plot_community_heatmap(community_df, row_label_size = 9)
Here’s how to order rows based on the symptom clusters identified
during network construction:
plot_community_heatmap(community_df, network = ptsd_network, row_label_size = 9)
We can also visualize the symptom severities of each community using
line plots. The community_line_plot function has two parameters:
data: output of get_communities
communities: ector specifying which communities to plot. If you
include 0 in the vector, it will plot the severity for all patients
combined. (Default: c(0))
community_line_plot(community_df, c(0,1,2))
We can also get a summary of each symptom in each community using the
get_community_summary function, which returns a dataframe with the mean,
median, upper bound, and lower bound of each symptom in each
community.
get_community_summary(community_df)
## # A tibble: 45 × 6
## # Groups: community [3]
## community variable mean median upper lower
## <dbl> <fct> <dbl> <dbl> <dbl> <dbl>
## 1 1 AVOID.REMINDERS.OF.THE.TRAUMA 0.551 0 0.716 0.387
## 2 1 BAD.DREAMS.ABOUT.THE.TRAUMA 0.309 0 0.431 0.187
## 3 1 DISTANT.OR.CUT.OFF.FROM.PEOPLE 0.838 0 1.02 0.657
## 4 1 FEELING.EMOTIONALLY.NUMB 0.640 0 0.802 0.478
## 5 1 FEELING.IRRITABLE 0.853 1 1.01 0.699
## 6 1 FEELING.PLANS.WONT.COME.TRUE 0.625 0 0.789 0.461
## 7 1 HAVING.TROUBLE.CONCENTRATING 1.37 1 1.55 1.18
## 8 1 HAVING.TROUBLE.SLEEPING 1.60 2 1.80 1.40
## 9 1 LESS.INTEREST.IN.ACTIVITIES 0.309 0 0.445 0.172
## 10 1 NOT.ABLE.TO.REMEMBER 0.397 0 0.535 0.259
## # ℹ 35 more rows
3. Construct networks for individual communities
We can construct GGM networks and conduct statistical assessment for
individual communities using the same pipeline as described in the
introduction and statistical assessment vignettes.
# Get symptom data for individual communities
library(dplyr)
## Warning: package 'dplyr' was built under R version 4.2.3
##
## Attaching package: 'dplyr'
## The following objects are masked from 'package:stats':
##
## filter, lag
## The following objects are masked from 'package:base':
##
## intersect, setdiff, setequal, union
comm1_df <- community_df[community_df$community==1,] %>% select(-community)
comm2_df <- community_df[community_df$community==2,] %>% select(-community)
comm1_ega <- construct_ggm(comm1_df, normal = FALSE)
## Conducting nonparanormal (npn) transformation via skeptic....done.
## Warning: The network input contains unconnected nodes:
## AVOID.REMINDERS.OF.THE.TRAUMA, BAD.DREAMS.ABOUT.THE.TRAUMA, FEELING.PLANS.WONT.COME.TRUE, HAVING.TROUBLE.CONCENTRATING, HAVING.TROUBLE.SLEEPING, NOT.THINKING.ABOUT.TRAUMA, PHYSICAL.REACTIONS, RELIVING.THE.TRAUMA, JUMPY.STARTLED.OVER.ALERT
## Warning: Some variables did not belong to a dimension: AVOID.REMINDERS.OF.THE.TRAUMA, BAD.DREAMS.ABOUT.THE.TRAUMA, FEELING.PLANS.WONT.COME.TRUE, HAVING.TROUBLE.CONCENTRATING, HAVING.TROUBLE.SLEEPING, NOT.THINKING.ABOUT.TRAUMA, PHYSICAL.REACTIONS, RELIVING.THE.TRAUMA, JUMPY.STARTLED.OVER.ALERT
##
## Use caution: These variables have been removed from the TEFI calculation
## Model: GLASSO (EBIC with gamma = 0)
## Correlations: auto
## Lambda: 0.301079593387959 (n = 100, ratio = 0.1)
##
## Number of nodes: 15
## Number of edges: 4
## Edge density: 0.038
##
## Non-zero edge weights:
## M SD Min Max
## 0.087 0.063 0.038 0.178
##
## ----
##
## Algorithm: Walktrap
##
## Number of communities: 2
##
## AVOID.REMINDERS.OF.THE.TRAUMA BAD.DREAMS.ABOUT.THE.TRAUMA
## NA NA
## DISTANT.OR.CUT.OFF.FROM.PEOPLE FEELING.EMOTIONALLY.NUMB
## 1 1
## FEELING.IRRITABLE FEELING.PLANS.WONT.COME.TRUE
## 1 NA
## HAVING.TROUBLE.CONCENTRATING HAVING.TROUBLE.SLEEPING
## NA NA
## LESS.INTEREST.IN.ACTIVITIES NOT.ABLE.TO.REMEMBER
## 1 2
## NOT.THINKING.ABOUT.TRAUMA PHYSICAL.REACTIONS
## NA NA
## RELIVING.THE.TRAUMA JUMPY.STARTLED.OVER.ALERT
## NA NA
## UPSETTING.REMINDERS
## 2
##
## ----
##
## Unidimensional Method: Louvain
## Unidimensional: No
##
## ----
##
## TEFI: -2.513
plot_ggm(comm1_ega, label.size = 2.5)
comm2_ega <- construct_ggm(comm2_df, normal = FALSE)
## Conducting nonparanormal (npn) transformation via skeptic....done.
## Model: GLASSO (EBIC with gamma = 0.25)
## Correlations: auto
## Lambda: 0.129502975706062 (n = 100, ratio = 0.1)
##
## Number of nodes: 15
## Number of edges: 47
## Edge density: 0.448
##
## Non-zero edge weights:
## M SD Min Max
## 0.081 0.068 0.001 0.343
##
## ----
##
## Algorithm: Walktrap
##
## Number of communities: 2
##
## AVOID.REMINDERS.OF.THE.TRAUMA BAD.DREAMS.ABOUT.THE.TRAUMA
## 1 1
## DISTANT.OR.CUT.OFF.FROM.PEOPLE FEELING.EMOTIONALLY.NUMB
## 2 2
## FEELING.IRRITABLE FEELING.PLANS.WONT.COME.TRUE
## 2 2
## HAVING.TROUBLE.CONCENTRATING HAVING.TROUBLE.SLEEPING
## 2 2
## LESS.INTEREST.IN.ACTIVITIES NOT.ABLE.TO.REMEMBER
## 2 2
## NOT.THINKING.ABOUT.TRAUMA PHYSICAL.REACTIONS
## 1 1
## RELIVING.THE.TRAUMA JUMPY.STARTLED.OVER.ALERT
## 1 1
## UPSETTING.REMINDERS
## 1
##
## ----
##
## Unidimensional Method: Louvain
## Unidimensional: No
##
## ----
##
## TEFI: -7.476
plot_ggm(comm2_ega, label.size = 2.5)
Because of its relatively low number of patients and low overall
symptom severity, the network for PC1 is sparse.
5. Assess the stability of patient communities using sub-sampling
analysis
To assess the stability of patient communities (aka how reliable they
are), we can use sub-sampling analysis. This technique takes random
subsamples of patients (default is subsamples of 100%, 99%, 95%, 90%,
80%, 70%, 60%, 50%, 40%, 30%, 20%, and 10% of the original dataset) and
reruns concordance network-based patient clustering on each of the
subsamples.
To perform this, we can use the cluster_subsamples function, which
takes 5 parameters:
data: Original dataframe of symptom severities (formatted as
required by PRONA)
detectAlgo: The community detection algorithm to use. It should
be the same as what you used when you ran get_communities. (Default:
WT)
simil_measure: The similarity measure to use. It should be the
same as what you used when you ran get_communities. (Default:
ARI)
simplify_graphs: Boolean that indicates whether to simplify the
graph by removing multi-edges and loops. It should be the same as what
you used when you ran get_communities. (Default: TRUE)
sampling_rates: A vector containing the percentages to subset the
data by. Default is c(1, 0.99, 0.95, 0.9, 0.8, 0.7, 0.6, 0.5, 0.4, 0.3,
0.2, 0.1). If you change this, the vector must start with 1 (meaning the
first subsample is 100% of the dataset) for downstream analysis to
work.
subsample_communities_df <- cluster_subsamples(reduced_df)
The output is a dataframe that contains the community designations
for each patient in each of the subsampling runs:
head(subsample_communities_df, 5)
## ID community_1 community_0.99 community_0.95 community_0.9
## 1 15_000651 2 2 1 1
## 2 15_000786 1 1 2 2
## 3 15_000969 2 2 3 3
## 4 15_000989 2 2 3 3
## 5 15_001332 2 3 1 1
## community_0.8 community_0.7 community_0.6 community_0.5 community_0.4
## 1 4 2 1 1 3
## 2 4 2 3 3 NA
## 3 4 2 1 1 4
## 4 1 3 4 4 4
## 5 3 2 1 1 3
## community_0.3 community_0.2 community_0.1
## 1 3 NA NA
## 2 NA NA NA
## 3 2 1 NA
## 4 2 1 NA
## 5 3 3 1
To see the distribution of patients across different communities in
each subsampling run, we can use the following code. It is important to
note that cluster labels do not necessarily overlap across different
runs (for example, community #2 is
“subsample_communities_df$community_1” does not necessarily correspond
to community #2 in “subsample_communities_df$community_0.99”). Thus, you
should compare the overall distribution of patients across different
communities in different runs, but not the number of patients under each
specific label:
table(subsample_communities_df$community_1)
##
## 1 2 3
## 136 218 5
table(subsample_communities_df$community_0.99)
##
## 1 2 3 4
## 99 124 128 4
table(subsample_communities_df$community_0.95)
##
## 1 2 3 4 5
## 159 79 70 30 3
table(subsample_communities_df$community_0.9)
##
## 1 2 3 4 5
## 155 80 62 24 2
table(subsample_communities_df$community_0.8)
##
## 1 2 3 4 5
## 78 33 125 49 2
table(subsample_communities_df$community_0.7)
##
## 1 2 3 4 5
## 49 71 56 73 2
table(subsample_communities_df$community_0.6)
##
## 1 2 3 4 5
## 110 32 36 35 2
table(subsample_communities_df$community_0.5)
##
## 1 2 3 4 5
## 92 29 26 30 2
table(subsample_communities_df$community_0.4)
##
## 1 2 3 4 5
## 13 33 70 25 2
table(subsample_communities_df$community_0.3)
##
## 1 2 3 4
## 9 36 60 2
table(subsample_communities_df$community_0.2)
##
## 1 2 3 4 5 6
## 17 11 33 7 2 1
table(subsample_communities_df$community_0.1)
##
## 1 2 3 4 5
## 12 4 2 15 2
Next, to compute the stability of each of the original clusters
across all the subsampling runs, we can use the
compute_cluster_stabilities and plot_cluster_stabilities functions.
These functions define “stability” as how often data points that belong
to the same cluster in the original data are still clustered together in
the subsamples (this is also known as Prediction Strength). It is
calculated by taking all the patients that belong to one of the
individual clusters, finding which communities those patients have been
grouped into in the new clustering iteration, and dividing the size of
the largest community in the new iteration by the size of the original
community.
For example, say clustering on a dataset of 1000 patients yields
three communities, one of 500 patients (community 1), one of 400
(community 2), and one of 100 (community 3). We rerun clustering on a
subsample of the data, and out of the 500 patients that originally
belonged to community 1, 400 of them still belong to the same community,
50 belong to a different community, and 50 were randomly removed during
subsampling. The “stability” of the original community 1 for this
subsampling iteration is 400/500 = 0.8 (aka 80%).
To perform this analysis on each of the original clusters for each of
the subsampling iterations, use the following code:
stabilities_df <- compute_cluster_stabilities(subsample_communities_df)
## Warning in max(freq_clusters[names(freq_clusters) != ""]): no non-missing
## arguments to max; returning -Inf
## SamplingRate Cluster 1 Stability Cluster 2 Stability Cluster 3 Stability
## 1 1 1.00000000 1.00000000 1.0
## 2 0.99 0.58715596 0.72794118 0.8
## 3 0.95 0.72935780 0.52941176 0.6
## 4 0.9 0.68348624 0.58823529 0.4
## 5 0.8 0.57339450 0.30147059 0.4
## 6 0.7 0.33486239 0.36029412 0.4
## 7 0.6 0.49082569 0.25735294 0.4
## 8 0.5 0.39908257 0.21323529 0.4
## 9 0.4 0.32110092 0.23529412 0.4
## 10 0.3 0.27522936 0.20588235 0.4
## 11 0.2 0.15137615 0.07352941 0.2
## 12 0.1 0.06880734 0.02941176 -Inf
To plot these results, use the plot_cluster_stabilities function:
plot_cluster_stabilities(stabilities_df)
## `geom_smooth()` using formula = 'y ~ x'
Finally, we can use the Adjusted Rand Index (ARI) to assess the
global stability of the community groupings across the different
subsampling iterations. This helps us understand if the overall
structure of clusters is consistent across different subsamples. We can
calculate the ARI between the original clusters (represented by
subsample_communities_df$community_1) and each of the subsamples using
the following code:
## Warning: package 'mclust' was built under R version 4.2.3
## Package 'mclust' version 6.1.1
## Type 'citation("mclust")' for citing this R package in publications.
adjustedRandIndex(subsample_communities_df$community_1, subsample_communities_df$community_0.99)
## [1] 0.4325011
adjustedRandIndex(subsample_communities_df$community_1, subsample_communities_df$community_0.95)
## [1] 0.4971679
adjustedRandIndex(subsample_communities_df$community_1, subsample_communities_df$community_0.9)
## [1] 0.5135012
adjustedRandIndex(subsample_communities_df$community_1, subsample_communities_df$community_0.8)
## [1] 0.3745413
adjustedRandIndex(subsample_communities_df$community_1, subsample_communities_df$community_0.7)
## [1] 0.2428197
adjustedRandIndex(subsample_communities_df$community_1, subsample_communities_df$community_0.6)
## [1] 0.5322329
adjustedRandIndex(subsample_communities_df$community_1, subsample_communities_df$community_0.5)
## [1] 0.4751524
adjustedRandIndex(subsample_communities_df$community_1, subsample_communities_df$community_0.4)
## [1] 0.5615193
adjustedRandIndex(subsample_communities_df$community_1, subsample_communities_df$community_0.3)
## [1] 0.664666
adjustedRandIndex(subsample_communities_df$community_1, subsample_communities_df$community_0.2)
## [1] 0.4221146
adjustedRandIndex(subsample_communities_df$community_1, subsample_communities_df$community_0.1)
## [1] 0.3372875
Although ARI is a helpful metric, it can be disproportionately skewed
if there are specific clusters that are particularly unstable, even if
the other clusters are generally stable. For this reason, we generally
prefer the above analysis of individual cluster stability over ARI
calculation because it gives you more specific information on which
clusters are reliable and which are not.
