avoid truncation effects. Though somewhat unrealistic, specification A serves
to separate the impact of modeling decisions or situations with more
(realistic) uncertainty for the user from the overall performance of the
multiFAMM.

For model specification B (truncation via total variation (TV)), the
truncation orders of the FPCs are chosen so that $95\%$ of the total variation
(6) are explained. In scenario C (truncation via univariate variation (UV)),
we choose the number of FPCs so that on every dimension $95\%$ of univariate
variance (7) are explained. Specifications D (TV with alternate scalar
product) and E (UV with alternate scalar product) use the cut-off criterions
analogous to B and C but we alternate the scalar product on which the MFPCA is
based: For data generated from a model based on an unweighted MFPCA, the
scalar product used in these scenarios is weighted by
$\frac{1}{\hat{\sigma}^{2}_{d}}$, and vice versa. Model specification F
(scedastic misspecification) evaluates misspecifying the multiFAMM with a
homoscedasticity assumption. We also contrast the multiFAMM with the
univariate approach of fitting independent univariate models. In modeling
scenarios denoted with U, we fit an independent FAMM on each dimension. We use
the FAMMs proposed by Cederbaum et al. (2016) so that we can apply the same
model specifications as for the multiFAMM (e.g. basis functions, penalties,
etc.). The number of FPCs in the model is then chosen so that $95\%$ of
univariate variation is explained.

Table 8: Description of model specifications.

Model Specification |  Description
---|---
A |  Univariate mean, multivariate mean, univariate auto-covariance estimation based on cubic P-splines with 8, 8, 5 (marginal) basis functions and choice of penalty as in the model used for generating the data setting (third order for 1-5, first order for 6). No univariate truncation for the MFPCA except negative eigenvalues. MFPCA based on scalar product as in the model used for generating the data setting (weights of one except for 5). Number of multivariate FPCs is fixed according to the data setting (1-5: ($\bm{B}$:3, $\bm{E}$:5), 6: ($\bm{B}$:6, $\bm{C}$:5, $\bm{E}$:5)). Weighted regression approach with weights obtained from univariate variance decompositions for heteroscedastic data settings (1-5).
B |  Number of FPCs chosen so that $95\%$ of the (weighted) sum of total variation in the data is explained. Rest as in A.
C |  Number of FPCs chosen so that on each dimension at least $95\%$ of the univariate variation in the data is explained. Rest as in A.
D |  Weighted scalar product for MFPCA with weights based on dimension specific measurement error variances if data setting is based on model using a scalar product with weights of one (all except 5) and vice versa (5). Number of FPCs chosen so that $95\%$ of the (weighted) sum of total variation in the data is explained. Rest as in A.
E |  Weighted scalar product for MFPCA with weights based on dimension specific measurement error variances if data setting is based on model using a scalar product with weights of one (all except 5) and vice versa (5). Number of FPCs chosen so that on each dimension at least $95\%$ of the univariate variation in the data is explained. Rest as in A.
F |  Homoscedasticity is assumed, no regression weights. Rest as in A.
U |  Univariate modeling approach. Independent sparse FLMM models for each dimension. Number of univariate FPCs chosen so that $95\%$ of univariate dimension is explained. Rest as in A.

#### Modeling Scenarios

For each of the combinations of data settings and modeling scenarios indicated
in Table 9, analyses were performed giving $500$ fitted models per
combination. The results presented in the main part correspond to modeling
scenarios 1A (Benchmark), 1B (Cut-Off Multi), 1C (Cut-Off Uni), 3A (Sparse
Data), and 4A (Uncentred Scores).

Table 9: Overview of the analysed modeling scenarios. |  | Data Setting
---|---|---
Model Specification | 1 | 2 | 3 | 4 | 5 | 6
| A | X | X | X | X | X | X
| B | X | X |  |  | X |
| C | X | X |  |  | X |
| D | X | X |  |  | X |
| E | X | X |  |  | X |
| F | X | X |  |  |  |
| U | X | X |  |  | X |

#### Evaluation Criteria

The accuracy of the estimated model components is measured using the rrMSE. We
analyse the accuracy of the covariance estimation, i.e. eigenfunctions,
eigenvalues, and measurement error variances as well as the mean estimation.
For evaluating the overall accuracy of the estimate $\hat{\bm{\zeta}}$ of the
multivariate functional component
$\bm{\zeta}=(\bm{\zeta}_{1},...,\bm{\zeta}_{S})^{\top}$, we define the mrrMSE
as

$\displaystyle\mathrm{mrrMSE}(\bm{\zeta},\hat{\bm{\zeta}})=\sqrt{\frac{\frac{1}{S}\sum_{s=1}^{S}|||\bm{\zeta}_{s}-\hat{\bm{\zeta}}_{s}|||^{2}}{\frac{1}{S}\sum_{s=1}^{S}|||\bm{\zeta}_{s}|||^{2}}}$

with the multivariate norm based on the unweighted scalar product. Note that
$S$ depends on the model component to be evaluated, e.g. for the fitted values
$\bm{y}_{ijh}(t)$, $S$ equals $N$ but for the subject-specific random effect
$\bm{B}_{i}(t)$, $S$ equals the number of individuals. The mrrMSE corresponds
to the evaluation criterion used in section 5 of the main part. We also define
a urrMSE

$\displaystyle\mathrm{urrMSE}\big{(}\bm{\zeta}^{(d)},\hat{\bm{\zeta}}^{(d)}\big{)}=\sqrt{\frac{\frac{1}{S}\sum_{s=1}^{S}||\zeta_{s}^{(d)}-\hat{\zeta}_{s}^{(d)}||^{2}}{\frac{1}{S}\sum_{s=1}^{S}||\zeta_{s}^{(d)}||^{2}}}$

for $\bm{\zeta}^{(d)}=(\zeta_{1}^{(d)},...,\zeta_{S}^{(d)})^{\top}$, which
allows to evaluate the dimension specific estimation accuracy as well as a
straightforward comparison to the univariate modeling approach. For scalar
estimates such as eigenvalues and error variances, we define the rrMSE as

$\mathrm{rrMSE}(\zeta,\hat{\zeta})=\sqrt{\frac{(\zeta-\hat{\zeta})^{2}}{\zeta^{2}}}$

with estimate $\hat{\zeta}$ of the scalar value $\zeta$. The rrMSE takes on
(unbounded) positive values with smaller values indicating a better fit. As
the rrMSE is a relative measure, small differences between estimate and true
component can result in large rrMSE values for true component norms close to
zero. Note that eigenfunctions are defined only up to a sign change. We thus
flip the estimated eigenfunction (multiply it by $(-1)$) if this results in a
smaller norm for the difference between the true function and its estimate.

Additionally, we evaluate the coverage of the point-wise CBs of the estimated
fixed effects.

###  D.2 Results of the Simulation Study

#### Impact of Model Specifications

The results for setting 1 (consonant assimilation data) demonstrate the
importance of the number of FPCs in the accuracy of the estimation (see Figure
22 and Table 10). With specifications A (true model) and F (scedastic
misspecification), the sets of FPCs are fixed giving overall low values for
the mrrMSE (the misspecification in F yields only slightly worse results).
Similarly, choosing the truncation order via the proportion of univariate
variance explained in scenario C gives models with roughly the same number of
FPCs as is used for the data generation. The cut-off criterion based on the
total amount of variance in scenario B results in more parsimonious models and
thus considerably higher mrrMSE values. The number of selected FPCs also
explains the wider boxplots of the scenarios B-E compared to A: rather than a
larger overall variance in estimation, we find separate clusters based on the
included FPC sets (see also Figure 23). For the modeling scenarios based on a
weighted scalar product (1D and 1E), the number of chosen FPCs is quite
similar regardless of the cut-off criterion but the overall mrrMSE values of
the fitted curves ($\bm{y}_{ijh}(t)$) are higher than for the unweighted
approach. The estimate of the mean $\bm{\mu}(\bm{x}_{j},t)$, however, is
comparatively stable over the different model specifications.

Figure 22: mrrMSE values of the fitted curves $\bm{y}_{ijh}(t)$, the mean
$\bm{\mu}(\bm{x}_{j},t)$, and the random effects $\bm{B}_{i}(t)$ and
$\bm{E}_{ijh}(t)$ for all modeling scenarios in data setting 1 (blue shades)
and the standard modeling scenario A in data settings 2 to 5.

| 1A | 1B | 1C | 1D | 1E | 1U | 2B | 2C | 2D | 2E | 2U
---|---|---|---|---|---|---|---|---|---|---|---
| dim1 | dim2 | dim1 | dim2
$\bm{B}_{i}(t)$ | 3.00 | 2.15 | 2.80 | 2.62 | 2.88 | 3.09 | 2.02 | 2.87 | 2.00 | 3.00 | 3.06
2.00 | 1.09 | 2.00 | 1.06
$\bm{E}_{ijh}(t)$ | 5.00 | 4.00 | 5.00 | 4.01 | 4.12 | 4.99 | 3.77 | 4.44 | 3.00 | 4.27 | 4.72
2.00 | 2.99 | 2.00 | 2.72

Table 10: Average number of eigenfunctions selected (scenario A fixed to
underlying truth) in the 500 simulation iterations per modeling scenario
(column) and random effect (row). For the univariate modeling approach
(scenario U) we report the total amount (top) and the number per independent
model (bottom).

Figure 23 shows the urrMSE values of the fitted curves $\bm{y}_{ijh}(t)$ for
different modeling scenarios (1U, 1B, 1C, 1D, 1E) depending on the number of
FPCs included in the model. For example in the univariate modeling scenario
1U, all 500 models choose two FPCs for each random effect (B2-E2) on dimension
$dim1$. On $dim2$, however, different combinations of FPCs lead to
considerably different urrMSE values. For the multivariate modeling scenarios,
we also find that the reduction in rrMSE values depends on which additional
FPC is included: $\bm{\psi}_{E5}$ reduces the values more than
$\bm{\psi}_{B3}$.

Figure 23: urrMSE values of the fitted values $\bm{y}_{ijh}(t)$ for different
scenarios. The urrMSE values are grouped by the number of FPCs included in the
models. B$x$-E$y$ means the model contains $x$ FPCs for $\bm{B}_{i}(t)$ and
$y$ for $\bm{E}_{ijh}(t)$.

Figure 24 shows the mrrMSE values of the fitted curves $\bm{y}_{ijh}(t)$, the
mean $\bm{\mu}(\bm{x}_{i},t)$, and the random effects $\bm{B}_{i}(t)$ and
$\bm{E}_{ijh}(t)$ for the different modeling scenarios of data settings 1, 2,
and 5. Again, we find considerable differences for different numbers of FPCs
in the models when there is strong heteroscedasticity between the different
dimensions (setting 2). The overall model fit seems to be better for applying
an unweighted scalar product for the MFPCA (2B, 2C compared to 2D, 2E). Later,
we will see that the model fit on single dimensions can be improved by
weighting the scalar product. Note that for setting 2F, misspecifying the
model assumption now has a larger negative impact on the fitted curves than in
setting 1F. In data setting 5, the data are generated based on the model using
a weighted scalar product. Then, the number of FPCs in modeling scenarios 5B
and 5C is very similar. Interestingly, basing the MFPCA on a scalar product
using weights of one can lead to a similar estimation accuracy than the
standard modeling scenario (5A compared to 5E).

Figure 24: mrrMSE values of the fitted curves $\bm{y}_{ijh}(t)$, the mean
$\bm{\mu}(\bm{x}_{i},t)$, and the random effects $\bm{B}_{i}(t)$ and
$\bm{E}_{ijh}(t)$ for different scenarios.

Our simulation study thus suggests that basing the truncation orders on the
proportion of explained variation on each dimension gives parsimonious and
well fitting models. If interest lies mainly on the estimation of fixed
effects, the alternative cut-off criterion based on the total variation in the
data allows even more parsimonious models. Furthermore, an unweighted scalar
product is a reasonable starting point for the multiFAMM.

#### Model Performance on Different Data Settings

To compare different data settings, we focus on model specification A (true
model) in Figure 22. Note that the mrrMSE values of the other data settings
cannot be directly compared as the denominator of the mrrMSE (slightly)
changes (except for the mean of 2A, 3A, and 4A). We find that strong
heteroscedasticity (2A) mainly negatively affects the fitted curves and the
smooth residual. Unsurprisingly, the results for scenario 3A suggest that the
estimation accuracy is lower for sparse functional data. However, the mean
estimation is comparable to more densely observed data. On the other hand, the
mean estimation is susceptible to violations of the assumption of uncorrelated
and centred realizations of the random effects (4A). The model then has
difficulties to separate intercept and random effects (see also Figure 31),
which does not necessarily translate to a worse overall fit to the data. The
results of scenario 5A (weighted scalar product) suggest that the accuracy of
the multiFAMM does not depend on the definition of the scalar product used for
the MFPCA.

Figure 25 shows the urrMSE values for the aforementioned model components of
modeling scenario 6A (snooker training data) where we compare the estimation
accuracy across the dimensions. The fit of the functional curves
$\bm{y}_{ijh}(t)$ suggests that there might be pronounced differences between
the estimation accuracy of the dimensions, the dimensions $dim1$
(corresponding to $elbow.x$) and $dim4$ (corresponding to $hand.y$) giving
high rrMSE values.

Figure 25: urrMSE values of the fitted curves $\bm{y}_{ijh}(t)$, the mean
$\bm{\mu}(\bm{x}_{i},t)$, and the random effects
$\bm{B}_{i}(t),\bm{C}_{ij}(t)$ and $\bm{E}_{ijh}(t)$ for modeling scenario 6A.
One outlying observation with high urrMSE value was removed from the boxplot
for $\bm{E}_{ijh}(t)$ and $dim1$.

We conclude that the proposed multiFAMM performs well even in more challenging
data settings such as sparse functional data or data with few grouping levels
for the random effects. Especially the estimation of the fixed effects seems
to be stable over the different analysed settings.

#### Comparison to Univariate Approach

We compare the univariate modeling approach U to multivariate modeling
scenarios C (UV) and E (UV with alternate scalar product), i.e. we make sure
that the proportion of explained univariate variation is also at least $95\%$
for each model. Table 10 shows that in data settings 1 (consonant assimilation
data) and 2 (strong heteroscedasticity) the number of included FPCs tends to
be higher for scenario U. Yet Figure 26 indicates that the multiFAMM yields
consistently lower urrMSE values on $dim1$ (smaller measurement error
variance). For $dim2$, the random effects seem to be estimated more accurately
and the fixed effects similarly well, but especially with the weighted scalar
product (scenarios E), the overall fit for $y$ can give higher urrMSE values.
Overall, the results suggest that the unweighted scalar product is to be
preferred in data situations with similar measurement error variances of the
dimensions such as setting 1, where it gives reliably good results across all
model components. However, downweighting the dimension $dim2$ with larger
error seems to be a reasonable modeling decision in setting 2 if interest lies
primarily on $dim1$ (lowest urrMSE values for 2E). Again, we point out that
the estimation of the fixed effects is relatively stable across approaches and
slightly better than for the univariate modeling approach.

Figure 26: urrMSE values of the fitted curves $\bm{y}_{ijh}(t)$, the mean
$\bm{\mu}(\bm{x}_{i},t)$, and the random effects $\bm{B}_{i}(t)$ and
$\bm{E}_{ijh}(t)$ for the univariate modeling scenario U (grey) and the
multivariate modeling scenarios C and E in data settings 1 (blue) and 2
(green).

Table 11 compares the standard errors of the fixed effects estimation for the
multiFAMM in scenario 1C with the univariate modeling approach of scenario 1U.
We look at the ratio of standard errors $\frac{se_{m}}{se_{u}}$, where
$se_{m}$ denotes the standard error of the multiFAMM and $se_{u}$ the standard
error of the corresponding univariate FAMM. For each dimension and partial
predictor, we calculate the proportion of ratios smaller and larger than one
over an equidistant grid of 100 time points and all simulation runs.
Especially on $dim1$, we find that the multivariate approach gives smaller
standard errors. On $dim2$, the proportions are more similar with slightly
more proportions larger than one. Overall, we find more ratios smaller than
one in our simulation study thus indicating that the multivariate approach can
yield smaller standard errors. Table 12 further shows that the coverage of
these two scenarios is very similar (see subsection on the estimation of fixed
effects).

Table 11: Proportion of standard error ratios of the multiFAMM ($se_{m}$) of
scenario 1C and the univariate modeling approach ($se_{u}$) of scenario 1U
over all simulation runs and all evaluated time points.

$d$ | Ratio | $\beta_{0}^{(d)}$ | $f_{0}^{(d)}(t)$ | $f_{1}^{(d)}(t)$ | $f_{2}^{(d)}(t)$ | $f_{3}^{(d)}(t)$ | $f_{4}^{(d)}(t)$ | $f_{5}^{(d)}(t)$ | $f_{6}^{(d)}(t)$ | $f_{7}^{(d)}(t)$
---|---|---|---|---|---|---|---|---|---|---
$dim1$ | $\frac{se_{m}}{se_{u}}<1$ | 0.67 | 0.93 | 0.65 | 0.74 | 0.76 | 0.69 | 0.85 | 0.56 | 0.42
$\frac{se_{m}}{se_{u}}>1$ | 0.33 | 0.07 | 0.35 | 0.26 | 0.24 | 0.31 | 0.15 | 0.44 | 0.58
$dim2$ | $\frac{se_{m}}{se_{u}}<1$ | 0.28 | 0.48 | 0.45 | 0.36 | 0.48 | 0.41 | 0.31 | 0.54 | 0.62
$\frac{se_{m}}{se_{u}}>1$ | 0.72 | 0.52 | 0.55 | 0.64 | 0.52 | 0.59 | 0.69 | 0.46 | 0.38

We conclude that the multivariate modeling approach can improve the mean
estimation but is especially beneficial for the prediction of the random
effects. In some cases, including weights in the multivariate scalar product
might further improve the modeling.

#### Covariance Estimation

Figure 27 shows the mrrMSE values for the estimated eigenfunctions of the two
random effects in the standard modeling scenario A of data settings 1 to 5. In
general, we find that the leading eigenfunctions tend to have lower mrrMSE
values. Especially the leading eigenfunctions of the smooth residual show a
high accuracy. There seems to be more variance in the estimation accuracy for
the subject-specific random effect. This can be confirmed with Figure 28,
which contains the estimated eigenfunctions of all 500 simulation iterations
for modeling scenario 1A (grey curves) and compares them to the data
generating function (red curve). Overall we find that the modes can be
reconstructed sufficiently well, albeit considerably better for the random
smooth residual.

Figure 27: mrrMSE values of the eigenfunction estimation for the random
effects $\bm{B}_{i}(t)$ and $\bm{E}_{ijh}(t)$ for different scenarios.

Figure 27 also allows to compare the eigenfunction estimation across the
different data settings. Keep in mind that scenario 5A is based on a different
model so the mrrMSE values are not directly comparable. We find that the
eigenfunction estimation of the smooth residual suffers from strong
heteroscedasticity (2A), whereas the subject-specific random effect seems
unaffected. The same effect (even more pronounced) can be observed for sparse
data (3A), where fewer information about the correlation within functions is
available. On the other hand, $\bm{B}_{i}(t)$ suffers from few different
individuals as this can lead to a considerable departure from the modeling
assumption when the independent draws of scores are not centred and
decorrelated (4A). With 720 different grouping levels of the smooth residual,
we typically get an empirical covariance for the random scores that is closer
to its assumption and thus smaller differences in mrrMSE values compared to
1A. The estimation of the eigenfunctions is somewhat less accurate in scenario
5A as the weights of the scalar product have to be estimated as well. This
additional uncertainty leads to a larger variance for the mrrMSE values and
makes it harder to correctly identify the data generating modes of variation.
Note that this does not affect the overall estimation of the random effects as
discussed above.

Figure 28: Estimated eigenfunctions in scenario 1A for all 500 simulation
iterations (grey curves). The red curves show the data generating
eigenfunctions.

We observe similar trends for the mrrMSE values of the estimated
eigenfunctions in scenario 6A as presented in Figure 29. Leading
eigenfunctions tend to be more accurately estimated and increasing the number
of grouping levels ($\bm{B}_{i}(t):25$, $\bm{C}_{ij}(t):50$,
$\bm{E}_{ijh}(t):300$) seems to have a diminishing effect on the variance of
the mrrMSE values.

Figure 29: mrrMSE values of the eigenfunction estimation for the random
effects $\bm{B}_{i}(t),\bm{C}_{ij}(t)$, and $\bm{E}_{ijh}(t)$ for scenario 6A.

Figure 30 shows the rrMSE values of the estimated multivariate eigenvalues and
measurement error variances. Compared to scenario 1A, the rrMSE values of the
smooth residual are higher for scenarios 2A (strong heteroscedasticity) and 3A
(sparse data), whereas the rrMSE values of the subject-specific random effect
are higher for scenario 4A (correlated scores). This is along the lines of the
findings for the eigenfunctions and the random effects. Scenario 5A is again
not directly comparable and the overall higher rrMSE values suggest that the
uncertainty in the estimation of the weights of the scalar product influences
the accuracy of the estimation of the eigenvalues. With regards to the error
variance, the rrMSE values are comparable over the different scenarios except
for the sparse data setting, where we find higher values.

Figure 30: rrMSE values of the eigenvalue and error variance estimation for
the standard scenario A of data settings 1-5.

We conclude that modes of variation can be recovered well in most of the model
scenarios. The nested random effect and its leading modes of variation can be
well captured by the multiFAMM.

#### Fixed Effects Estimation

Figure 31 shows the urrMSE values of the estimated effect functions in
scenario 1A (blue). We find that the estimation of the functional intercept
$\bm{f}_{0}(t)$ and the covariate effect of order $\bm{f}_{1}(t)$ yield low
urrMSE values on both dimensions. The estimation of the other covariates
($\bm{f}_{2}(t)$ to $\bm{f}_{4}(t)$) and especially the interactions of order
with the other covariates ($\bm{f}_{5}(t)$ to $\bm{f}_{7}(t)$) give larger
urrMSE values and a higher variance of these values. The yellow boxplots show
the corresponding values in scenario 4A, thus indicating that only the
estimation of the intercept is affected by correlated and uncentred scores (as
an empirical score mean different from zero times the corresponding
eigenfunctions is captured by the intercept). Figure 32 plots all estimated
effect functions against the data generating effect functions in scenario 1A.
This suggests that the multiFAMM can overall capture characteristics of the
true underlying effect functions.

Figure 31: urrMSE values of the estimated effect functions in scenarios 1A and
4A. Figure 32: Estimated effect functions in scenario 1A for all 500
simulation iterations (grey curves). The red curves show the data generating
effect functions.

Figure 33 shows the urrMSE values for each effect function in scenario 6A.
Overall, the functional intercept $\bm{f}_{0}(t)$ and the covariate effect of
skill $\bm{f}_{1}(t)$ have a higher estimation accuracy as implied by lower
urrMSE values. The other effect functions are also somewhat less pronounced,
which is why the scaling with the inverse mean norm can give large urrMSE
values. Figure 34 shows all estimated effect functions (grey curves) and the
corresponding data generating functions in red. This plot suggests that
especially on the dimensions $dim1,dim5$, and $dim6$ the effect estimation
shows a larger variance across the simulation runs. These dimensions
correspond to those dimensions in the data set ($elbox.x,shoulder.x,$ and
$shoulder.y$ in Figure 5) where the response functions are relatively constant
over $t$ compared to the other dimensions.

Figure 33: urrMSE values of the estimated effect functions in scenario 6A.
Figure 34: Estimated effect functions in scenario 6A for all 500 simulation
iterations (grey curves). The red curves show the data generating effect
functions.

For all modeling scenarios, we use the average point-wise coverage to evaluate
the $95\%$ CBs of the estimated fixed effects (see Table 12). For scenario 1A,
we find the CBs to cover the true effect $88-95\%$ of the time. Additional
uncertainty e.g. about the number of FPCs further reduces the coverage (for
example 1B). Overall, the coverage of fixed effects by the corresponding CBs
is comparable over the different data settings (see Appendix Table 12). In
data settings 2 and 3 the uncertainty in the data is increased which allows
for wider CBs and thus a slightly better coverage. The averaged pointwise
coverage in 4A is comparable to 1A except for the functional intercept. Here,
the averaged point-wise coverage of the scalar and functional intercept lie
well below $70\%$ as the data setting makes it particularly hard to identify
the true underlying mean (the scores of the random effects are not centred or
decorrelated and any mean difference to zero is absorbed by the intercept). In
scenario 5A, the coverage tends to be lower than for models based on an
unweighted scalar product, possibly due to the added uncertainty from
estimation of the weight in the scalar problem. We find averaged point-wise
coverage to be considerably lower than its nominal value in scenario 6A. In
this scenario, we find that in many of the simulation runs the effects are
wrongly estimated as constants (cf. also Figure 34), which can lead to
undercoverage also in scalar additive mixed models (Greven and Scheipl; 2016).

Table 12: Averaged point-wise coverage of the point-wise CBs for the estimated
effect functions and the averaged coverage of the scalar intercept
$\beta_{0}^{(d)}$ for different scenarios. The coverage is averaged over the
500 simulation iterations and over 100 evaluation points along $t$.

$d$ | Scenario | $\beta_{0}^{(d)}$ | $f_{0}^{(d)}(t)$ | $f_{1}^{(d)}(t)$ | $f_{2}^{(d)}(t)$ | $f_{3}^{(d)}(t)$ | $f_{4}^{(d)}(t)$ | $f_{5}^{(d)}(t)$ | $f_{6}^{(d)}(t)$ | $f_{7}^{(d)}(t)$
---|---|---|---|---|---|---|---|---|---|---
$dim1$ | 1A | $93.6$ | $90.6$ | $92.6$ | $89.5$ | $93.1$ | $93.9$ | $91.3$ | $89.0$ | $90.4$
1B | $93.8$ | $83.5$ | $90.6$ | $86.6$ | $89.7$ | $90.2$ | $87.9$ | $86.6$ | $88.0$
1C | $93.8$ | $90.7$ | $92.5$ | $89.5$ | $93.1$ | $94.1$ | $91.3$ | $88.9$ | $90.8$
1D | $93.4$ | $87.7$ | $91.3$ | $86.2$ | $89.5$ | $91.3$ | $87.9$ | $87.1$ | $88.0$
1E | $92.8$ | $88.7$ | $90.7$ | $86.9$ | $90.2$ | $91.3$ | $88.0$ | $86.9$ | $87.8$
1F | $94.6$ | $92.2$ | $93.7$ | $90.1$ | $94.1$ | $94.7$ | $91.5$ | $90.3$ | $93.2$
1U | $92.6$ | $90.8$ | $92.5$ | $90.2$ | $93$ | $92.6$ | $91.7$ | $89.5$ | $88.2$
2A | $95.4$ | $92.8$ | $94.2$ | $90.6$ | $93.9$ | $94.8$ | $92.9$ | $90.5$ | $92.0$
3A | $94.4$ | $89.2$ | $93.0$ | $91.2$ | $93.0$ | $91.7$ | $90.8$ | $92.4$ | $93.5$
4A | $68.4$ | $48.7$ | $92.7$ | $89.7$ | $92.7$ | $93.4$ | $91.4$ | $89.4$ | $90.2$
5A | $94.0$ | $88.4$ | $91.4$ | $89.8$ | $92.8$ | $92.7$ | $89.7$ | $91.4$ | $94.2$
$dim2$ | 1A | $93.4$ | $87.6$ | $93.3$ | $91.8$ | $94.6$ | $92.9$ | $93.3$ | $93.5$ | $94.1$
1B | $95.0$ | $86.1$ | $93.1$ | $91.9$ | $94.0$ | $93.3$ | $93.6$ | $92.2$ | $92.5$
1C | $93.8$ | $87.5$ | $93.2$ | $91.9$ | $94.9$ | $93.0$ | $93.4$ | $93.5$ | $94.0$
1D | $91.2$ | $86.5$ | $90.7$ | $89.4$ | $89.4$ | $89.5$ | $90.0$ | $88.2$ | $90.0$
1E | $91.0$ | $86.6$ | $90.8$ | $90.0$ | $90.0$ | $89.3$ | $90.2$ | $88.5$ | $90.0$
1F | $93.4$ | $87.8$ | $92.7$ | $92.3$ | $93.7$ | $92.5$ | $93.2$ | $92.9$ | $94.0$
1U | $93.2$ | $87.1$ | $92.3$ | $93.6$ | $94.7$ | $92.0$ | $94.0$ | $93.9$ | $93.6$
2A | $94.8$ | $89.4$ | $92.5$ | $93.2$ | $94.5$ | $93.4$ | $93.8$ | $93.4$ | $92.7$
3A | $94.6$ | $90.1$ | $92.7$ | $92.8$ | $95.8$ | $94.1$ | $94.0$ | $95.2$ | $91.9$
4A | $68.6$ | $54.7$ | $93.2$ | $92.3$ | $94.3$ | $92.7$ | $93.3$ | $93.7$ | $93.8$
5A | $90.4$ | $83.5$ | $88.1$ | $88.6$ | $93.9$ | $87.5$ | $86.7$ | $91.9$ | $91.4$

$d$ | Scenario | $\beta_{0}^{(d)}$ | $f_{0}^{(d)}(t)$ | $f_{1}^{(d)}(t)$ | $f_{2}^{(d)}(t)$ | $f_{3}^{(d)}(t)$ | $f_{4}^{(d)}(t)$
---|---|---|---|---|---|---|---
$dim1$ | 6A | $87.4$ | $84.9$ | $90.8$ | $80.7$ | $74.4$ | $76.0$
$dim2$ | $87.6$ | $88.3$ | $87.9$ | $86.9$ | $75.2$ | $81.4$
$dim3$ | $82.4$ | $85.3$ | $87.1$ | $84.2$ | $81.6$ | $71.7$
$dim4$ | $82.4$ | $82.2$ | $83.6$ | $78.5$ | $78.3$ | $74.1$
$dim5$ | $85.4$ | $82.1$ | $84.3$ | $84.3$ | $78.0$ | $76.3$
$dim6$ | $86.2$ | $74.1$ | $83.0$ | $84.2$ | $79.3$ | $65.7$

Figure 35 shows the point-wise coverage average over the 500 simulation
iterations in scenario 1A. It suggests that coverage tends to lie below the
nominal value in areas close to the border of the functions.

Figure 35: Point-wise averaged coverage over the functional index for the
estimated effect functions in scenario 1A. The nominal value is $95\%$ (red
dotted line).


$\operatorname{rank}\nolimits_{\text{full}}(\xi^{\prime})=(r_{d}^{\prime},\ldots,r_{0}^{\prime})$,
then there exists $3\leq i\leq d$ with $r_{i}^{\prime}<r_{i}$ and
$r_{j}^{\prime}=r_{j}$ for all $i<j\leq d$. This implies
$\operatorname{rank}\nolimits(\xi^{\prime})<_{\text{lex}}\operatorname{rank}\nolimits(\xi)$.

###  C.7 Proof of Lemma 5.17

See 5.17

######  Proof C.11.

If a clean linear KLM sequence $\xi$ is not normal, it is either bounded,
unpumpable, or not rigid. We apply one of Lemmas 5.14, 5.15, and 5.16 to
decompose $\xi$ into a set $\Xi$ of linear KLM sequences such that
$\bigcup_{\xi^{\prime}\in\Xi}L_{\xi^{\prime}}=L_{\xi}$, and
$\operatorname{rank}\nolimits(\xi^{\prime})<_{\text{lex}}\operatorname{rank}\nolimits(\xi)$
and $|\xi^{\prime}|\leq|\xi|^{O(|\xi|)}$ for every $\xi^{\prime}\in\Xi$. Then
apply Lemma 5.12 to clean each sequence $\xi^{\prime}\in\Xi$. We are done by
setting
${\textrm{dec}({\xi})}:=\bigcup_{\xi^{\prime}\in\Xi}{\textrm{clean}({\xi^{\prime}})}$.

###  C.8 Proof of Lemma 5.18

See 5.18

The proof of 5.18 is basically just a repetition of that of [15, Lem. 4.19],
except that we need to take additional care of the part of linear path
schemes. Fix a normal linear KLM sequence
$\xi=\left\langle{p_{0}(\bm{x}_{0})G_{0}q_{0}(\bm{y}_{0})}\right\rangle\Lambda_{1}\cdots\Lambda_{k}\left\langle{p_{k}(\bm{x}_{k})G_{k}q_{k}(\bm{y}_{k})}\right\rangle$
throughout this section, where $G_{i}=(Q_{i},T_{i})$ for each $i=0,\ldots,k$.
First we need a bound on the models of $\mathcal{E}({\xi})$ and
$\mathcal{E}^{0}({\xi})$.

######  Claim 1.

There is a model $\bm{h}_{0}\models\mathcal{E}^{0}({\xi})$ with
$\left\|{\bm{h}_{0}}\right\|\leq|\xi|^{O(|\xi|)}$ such that for every
$i=0,\ldots,k$, every $j\in[d]$ and every $t\in T_{i}$,

  * •

$\bm{m}_{i}^{\bm{h}_{0}}(j)>0$ whenever $\bm{x}_{i}(j)=\omega$,

  * •

$\bm{n}_{i}^{\bm{h}_{0}}(j)>0$ whenever $\bm{y}_{i}(j)=\omega$,

  * •

$\phi_{i}^{\bm{h}_{0}}(t)>0$.

Since $\xi$ is saturated and unbounded, we are done by taking $\bm{h}_{0}$ as
a sum of at most $2d(k+1)+\sum_{i=0}^{k}|T_{i}|\leq 2|\xi|$ models of
$\mathcal{E}^{0}({\xi})$ each of whose norm is bounded by $|\xi|^{O(|\xi|)}$
by Lemma 5.7.

######  Claim 2.

There is a model $\bm{h}\models\mathcal{E}({\xi})$ with
$\left\|{\bm{h}}\right\|\leq|\xi|^{O(|\xi|)}$ such that
$\phi_{i}^{\bm{h}}(t)>0$ for every $i=0,\ldots,k$ and every $t\in T_{i}$.

By Lemma C.5 there is a model $\bm{h}^{\prime}\models\mathcal{E}({\xi})$ with
$\left\|{\bm{h}^{\prime}}\right\|\leq|\xi|^{O(|\xi|)}$. We are done by taking
$\bm{h}=\bm{h}^{\prime}+\bm{h}_{0}$ where $\bm{h}_{0}$ is given by Claim 1.

We also need to borrow the following result from [15], which bounds the
lengths of the runs that “pump” each input and output configurations in $\xi$.

######  Claim 3 ([15, Claim C.4]).

For each $i=0,\ldots,k$ there are paths $u_{i},v_{i}$ of $G_{i}$ with
$|u_{i}|,|v_{i}|\leq|\xi|^{O(1)}$ and vectors
$\bm{x}_{i}^{\prime},\bm{y}_{i}^{\prime}\in\mathbb{N}_{\omega}^{d}$ with
$\bm{x}_{i}^{\prime}\geq\bm{x}_{i},\bm{y}_{i}^{\prime}\geq\bm{y}_{i}$ such
that

  * •

$p_{i}(\bm{x}_{i})\xrightarrow{u_{i}}p_{i}(\bm{x}_{i}^{\prime})$ and
$\bm{x}_{i}^{\prime}(j)>\bm{x}_{i}(j)$ whenever
$\textsc{Facc}_{G_{i},p_{i}}(\bm{x}_{i})(j)=\omega\neq\bm{x}_{i}(j)$ for every
$j\in[d]$,

  * •

$q_{i}(\bm{y}_{i}^{\prime})\xrightarrow{v_{i}}q_{i}(\bm{y}_{i})$ and
$\bm{y}_{i}^{\prime}(j)>\bm{y}_{i}(j)$ whenever
$\textsc{Bacc}_{G_{i},q_{i}}(\bm{y}_{i})(j)=\omega\neq\bm{y}_{i}(j)$ for every
$j\in[d]$.

Now we can prove Lemma 5.18.

######  Proof C.12 (Proof of Lemma 5.18).

Let $\bm{h}$ and $\bm{h}_{0}$ be the models of $\mathcal{E}({\xi})$ and
$\mathcal{E}^{0}({\xi})$ given by Claim 2 and 1. For each $i=0,\ldots,k$ let
the paths $u_{i},v_{i}$ and the vectors
$\bm{x}_{i}^{\prime},\bm{y}_{i}^{\prime}$ be given by Claim 3. Denote by
$\psi_{u_{i}}$ and $\psi_{v_{i}}$ the Parikh image of $u_{i}$ and $v_{i}$,
respectively.

Let $\phi_{i}:=r_{i}\phi_{i}^{\bm{h}_{0}}-(\psi_{u_{i}}+\psi_{v_{i}})$ where
we choose $r_{i}:=\left\|{T_{i}}\right\|\cdot(|u_{i}|+|v_{i}|)+1$. By Claim 3
we have $r_{i}\leq|\xi|^{O(|\xi|)}$. Note that $\phi_{i}(t)>0$ for every $t\in
T_{i}$. Now observe that $\phi_{i}\models K_{G_{i},p_{i},q_{i}}^{0}$. Since
$G_{i}$ is strongly connected, by Euler’s Lemma there is a path $w_{i}$ in
$G_{i}$ whose Parkih image is $\phi_{i}$. We have $|w_{i}|=\sum_{t\in
T_{i}}\phi_{i}(t)\leq r_{i}|\xi|^{O(|\xi|)}\leq|\xi|^{O(|\xi|)}$. Also note
that

$\Delta(w_{i})=\Delta(\phi_{i})=r_{i}\Delta(\phi_{i}^{\bm{h}})-(\Delta(u_{i})+\Delta(v_{i}))=r_{i}(\bm{n}_{i}^{\bm{h}_{0}}-\bm{m}_{i}^{\bm{h}_{0}})-(\Delta(u_{i})+\Delta(v_{i})).$
(80)

Let $\bm{m}_{i}^{0}:=r_{i}\bm{m}_{i}^{\bm{h}_{0}}+\Delta(u_{i})$ and
$\bm{n}_{i}^{0}:=r_{i}\bm{n}_{i}^{\bm{h}_{0}}-\Delta(v_{i})$, then
$\Delta(w_{i})=\bm{n}_{i}^{0}-\bm{m}_{i}^{0}$. Now we prove that
$\bm{m}_{i}^{0},\bm{n}_{i}^{0}\geq\bm{0}$ and
$\bm{m}_{i}^{0}(j),\bm{n}_{i}^{0}(j)>0$ for every coordinate $j$ not fixed by
$G_{i}$. For coordinates $j$ fixed by $G_{i}$ we know that
$\Delta(u_{i})(j)=\Delta(v_{i})(j)=0$ since $u_{i},v_{i}$ are cycles, so
$\bm{m}_{i}^{0}(j)=r_{i}\bm{m}_{i}^{\bm{h}_{0}}(j)\geq 0$ and
$\bm{n}_{i}^{0}(j)=r_{i}\bm{n}_{i}^{\bm{h}_{0}}(j)\geq 0$. For coordinates $j$
not fixed by $G_{i}$, if $\bm{x}_{i}(j)\in\mathbb{N}$ then
$\textsc{Facc}_{G_{i},p_{i}}(\bm{x}_{i})(j)=\omega$ and $\Delta(u_{i})(j)>0$,
so $\bm{m}_{i}^{0}(j)\geq\Delta(u_{i})(j)>0$; if $\bm{x}_{i}(j)=\omega$ then
$\bm{m}_{i}^{\bm{h}_{0}}(j)>0$ and we have $\bm{m}_{i}^{0}(j)\geq
r-\Delta(u_{i})(j)>0$ since
$r>\left\|{T_{i}}\right\|\cdot|u_{i}|\geq\Delta(u_{i})(j)$. Symmetrically we
can prove $\bm{n}_{i}^{0}(j)>0$ for $j$ not fixed by $G_{i}$.

Since $u_{i}$ is a run from
$\bm{x}_{i}=\bm{m}_{i}^{\bm{h}}+\omega\bm{m}_{i}^{\bm{h}_{0}}$, and
$\Delta(u)_{i}(j)\leq\left\|{T_{i}}\right\|\cdot|u_{i}|<r_{i}$ for every
$j\in[d]$, we deduce that $u_{i}$ can be fired from the configuration
$p_{i}(\bm{m}_{i}^{\bm{h}}+r_{i}\bm{m}_{i}^{\bm{h}_{0}})$. Then we have

$p_{i}(\bm{m}_{i}^{\bm{h}}+r_{i}\bm{m}_{i}^{\bm{h}_{0}})\xrightarrow{u_{i}}p_{i}(\bm{m}_{i}^{\bm{h}}+\bm{m}_{i}^{0}).$
(81)

Let
$s_{i}:=1+\left\|{T_{i}}\right\|\cdot\max\\{|w_{i}|,|T_{i}|\cdot\left\|{\bm{h}}\right\|\\}\leq|\xi|^{O(|\xi|)}$,
by monotony we have

$p_{i}(\bm{m}_{i}^{\bm{h}}+s_{i}r_{i}\bm{m}_{i}^{\bm{h}_{0}})\xrightarrow{u_{i}^{s_{i}}}p_{i}(\bm{m}_{i}^{\bm{h}}+s_{i}\bm{m}_{i}^{0}).$
(82)

Symmetrically, we have

$q_{i}(\bm{n}_{i}^{\bm{h}}+s_{i}\bm{n}_{i}^{0})\xrightarrow{u_{i}^{s_{i}}}q_{i}(\bm{n}_{i}^{\bm{h}}+s_{i}r_{i}\bm{n}_{i}^{\bm{h}_{0}}).$
(83)

Note that $s_{i}>|w_{i}|\cdot\left\|{T_{i}}\right\|\geq\Delta(w_{i})(j)$ for
all $j\in[d]$ and that $\bm{m}_{i}^{0}(j),\bm{n}_{i}^{0}(j)>0$ for every
coordinate $j$ not fixed by $G_{i}$, we deduce

$p_{i}(\bm{m}_{i}^{\bm{h}}+s_{i}\bm{m}_{i}^{0})\xrightarrow{w_{i}^{s_{i}}}p_{i}(\bm{m}_{i}^{\bm{h}}+s_{i}\bm{n}_{i}^{0}).$
(84)

Since $\phi_{i}^{\bm{h}}\models K_{G_{i},p_{i},q_{i}}$ and
$\phi_{i}^{\bm{h}}(t)>0$ for every $t\in T_{i}$, again by Euler’s Lemma we
know there exists a path $\sigma_{i}$ from $p$ to $q$ whose Parikh image is
$\phi_{i}^{\bm{h}}$. Observe that

$p_{i}(\textsc{Facc}_{G_{i},p_{i}}(\bm{x}_{i}))\xrightarrow{\sigma_{i}}q_{i}(\textsc{Bacc}_{G_{i},q_{i}}(\bm{y}_{i})).$
(85)

Also note that
$\textsc{Facc}_{G_{i},p_{i}}(\bm{x}_{i})=\bm{m}_{i}^{\bm{h}}+\omega\bm{n}_{i}^{0}$
and
$\textsc{Bacc}_{G_{i},q_{i}}(\bm{y}_{i})=\bm{n}_{i}^{\bm{h}}+\omega\bm{n}_{i}^{0}$.
Since $|\sigma_{i}|=\sum_{t\in
T_{i}}\phi_{i}^{\bm{h}}(t)\leq|T_{i}|\left\|{\bm{h}}\right\|$, we have
$s_{i}>|\sigma|\cdot\left\|{T_{i}}\right\|\geq\Delta(\sigma_{i})(j)$ for all
$j\in[d]$. Thus

$p_{i}(\bm{m}_{i}^{\bm{h}}+s_{i}\bm{n}_{i}^{0})\xrightarrow{\sigma_{i}}q_{i}(\bm{n}_{i}^{\bm{h}}+s_{i}\bm{n}_{i}^{0}).$
(86)

Let $s:=\max_{0\leq i\leq k}s_{i}$ and $r:=\max_{0\leq i\leq k}r_{i}$, one can
verify that the above discussion is still valid with $r_{i},s_{i}$ replaced by
$r,s$. Combining (82), (84), (86) and (83), let
$\pi_{i}:=u_{i}^{s}w_{i}^{s}\sigma_{i}v_{i}^{s}$, we have

$p_{i}(\bm{m}_{i}^{\bm{h}}+sr\bm{m}_{i}^{\bm{h}_{0}})\xrightarrow{\pi_{i}}q_{i}(\bm{n}_{i}^{\bm{h}}+sr\bm{n}_{i}^{\bm{h}_{0}})$
(87)

holds for all $i=0,\ldots,k$. Observe that $|\pi_{i}|\leq|\xi|^{O(|\xi|)}$.

Let $\bm{h}^{\prime}:=\bm{h}+sr\bm{h}_{0}$, then
$\bm{h^{\prime}}\models\mathcal{E}({\xi})$, especially
$(\bm{n}_{i-1}^{\bm{h}^{\prime}},\bm{e}_{i}^{\bm{h}^{\prime}},\bm{m}_{i}^{\bm{h}^{\prime}})\models\mathcal{E}_{\text{LPS}}({\Lambda_{i}})$
for each $i=1,\ldots,k$. By Lemma 4.2 there exists a path $\rho_{i}$ admitted
by $\Lambda_{i}$ such that

$q_{i-1}(\bm{n}_{i-1}^{\bm{h}}+sr\bm{n}_{i-1}^{\bm{h}_{0}})\xrightarrow{\rho_{i}}p_{i}(\bm{m}_{i}^{\bm{h}}+sr\bm{m}_{i}^{\bm{h}_{0}}).$
(88)

Note that
$|\rho_{i}|\leq|\Lambda_{i}|+|\Lambda_{i}|_{*}\cdot|\Lambda_{i}|\cdot\left\|{\bm{e}_{i}^{\bm{h}^{\prime}}}\right\|\leq|\xi|^{O(|\xi|)}$.

Combining (87) and (88), we know that the run
$\pi:=\pi_{0}\rho_{1}\pi_{1}\ldots\rho_{k}\pi_{k}$ is admitted by $\xi$, thus
$\llbracket{\pi}\rrbracket\in L_{\xi}\neq\emptyset$. Observe that
$|\pi|\leq(2k+1)\cdot|\xi|^{O(|\xi|)}\leq|\xi|^{O(|\xi|)}$.

##  Appendix D Proof of Lemma 6.4

See 6.4

######  Proof D.1.

We need to use the fact that if $h$ is a monotone inflationary function, then
so are the functions $h^{\alpha}$ in the Hardy hierarchy, which can be proved
by induction on $\alpha$. Let $G:=H^{\omega^{b}\cdot c}$, first we show the
following inequality holds for all $x\geq 2c$:

$H^{\omega^{b+a}}((c+1)x)\geq(c+1)G^{\omega^{a}}(x).$ (89)

We proceed by induction on $a$. For $a=1$:

$\displaystyle(c+1)G^{\omega}(x)$ $\displaystyle=(c+1)G^{x+1}(x)$
$\displaystyle=(c+1)H^{\omega^{b}\cdot c(x+1)}(x)$ $\displaystyle\leq
H^{\omega^{b}\cdot c}(H^{\omega^{b}\cdot c(x+1)}(x))$
$\displaystyle\text{since }H^{\omega^{b}\cdot c}(x)\geq H^{\omega\cdot
c}(x)\geq 2^{c}\cdot x$ $\displaystyle=H^{\omega^{b}\cdot c(x+2)}(x)$
$\displaystyle\leq H^{\omega^{b}\cdot((c+1)x+1)}(x)$ $\displaystyle\text{by
}x\geq 2c$ $\displaystyle\leq H^{\omega^{b}\cdot((c+1)x+1)}((c+1)x)$ by
monotony $\displaystyle=H^{\omega^{b+1}}((c+1)x).$

Next let’s consider the case $a+1$, assuming (89) holds for $a$. We first
prove the following inequality holds for all $y\geq 2c$ by induction on $j$:

$H^{\omega^{b+a}\cdot j}((c+1)y)\geq(c+1)G^{\omega^{a}\cdot j}(y).$ (90)

The base case $j=0$ holds trivially. For the induction step, note that

$\displaystyle H^{\omega^{b+a}\cdot(j+1)}((c+1)y)$
$\displaystyle=H^{\omega^{b+a}}(H^{\omega^{b+a}\cdot j}((c+1)y))$
$\displaystyle\geq H^{\omega^{b+a}}((c+1)G^{\omega^{a}\cdot j}(y))$ by
induction hypothesis of (90)
$\displaystyle\geq(c+1)G^{\omega^{a}}(G^{\omega^{a}\cdot j}(y))$ by induction
hypothesis of (89) $\displaystyle=(c+1)G^{\omega^{a}\cdot(j+1)}(y).$

Now we can prove (89) for $a+1$:

$\displaystyle H^{\omega^{b+a+1}}((c+1)x)$
$\displaystyle=H^{\omega^{b+a}\cdot((c+1)x+1)}((c+1)x)$ $\displaystyle\geq
H^{\omega^{b+a}\cdot(x+1)}((c+1)x)$
$\displaystyle\geq(c+1)G^{\omega^{a}\cdot(x+1)}(x)$ by (90)
$\displaystyle=(c+1)G^{\omega^{a+1}}(x).$

This completes the proof of (89).

Now we return to the proof of the lemma. For all $x\geq\max\\{2c,x_{0}\\}$ we
have

$\displaystyle h^{\omega^{a}}(x)$ $\displaystyle\leq G^{\omega^{a}}(x)$ since
$h(x)\leq G(x)$ for $x\geq x_{0}$ $\displaystyle\leq(c+1)G^{\omega^{a}}(x)$
$\displaystyle\leq H^{\omega^{b+a}}((c+1)x)$ $\displaystyle\text{by
(\ref{eq:fast-grow-func-claim})}.$

This completes the proof.