Unlike our prior experimental design that relied on explicit geometric registration, this study utilized a Hassani-style background pattern [15] to enable direct comparison with previous paradigms. The background pattern was printed on a paper sheet measuring 54 × 60 cm (all dimensions specified as width × height throughout this paper). Its central 43 × 30 cm region formed a 3 × 3 matrix of nine rectangular color patches, each subtending a visual angle of approximately 8.1∘×5.7∘. The two outer columns contained chromatic patches that served as the target backgrounds, whereas the central column contained three uniform black patches that merged into a seamless black strip, serving as the test backgrounds for the matching task. From top to bottom, the left column comprised blue, gray, and white patches while the right column comprised green, brown, and red patches. The peripheral region surrounding the matrix was filled with a neutral gray background though its visibility was heavily restricted by the aperture.

The background booth accommodated two ambient illumination levels: a high-luminance condition and a low-luminance condition, providing ambient luminances of approximately 31.8 cd/m2 and 7.06 cd/m2 (corresponding to illuminances of 100 lx and 22.2 lx), respectively. The illumination reproduced Hassani’s “cool light”, where the nominal CCT was 4334 K. The spectral power distributions of the light sources were optimized along with the background pattern to reproduce the colorimetric profiles of the background patches used in Hassani’s study although minor chromatic errors remained. The detailed colorimetric coordinates of the background patches (along with their corresponding target mixtures) under both illuminance conditions are reported in Table I.

The foreground display operated at a 10-bit color depth with a native resolution of 2560 × 1440 pixels and a pixel pitch of 0.2331 mm, subtending approximately 0.8 arcmin per pixel. The display was colorimetrically characterized to enable rapid and accurate transformations between CIEXYZ tristimulus values and digital RGB values with minimal errors. To align with the two illuminance conditions, the peak display luminance was manually toggled via the monitor’s physical hardware controls to either 300 cd/m2 or 75 cd/m2, with the display characterization model updated to reflect the corresponding luminance setting. This adjustment maximized the display’s usable dynamic range under both illuminance conditions, preventing quantization errors and ensuring finer control during color matching.

Across all experiments, foreground targets and background patches were generally paired one to one using coordinates from the previous literature [15]. In Experiment 1, two backgrounds were additionally paired with a second foreground to investigate background chromatic induction. The target colors were defined by the tristimulus values of the additive mixture, representing the physical sum of the background and foreground (XY Zmixture = XY ZFG + XY ZBG). During the trials, the control program calculated the required digital RGB values based on the measured background tristimulus values and the characterization profile scaled to the selected luminance level. This generated a virtual foreground target stimulus that was spatially nested within the boundary of the background patch. Concurrently, an identical-sized adjustable matching stimulus was displayed on the adjacent black patch in the central column.

The psychophysical matching program was implemented in MATLAB using the Psychophysics Toolbox Version 3 [20]. Observers independently manipulated chroma, hue, and lightness within the CIELCH color space. The program provided interactive functionalities allowing observers to toggle between the current matching state and the initial starting color coordinates, to save and recall up to three intermediate matching states, and to confirm the final match to advance the trial.

The white patch located in the lower-left corner of the matrix was initially selected as the reference white point for the CIELCH transformation, but it was noticeably yellowish. This caused the color trajectory during chroma adjustments to deviate from subjective neutral white, introducing unintended hue shifts and increasing task difficulty. To address this, a white-point adjustment mode was introduced to customize the color space transformation for each individual and to quantify the observer’s actual chromatic adaptation state under the current ambient conditions. Subsequent statistical analysis confirmed that individual white-point adjustment did not exert a significant effect on the final matching coordinates (see Appendix A.5, Supplementary Material). Full procedural details and session coverage are provided in Appendix A.1 and Table B1 (Supplementary Material).

Matching precision and overall perceptual bias were quantified using mean Euclidean distance from the mean (MEDM) and mean Euclidean distance from the target (MEDT) within each subgroup in CIECAM16-UCS: where xi represents the matched coordinate vector of trial i, x¯ is the condition-wise arithmetic mean, and xtarget denotes the physical target coordinate vector.

The cleaned dataset was analyzed in R (version 4.5.2) [28] using the glmmTMB package [4]. A linear mixed-effects model [25] provides a regression-based framework capable of simultaneous modeling of fixed and random effects, making it highly effective for analyzing unbalanced data structures containing repeated measures while controlling for random variations across individual observers. Extending this framework, GLMM [3] directly accommodates high-dimensional datasets and permits conditional response distributions to follow non-normal profiles while explicitly supporting the configuration of heteroscedastic structures. Although less conventional in standard color science, this statistical approach is well suited for the complexities of our psychophysics tasks. Specifically, it allows us to systematically control for interobserver adjustment preferences as random effects while analyzing repeated-measures data. Furthermore, unlike the traditional analysis of variance, the GLMM framework gracefully handles the unbalanced data structures resulting from our multifactorial design and enables the configuration of heteroscedastic error distributions directly aligned with human color-matching volatility.

The statistical analysis focused on seven experimental factors—viewing mode, stimulus size, aperture presence, foreground spatial offset, starting point, ambient luminance, and presentation medium—along with two specific interactions: between the viewing mode and the stimulus size, and between the starting point and the presentation medium. The response variable was defined as the colorimetric difference vector (Diff), calculated as the matched color minus the intended target color within the CIECAM16-UCS space. To model the three-dimensional color coordinates, the data was converted into a long format, introducing a new fixed factor, coordinate axis (axis). All explanatory variables, including both fixed and random effects, were treated as categorical factors. The random-effects structure with observers as the grouping variable was specified as (0 + axis | observer) to account for interobserver variability across the three perceptual dimensions. While structurally specified as random slopes, this configuration functionally assigns distinct but correlated random intercepts to each axis.

The dispersion model was specified as ∼axis * target_id, which allowed conditional residual variances to be heteroscedastically distributed across different target colors and coordinate axes. The conditional response distribution was specified using t_family() to minimize the influence of extreme values and enhance the structural robustness of the model estimates against influential outliers and heavy-tailed residual distributions.

To establish a unified and directly comparable hierarchy of all experimental factors, our modeling strategy adopted a two-stage approach: we first screened for significant predictors within individual experiments and subsequently consolidated these key factors into a single pooled model.

In the initial screening stage, model selection followed a top-down refinement process starting from a full multivariate GLMM within each of the four experiments. Candidate models incorporating experimental factors of interest, their mutual interactions, and their interactions with the control variables (coordinate axis and target combination) were evaluated using standard maximum likelihood estimation, and the model yielding the minimum corrected Akaike Information Criterion (AICc) [29] was selected as the optimal structure, with lower values indicating a superior balance between model fit and parsimony. During this stage, candidate dispersion formulas were also evaluated, with ∼axis * target_id selected as it maximized the reproduction of empirical variance characteristics while maintaining strict model convergence.

In the second stage, the significant predictors identified from individual experiments were consolidated into a unified formula and applied to the pooled dataset. Top-down AICc comparisons supported retaining (0 + axis | observer) as the optimal random-effects structure without an interexperiment random intercept. Secondary procedural factors—specifically white-point adjustment status (wp_adjusted) and observer awareness of the display medium (disp_aware)—were evaluated and excluded due to nonsignificant contributions, confirming that perceptual matching was robust against white-point adjustment and medium awareness (see Appendix A.5, Supplementary Material for selection details).

Final parameter estimates for the consolidated pooled model were computed using restricted maximum likelihood to ensure unbiasedness and robustness [25]. Pooling the data across experiments provides two methodological advantages: first, integrating conditions with shared factor levels maximizes the available control data, thereby enhancing statistical power and estimation stability; second, a single consolidated model enables a direct, quantitative comparison of the relative effect sizes across all significant factors within a unified framework.

Model adequacy and distributional assumptions were first evaluated through standard diagnostic procedures. Residual distribution adequacy, conditional dispersion, and potential outliers were verified using DHARMa-based diagnostics. Multicollinearity among predictors was assessed using variance inflation factors. Model precision and variance partitioning were additionally summarized using mean absolute errors (MAEs) across coordinate axes, observer random-effects standard deviations, and the interclass correlation coefficient (ICC).

Upon confirming model validity, the statistical contribution of experimental factors in the final pooled model was evaluated using Type III Wald χ2 tests. To interpret significant main effects, estimated marginal means (EMMs) [21] were computed for all unique combinations of the significant factors. To quantify the direction and the magnitude of each factor’s effect, pairwise contrasts were additionally computed under standardized reference conditions: for each factor, the levels being compared were evaluated while holding the remaining factors at their default baseline levels.

Among all primary predictors, stimulus size exerted the most pronounced influence on matching performance, with every model term involving stimulus size attaining statistical significance (all p < 0.001, Table III). This spanned main effects, lower-order interactions, and notably a highly significant three-way interaction with coordinate axis and target combinations (axis × stim_size × target_id, χ2 = 227.06, p < 0.001). This effect is visually evident in Fig. 3, where the difference between large (Fig. 3a) and small (Fig. 3b) stimulus sizes is substantial and varies across targets.

Across the four evaluated target combinations (Targets 3, 7, 13, and 19), pairwise contrast analysis between the small and large conditions yielded overall Euclidean differences ranging from 2.22 to 6.80 (Table V and Appendix Table B5.2, Supplementary Material). As is also visually apparent in Fig. 3, pairwise contrast analysis along the chromatic plane (a ′–b ′) revealed a systematic spatial rule: for all coordinate shifts that attained statistical significance (p < 0.05), the directional sign of the “small minus large” contrast was consistently opposite to the sign of the corresponding background chromaticity coordinate. For Target 3 (BG: green, FG: green; a ′ = −24.56, b ′ = 14.41 under the reference high luminance), reducing stimulus size drove significant shifts toward positive a ′ (+2.04) and negative b ′ (−3.43). Similarly, Target 7 (BG: white, FG: blue; a ′ = 0.35, b ′ = 10.01) exhibited a significant negative shift along b ′ (−2.19) opposing its positive background b ′ coordinate, accompanied by a sharp drop in lightness (ΔJ  ′ = −6.40), whereas Target 19 (BG: red, FG: green; a ′ = 16.26, b ′ = 12.56) shifted toward negative a ′ (−1.17). Nonsignificant chromatic shifts, such as those in Target 13 (BG: green, FG: red; p > 0.05), remained within the margin of statistical uncertainty.

The effect of ambient luminance was evaluated across six target combinations (Targets 1, 3, 4, 6, 7, and 9) by comparing high (31.8 cd/m2) and low (7.06 cd/m2) illumination environments. In the pooled model, although the main effect ambient_lum alone was not statistically significant (p = 0.133), it demonstrated highly significant interactions with target combination (ambient_lum × target_id, χ2 = 17.64, p = 0.003) and coordinate axis (axis × ambient_lum × target_id, χ2 = 121.65, p < 0.001; Table III), confirming that its effect is strongly target- and axis-dependent.

Pairwise contrast analysis between high and low ambient luminance conditions—as visualized in Figure 4—revealed that statistically significant perceptual shifts were concentrated in Targets 4 (represented in gray) and 7 (represented in black) (Table V and Appendix Table B5.2, Supplementary Material). Substantial Euclidean distances occurred for Target 4 (BG: gray, FG: purple; ΔE  ′ = 3.61) and Target 7 (BG: white, FG: blue; ΔE  ′ = 3.21). For Target 4, high ambient luminance drove significant shifts along a ′ (−0.99, p = 0.028) and b ′ (+2.70, p = 0.001). Similarly, Target 7 exhibited significant chromatic shifts along both a ′ (−0.86, p < 0.001) and b ′ (+3.09, p < 0.001). This statistical difference is visually mirrored in Fig. 4, where the confidence intervals for Target 7 along the b ′ axis are widely separated with no overlap. Target 9 (BG: red, FG: orange) showed a statistically significant shift along a ′ (+0.79, p = 0.023) though its overall Euclidean shift remained small (ΔE  ′ = 0.85). Other targets showed no significant shifts (p > 0.05). Across all six evaluated targets, lightness contrasts (ΔJ  ′) failed to attain statistical significance (p > 0.05; e.g., ΔJ  ′ = 2.18, p = 0.142 for Target 4). Similar factor-wise comparisons for other experimental factors, as well as forest plots across the 11 unique conditions, are provided in Appendix Figures C2 and C3 (Supplementary Material).

The remaining primary predictors—starting point, viewing mode, and presentation medium—exhibited statistically significant but numerically small effects on perceived color matching. Interactions between these three factors and target combination were excluded during model selection as they did not improve the model fit, leaving their contributions captured via main effects and axis interactions (axis × start_point, χ2 = 51.02, p < 0.001; axis × view_mode, χ2 = 14.17, p < 0.001; axis × medium, χ2 = 8.99, p = 0.011; Table III and Appendix Table B4, Supplementary Material).

Pairwise contrast analysis confirmed that overall Euclidean color shifts for all three factors remained small (ΔE  ′ < 1.00; Table V and Appendix Table B5.2, Supplementary Material). For the starting point, evaluating matches between mid and low initial color coordinates yielded a small shift (ΔE  ′ = 0.996) that attained statistical significance across all three coordinate axes (Table V). In contrast, the effects of the viewing mode and presentation medium were smaller (ΔE  ′ = 0.69 and 0.56, respectively) and specifically confined to the blue–yellow axis (Δb ′ = −0.41, p < 0.001 for binocular vs monocular viewing; Δb ′ = +0.41, p < 0.001 for AR vs display-only presentation; Table V).