and for k = 1 , . . . , r and l = 1 , . . . , r . As before, ˆ G is obtained replacing expectations by averages. Since the random-effect estimators ̂ u 2 $_{ik}$, ̂ v 2 $_{ijk}$, ̂ η ⊗ 2 i and ̂ ξ ⊗ 2 ij are by-products of the EM algorithm, no extra computational costs are incurred in computing ˆ F and ˆ G . Finally, we note that since ˆ h$_{z}$ and ˆ h$_{w}$ live in the interval [0 , 1] , a transformation like $_{arcsin}$√$_{h}$ usually provides a better Normal approximation when ˆ h$_{z}$ or ˆ h$_{w}$ are close to the boundaries. The asymptotic variance of arcsin √ ˆ h is given by avar( ˆ h ) / { 4 ˆ h (1 − ˆ h ) } . The simplest procedure to derive a confidence interval for ˆ h in that case is to construct a standard confidence interval for arcsin √ ˆ h and then back-transform the endpoints. In this section we study the finite-sample behavior of the newestimators by simulation. The main goals are to determine if the new method (i) represents a substantial improvement over common functional ANOVA in presence of time variability, (ii) is at least comparable to the naive approach of pre-warping the data using an existing warping method, (iii) is robust to mild departures from the normality assumptions, and (iv) does not overfit, i.e. is not worse than common functional ANOVA in absence of time variability. To this end we generated data from ten different models, all balanced, with I = 10 groups and J = 5 observations per group. The raw data (5) was sampled on an equally- spaced time grid of ν = 20 points in [0 , 1] , and the noise variance was σ 2 = . 1 2 in all cases. The mean function was µ ( t ) = . 6 ϕ ( t, . 3 , . 1) + . 4 ϕ ( t, . 6 , . 1) in all cases, where ϕ ( t, a, b ) 15 5 Simulations G$_{r}$ + k,r + l = − J 2 4 ( Ω − 1 ) kk ( Ω − 1 ) ll + 1 4 { Ω − 1 T · k ⊗ ( Ω − 1 ) T · k } E ( J ∑ j =1 ̂ ξ ⊗ 2 ij J ∑ j =1 ̂ ξ ⊗ 2 ij T ) {( Ω − 1 ) · l ⊗ ( Ω − 1 ) · l } ,