Every year, more than 50 million people worldwide experience a traumatic brain injury (TBI). It is estimated that about half the world’s population will have one or more TBIs in their lifetime. TBI is the leading cause of mortality in young adults and a major cause of death and disability across all ages worldwide, as recently reported in The Lancet Neurology [1]. Cognitive impairments due to TBI are the significant sources of morbidity in the affected individuals, their family members, and in the society. Disturbances in attention, memory, and executive functioning are the most common cognitive consequences of TBI at all levels of severity [2,3]. The clinical picture of TBI is characterized by a wide heterogeneity because of the nature and location of the injury [4]. Patients with TBI can show various combinations of motor, cognitive, behavioral, psychosocial, and environmental issues that have a huge impact on everyday activities [5], and these issues can greatly interfere with the effectiveness of rehabilitation interventions. It has been proposed that the efficacy of the rehabilitation would increase if programs moved from disease-centered to person-centered issues such that the rehabilitation is tailored to individual needs [6,7]. A number of studies have suggested that brain injury does not have any prototypical pattern of cognitive performance and outcome but may be best characterized by heterogeneity, both in regard to cognitive deficit and ultimate level of functioning [8]. TBI is an extremely heterogeneous disorder ranging from mild reversible conditions, often characterized as concussion, to severe massively destructive trauma, sometimes resulting in death. Saatman et al [9] highlighted the problem as follows: “The heterogeneity of TBI is considered as one of the most significant barriers to finding effective therapeutic interventions.”

TBI is a heterogeneous disease, and the mechanism/location of injury, premorbid functioning, secondary complications, and numerous other factors can influence cognitive performance [10]. As cognitive performance is a robust indicator of the current functioning and the prognostic outcome [11], it is critical to identify subgroups of patients who have distinct cognitive profiles that, in turn, can assist in treatment planning and patient care [12]. This can be empirically accomplished using cluster analysis, which is a multivariate classification technique that allows for statistical grouping of like cases into homogeneous subsets (or clusters) based on their similarity across one or more characteristics. Cluster analysis allows for the identification of homogeneous subgroups wherein cognitive heterogeneity is present based on the similarities in performance on neuropsychological tests.

Cluster analysis has been extensively applied in the study of TBI in the last 30 years [13-31]. Nevertheless, we have identified several common limitations such as the number of TBI patients that were clustered (<100 in many studies), the clustering approaches (only hierarchical clustering and k-means and not discussing other possible techniques), the specific implementation of such techniques (most of them restricted to only commercial products), as well as the lack of relation between the obtained clusters and rehabilitation tasks. The details are presented in Supplementary Material Table A1 (see Multimedia Appendix 1).

Cognitive rehabilitation has been playing an ever-increasing role in the treatment of patients with TBI who have cognitive deficits. The data gathered support the idea that improvements attributed to rehabilitation may generalize beyond task-specific skills [32]. Since the number of patients that could be eligible for this type of treatment is ever increasing, it is essential to develop new strategies that may improve access without elevating the costs to deliver such care [33]. The incorporation of computers and information technology-based systems in current clinical practice contributes to optimizing cognitive interventions, that is, their intensity, personalization, patient adherence, and quality of professional monitoring [34,35]. The types of cognitive rehabilitation programs that are the most effective in improving cognitive skills are still unclear [36]. Approaches that are designed to accommodate each individual’s cognitive strengths and weaknesses, offer instant item-specific feedback, and dynamically adapt the rehabilitation program accordingly appear to be the most effective, especially in populations with particular cognitive needs [37]. The objective of this study was to contribute to the personalization of web-based cognitive rehabilitation and to identify and characterize subgroups of patients who have distinctive profiles obtained from standard neuropsychological tests administered to patients before starting the rehabilitation.

In the following subsections, we describe the main characteristics and specific objectives of this study.

We proposed to identify and characterize cognitive profiles in a web-based cognitive rehabilitation platform by using cluster analysis with the following specific aims:

Our study consisted of patients with TBI who were admitted in the Rehabilitation Unit of the Acquired Brain Injury Department of a tertiary institution (Institut Guttmann, Spain). The period of the study was from August 2008 to July 2019.

This study was performed in accordance with the Declaration of Helsinki of the World Medical Association and approved by the ethics committee of the Clinical Research of this institution. Signed informed consent was obtained from every patient or their relatives after full explanation of the procedures. The inclusion criteria for the study were as follows: adult patients with the diagnosis of TBI and without any previous comorbidities leading to disability. Participants were excluded for illiteracy and inability to undergo formal cognitive evaluation for clinical reasons (eg, excessive sleepiness, bedridden patients, or uncontrolled sharp pain).

Initial cognition assessments used as input to cluster analysis were obtained through standardized administration of neuropsychological tests on admission; most of them were also applied to the state-of-the-art cluster analysis, as shown in Table A1 (Multimedia Appendix 1): Wisconsin Card Sorting Test, Barcelona Test, Rey Auditory Verbal Learning test, Wechsler Adult Scale III (digit span forward and backward), and Trial Making Test (Part A and Part B). All direct punctuations obtained by patients in each test were then mapped to the 0.4 ICF values. Details on the mapping of assessment instruments to ICF are presented in a previous study [44].

In this study, we took the 12 cognitive functions assessments (each one ranging from 0 to 4) as input to clustering techniques. For agglomerative hierarchical clustering, we applied the hclust function of the stats R package [45] and the AGNES function of the cluster [46] R package. For divisive hierarchical clustering, we applied the DIANA function of the cluster R package. The eclust function of the factoextra [47] R package was applied for the classic k-means implementation. The PAM function of the cluster R package was applied for PAM clustering, and similarly, the CLARA function of the same package was applied. For model-based clustering, the MClust [48] R package was applied.

We proceeded using the following steps [43]:

We then considered an alternative approach, which combined dimensionality reduction and clustering: the hierarchical clustering on principal components (HCPC) function of the FactoMineR [50] R package. It involves the following steps:

We then proposed to compare the internal validity (based only on the clustered data) of the resulting clusters based on the CVIs. These include average silhouette width [51], average Pearson gamma [52], entropy [53], Dunn index [52], and within-between cluster ratio (a higher metric of the former 3 statistics and a smaller within-between cluster ratio indicating a better fitting; eg, Clinical Cancer Research [54]). We focused especially on average silhouette width based on the conclusions in a recent review [55]. We applied the cluster.stats function of the fpc R package [56] to each of the proposed techniques for different number k of clusters, in order to obtain the CVIs. We focused on the average silhouette width by considering the following criteria [51]: 0.71-1.0, a strong structure has been found; 0.51-0.70, a reasonable structure has been found; 0.26-0.50, a weak structure has been found and could be artificial; and <0.25, no substantial structure has been found. In order to assess if the cluster holds up under plausible variations in the dataset (stability), our approach was to perform bootstrap resampling to evaluate the stability of a given cluster [57]. The cluster stability of each cluster in the original clustering is the mean value of its Jaccard coefficient over all the bootstrap iterations.

As in previous publications presented in Table A1 (Multimedia Appendix 1), in order to validate any cluster solution, it is important to compare the resulting clusters on variables that were not included in the original clustering process [25]. Various demographic variables were examined for this purpose. Regarding statistical analysis, first, analysis of the homogeneity of variance by Levene’s test and normality of distribution by the Kolmogorov-Smirnov test were conducted. Chi-square tests were conducted for most of these variables because of their ordinal nature (eg, gender), whereas analyses of variance were performed with interval variables such as age. P<.05 was considered statistically significant. We included external variables that were described in previous studies such as gender, age, age ranges, education level, FIM [58], and severity at admission measured using the GCS. In Table A2 (Multimedia Appendix 1), we have included a detailed description of FIM and GCS.

A standard cognitive rehabilitation treatment in GNPT takes 2-5 months, which is distributed in 2-5 sessions a week, and each session is composed of 4-10 cognitive training tasks. GNPT integrates a set of about 100 web-based cognitive tasks, each of which mainly addresses one of the 12 functions described above. Typically, each patient executes a different number of tasks along with treatment and in a different order. For each execution, the patient obtains an immediate result (ranging from 0 to 100, as the percentage of compliance) [59].

A final sample of 574 adult patients with TBI who performed web-based cognitive rehabilitation training in the GNPT platform were included in this study. The study was performed from August 1, 2008 to July 1, 2019. Of the 574 patients, 105 (18.3%) were women and 469 (81.7%) were men. Their distribution in the age ranges was as follows: 241 (42.0%) in the 17-30 years range, 259 (45.1%) in the 31-55 years range, and 74 (12.9%) in the >56 years range. With respect to the education level, of the 574 patients, 9 (1.6%) patients had completed primary education, 259 (45.1%) had completed secondary education, 205 (35.7%) completed tertiary education, and 101 (17.6%) completed post-tertiary education. The data of the severity of TBI at admission was available for 455 of the 574 patients (79.3%) by using the GCS, and the data were as follows: 44 (9.6%) had mild head injury, 57 (12.5%) had moderate head injury, and 354 (77.8%) had severe head injury.

In order to run the implementations of the different algorithms presented in the Methods section, input parameters were selected as mentioned in previous state-of-the-art publications presented in Table A1 (Euclidean distance and Ward criteria). As the initial preprocessing phase, we performed Spearman correlation analysis by using the corrplot [60] R package in order to identify highly correlated variables. Figure 1 shows the correlation matrix among the 12 initial variables, which is colored according to the correlation coefficient. We observed the following 3 variables with r>0.80 and P<.001: flexibility, sequencing, and working memory. Therefore, we removed them for clustering.

Table 1 shows the internal validation results for different k values and for the 6 proposed clustering techniques.

We then calculated random forest classification with 2000 trees as input parameters, and we obtained the following overall out-of-bag errors for the different k values: 1.05% (k=3), 3.83% (k=4), and 5.23% (k=5). In Supplementary Material Table A3 (Multimedia Appendix 1), we present the confusion matrix for the different k values. When calculating variable importance, there was a loss of 20% in accuracy when removing the less important variable (visual gnosis) and 25% loss when removing inhibition, as shown in Supplementary Material Figure A2 (Multimedia Appendix 1). Therefore, no variable was removed, and we did not proceed to steps 4 and 5 of the methodology.

Since FactoMineR uses a singular value decomposition algorithm, the PCA is calculated over the standardized correlation matrix, wherein a matrix of 40 uncorrelated components is obtained. Table S1 in Supplementary Material (Multimedia Appendix 1) shows the percentage of variance and the eigenvalues for the first 9 components of this matrix. The remaining components (31) correspond to a residual amount of variance. By selecting only the first 3 principal components, we reduced the dimensionality of the multivariate description so that the graphical representation and its subsequent interpretation were simplified. The first 3 principal components described 75.53% of the total variance. The first component described 55.04% of the variance, the second one described 13.42%, and the third component described 7.06%. In the case of the goodness of fit, we relied on the following metrics to verify the choice of the first 3 components: the root mean square of the residuals is 0.05 and the fit based upon off-diagonal values is 0.99.

We then ran the HCPC function with the following parameters: min=2, max=10, distance=Euclidean, criteria=Ward, and agglomerative hierarchical clustering.

When specifying min=2 and max=10 as parameters, HCPC identified the optimal k value maximizing the inertia gain. As shown in Supplementary Material Figure A3 (Multimedia Appendix 1), inertia gain dramatically decreased after the third class; therefore k=3 is the optimal partition proposed by HCPC.

When testing HCPC internal validation with the same indicators as presented in Table 1, we obtained the following CVIs: within-between ratio, 0.3706104; entropy, 0.9873104; Dunn index, 1.849996; Pearson gamma, 0.6511913; and average silhouette width, 0.515794. These CVIs clearly outperformed the CVIs presented in Table 1. For the individual approaches, the best average silhouette width was obtained by PAM for k=2 (0.395) and by k-means for k=3 (0.358). When the average silhouette width ranges from 0.26 to 0.50, the identified structure is weak and can be artificial. We focused especially on the average silhouette width, based on the conclusions in a recent CVI review [55], where 30 different indices with 720 synthetic and 20 real datasets were compared. A group of 10 indices was found to be the most recommended, with silhouette at the top in both synthetic and real datasets. Nevertheless, when considering the other CVIs in Table 1, the within-between ratio (the lower the better) HCPC was also the lowest, and Pearson gamma (the higher the better) was also higher for HCPC than any other in Table 1.

In relation to the random forest approach, when calculating variable importance, there was a loss of 20% in accuracy when removing the less important variable (visual gnosis) and 25% loss when removing inhibition. A previous study [43] removed variables leading to less than 5% loss in accuracy. In our case, no variable was removed, and therefore, we did not proceed to steps 4 and 5 of the methodology.

As presented in Table 2, the following clusters were found: cluster 1 (n=170), cluster 2 (n=311), and cluster 3 (n=93).

Table 2 shows statistically significant results for the education level of the participants as well as for all the involved cognitive functions. Analysis of cluster rationale indicated that cluster 1 is characterized by the highest level of education with almost 66% (66/170, 38.8% + 45/170, 26.5%) of its participants having tertiary or post-tertiary education. Meanwhile less than half of the participants in the other two clusters reach such educational levels: 49.2% (42/311, 13.5% + 111/311, 35.7%) of cluster 2 participants and 45.2% (14/93, 15.1% + 28/93, 30.1%) of cluster 3 participants. Furthermore, cluster 3 was characterized as complete impairment in all cognitive functions. Therefore, this cluster was characterized as very severe cognitive impairment. Meanwhile, cluster 1 presented mild impairment in working memory, visual gnosis, spatiotemporal perception, and inhibition and moderate impairment in categorization, divided attention, flexibility, planning, and sequencing. We characterized this cluster as highly educated, moderate attentional impairment, and moderate dysexecutive syndrome with mild memory impairment, and good spatiotemporal perception. Cluster 2 presented severe impairment in executive functioning (flexibility, categorization, and planning) and presented the highest degree of impairment in divided attention, as well as severe impairment in selective attention. Therefore, this cluster was characterized by severe dysexecutive syndrome with severe attentional and memory impairment and good spatiotemporal perception.

We performed twofold external validation: (1) by using demographic and clinical variables (age, gender, education level, age ranges) and then by using FIM and GCS evaluations at admission and (2) considering all cognitive tasks executed by the patients in GNPT during the period under study. We found no statistically significant differences when considering age, gender, or age ranges. The total number of available FIM assessments at admission was 439 of the original 574 participants (76.5%). Table 3 shows the number of participants, the mean, median, and IQRs for total FIM as well as the motor and cognitive subtotals for each cluster.

Regarding total FIM, patients in the 3 clusters required assistance for up to 25% of the tasks but cluster 3 was quite close to requiring assistance for 50% of the tasks. When considering the motor subtotal score with a maximum possible score of 91, patients in cluster 1 obtained 60.91, while cluster 2 obtained less than 50 and cluster 3 obtained 47.58. Regarding the cognition subtotal score (maximum score 35), cluster 1 was almost 30 while clusters 2 and 3 were close to 20.

In relation to GCS, the total number of available GCS assessments at admission was 455 (79.3%) of the original 574 participants. Table 4 shows the number of participants, mean, median, and IQRs for each cluster, and it shows the highest values for cluster 1, followed by cluster 2, and the lowest for cluster 3. Further, the IQR for cluster 3 ranged from 3 to 7, which was lower than that in clusters 1 and 2.

Regarding the second external validation, in GNPT, each task addresses a specific cognitive function. Table 5 shows the number of tasks for each function executed by cluster, with a total of 151,763 executions during the whole period under study.

Figure 2 shows the tasks result boxplots for 5 representative functions. Cluster 1 (at the left of each subplot) shows higher performance (punctuations closer to 100) than cluster 2, with cluster 3 showing lower punctuations. As shown in Table 2, for example, for the categorization function, the respective mean values for clusters 1, 2, and 3 were as follows: 2.14 (1.20), 3.72 (0.64), and 4.00 (0.00). The Figure 2 boxplots for the categorization function somehow reflect such different levels. Figure 3 represents the obtained results in every task execution for 2 functions: verbal memory and working memory. Verbal memory was the function with the largest number of executions, as shown in Table 5: 19.2% (29,148 of the total 151,763 task executions). In Figure 3, we present only cluster 1 (blue) and cluster 2 (red) in order to visually show their results, summarized weekly and plotted yearly during the whole period under study. Figure 3 shows that the working memory tasks have been integrated to the system in 2010, whereas verbal memory task executions started in 2008. For verbal tasks, cluster 1 patients outperformed cluster 2 during almost the whole period under study. Working memory tasks behave similarly, with a higher performance of cluster 2 patients.

Values between 0.60 and 0.75 indicate that the cluster is measuring a pattern in the data, but there is no high certainty about which points should be clustered together. Clusters with stability values above 0.85 can be considered highly stable (they are likely to be real clusters). The obtained values by cluster were 0.7524206, 0.6647378, and 0.9910572. Therefore, there were 2 clusters with stability >0.75. As a rule of thumb, clusters with a stability value less than 0.60 should be considered unstable, which is not our case. Therefore, meaningful valid clusters as the ones identified in our study should not disappear if the data set is changed in a nonessential way. Nevertheless, it could also be of interest whether clusters remain stable under the addition of outliers; such cases should be individually considered by clinicians (eg, in case of the lowest GCS assessment values).

In this study, we proposed the application of cluster analysis to a chronic health condition in a GNU framework by using a set of publicly available R libraries (R-3.5.1) in the context of a web-based cognitive platform. We proposed 6 specific clustering techniques (ie, PAM, CLARA, AGNES, DIANA, k-means, and MClust) and 2 combined approaches (HCPC=PCA+AGNES and random forest+PAM) and evaluated them by using state-of-the-art CVIs. It is straightforward to apply both the individual techniques and the combined approaches to other acquired brain injury populations in the same web-based platform (GNPT) or in others. For example, in the Multimedia Appendix 1, we present an initial correlation analysis for patients who had an ischemic stroke that we will address in future work. We obtained the best CVIs with the combined HCPC=PCA+AGNES hierarchical clustering, with average silhouette over 52%; therefore, a reasonable structure has been found. We performed stability analysis, and clusters 1 and 2 were rated as 0.60 and 0.75, indicating that the clusters are measuring a pattern, and cluster 3 was rated as highly stable. We identified 3 clearly different profiles. Cluster 1 was characterized as highly educated, moderately distracted, with dysexecutive syndrome and good working memory. Cluster 2 was characterized as severe dysexecutive syndrome and severely distracted. Cluster 3 identified a group of patients with severe symptoms in all the involved functions. External validity in functional independence confirms this characterization by means of severity using GCS and functionality in the activities of daily living, especially when considering the motor FIM subtotal. When considering the performance in the cognitive tasks executed during the whole period, task results confirmed the identified profiles, with cluster 1 visual representation showing higher values during the whole period than cluster 2. Similar results were obtained when visualizing cluster 3.

The actual GNPT implementation integrates an automatic therapy planning functionality, the intelligent therapy assistant (ITA) [61]. The ITA provides therapists with a recommended schedule of cognitive tasks to be executed by each patient during a given period of time. The recommendations provided by the ITA can always be manually modified by therapists according to their own clinical criteria. The ITA takes a predefined set of patient’s cognitive profiles as the starting point, which have been obtained using the baseline cognitive evaluation (mapped to ICF as described in the Methods section) as input to CA. When a new patient starts cognitive training in GNPT, the ITA dynamically assigns the patient to the appropriate cluster. The ITA then schedules different cognitive tasks during a user-defined rehabilitation period to the new patient, according to several criteria (eg, usage score, improvement score, clinical score) as described in previous studies. Therefore, the first clinical implication involves the ITA starting point to configure patients’ treatments. During therapy, when the patient executes a task (and obtains the result ranging from 0 to 100), GNPT automatically generates another version of the task with a higher or lower difficulty level—increasing the difficulty if the result was “too high” or decreasing the difficulty if the result was “too low” [62]. A second clinical implication involves linking cognitive profiles with performance in task execution. As shown in Figure 3, this allows therapists to identify patterns in performance, for example, results seem to be too close to 50 for cluster 2 in verbal memory tasks during the 2013-2016 period. The current clinical working hypothesis in relation to patient’s performance in GNPT tasks is that the optimal range of results is 65-85 [63]. Therefore, Figure 3 (top, verbal memory) suggests that difficulty levels in such tasks might have been too high for patients in cluster 2 during the 2013-2016 period. A more appropriate approach regarding the optimal range of results could be to consider such ranges to vary in relation to clusters. Therefore, a patient in cluster 1 would have a different optimal range than a patient in cluster 2. The next step is to consider the optimal range of the results depending on the cognitive profiles identified by cluster analysis (instead of considering a fixed optimal range as it is now). Future work should also include comparing ITA current cluster analysis results [61] with clusters 1, 2, and 3 obtained in this work for patients with TBI. The integration of cluster analysis as the initial phase of an ITA process also allows for a straightforward extension of a similar approach to other medical conditions, for example, patients who had a stroke, as we present in the Supplementary Material (Multimedia Appendix 1).

First, we conducted a single-center study; an advantage of this is that data were obtained and included by clinicians trained in neurological rehabilitation, and all patients were managed under the same TBI rehabilitation protocols. The GNPT platform is already integrated into the clinical practice of several acquired brain injury centers; nevertheless, their patients were not included in this analysis. A multicenter TBI study may include an initial preprocessing phase, wherein patients are grouped according to their initial GCS severity in order to avoid additional heterogeneity. Thereafter, cluster analysis techniques, as those proposed in this study, may be applied within such groups. External validation assessments, common to all participating centers, is also an important aspect to be addressed in this future multicenter study. Second, the health area studied belongs mainly to the urban population, with a small rural population or populations from other regions.

Third, our analysis lacked computerized tomography or magnetic resonance imaging examinations that describe the presence of contusion, hematoma, hemorrhage, ischemia, or other signs of parenchymal lesion on frontal, temporal, parietal, occipital, and cerebellar lobes or diffuse axonal injury. Fourth, our sample did not include any patient with missing data. All data used as input to cluster analysis are complete. Although there are several R packages addressing the subject (MICE, MissForest, HMISC), we decided to address the problem of missing data in a separate future analysis in order to consider not only the possible imputation strategies but also the reasons for missing data and include such reasons when characterizing the clusters. Fifth, our analysis did not include indicators of mental health or other comorbidities. Persons who experience TBI may have 1 or more preexisting medical comorbidities at the time of injury (eg, alcohol use and depression). Other medical conditions may occur simultaneously with TBI, such as orthopedic trauma, or these conditions may develop afterward as a direct consequence of the TBI such as epilepsy. Still, other medical comorbidities may begin months or years following injury in comparison to uninjured control groups. Studies have suggested that individuals with TBI have more than twice the rates of pain, growth hormone deficiency, insomnia, fatigue, new-onset stroke, urinary incontinence, and epilepsy [64]. Therefore, we aim to include comorbidity analysis in future research studies.

We have worked with public GNU libraries, as opposed to the state-of-the-art publications presented in Table A1, wherein most techniques were implemented using commercial packages [15-18,20-23,25-27,29-31]. Previous research presented in Table A1 applied clustering techniques in a batch mode as desktop applications. In our case, the work was integrated in the context of a web-based cognitive training platform. Our baseline assessment consisted of 12 cognitive functions, thereby allowing for a comprehensive description of the patient’s profiles, involving cognitive aspects addressed by such different functions, ranging from visual attention to gnosis. Meanwhile, previous clustering research presented in Table A1 addresses specific functions—only one of them in most cases: memory [14,16,18-21,24-26,30], executive functions [17,21,31], or attention [22]. We have proposed different clustering techniques and applied state-of-the-art CVIs to all of them. We have taken advantage of the web-based platform by increasing the number of participants, whereas in only 3 of the 20 studies in Table A1, n is larger than 300 [20,25,30]. We have included the whole set of cognitive tasks performed by all participants as part of the external validation during the whole period under study (more than 150,000 task executions). We have visually mapped such executions to the obtained clusters along time. To the best of our knowledge, the linking of specific rehabilitation tasks to the obtained clusters has not been yet performed in the state-of-the-art publications presented in Table A1.

Cluster analysis in web-based cognitive rehabilitation treatments allows for identifying and characterizing strong patterns of response to neuropsychological tests, externally validating the obtained clusters by using important aspects of TBI rehabilitation such as severity or functional independence in activities of daily life, tailoring cognitive web-based tasks available in the web platform to the identified profiles by providing clinicians a tool for treatment personalization, which were not addressed in previous traditional cluster analyses, and straightforward extension of a similar approach to patients with other medical conditions, for example, for patients who have had a stroke.