{'fraction_non_alphanumeric': 0.035215850748427216, 'fraction_numerical': 0.01590859787403283, 'mean_word_length': 4.582038345105953, 'pattern_counts': {'":': 0, '<': 3, '<?xml version=': 0, '>': 1, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 11, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Defying the Circadian Rhythm: Clustering Participant Telemetry in the UK Biobank Data

The UK Biobank dataset follows over 500,000 volunteers and contains a diverse set of information related to societal outcomes. Among this vast collection, a large quantity of telemetry collected from wrist-worn accelerometers provides a snapshot of participant activity. Using this data, a population of shift workers, subjected to disrupted circadian rhythms, is analysed using a mixture model-based approach to yield protective effects from physical activity on survival outcomes. In this paper, we develop a scalable, standardized, and unique methodology that efficiently clusters a vast quantity of participant telemetry. By building upon the work of Doherty et al. (2017), we introduce a standardized, low-dimensional feature for clustering purposes. Participants are clustered using a matrix variate mixture model-based approach. Once clustered, survival analysis is performed to demonstrate distinct lifetime outcomes for individuals within each cluster. In summary, we process, cluster, and analyse a subset of UK Biobank participants to show the protective effects from physical activity on circadian disrupted individuals.


Introduction
The UK BioBank dataset is a large collection of participant data with a broad selection of information that is relevant for statistical investigation (Biobank, 2014). The use of accelerometer data for predicting lifetime outcomes has become a key research focus in the last decade. There is an overwhelming consensus that physical inactivity has detrimental effects on lifetime outcomes (Lee et al., 2012). Those include inter alia; breast cancer, coronary heart disease, and type 2 diabetes (La Vecchia et al., 2012). As a result, physical inactivity directly affects life expectancy (Cunningham et al., 2020); leading to early death and onset of disease.
Additional factors, often not considered, are circadian rhythm defying lifestyle choices.
Many workers who undertake their job in some form of shift work or extended hours have implications on respective sleeping patterns. Studies have shown that these lifestyles can have negative health effects (Harrington, 2001;Dembe et al., 2005;Fan et al., 2020). However, Roveda et al. (2017) showed a protective effect where physical activity may be beneficial against the detrimental health effects typically associated with sleep disruption. Such studies are not uncommon as found in Yamanaka et al. (2006), Montaruli et al. (2017) and Nohara et al. (2015).
As highlighted by Doherty et al. (2017), a major short coming of studies in this area have been based on self-reported evidence. As a consequence, such studies are subject to inaccuracies as self-reported evidence is heavily biased; particularly when concerning physical activity as outlined in Brenner and DeLamater (2016). Again, Doherty et al. (2017) emphasizes the growing use of objective measures for physical activity in the form of telemetry collected from wrist-worn accelerometers. However, the use of accelerometer telemetry beckons another issue. The sheer overwhelming quantity and vastness of telemetric data.
At the time of analysis, the size of accelerometer data totalled to 26 terabytes of memory.
Each participant generates approximately 217 megabytes of data or, the equivalent of two and a half hours of social media consumption (Cisco, 2016). For specifics regarding wear time compliance, Doherty et al. (2017) provides a thorough analysis across age and sex.
From a computational perspective, the size of data poses massive challenges to process, manage, and analyse lifetime outcomes. Our paper builds upon Doherty et al. (2017) and Willetts et al. (2018), by developing a standardized, low-dimensional feature that captures physical activity behaviour from accelerometer telemetry. We use a mixture model-based approach to cluster and identify participant sub-groups. Finally, we perform survival analysis to demonstrate how physical inactivity directly affects lifetime outcomes.

Methodology
This section is broken down in chronological order on how the analysis is performed. Firstly, a cohort is selected where lifetime outcomes are well known from an occupational health perspective. Second, detail is provided on the use of the activity extraction tool outlined in Doherty et al. (2017). Thirdly, a low-dimensional feature that standardizes and captures physical activity is developed. Furthermore, a mixture model-based approach is used to cluster said feature. Finally, we summarize two popularized methods used for survival analysis.

Study Population Selection
There is a well known phenomenon that life expectancy is lower among workers that defy their natural circadian rhythm through lifestyle factors such as shift work or extended hours (Gu et al., 2015). Consequently this results in unusual sleeping patterns. Several studies from an occupational health perspective, highlight cases of early onset non-communicable diseases, thus lowering their life expectancy (Ijaz et al., 2013). Within the UK Biobank dataset, a large subset of participants reported their employment history regarding shift work. We investigate this phenomenon within the aformentioned population to achieve a consesus with previous literature.
At the time of analysis the UK Biobank dataset contained 502, 536 individuals. First, a subset of participants is selected based on their response to a questionnaire regarding employment history. The subset is split between two groups, regular shift, and late shift.  Tables 1 and 2 show a fairly even split between regular and late shift workers. There is a clear majority of males among late shift workers. As a consequence, known survival outcomes based on sex differences may play a role in our analysis. Typically, literature has shown to favour females over males (Topjian et al., 2010;Antero et al., 2020;Di Lego et al., 2020). These differences are taken into account during our survival analysis. However, as a simplification, we proceed without separating based on sex during our clustering procedures.

Data Processing
Recent developments on extracting behavioural patterns from the UK BioBank dataset has yielded an opportunity for prediction of lifetime outcomes. Doherty et al. (2017) (Gramacki, 2019). This type of feature has been well studied from an engineering and actuarial perspective (Kamble et al., 2009;Wüthrich, 2017). Our force map is constructed as follows. Let f t denote force measurement at time interval t.
Force is measured in milli-gravity units (mg), where it is averaged over a 5 second time interval. Let f t denote the numerical first derivative calculated as Note the denominator, as each interval differs by 5 seconds. f t is considered to be average change in force over some time interval t. Wüthrich (2017) elaborates that this calculation is determined by the "average" change over the time interval, and not reminiscent of the true instantaneous measurement. To construct a force map, let R be some rectangle where is considered to be values of a two dimensional coordinate system. Next, construct a partitioning of R into M equally sized rectangles as Finally, consider a probability distribution of F ∈ P(R) having probability weight The selection of F can allow for smooth approximations to the natural behaviour. As Pocuca (2019) demonstrates, the selection of a two dimensional Gaussian kernel is superior to other methods used in Wüthrich (2017) and Gao et al. (2019).
, corresponding to the bivariate collection of all force measurements for a particular participant i. Next, consider the use of KDE for the problem of estimating an unknown joint probability density p(f ) on the space R ⊆ R 2 . Let K be a bivariate function defined on the space R 2 . Furthermore, let H be a constant, positive definite, symmetric matrix defined as the smoothness parameter. The KDE of p(f ) is written aŝ Due to popularity and radial symmetrical properties, the standard bivariate Gaussian distribution is selected as the kernel function K. As a result, K is formulated as This selection of K allows the kernel estimator to be the weighted sum of normal densities centred at force measurements f it . In summary, the true density p(f ) is approximated using the kernel estimatorp Here, φ is the Gaussian density with covariance matrix H. The selection of the smoothness matrix H affects both the shape and orientation of the kernels on the two dimensional space. The bandwidth is a key parameter for optimizing performance of KDE. For the purposes of simplicity, the normal scale selector (Chacón et al., 2011) is selected as whereΣ is considered to be the sample covariance matrix. The use of this method for estimating the probability density function is superior to that of a histogram or an empirical kernel. Furthermore, the KDE method captures the continuous nature of data even when latency is large. For specifics, see Pocuca (2019). For implementation, see the MASS package regarding the use of two dimensional KDE (Venables and Ripley, 2002).
Finally, construct the matrix variate object X i whose matrix entries directly correspond to the probability weights x im over region R m . Probability weights are calculated as This matrix variate object is defined as the force map for participant i. Note that this extension differs from the original construction by Wüthrich (2017) in one key distinction.
That distinction being the use of KDE to estimate the joint density on R. The choice of a smooth kernel allows the matrix X i to have non-zero entries which more accurately resemble a participant's behaviour. to have another localized area of force measurements centred at 175 units. This is again reflected in the force map in Figure 1b. In summary, for this participant, we capture their walking behaviour characterized by this force map. Repeating the same process for each behaviour, we visualize force maps as shown in Figure 2. We

Model-based Clustering
The framework of model-based clustering utilizes the underlying assumption that a finite mixture model embodies the representation of heterogeneous data. Consider a random variable X from a G-component finite mixture model with probability density function of the form where ϑ = {π 1 , . . . , π G , θ 1 , . . . , θ G }, x is a realization of X , π g is a mixing proportion where π g > 0, G g=1 π g = 1, and p g is a probability density function parametrized by θ g . The distribution for each cluster is usually taken to be the same density in (1) and is simplified as p g (x|θ g ) = p(x|θ g ) ∀g.
For the purposes of clustering high dimensional data, issues arise due to the curse of dimensionality (Di Marzio and Lafratta, 1999). A standard approach is to reduce the number of dimensions by considering a series of underlying factors with a lesser dimension (Spearman et al., 1950). Let X i represent an r dimensional random vector, with x i as a realization. The factor analyzers model for X 1 , . . . , X N , is given by where µ is a mean location vector, Λ is a r × s matrix of factor loadings, with s < r, U i ∼ N r (0, I) denoting the latent factors, and i ∼ N r (0, Ψ) where Ψ = diag(ψ 1 , ψ 2 , ..., ψ r ).
Here, N r denotes the r-dimensional multivariate normal. Furthermore, the latent factors U i and noise i are independent of each other. It is noted that the probabilistic principal component analysis (PPCA) is a special case of the factor analysis model with a specific isotropic constraint on Ψ (Tipping and Bishop, 1999). The factor analyzers model is a flexible extension of PPCA. By removing the isotropic constraint on Ψ, this allows parameters to vary. The factor analyzers model is considered to be the best choice for dealing with telemetric data as it is highly efficient in reducing dimensionality (Inui et al., 2009).
Suppose a matrix is considered to be an observation sampled from a distribution. Naturally, an appropriately sized matrix variate distribution should be considered to model randomness. Consider the matrix variate normal distribution (Gupta and Nagar, 1999).
Let X be a random variable with an r × c matrix X as a realization. As a result, X is distributed according to a matrix variate distribution. The random matrix X(r × c) is said to have a matrix variate normal distribution with mean matrix M (r × c) and covariance matrix Ψ ⊗ Σ. Each matrix is appropriately sized as Σ(r × r), Ψ(c × c), where vec(X) ∼ N rc (vec(M ), Ψ ⊗ Σ). Here, ⊗ refers to the Kronecker Product and vec denotes the vectorization of a matrix. Given this specification, the density is formulated as The matrix variate normal distribution is defined through a vectorization of a multivariate normal. Note that the covariance matrices of row and column are non-unique as they are defined through a Kronecker product (Dutilleul, 1999). As a result, both densities are parametrized by the product and not individual co-variance matrices (Gupta and Nagar, 1999). There are benefits for using a matrix variate representation. The main benefit of the matrix variate approach is the speed at which model parameters are estimated in high dimensional settings.
Due to the issues of high dimensionality, an analogous extension of the factor analyzers model for matrix variate data is implemented. The mixtures of matrix variate bilinear factor analyzers model (MBI) is a powerful approach for dealing with both high-dimensional data, and the presence of a mixture of populations (Gallaugher and McNicholas, 2018). Suppose latent factors of size s < r, v < c, for a matrix variate random variable constitute the data with probability π g of occurring as is a matrix random variate of latent factors, A g (r × s) are column factor loadings, and B g (c × v) are row factor loadings respectively. Finally, the noise is distributed according to For applications in accelerometer telemetry, the MBI model is used to cluster matrix variate objects pertaining to the heterogeneous population of participants.
The estimation procedure for MBI is based on local maximum likelihood estimation.
The most common approach for estimating finite mixture models is with the expectation maximization (EM) algorithm (Dempster et al., 1977). However, when dealing with latent factor models, McNicholas and Murphy (2008)  For example, the component membership for observation 1 is given as z 1 := (z 11 , . . . , z 1G ).
Suppose observation i is in group g, the formulation of factor analysers in the matrix variate case has the density as When written in this formulation, the complete data-likelihood is taken to be  Gallaugher and McNicholas (2018). The AECM algorithm consists of three stages. Within the first stage, the complete-data is taken to be the observed matrices X 1 , ..., X N , and the component memberships z = (z 1 , ..., z N ). The updates for π g , and M g , for some iteration t are calculated aŝ In the second stage, the complete-data is taken to be the observed X 1 , ..., X N , the component memberships z, and the r × s latent matrices for column factors. In addition, N ig . The updates for U are taken to bê In the third stage, the complete-data is taken to be the observed X 1 , ..., X N , the component memberships z, and the c×v latent matrices for column factors allowing to compute estimate for V asV Convergence of the AECM algorithm is based on the Aitken acceleration criterion (Aitken, 1926) defined as

1−a (t)
be the observed estimate after many iterations at t + 1. Termination of the algorithm occurs when l (t+1) ∞ − l (t) ∈ (0, ε) for some pre-specified ε. Model selection is based on the Bayesian Information Criterion (BIC, Schwarz et al., 1978). The BIC is a measure to assess the performance of the model fit, while penalizing for the number of parameters used. For interpretability, the positive scale BIC is used for assessing model performance (larger is better). Let ρ be the number of parameters used. The positive scale BIC is then formulated as BIC = 2l(θ) − ρ log N. With all methods relating to MBI introduced, the clustering problem is formulated as follows. Assuming there exists a heterogeneous population of participants of up to G types. Let X i be a force map of participant i. Formally, X ∼ N G r,c (θ) with probability density function where θ = (θ 1 = (π 1 , M 1 , U 1 , V 1 , A 1 , B 1 ), . . . , θ G ). Estimation of the model is done in accordance with the AECM algorithm as mentioned previously. Once the model has been estimated, classification of participants into one of G types is done in accordance with taking the maximum aposteriori of component memberships. The classification C i for participant i is given by C i = argmax(ẑ i1 , . . .ẑ iG ). With each participant classified into a group, we pursue survival analysis to determine differences in lifetime outcomes.

Survival Analysis
As a standard practice, survival analysis in this work is performed by the use of the Cox proportional hazards model (Cox, 1972). In addition, from a non-parametric perspective, we also fit Kaplan-Meier (KM) estimators to the data (Kaplan and Meier, 1958). For both models, we test for significance using the log-rank test (Kleinbaum and Klein, 2012). For implementation, see survival and survminer packages (Therneau, 2020;Kassambara et al., 2017). In summary, we use these methods to determine whether or not survival outcomes differ between groups.
We report our results in chronological order of analysis. We first report our survival analysis on both SW-0 and SW-P to establish a baseline understanding of our shift work population.
All survival analysis is performed with the R programming language (R Core Team, 2020).
Next, we report our clustering results of SW-P's force maps. All clustering and data manipulation is performed using the Julia programming language (Bezanson et al., 2017). For an in-depth look on how to interoperate between R and Julia see McNicholas and Tait (2019).
Finally, we report our survival analysis results by cluster and sex.

Work Shift Survival Analysis
Our goal for this portion of our analysis is twofold. First is to reach consensus with the standard occupational health viewpoint that life expectancy suffers when participants defy their natural circadian rhythm. Our two populations, SW-0, and SW-P, contain a sufficient number of individuals for analysis to be valid. We fit KM estimators on both SW-0 and SW-P. In addition, we performed a log-rank test. Both populations SW-0, and SW-P yield statistical significant results indicating the following. Late shift workers have a lower survival rate than regular shift workers. Figures 3a, and 3b illustrate the KM estimates over time, where time is measured in years. The curves between both of the populations are virtually the same. In both plots, both types of shift workers have the same survival estimate until the age of 57. At this point, the two shift worker groups start to diverge, until reaching the age of 80. The log-rank tests for both SW-0, and SW-P yield a p-value of less than 0.001. Despite there being a loss in the number of participants between SW-0 and SW-P (see section 2.2), SW-P still retains the relevant survival outcomes of its predecessor SW-0.
The results of the Cox regression yields similar results. Figures 3c, and 3d report the hazard ratios in the SW-0, and SW-P populations respectively. Again we find that the regular shift workers have a better survival outcome than late shift workers. As expected, there are ostensibly no differences between the results of SW-0, and SW-P. However, the range of SW-P is slightly higher with a right bound of 0.89, than SW-0's 0.88. In summary,  both analyses yield statistically significant results, and conclude the same survival outcomes.
Regular shift workers have a hazard rate that is less than late shift workers. This is in-line with the general consensus of Gu et al. (2015), and Ijaz et al. (2013).

Force Map Clusters
As a preparation step, we perform a logit transformation on each entry across all force maps to map the domain from R[0, 1] to R. This allows the MBI model to have a better fit for the data while retaining the force map structure across each participant. Another preparation step is regarding initialization of cluster memberships. There are a number of initializations possible, however for simplicity we use k-means, and random soft (k meaning the number of groups in this context). For specifics regarding these two methods, see McNicholas (2015).
There are a number of settings to consider for the MBI model. We search for the best model with the following settings: q = {1, . . . , 7}, s = {1, . . . , 7}, G = {1, . . . , 6}. Table 3 shows the results of our model search. According to BIC, we report the top 5 performing MBI models.
The best model is selected to be G = 3 for the number of groups, a BIC of −9527009, and q = 5, s = 4 for the dimension of latent factors.
By clustering each participant of SW-P into one of three groups, we define three new populations of shift workers as SW-1, SW-2, and SW-3. Table 4 displays the cluster results for each group. There are 407 participants clustered in SW-1 which account for smallest portion of the total population of SW-P. In contrast, there are 1, 708 participants clustered into SW-2. These participants account for 33.85% of SW-P. Finally, SW-3 accounts for the majority with 2, 911 participants. We draw attention to the mean force map of each cluster in Figure 4. The mean force map for SW-1 seen in Figure 4a, has most of the readings around 25 force units and is quite concentrated. In contrast, the mean force map for SW-2 has most of the readings centred around 50 force units but is quite spread. Finally, the mean force map of SW-3 has most of the readings centred around 73 force units, and is the most spread by far. The interpretation of each force map is as follows. Participants who exhibit stronger physical activity during moderate tasks have higher force readings. As a result, their force map becomes variable and less concentrated. Each cluster from SW-1 to SW-3 embodies a gradual increase in physical exertions. In SW-1 we have participants who have the lowest level of physical exertion, while in SW-3 we have the highest. In summary, our clusters embody the heterogeneous levels of physical behaviour within the SW-P population. Now that we have established our populations of interest, we further investigate our clustering results using survival analysis.

Cluster Survival Analysis
Using our previously defined populations SW-1 to SW-3, we proceed with survival analysis by cluster. Figure 5 shows three decreasing levels of risk. We report a statistically significant global log-rank a p-value of less than 0.001. Cluster 1 is used as a reference for the cox proportional hazards model. Both cluster 2 and cluster 3 have a lower risk than the reference cluster. We further fit a Cox regression between cluster 2 and 3, this time using cluster 2 as a reference. Here, we see that cluster 3 has indeed a lower risk than cluster 2. Again, the log-rank test shows a statistically significant result. z Figure 6: Cox regression results between SW-2 and SW-3.
To summarize, we describe three decreasing levels of risk where each cluster has a statistically significant difference in hazard rates. These survival outcomes are consistent with our force map analysis in Section 3.2. Cluster 1 force maps have the lowest level of physical exertion which increases between clusters. Similarly in our survival analysis, cluster 1 participants have the highest risk that decreases from cluster 2 to 3. We also fit a cox regression to analyse a mixed interaction between clusters and type of shift work. Particularly to individuals belonging to cluster 1, and engaging in regular work, as well as individuals who belong to cluster 2 and 3, that engage in late shift work. This enables us investigate if individuals who engage in late shift work and have higher physical activity, may offset their risk when compared to regular shift workers that do not engage in physical activity. Figure   7 shows the hazard ratios of these populations. A statistically significant difference between cluster 1 regular shift workers, and cluster 3 late shift workers is reported. We further fit the same cox regression only on these two populations. As expected, Figure 8 shows a statistically significant difference between rates of these populations. Consistently, the KM estimators capture the same decreasing survival outcome as shown in Figure 9. We interpret these results as follows. Individuals who are defined as late shift workers, and exert the highest level of physical activity, have a better survival outcome than individuals who undergo regular shift work, and exert the lowest physical activity. Figure 10 shows separated KM curves in consensus with our cox regression results in Figure 8. Generally, across all KM curves, there is no evidence of censoring within each population. Each cluster or sub-group contains an ample amount of individuals to reliably estimate survival probabilities using a KM estimator.

Cluster Survival Analysis by Sex
When taking into consideration sex differences, consider the proportion of males and females within each cluster. Table 5 shows a clear majority of males among clusters 1 and 2, while a female majority among cluster 3. As mentioned previously, there are more favourable survival outcomes for females when compared to males. However, when performing Cox regression by sex, we have the same three decreasing levels of risk as in Section 3.3. Figures   11 and 12 show statistically significant differences between each of the clusters. However, females have stronger differences in their hazard ratios when compared to males.  Figure 14: Hazard ratios for males by mixed groups between type of shift work, and cluster.

Conclusion
For individuals that defy their circadian rhythm, through for example, working shift patterns, our results indicate that an increase in physical activity may be enough to offset the additional risk associated with late shift work. Although our results indicate that there are no statistically significant differences for males, it beckons the question of why this protective effect is only present for SW-P females. During our clustering procedure, male and female force maps were not separated, and therefore, the stronger signal from the female population may have caused this phenomenon. Nevertheless, a segmentation by sex approach can be easily incorporated into our methodology.
By considering the sheer size and computational cost of managing this dataset, our approach is feasible and scalable for larger cohorts. The methodology proposes a novel set of statistical tools for capturing participant behaviour within accelerometer data. The force maps characterize participant behaviour and allow for visually interpretable results. Furthermore, by introducing this matrix variate object, we standardize an otherwise incomparable time series; allowing for comparisons to be made between participants.
In conclusion, we utilize accelerometer data in a manageable and efficient process to cluster participants accordingly. The results of these clusters indicate valid survival outcomes with great potential for improving health from an occupational perspective. We developed a low-dimensional feature that is interpretable, and captures participant behaviour in a unique way. Future work, from a clinical perspective, may consider the possibility of investigating the protective nature of physical activity for those who defy their circadian rhythm.