Decomposing Juggling Skill into Sequencing, Prediction, and Accuracy: A Computational Model with Low-Gravity VR Training

TL;DR

This study breaks down juggling into key components and uses VR to improve learning, showing how different skills contribute to overall performance.

Contribution

The first computational decomposition of juggling into sequencing, prediction, and accuracy using a multimodal evaluation system.

Findings

Sequencing is the dominant factor in early juggling skill acquisition.

Prediction and accuracy become increasingly important as skill develops.

Reduced gravity VR training enhances early-stage motor learning.

Abstract

Juggling is a complex motor skill that requires multiple sub-skills and cannot be mastered without extensive practice. Although prior studies have quantified performance differences between novice and expert jugglers, none have attempted to quantitatively decompose these components or model their contribution to juggling performance. This longitudinal study presents a multimodal evaluation system that integrates computer vision, motion capture, and biosensing to quantify three key elements of juggling ability: Sequencing, Prediction, and Accuracy. Twenty beginners completed a 10-day, three-ball juggling experiment combining visuo-haptic virtual reality (VR) and real-world practice, with half training in reduced gravity, previously shown to enhance early-stage motor learning. The fitted Gamma-Log generalized linear model (GLM) indicated that Sequencing is the dominant factor of early skill acquisition, followed by Prediction and Accuracy. This study provides the first computational decomposition of juggling, demonstrates how multiple elements jointly contribute to performance, and results in a principled approach to characterizing motor learning in complex real-world tasks.

Figures

Overview of processing pipeline for constructing Sequencing, Prediction, and Accuracy. (a) Display of trimming the RGB video into proper trials. The red horizontal line is a shoulder threshold, while the blue line is a knee threshold. LED glow balls’ positions are estimated and marked as green circles. (b) Display of detecting balls and hands positions of each frame. The balls are labeled as green circles, while the hands are labeled as red circles. (c) A graphical example of reconstructed tracks of balls and hands by the tracker. The x-axis is time (s), and the y-axis is height (pixel). The b

Full text

1. Introduction

The acquisition of motor skills refers to long-lasting improvements in movement ability through practice, which are shown through their retention and transfer to new situations [1]. Motor learning occurs on multiple time scales: rapid gains can occur within a single practice session, while slower consolidation occurs over days or weeks [2]. Classic stage models describe skill acquisition as a progression from cognitively demanding, error-prone performance to more automatic, stable execution [3,4]. In parallel, dynamical systems and computational motor control theories emphasize how coordinated movement patterns emerge and how internal models and optimal feedback control are used to reduce variability and improve efficiency [5,6,7]. Together, these views suggest that complex skills are not supported by a single factor but by multiple interacting sub-processes (e.g., timing, prediction, spatial precision) that may change at different rates during learning. Although motor learning has been widely studied at the behavioral and training levels, building computational models that map these theoretical ideas onto concrete, measurable components remains challenging because human movement involves many redundant degrees of freedom [5], individuals differ in their strategies and constraints, and performance is strongly affected by task context [7]. In general, these factors make it difficult to isolate the contribution of specific elements and to connect them to overall skill.

Three-ball cascade juggling provides a useful model system for studying complex motor learning. Performance can be quantified in a relatively objective way, for example, by measuring how long a participant can maintain the pattern or by counting successful throws and catches [8]. At the same time, juggling is a complex skill because it requires several physical and cognitive processes to be coordinated simultaneously. Learners must perform consistent throws (e.g., stable speed and release angle), track multiple balls visually, and maintain rhythmic timing for continuous catch-throw cycles [9]. Because these components must be tuned together, most people cannot acquire juggling skill without sufficient practice.

For these reasons, juggling has been used in many physiology and neuroscience studies as a model of complex motor skill acquisition. Training-related changes in white matter architecture, including regions near the intraparietal sulcus, have been reported using diffusion tensor imaging (DTI) [10], and increases in gray matter volume in visual-processing areas have been observed using voxel-based morphometry (VBM) [11]. Behavioral studies also show that expert jugglers use more efficient gaze strategies, with tighter gaze–ball coupling and reduced horizontal gaze shifts, consistent with predictive perception–action control [12]. In addition, spatiotemporal variability in cascade juggling shows regular patterns rather than random noise, so rhythm stability and coordination can be quantified [13].

However, most prior studies did not systematically decompose juggling into multiple interacting motor and cognitive elements. Some work has begun to extract lower-dimensional, interpretable components. Post et al. used principal component analysis (PCA) to reveal low-dimensional coordination structures underlying three-ball cascade juggling [14]. Müller and Sternad showed that skill gains reflect structured changes in execution variability, namely noise reduction, tolerance exploitation, and covariation (N-T-C), rather than a simple decrease in variance [15]. Ramsay et al. applied functional data analysis (FDA) to juggling trajectories and introduced landmark and continuous registration methods to separate phase and amplitude, allowing cycle-wise alignment across participants and sessions [16]. Finally, Yamamoto and colleagues found that learning does not change smoothly. Instead, it advances in stages constrained by the task, with stable dwell-time ratio levels (related to the Farey sequence), and adaptability differs between coordination patterns, which indicates pattern-specific learning dynamics [17,18]. Despite these advances, two key gaps remain in a computational explanation of how multiple elements work together to determine overall performance. First, most studies focus on a single facet (e.g., PCA [14], FDA [16], or N-T-C [15]), so there is no integrated framework that can estimate multiple candidate elements at the same time. Second, decomposed measures are rarely linked to a single predictive model that explains overall juggling performance with formal validation.

Outside the juggling domain, recent multimodal studies in transportation and human–machine interactions show how synchronized behavioral and sensor data can be used to estimate latent cognitive or motor states, for example, by combining gaze and control inputs to model pilots’ situation awareness [19] and integrating inertial, physiological, and environmental signals for robust human activity recognition [20]. More general data-driven frameworks for motor learning also emphasize the fusion of kinematic, physiological, and contextual information to describe how complex skills change over time [21]. In rehabilitation engineering, adaptive exoskeleton controllers based on admittance models adjust assistance according to the user’s evolving motor state [19]. These studies highlight the value of multimodal sensing and component-level descriptions of motor ability.

Building on these insights, we develop a component-based computational model of early juggling learning. In this model, juggling skill is described by three complementary facets suggested by motor control theory: spatiotemporal coordination (Sequencing), predictive state estimation (Prediction), and spatial precision (Accuracy). Sequencing (S) captures the rhythmic, phase-locked organization of throws and catches across cycles, corresponding to the stable coordination patterns and movement synergies in dynamical system frameworks [4,5]. Prediction (P) reflects the ability to estimate ball states during visual occlusion and to perform appropriate catch actions based on internal models of object dynamics, in line with optimal feedback control and internal model theories [6,7]. Accuracy (A) quantifies the spatial consistency of second-catch locations in a hand-centered frame, indexing endpoint precision and control of execution variability, as emphasized in variability-based accounts of skill [15]. These three sub-skills form a low-dimensional, interpretable set that can be directly derived from our kinematic and occlusion data and linked to overall juggling performance in a single computational model.

We also considered other possible decompositions. One alternative is to treat juggling as mainly feedforward once the cascade becomes stable, where performance could be explained mostly by Sequencing and Accuracy without an explicit Prediction component. However, in early learning, the ball trajectory is not consistent across cycles, so novices likely rely on prediction-based state estimation to support rapid online corrections. This view aligns with optimal feedback control and internal model accounts, which emphasize predicting movement outcomes and using feedback to reduce the impact of motor variability [6,7,22]. For this reason, we retain Prediction as a distinct component and probe it using the visual-occlusion catching task.

Also, among these components, Sequencing reflects the ordered and rhythmic structure of the cascade, which is often the main difficulty when beginners start learning. Based on this idea, our visuo-haptic VR training system was designed to slow down ball movement by reducing gravity so that novices can perceive and practice the juggling rhythm with less time pressure. Because the only experimental intervention between groups was this low-gravity slow-tempo training, we also explore whether group differences appear mainly in Sequencing as a behavioral check of the intended training effect.

Based on this framework, we tested three hypotheses regarding the computational structure of juggling skill:

Hypothesis 1(H1). Juggling can be decomposed into three latent components: Sequencing (S), Prediction (P), and Accuracy (A) (Figure 1). Following our previous definitions [23], S represents spatiotemporal coordination among balls across cycles, P represents the ability to predict the landing position of a visually occluded ball, and A represents the spatial consistency of catch locations.

Hypothesis 2(H2). Day-level juggling performance ( $[eqn]$ ) can be modeled as a function of S, P, and A, i.e., $[eqn]$ , with independent and complementary contributions from each component.

Hypothesis 3(H3). During early learning, S, P, and A do not increase at the same rate, but show lead–lag relationships across days (asymmetric development). Because Sequencing is a main bottleneck for beginners, we expect slow-tempo low-gravity VR training to mainly promote early Sequencing gains.

The contribution of this study is four-fold:

1. Building on our previous VR study, we report, to our knowledge, the first 10-day juggling training multimodal evaluation system that combines computer vision (ball and hand kinematics), motion capture (3D wrist and shoulder trajectories), and biosignals (EEG and EMG), all synchronized on a common time base.
2. We develop a fully automatic analysis pipeline that replaces manual video scoring and extracts trial duration, number of successful catches, and ball trajectories at scale from the multimodal recordings.
3. Using these extracted measures, we innovatively define three components—Sequencing, Prediction, and Accuracy—and connect them in a computational model that jointly explains day-level juggling performance.
4. Leveraging the longitudinal design with repeated measurements in novices, we quantify how these three components change across days and demonstrate that they follow distinct learning rates during early acquisition.

The remainder of this paper is structured as follows: Section 2 outlines the experiment setup, including the research hypotheses; experimental design; sensors and data collection; the construction of Sequencing, Prediction, and Accuracy indices; and statistical modeling. Section 3 presents the main findings, followed by further interpretation and discussion in Section 4, together with limitations and future works. Section 5 concludes this study.

2. Materials and Methods

This study used the same ten-day juggling experiment as our previous work on VR-based motor learning, which introduced the behavioral effects of VR training on real-world juggling. We introduced different gravity-level settings in the VR training in two groups: a progressive low-to-normal gravity group (hereafter, Low-g group) or a constant normal-gravity group (hereafter, Normal-g group). This paper focuses on building three skill indices, Sequencing (S), Prediction (P), and Accuracy (A), from kinematic data and modeling day-level juggling performance as a function of these indices. Detailed descriptions of the visuo-haptic VR setup, gravity intervention, and biosensor recordings (EEG, EMG, and functional magnetic resonance imaging (fMRI)) are provided in the work of Cho et al. [23].

2.1. Participants

Twenty healthy adults with no prior juggling experience joined the experiment (14 males, 6 females; age $[eqn]$ years). Participants were randomly assigned, with matched sex ratios, to either a Low-g group (n = 10; 7 males, 3 females; age $[eqn]$ years) or a Normal-g group (n = 10; 7 males, 3 females; age $[eqn]$ years). All participants provided written informed consent. The study protocol was approved by the ethics committee of the Tokyo Institute of Technology (currently, Institute of Science Tokyo) (approval no. 2020107).

2.2. Experimental Design

Each participant completed 10 consecutive days of training. On each day, the following four sessions were performed in a fixed order:

• Session 1:Real-world three-ball juggling (pre).
• Session 2:VR three-ball juggling.
• Session 3:Visual-occlusion juggling with shutter glasses.
• Session 4:Real-world three-ball juggling (post).

The visuo-haptic VR training in Session 2 includes low-gravity intervention and group allocation (progressive Low-to-Normal gravity vs. constant gravity), as described in detail by Cho et al. [23]. The training protocol used a slow-tempo, low-gravity visuo-haptic VR system; we hope that by reducing gravitational acceleration and varying the balls’ flight time, the learning of cascade rhythms may be promoted, potentially improving the development of Sequencing during early practice. Briefly, the Low-g group practiced with progressively increasing gravity (0.2 g, 0.4 g, and 0.6 g on Days 1–3, 4–6, and 7–10, respectively), whereas the Normal-g group trained at a constant 0.6 g, as shown in Figure 2. This study focuses on kinematic indices and their relationship to juggling performance, while group-wise patterns are presented descriptively to explore the behavioral trends. The full experiment was designed as a multimodal platform combining kinematic recordings, VR logs, EEG, EMG, and pre-/post-training fMRI.

2.3. Sensors and Data Collection

2.3.1. Overview of Sensors

The experiment collected multimodal data from kinematic, VR, and neurophysiological systems. All devices were time-synchronized using Lab Streaming Layer (LSL; liblsl v1.13.0) [24]. In the present analysis, we only used the kinematic streams from motion capture, RGB-D video, and shutter timing to compute the Sequencing (S), Prediction (P), and Accuracy (A) indices. The remaining modalities are listed briefly to describe the overall study design.

Optical motion capture (used for P, A): An OptiTrack (NaturalPoint Inc., Corvallis, OR, USA) system recorded 3D positions of four reflective markers placed on both shoulders and both wrists at 200 Hz.RGB-D camera (used for S, P, A): A front-facing Intel RealSense D455 (Intel Corporation, Santa Clara, CA, USA) recorded color video and depth at 60 Hz during real-world juggling and the occlusion task, capturing ball trajectories and hand movements (librealsense2 v2.49.0).Shutter-glasses trigger (used for P): signals indicating the opening and closing of liquid-crystal shutter glasses (PLATO Visual Occlusion Spectacles; Translucent Technologies Inc., Toronto, ON, Canada [25]) were recorded as an LSL stream and aligned with motion-capture and video data.

As part of the broader project [23], additional signals were also recorded:

• VR system logs: Positions of virtual hands and balls in the visuo-haptic VR environment.
• EEG: A 64-channel EEG (ActiveTwo; BioSemi B.V., Amsterdam, The Netherlands; 2048 Hz) on selected days.
• EMG: Wireless EMG (Trigno; Delsys Inc., Natick, MA, USA; 12 sensors) from upper-limb muscles on selected days.
• fMRI: Functional MRI scans before and after the entire training experiment.

Analyses of EEG, EMG, and fMRI are outside the scope of this study and will be reported separately.

2.3.2. Real-World Three-Ball Juggling

In Sessions 1 and 4, participants performed real-world three-ball cascade juggling using LED glow balls under controlled lighting to improve tracking quality. Each daily session lasted about 5 min (15 min on Days 1, 5, and 10). Motion capture and RGB-D video were recorded continuously. These sessions were used to evaluate overall real-world juggling performance ( $[eqn]$ ) and to provide the data for computing the Sequencing index (S).

2.3.3. Visual-Occlusion Task

In Session 3, participants performed a two-ball juggling task while wearing liquid-crystal shutter glasses (PLATO Visual Occlusion Spectacles [25]). They alternately threw and caught two balls with opposite hands according to verbal instructions. An RGB-D camera recorded ball trajectories and hand movements, and a motion-capture system tracked markers placed on both shoulders and wrists.

Visual occlusion was controlled by an apex-locked system. Ball centroids were tracked online from a second RGB-D camera (Intel RealSense D435; Intel Corporation, Santa Clara, CA, USA; librealsense2 v2.49.0). The vertical position of each ball was temporally smoothed, and apex events were detected when the vertical velocity crossed zero while the vertical acceleration was negative. When the second ball reached its apex, both lenses of the shutter glasses were closed for 0.5 s, transiently blocking vision of the descending ball.

On most training days, participants completed 60 throw–catch repetitions. On Days 1, 5, and 10, the task consisted of four blocks of 75 attempts (300 trials in total), evenly divided among three conditions: (i) apex-locked occlusion, (ii) no occlusion, and (iii) randomly timed occlusion, with condition order randomized within each block. These trials provided the data to be used to compute the Prediction (P) and Accuracy (A) measures in the subsequent analyses.

2.4. Construction of Sequencing, Prediction, and Accuracy Indices

2.4.1. Sequencing (S) from Real-World Three-Ball Juggling

Sequencing (S) quantifies how well the three balls maintain a stable 120° phase relationship during real-world three-ball cascade juggling. To avoid extensive manual labeling, we designed an automatic Python (v3.12.7) pipeline to compute Sequencing with four stages (see Figure 3a–d): (i) video trimming and trial timing, (ii) frame-wise ball/hand detection, (iii) ball tracking and event labeling, and (iv) phase estimation and Sequencing score computation.

Video trimming and trial timing. The RGB video from Sessions 1 and 4 was segmented into individual trials using Python. Mediapipe Pose (v0.10.18) [26] was used to track both hands, and OpenCV (v4.12.0)-based blob detection [27] identified candidate ball positions. Based on the relative positions and movements of the hands and balls, the program automatically determined when each juggling trial began and ended (Figure 3a).Frame-wise ball and hand detection. Within each trimmed juggling trial, videos were processed frame by frame in the second stage. For each frame, it recorded the detected ball candidates and the hand positions estimated by OpenCV [27] and Mediapipe [26]. These frame-wise detections were stored in a compact columnar format (parquet) and served as the input for the tracking stage (Figure 3b).Ball tracking and event labeling. Using the frame-wise detections, the third stage reconstructed three continuous ball trajectories for each trial. Throws were labeled when a ball left the vicinity of a hand, and catches were labeled when a ball entered and stayed near the opposite hand. By connecting detections across frames, the algorithm followed each ball over time and produced smooth trajectories with associated throw and catch events (Figure 3c).Phase estimation and Sequencing score. In the final stage, we used the reconstructed trajectories to compute phases and derive a Sequencing score (Figure 3d). For each ball $[eqn]$ , the smoothed and mean-subtracted vertical position $[eqn]$ was converted to an instantaneous phase $[eqn]$ using the analytic signal from the Hilbert transform [28],

$$[eqn]$$

where $[eqn]$ denotes the Hilbert transform.For windows containing three consecutive apexes from three different balls, we quantified how closely the three phases maintained the ideal $[eqn]$ spacing. Let $[eqn]$ be the phase of ball i in degrees and $[eqn]$ the usable time samples within window w ( $[eqn]$ ). The window-wise Sequencing score was defined as

$$[eqn]$$

where $[eqn]$ is the circular difference (modulo $[eqn]$ ) between $[eqn]$ and the ideal $[eqn]$ spacing. Thus, $[eqn]$ takes values in $[eqn]$ , with $[eqn]$ indicating perfectly maintained $[eqn]$ phase offsets. Figure 3. Overview of processing pipeline for constructing Sequencing, Prediction, and Accuracy. (a) Display of trimming the RGB video into proper trials. The red horizontal line is a shoulder threshold, while the blue line is a knee threshold. LED glow balls’ positions are estimated and marked as green circles. (b) Display of detecting balls and hands positions of each frame. The balls are labeled as green circles, while the hands are labeled as red circles. (c) A graphical example of reconstructed tracks of balls and hands by the tracker. The x-axis is time (s), and the y-axis is height (pixel). The blue, orange, and green lines describe trajectories of 3 different balls, while the pink and light blue lines show the vertical movements of both hands. Sequencing is computed based on these reconstructed balls’ phases. (d) Phase estimation in the Phasor stage. The top panel shows vertical trajectories of the three balls, with apex times used to define windows. The middle panel shows the instantaneous phase of each ball (degrees), and the bottom panel shows pairwise phase differences compared with the ideal 120° phase difference. (e) Example of success and failure of the visual-occlusion task for computing Prediction (P). (f) Illustration of probability ellipsoid for computing Accuracy (A). The x-axis is the lateral distance (m) from the center of a participant’s body, while the y-axis is the front–back distance (m) [29]. Each black dot connected with a blue line means the position of a first catch with the right hand (in this case), and each black dot connected with a red line means the position of a second catch with the left hand. Pink ovals represent the 95% probability ellipsoid of second catch positions of each day. Most participants show a decrease in the ellipsoid volume throughout the training period (left → right). To quantify the sequencing ability for each day, we pooled all valid sequencing windows ( $[eqn]$ ) obtained from trials yielding at least one valid window across the analyzed sessions (Sessions 1 and 4). The day-level Sequencing index (S) was computed as the arithmetic mean of these pooled window values.

2.4.2. Prediction (P) and Accuracy (A) from the Visual-Occlusion Task

Prediction (P) and Accuracy (A) were derived from the visual-occlusion task in Session 3 using synchronized motion-capture and video recordings (Figure 3e,f).

Trial segmentation and labeling. We used Lab Streaming Layer (LSL) XDF files containing video timestamps, motion-capture markers for both wrists and shoulders, and shutter timing signals. In MATLAB (R2021b), continuous recordings were segmented into individual trials based on changes in hand velocity and acceleration, and each trial was saved as a separate .mat file. For each trial, we annotated both the timing and success of the catches for each hand (including the second catch). Trials with shutter malfunction, severe marker loss, or unclear outcomes were excluded from subsequent analyses.Prediction (P): Catch success in the visual-occlusion task. Prediction (P) summarizes how reliably both hands complete the catches during the two-ball visual-occlusion task (Figure 3e). For each participant p and day d, let $[eqn]$ denote the number of valid trials in Session 3 in which both hands successfully completed the second catch, and let $[eqn]$ be the total number of valid trials on that day. The day-level Prediction index was defined as

$$[eqn]$$

with higher values indicating more reliable catch performance in the visual-occlusion task.Accuracy (A): Spatial dispersion of catch locations. Accuracy (A) quantifies how tightly the second-catch positions cluster in 3D space (Figure 3f). This two-ball shutter-glasses paradigm was chosen because it has many controlled catch samples and allows us to quantify hand-centered spatial precision under reduced visual feedback, which is harder to measure reliably in short and unstable three-ball cascades during early learning. For each participant p and day d, we expressed the 3D position of each second catch in a hand-centered coordinate system defined by the corresponding wrist marker. In Session 3 on Days 1, 5, and 10, two-ball juggling was performed in four blocks; for each valid block $[eqn]$ and each hand $[eqn]$ , all second-catch positions were pooled to form a 3D point cloud, and a 95% probability ellipsoid was fitted to this cloud, yielding an ellipsoid volume $[eqn]$ .Day-level Accuracy was then defined as the geometric mean of these volumes,

$$[eqn]$$

so that smaller values of $[eqn]$ indicate tighter clustering and, therefore, higher spatial Accuracy of the second catches. Blocks in which the ellipsoid fit was unstable (e.g., nearly singular covariance) were excluded before computing the geometric mean.Aggregation and data merging. Day-level indices P, A, and S were summarized for Days 1, 5, and 10 and merged into a single dataset with participant and group information for the subsequent statistical analyses.

2.5. Statistical Modeling of Performance

2.5.1. Data Preprocessing

All statistical analyses were conducted in MATLAB. From the preceding processing steps, we obtained a merged day-level dataset for each participant p and training day $[eqn]$ : (i) the three skill indices $[eqn]$ , $[eqn]$ , and $[eqn]$ ; (ii) a behavioral performance measure $[eqn]$ ; and (iii) group labels (Low-g/Normal-g) and participant IDs. To ensure comparability and avoid bias from different practice amounts, the main analysis focused on three standardized test days (Days 1, 5, and 10), while the remaining days were used to describe trends and test robustness.

Day-level juggling performance $[eqn]$ was defined from the real-world three-ball sessions (Sessions 1 and 4). For each participant and day, we pooled all trials from these two sessions, counted the number of successful catches in each trial, and identified the 10 trials with the longest successful catches. The performance measure $[eqn]$ was then defined as the mean of successful catches across these top 10 trials.

Sequencing values $[eqn]$ were taken from the phase-based analysis described above. To ensure statistical reliability, we identified and excluded outliers based on the distribution of valid window counts ( $[eqn]$ ). Given the positively skewed nature of the data, we applied a log transformation and established a conservative lower bound of three standard deviations below the mean ( $[eqn]$ ). This statistical analysis yielded a cutoff of approximately 4.5 windows. Accordingly, sessions falling below this threshold ( $[eqn]$ )—specifically, three participants on Day 1 who produced 0, 1, and 4 valid windows—were identified as statistical outliers. Since the Sequencing index is computed as a mean, such extremely low sample sizes inherently lead to inflated standard errors and unstable estimates. Therefore, these cases were set to missing for the day-level aggregation. We evaluated the robustness of excluding these cases in a sensitivity analysis reported in Section 3.3.

Since S, P, and A had different ranges and distribution shapes, we applied corresponding logit/log transforms and then z-transformation to them so that model coefficients reflect comparable effects.

Prediction values P are bounded between $[eqn]$ and left-skewed, so we applied logit transform to it as follows:

$$[eqn]$$

Since Accuracy values A are positive with a long right tail, we applied a negative log transform so that larger A correspond to better Accuracy while stabilizing the distribution,

$$[eqn]$$

Sequencing values were used directly as follows: $[eqn]$ .

For each transformed index $[eqn]$ , we then applied a robust z-transformation based on the median and median absolute deviation (MAD) to reduce sensitivity to outliers [30].

$$[eqn]$$

Figure 4 shows how preprocessing through monotonic transforms and z-transforms changed the distributions of the three indices. Raw Prediction (P) was bounded in $[eqn]$ and highly concentrated near 1.0, showing a strong ceiling effect [31]. After the logit transform, the upper range was expanded, and the distribution became more symmetric and more suitable for linear modeling. Raw Accuracy (A) was strictly positive and right-skewed, with a small number of large-volume trials. The negative log transform reduced this skewness and reversed the direction, making the distribution more symmetric. Raw Sequencing (S) was already relatively well-behaved, so only the z-transform was applied to make it on the same standardized scale as P and A. Overall, the distinct raw distributions indicate that S, P, and A capture different behavioral properties rather than a single common factor, and the distributions of three skill metrics were improved by transforms to be more stable and suitable for fitting generalized linear models.

2.5.2. Generalized Linear Models

Comparing with the Normal-identity generalized linear model (Normal-identity GLM) [32], this work presents a better fit using a Gamma-Log generalized linear model (Gamma-Log GLM) [33] to fit the pooled Day 1/5/10 data. In a Gamma-Log GLM, the expected performance $[eqn]$ is linked to a linear predictor through the log function, while the response itself is assumed to follow a Gamma distribution:

$$[eqn]$$

Gamma distribution is a better fit because the day-level juggling performance ( $[eqn]$ ) is strictly positive and right-skewed, where there are many low values with a tail of high values. Compared to Normal distribution, Gamma distribution is specifically designed for positive continuous outcomes with skewness, while it also naturally captures the pattern where variability grows with the mean. With a log link, the model ensures positive predictions, while each skill index changes performance by a constant factor regardless of whether someone starts low or high.

2.5.3. Model Evaluation and Comparison

In addition to reporting coefficient estimates and p-values, we evaluated model fitness using complementary fit criteria that are standard for GLMs and appropriate for our positive-valued performance outcome. Since the Gamma-Log GLM is estimated by maximum likelihood rather than least squares [33], likelihood-based indices (log-likelihood, Akaike information criterion (AIC) [34] and Cox–Snell pseudo- $[eqn]$ [35]) provide the most direct measures for illustrating goodness-of-fit for the model. To compare the predictive accuracy with Normal-Identity GLM, we also report the coefficient of determination ( $[eqn]$ ) and Root Mean Square Error (RMSE), together with mean absolute error (MAE), which yields an interpretable error magnitude that is less sensitive to extreme outliers than RMSE [36]. Finally, given the large inter-individual variability in early motor learning, we used leave-one-subject-out cross-validation (LOSO-CV) [37] to test whether the fitted model generalizes to unseen participants. Together, these metrics show a comprehensive assessment of both in-sample fit and cross-subject generalization.

3. Results

3.1. Skill Metrics over Days

Figure 5 and Table 1 summarize day-wise medians and interquartile ranges (IQR; 0.25–0.75) of juggling performance ( $[eqn]$ ) (Figure 5a) and the three standardized skill indices ( $[eqn]$ , $[eqn]$ , $[eqn]$ ) (Figure 5b–d) for each group. We present medians with IQR because they provide a robust summary of central tendency and variability and are less sensitive to outliers. While all metrics show increases from Day 1 to Day 10, the improvement differences in each group are different. Although group differences were not statistically significant, the Low-g group showed a trend of higher median values in juggling performance ( $[eqn]$ ) and Sequencing ( $[eqn]$ ) than the Normal-g group, while this trend was not observed in Prediction ( $[eqn]$ ) nor Accuracy ( $[eqn]$ ).

Juggling performance ( $[eqn]$ ) increased in both groups through days, with a higher median trajectory in the Low-g group (Figure 5a). The Low-g group improved from Days 1 to 10 (median of 3.5 → 17.8), whereas the Normal-g group also improved but reached a lower median (median of 3.0 → 9.75). Sequencing ( $[eqn]$ ) also increased over time in both groups, and the Low-g group consistently showed higher median scores than the Normal-g group throughout the whole period (Figure 5b). In the Low-g group, $[eqn]$ increased from $[eqn]$ (Day 1) to $[eqn]$ (Day 10), while in the Normal-g group, $[eqn]$ increased from $[eqn]$ (Day 1) to $[eqn]$ (Day 10). Prediction ( $[eqn]$ ) also increased across days in both groups and showed very limited group difference by Day 10 (Figure 5c). Low-g changed from $[eqn]$ (Day 1) to $[eqn]$ (Day 10), while Normal-g changed from $[eqn]$ (Day 1) to $[eqn]$ (Day 10). Accuracy ( $[eqn]$ ) also increased over days in both groups (Figure 5d). However, the Normal-g group showed higher $[eqn]$ than the Low-g group on all measured days. Low-g increased from $[eqn]$ (Day 1) to $[eqn]$ (Day 10), whereas Normal-g increased from $[eqn]$ (Day 1) to $[eqn]$ (Day 10).

While Low-g group showed higher $[eqn]$ and Sequencing at all measured days, group differences were not significant (all p > 0.35), with small-to-moderate effect sizes. Considering the small sample size and limited training days, power was limited by the sample amount and high inter-individual variability.

Apart from group-median performance, each individual participant’s highly variable trajectories of $[eqn]$ , $[eqn]$ , and $[eqn]$ across days are also presented in Figure 5e–j. Although the median robust z-scores increased from Days 1 to 10 in both groups, participants were highly different in baseline level and rate of change. Several individuals showed non-monotonic trajectories, including decreases across days, even within the Low-g group. This heterogeneity indicates that early juggling learning is strongly participant-specific, with a small sample number $[eqn]$ per group, a small number of fluctuating performers can substantially influence group-level summaries.

Before fitting the generalized linear models, we examined the relationships among the three standardized indices ( $[eqn]$ , $[eqn]$ , $[eqn]$ ). Pairwise Pearson correlations [38] were $[eqn]$ ( $[eqn]$ ), $[eqn]$ ( $[eqn]$ ), and $[eqn]$ ( $[eqn]$ ). To assess collinearity in the regression context, we also computed variance inflation factors (VIFs) [39] from an ordinary least squares model with $[eqn]$ as the outcome; VIF values were low (VI $[eqn]$ , VI $[eqn]$ , VI $[eqn]$ ), which is well below common thresholds for collinearity (e.g., 5–10), indicating no problematic multicollinearity among the three predictors.

3.2. Results of Computational Modeling

Table 2 reports the Gamma-Log GLM coefficients estimated from standardized predictors $[eqn]$ . All coefficients are positive, indicating that higher Prediction, Accuracy, and Sequencing are associated with better juggling performance. Sequencing shows the strongest effect and is highly significant ( $[eqn]$ , SE = 0.1036, $[eqn]$ ), followed by Prediction ( $[eqn]$ , SE = 0.1766, $[eqn]$ ). Accuracy has a smaller coefficient and is not significant in the model ( $[eqn]$ , SE = 0.1838, $[eqn]$ ), with a 95% CI that includes zero. Because the model uses a log link, $[eqn]$ gives the multiplicative change in expected performance in each predictor: $[eqn]$ and $[eqn]$ correspond to about 70% and 50% increases in expected performance, respectively, while $[eqn]$ corresponds to an 18% increase but without statistical support.

We compared two generalized linear models using the standardized skill indices $[eqn]$ pooled across Days 1, 5, and 10: (i) a Gamma GLM with a log link (Gamma-Log) and (ii) a Normal GLM with an identity link (Normal-Identity). Table 3 shows that the Gamma-Log model performed better across all evaluation criteria. Gamma-Log achieved higher explanatory power and lower prediction error ( $[eqn]$ , RMSE $[eqn]$ , MAE $[eqn]$ ) than Normal-Identity ( $[eqn]$ , RMSE $[eqn]$ , MAE $[eqn]$ ). Likelihood-based evidence also supported Gamma-Log GLM, with a higher log-likelihood ( $[eqn]$ vs. $[eqn]$ ), lower AIC (385.26 vs. 553.43), and higher Cox–Snell pseudo- $[eqn]$ (0.701 vs. 0.301). Leave-one-subject-out cross-validation confirmed better generalization of Gamma–Log (LOSO $[eqn]$ , RMSE $[eqn]$ , MAE $[eqn]$ ) compared with Normal-Identity (LOSO $[eqn]$ , RMSE $[eqn]$ , MAE $[eqn]$ ). These results illustrate that Gamma-log GLM is a better fit as the computational model of juggling performance.

Table 4 reports fit metrics computed separately for each day. On Day 1, both models report $[eqn]$ , indicating no explained variance at that early stage; still, the Gamma-Log model shows a smaller RMSE than the Normal-Identity model (5.17 vs. 15.10). On Day 5, the Gamma-Log model explains more variance ( $[eqn]$ ) and has lower error (RMSE = 17.90) than the Normal-Identity model ( $[eqn]$ , RMSE = 19.28). On Day 10, both models improve, with Gamma-Log again higher in $[eqn]$ (0.375 vs. 0.265) and lower in RMSE (39.36 vs. 42.69). Across days, $[eqn]$ increases from Days 1 to 10 for both models, while RMSE also increases, reflecting larger absolute performance values at later days. Overall, the Gamma-Log model provides better day-wise fit than the Normal-Identity model on Days 5 and 10.

Figure 6 summarizes the fit of the Gamma-Log model, with the Normal-Identity model shown as a reference check. Figure 6a shows the predicted versus observed performance for the Normal-Identity GLM. The wide spread around the identity line, especially at higher $[eqn]$ , indicates limited fit for this one-sided outcome. Figure 6b shows the same comparison for the Gamma-Log GLM. Predictions follow the identity line more closely within the dense low-to-mid performance region, which covers most observations, although several high-performance points still deviate from the model. Figure 6c plots residuals against predicted performance for both models. The non-outlier range was defined on observed $[eqn]$ by excluding MAD-based upper-tail outliers, i.e., trials with $[eqn]$ (cutoff $[eqn]$ ). Within this non-outlier range, Gamma-Log residuals show smaller robust bias and absolute deviation than the Normal model (Gamma-Log: $[eqn]$ , $[eqn]$ ; Normal: $[eqn]$ , $[eqn]$ ), consistent with the tighter, more zero-centered non-outlier ellipse in Sub-figure (c). Both models still show large positive residuals for a small number of high-performance cases, reflecting underestimation in the upper tail. Overall, the figure supports using the Gamma-Log GLM as a suitable computational description of early juggling performance, consistent with the pooled and by-day fit metrics in Table 3 and Table 4.

3.3. Sensitivity Analyses

As an exploratory extension of the main model, we added interaction terms between the standardized indices, as shown in Table 5. A Gamma-Log GLM including $[eqn]$ , $[eqn]$ , $[eqn]$ , and their pairwise interactions ( $[eqn]$ , $[eqn]$ , $[eqn]$ ), improved overall fit relative to the base model (AIC = 382.96 vs. 385.26; $[eqn]$ , RMSE = 20.21), while preserving the pattern that Sequencing and Prediction had positive coefficients and Accuracy was not significant. When we further included Group and Day as fixed effects together with the interaction terms, model fitness increased again (AIC = 382.25; $[eqn]$ , RMSE = 18.37). In this extended model, $[eqn]$ remained a robust positive variable, $[eqn]$ showed a positive trend ( $[eqn]$ ), and the $[eqn]$ interaction was significant ( $[eqn]$ ), suggesting that the effect of Prediction is stronger when Sequencing is high. Other interactions were not significant.

To partially account for participant-specific baselines, we also fitted an extended GLM that added fixed intercepts for each participant with the skill interactions, Group, and Day. This model achieved the best overall fit (AIC = 303.31; $[eqn]$ , RMSE = 10.18). Sequencing still remained a significant positive predictor ( $[eqn]$ , $[eqn]$ ), while the main-effect coefficients for Prediction and Accuracy became slightly negative ( $[eqn]$ , $[eqn]$ ; $[eqn]$ , $[eqn]$ ) despite significant positive interactions with Sequencing. Given the small sample size ( $[eqn]$ ) and large number of parameters in this model, we interpret these sign changes cautiously as likely reflecting overfitting and partialing among correlated terms. These sensitivity analyses support the robustness of the main Gamma-Log results: Sequencing shows the most stable association with performance, Prediction contributes in interaction with Sequencing, and Accuracy has only a weak independent effect in this dataset.

Finally, to assess whether excluding the three Day 1 statistical outliers ( $[eqn]$ ) introduced selection bias, we refitted the base Gamma-Log GLM on the full dataset ( $[eqn]$ ). For this robustness check, we included the actual calculated (but unstable) $[eqn]$ values for the two participants with one and four windows, and imputed the standardized median ( $[eqn]$ ) for the single participant with no valid windows. Coefficients for the core predictors changed only minimally (all shifts remained well within their respective standard errors) relative to the complete-case fit ( $[eqn]$ ): $[eqn]$ , $[eqn]$ , and $[eqn]$ . Consequently, the qualitative conclusions were unchanged (Sequencing and Prediction: positive; Accuracy: non-significant). Fit indices were also comparable (complete case: AIC = 385.26, RMSE = 25.77; full-dataset imputation: AIC = 396.99, RMSE = 25.39).

4. Discussion

This study built a component-based computational model for early juggling learning using three sub-skills: Sequencing (S), Prediction (P), and Accuracy (A). With pooled data from Days 1, 5, and 10, we explored our three hypotheses about how these components form juggling skill and relate to juggling performance. The main findings are as follows:

1. H1. The three indices showed different distributions and behaviors, supporting the idea that early juggling skill can be separated into Sequencing, Prediction, and Accuracy (Figure 4).
2. H2. Day-level juggling performance can be explained by the three components together. All components had positive effects on juggling proficiency $[eqn]$ , and Sequencing and Prediction showed the strongest independent contributions to it in the model (Table 2).
3. H3. The sub-skills did not increase at the same rate across days, showing lead–lag development. In line with our expectation, Sequencing rose steadily and remained higher in the Low-g group, suggesting that slow-tempo low-gravity VR training mainly supported early Sequencing gains (Figure 5).
4. Beyond these hypotheses, we observed large individual differences in learning trajectories, and a few very high-performance points were not well captured by the model (Figure 6).

4.1. Sub-Skill Contributions to Real-World Juggling Performance

In the Gamma-Log model, all three sub-skills had positive coefficients, meaning that higher S, P, and A were all linked to higher real-world juggling performance. Among them, Sequencing had the largest effect and was strongly significant (Table 2). This suggests that early performance mainly depends on building stable three-ball coordination over time, captured by rhythmic and phase-based measures. This agrees with the general idea that early jugglers first need to form a reliable cascade rhythm before they can sustain long runs.

Prediction was also significant. This suggests that the Prediction index, measured from apex-occlusion trials, captures a more general ability to estimate ball flight and landing without relying on continuous visual feedback. Participants who showed more reliable catching under occlusion tended to have higher three-ball performance even after controlling for S and A. Therefore, Prediction appears to be a real and measurable part of juggling skill that supports real-world performance rather than a task-specific artifact of the occlusion session.

Accuracy had a positive coefficient but was not significant in the computational model. This means that better spatial consistency of catch locations was related to higher $[eqn]$ , but its unique contribution was smaller than Sequencing and Prediction at this early stage. We interpret this primarily as a limitation of the current metric rather than evidence that spatial accuracy is unimportant for juggling. By design, A is estimated from second-catch clustering in a two-ball shutter-glasses task, which provides many controlled samples and isolates hand-centered spatial precision under reduced visual feedback, but does not capture all spatial demands of three-ball cascade juggling (e.g., full height and azimuth variability over long sequences). In the present dataset, the number of stable three-ball cycles with reliable 3D catch positions was too limited and heterogeneous across days to support a robust, parallel kinematic variability index, so a formal validation of A against three-ball dispersion was not feasible. Therefore, we treat A as a conservative occlusion-based proxy for spatial control in early learning.

Although Prediction and Accuracy were both derived from the visual-occlusion session, they did not collapse onto a single measure. Prediction ( $[eqn]$ ) and Accuracy ( $[eqn]$ ) showed a moderate correlation ( $[eqn]$ , $[eqn]$ ), but multicollinearity diagnostics remained low (all VIF $[eqn]$ ), suggesting that the two indices contribute separable information when modeled jointly. This aligns with their operational definitions: Prediction reflects the probability of successful catches under occlusion, whereas Accuracy captures the spatial clustering of second-catch locations.

As stated in our research Hypothesis H1, juggling skill can be decomposed into three components. The different distributions of S, P, and A and their positive contributions support this hypothesis, and that performance can be modeled as their joint function (Hypothesis H2).

4.2. Learning-Stage Dependence and Asymmetric Development

Day-wise model evaluation showed that fitness increased from Days 1 to 10 (Table 4). Both models explained almost no variance on Day 1, while fitness improved on Day 5 and further on Day 10, indicating that the three sub-skill indices relate to real-world performance more clearly as learners build a stable cascade. On Day 1, most trials were short and unstable, and some indices were close to baseline or missing, especially Sequencing when too few valid phase windows were available. With practice, the indices expanded in range and tracked $[eqn]$ more consistently.

The three components also changed at different rates across days. Sequencing increased steadily and remained higher in the Low-g group, suggesting that slow-tempo low-gravity VR training may preferentially support early Sequencing gains. Prediction improved strongly in both groups and became similar by Day 10, while Accuracy showed a different pattern and stayed higher in the Normal-g group (Figure 5). Given the modest sample size, these group patterns should be interpreted as trends rather than definitive group effects. These non-parallel trajectories support the lead–lag idea in H3, showing that early juggling learning involves multiple sub-skills that develop at different speeds, rather than a single skill improving uniformly.

4.3. Inter-Individual Variability and High-Performance Tail

Another clear outcome was strong variability between participants. Individuals differed in starting level and learning slope, and some showed non-monotonic changes. Baseline differences between participants were already large on Day 1, even though they were all beginners. Performance also fluctuated across days, which is common in motor learning, for example, due to fatigue, attention changes, or health conditions. With a small sample size (only $[eqn]$ ) per group and only three modeled days, these individual differences can largely influence group summaries, limit day-specific fit, and lower cross-subject prediction. The LOSO-CV result for Gamma-Log ( $[eqn]$ ) shows that the model can generalize to new participants, but also that a large part of performance variance is still person-specific.

Although the Gamma-Log GLM provided the appropriate description of the data overall, the residual plots (Figure 6c) still show systematic misfit for a few high-performance observations. One likely reason is the sparse upper tail of the performance distribution: only a small number of participants achieved very long cascades, so the GLM is dominated by the dense low-to-moderate range and shrinks extreme values toward the central trend. In addition, there may be practical ceiling effects due to the finite trial duration and experimental constraints. Finally, high performers may rely on qualitatively different control strategies (e.g., more automated rhythm control or different gaze patterns), which are not well captured by a single global log-linear relationship between $[eqn]$ and performance. Future work could explore alternative GLM families that allow more flexible variance structures, such as Tweedie [40] or inverse Gaussian distributions [32], but the current sample size is limited for reliably estimating and comparing such larger models. Therefore, we interpret the Gamma-Log model as a good description of typical early-stage performance while acknowledging reduced accuracy in the extreme high-performance tail.

4.4. Computational Model and Goodness-of-Fit

We used a Gamma-Log GLM as the main computational model because the performance outcome is strictly positive and right-skewed in early juggling. In this case, a Gamma distribution with a log link matches the scale of the data and naturally allows variance to increase with higher performance. The results showed that the Gamma-Log model captured the non-outlier performance range well and generalized reasonably across subjects in LOSO-CV (Table 3).

We also fitted a Normal-Identity GLM as a reference check. As expected for a one-sided outcome, its fit was weaker overall, and residual patterns were more biased than those of Gamma-Log (Figure 6). This study’s choice of the Gamma-Log model provides an appropriate statistical description of positive and heteroscedastic juggling performance, and this approach can be extended to other motor learning studies with similar outcome distributions.

Nonetheless, sensitivity analysis with pairwise interaction terms (Table 5) indicated a significant $[eqn]$ interaction when Group and Day were included, implying that the benefit of Prediction may be higher when Sequencing is already stable. Importantly, these extended models did not change the qualitative pattern of the main results: Sequencing and Prediction remained positive predictors of performance, whereas Accuracy showed only a small, non-significant contribution.

Because observations from the same participant were collected on multiple days, we also fitted an extended Gamma-Log GLM that included fixed intercepts for each participant in addition to Group, Day, and the skill interactions. In this model, Sequencing remained a significant positive predictor, while the main-effect coefficients for Prediction and Accuracy became slightly negative, despite positive interactions with Sequencing. Given the small sample size ( $[eqn]$ ) and the large number of parameters, we interpret these sign changes cautiously as likely reflecting overfitting and partialing among correlated terms, rather than a reversal of the underlying relationships. For this reason, we retain the simple Gamma-Log GLM as the primary model and use the interaction effects and participant-intercept models only as complementary sensitivity checks.

In addition to these model-form robustness checks, we also examined whether excluding the Day 1 statistical outliers ( $[eqn]$ ) introduced selection bias. We refitted the base Gamma-Log GLM on the full dataset ( $[eqn]$ ) by retaining the actual $[eqn]$ values for the cases with 1 and 4 windows despite their instability, and imputing the standardized median ( $[eqn]$ ) for the case with no valid windows. This produced only minimal changes in the core coefficients ( $[eqn]$ , $[eqn]$ , $[eqn]$ ) and similar fit indices, indicating that the main inferences are not driven by the Day 1 exclusions.

4.5. Limitations

This study has several limitations. First, the sample size ( $[eqn]$ ) was relatively small for evaluating the training intervention. Although the Low-g group exhibited a trend of higher performance and Sequencing scores, the observed effect sizes were small to moderate (Hedges’ g ≈ 0.26–0.40 for $[eqn]$ ). A post hoc power analysis indicated that a sample size of over 100 participants (approx. 139) would be required to definitively confirm the efficacy of Low-g training with sufficient statistical power. However, as the primary focus of this study was to decompose juggling skills into computational components rather than to validate the training effect itself, we interpret the group-level trends descriptively.

Second, the main Gamma-Log GLM treated Participant–Day samples as independent, although each participant contributed multiple days. This may induce correlated residuals; an exploratory mixed-effects model with participant-level random intercepts showed similar patterns but was not adopted as the primary model, and future work with larger samples should use mixed-effects Gamma models more systematically. Second, both tested GLMs were additive and still underestimated a few very-high-performance cases, suggesting that models with nonlinear or interaction terms, or alternative positive distributions, could better capture the upper tail. Third, Accuracy was measured only by ellipsoid volume from the two-ball occlusion task. Spatial accuracy in three-ball juggling also includes throw height and direction consistency, which may become more important later.

4.6. Future Work

Future studies should extend the model to all training days and increase the sample size so that more complex structures can be estimated reliably. In particular, mixed-effects Gamma models with participant-specific random intercepts could more reliably model the repeated-measures structure directly and test whether component weights change across learning stages (e.g., Day $[eqn]$ , Day $[eqn]$ , Day $[eqn]$ ) with sufficient data. In addition, the underestimation of very high performers suggests that nonlinear terms or interaction effects among the components (e.g., $[eqn]$ ) may capture residual variance that the current simple GLM cannot explain.

The current Accuracy index captures only part of the spatial precision required in full three-ball cascade juggling because it is based on two-ball occlusion catches. Future work should derive complementary Accuracy measures directly from three-ball kinematics (e.g., variability of catch height and coordinates, horizontal error distributions) and test how these measures relate to $[eqn]$ and to the present two-ball-based index A. This will allow a more complete assessment of spatial control as a component of juggling skill.

Beyond the behavioral and kinematic results presented here, the 10-day multimodal evaluation system also included EEG, EMG, and pre-/post-fMRI recordings [23]. In future work, these modalities could be integrated with the three-component model to relate Sequencing, Prediction, and Accuracy to patterns of brain activity, functional connectivity, and muscle activation. For example, EEG and fMRI could be used to examine how training-related changes in cortical dynamics covary with improvements in each sub-skill, and whether low-gravity VR training preferentially modifies networks related to temporal coordination, while EMG could quantify how motor efficiency and muscle usage evolve with practice. This multimodal extension would provide a more complete picture of how complex motor skills such as juggling are acquired through training in both brain and body.

Although the present work focuses on three-ball cascade juggling, the proposed three-component structure is likely relevant for other complex motor skills. In racket and ball sports (e.g., tennis, table tennis, basketball), athletes must learn stable action sequences and rhythms (Sequencing), predict ball trajectories and player movements under uncertainty (Prediction), and place the racket, hand, or foot at precise spatial targets (Accuracy). Similarly, in musical performance (e.g., piano and drumming), learners acquire ordered finger and limb patterns (Sequencing), anticipate upcoming notes and phrases (Prediction), and control spatial endpoints such as key locations or hand positions (Accuracy). These examples suggest that a Sequencing–Prediction–Accuracy decomposition may provide a general template for analyzing how multiple sub-skills support proficiency in real-world motor tasks beyond juggling. Our results show that these elements can be quantified separately and may develop at different speeds in early learning. Thus, the present framework may provide a general starting point for decomposing and modeling other complex motor skills in future work.

5. Conclusions

This study introduces a longitudinal learning study focusing on the behavioral and kinematic results from our larger multimodal juggling project. In conclusion, our computational model using Gamma-Log GLM shows that early juggling performance can be described as a combined function of Sequencing, Prediction, and Accuracy. Among them, Sequencing has the strongest positive associations with juggling performance, followed by Prediction, and then Accuracy. Across days, the three components followed different learning trajectories, and Sequencing increased most clearly in the low-gravity group, consistent with the expectation that slow-tempo low-gravity VR practice may potentially promote early rhythmic coordination. Overall, this study provides a structured, multi-component view of early juggling acquisition and offers a clear computational framework for interpreting motor learning in complex real-world tasks.