An Efficient Two-Dimensional Functional Mixed-Effect Model Framework for Repeatedly Measured Functional Data
Cheng Cao Jiguo Cao Hao Pan Yunting Zhang Fan Jiang Xinyue Li




#1 1




1
  
  
  An Efficient Two-Dimensional Functional Mixed-Effect Model Framework for Repeatedly Measured Functional Data
    CHENG CAO 

    Department of Data Science, City University of Hong Kong 

    and 

    JIGUO CAO 

    Department of Statistics and Actuarial Science, Simon Fraser University 

    and 

    HAO PAN,  YUNTING ZHANG, FAN JIANG 

    Department of Developmental and Behavioral Pediatrics, 
 Shanghai Children’s Medical Center, 
School of Medicine, Shanghai Jiao Tong University 
 
    and 

    XINYUE LI  Corresponding author: Xinyue Li, Email: xinyueli@cityu.edu.hk

    Department of Data Science, City University of Hong Kong
    September 9, 2024
=====================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================

 

0
  
  An Efficient Two-Dimensional Functional Mixed-Effect Model Framework for Repeatedly Measured Functional Data
    CHENG CAO 

    Department of Data Science, City University of Hong Kong 

    and 

    JIGUO CAO 

    Department of Statistics and Actuarial Science, Simon Fraser University 

    and 

    HAO PAN,  YUNTING ZHANG, FAN JIANG 

    Department of Developmental and Behavioral Pediatrics, 
 Shanghai Children’s Medical Center, 
School of Medicine, Shanghai Jiao Tong University 
 
    and 

    XINYUE LI  Corresponding author: Xinyue Li, Email: xinyueli@cityu.edu.hk

    Department of Data Science, City University of Hong Kong
    September 9, 2024
=====================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================

 

§ ABSTRACT 

With the rapid development of wearable device technologies, accelerometers can record minute-by-minute physical activity for consecutive days, which provides important insight into a dynamic association between the intensity of physical activity and mental health outcomes for large-scale population studies. Using Shanghai school adolescent cohort we estimate the effect of health assessment results on physical activity profiles recorded by accelerometers throughout a week, which is recognized as repeatedly measured functional data. To achieve this goal, we propose an innovative two-dimensional functional mixed-effect model (2dFMM) for the specialized data, which smoothly varies over longitudinal day observations with covariate-dependent mean and covariance functions. The modeling framework characterizes the longitudinal and functional structures while incorporating two-dimensional fixed effects for covariates of interest. We also develop a fast three-stage estimation procedure to provide accurate fixed-effect inference for model interpretability and improve computational efficiency when encountering large datasets. We find strong evidence of intraday and interday varying significant associations between physical activity and mental health assessments among our cohort population, which shed light on possible intervention strategies targeting daily physical activity patterns to improve school adolescent mental health. Our method is also used in environmental data to illustrate the wide applicability. Supplementary materials for this article are available online.




Keywords: Functional Mixed Effect Model, Wearable Device Data, Physical Activity Data, Mental Health


1.9 


§ INTRODUCTION



A growing amount of research suggests an essential relationship between adolescents’ physical activity and mental health. For instance, prolonged sedentary behaviors have been found to elevate the risk of depression (). Recent studies have further discovered that the context of when and where physical activity occurs are also influential factors (). Wearable devices can continuously record an individual's physical activity profiles over consecutive days, which offers unprecedented opportunities to analyze the varying associations with health outcomes for large-scale population studies (). However, the dynamicity of the association is little studied as it relies on both functional variations of repeatedly measured activity profiles characterizing temporal differences across different times of the day and longitudinal variations capturing different days of the week. The development of new functional methodological tools for repeatedly measured functional data is required to address this scientific problem.

Our motivating dataset comes from a Shanghai school adolescent study, which aims to examine whether a student's physical activity pattern is associated with mental health outcomes after adjusting for covariates such as demographic information. In this study, all adolescent participants wore ActiGraph for consecutive seven days to obtain activity count signals. A total of 2,313 students aged between 11 and 18 years from several schools in Shanghai participated in the study. The recorded accelerometer signals were summarized into minute-by-minute activity counts, resulting in a total of 1,440 functional grids over a day, and one example of the physical activity profile for one participant is shown in Figure <ref>. For each subject, we collected his/her demographic information and mental health assessment results. More details will be provided in Section <ref>.



The physical activity recorded by advanced accelerometers can be seen as repeatedly measured functional observations with covariate-dependent mean and covariance functions. Understanding the intraday and interday relationships between physical activity and key factors can provide insights into intervention strategies for school adolescent health. However, the specialized data structure makes it difficult to characterize using one-dimensional models alone, which implies the importance of two-dimensional modeling to adequately capture the data natures. Furthermore, many methods struggle with scalability even when encountering smaller sample sizes and functional sampling grids than our motivating data (). To improve the computational efficiency of model estimation and inference is also challenging for practical application. 

Currently, repeatedly measured physical activity data monitored by wearable devices are commonly studied as one special case of longitudinal functional data analysis using functional mixed effect model (FMEM) (). The method aims at fixed-effects estimation and inference under the mixed-effect framework by applying splines (), functional principal component analysis (FPCA) (), or Bayesian methods (). While longitudinal functional data analysis mainly focuses on sparse and irregular longitudinal-based design (), repeatedly measured functional data are often densely sampled and measured over regular longitudinal visits. Hence, the analysis of repeatedly measured functional data is of independent interest (), where the methodology can be refined in three aspects when taking into account the practical data structure. 

First, the variations of effect along both longitudinal and functional directions underscore the need for evaluation from two-dimensional rather than one-dimensional effects as in conventional longitudinal functional studies. While some research suggests bivariate functional models for symmetrical two-dimensional functional data, i.e. images, the modeling framework and estimation procedures rely essentially on their univariate cases and lead to extra computational burden (). Besides, they also inadequately account for complex four-dimensional correlation structures presented in longitudinal and functional contexts. Second, enhancing the correlation structure is crucial for improving model flexibility and universality. The random effect component of traditional one-dimensional FMEM incorporates longitudinal visits through a linear framework with additive assumptions (). In comparison, a nonparametric four-dimensional covariance function that considers continuity along two directions and enables a flexible model structure with minimal assumptions is a potential alternative. Third, there is a growing need to boost the computational efficiency of existing methods to apply to population studies. They are often computationally expensive due to challenges in modeling the complex four-dimensional correlation structure and conducting fixed-effect inference. Although some techniques have been proposed to enhance computational efficiency (), implementing FMEM in large-scale population studies remains challenging.

In this study, we aim to propose a novel two-dimensional functional mixed-effect model (2dFMM) framework for repeatedly measured functional data. Our model considers two-dimensional fixed effects for the specialized dense and regular longitudinal designs and leverages a nonparametric four-dimensional covariance structure. This flexible structure can effectively avoid covariance structure assumption violations that may often occur in longitudinal settings. We also develop a fast three-stage estimation procedure using pointwise and smoothing techniques for accurate bivariate coefficient function estimation, and it can effectively achieve good model interpretability and at the same time ensure fast computation. Moreover, the fixed-effect inference is built under the weak separability assumption of spatiotemporal covariance function to decompose marginal effects () and enable efficient estimation of the four-dimensional covariance function, leveraging FPCA and basis splines to relieve the computational burden. 

The rest of this paper is organized as follows. In Section <ref>, we propose 2dFMM with estimation and inference procedure. Asymptotic results of the proposed estimator are also provided. Extensive simulation studies are conducted in Section <ref> to evaluate the performance of 2dFMM and compare it with existing approaches. Section <ref> further applies 2dFMM to the Shanghai school adolescent study and Australia electricity demand data. Conclusion and discussion are in Section <ref>. 




§ METHODS





 §.§ Two-Dimensional Functional Mixed-Effect Model



The functional response is denoted by Y_i(s,t), ith subject's profile functional time t ∈𝒯, repeatedly measured at longitudinal time s ∈𝒮, where i=1,…,N. We assume that the covariate-dependent mean function μ(s,t,𝐗_i) is the linear combination of P-dimensional time-invariant or time-varying covariates of interest 𝐗_i = (1, x_i,1(s,t), …, x_i,P(s,t))^T at time s and t, which is also known as a standard linear concurrent model. The proposed two-dimensional functional mixed-effect model (2dFMM) is 

    Y_i(s,t)    = μ(s,t,𝐗_i) + η_i(s,t) + ϵ_i(s,t) 
       = β_0(s,t) + ∑_p=1^Px_i,p(s,t)β_p(s,t) + ∑_j=1^∞ξ_i,j(t)ψ_j(s) + ϵ_i(s,t),

where {β_0(s,t), β_1(s,t), …, β_P(s,t)} are corresponding coefficient functions, η_i(s,t) is a bivariate random process with mean zero and four-dimensional covariance function C(s,t;u,v), while random measurement error process ϵ_i(s,t) is mean zero with some covariance function ℰ(s,t;u,v) and independent of η_i(s,t). 

For dimension reduction,  η_i(s,t) is typically decomposed via well-developed two-dimensional FPCA and represented in Karhunen-Loève expansion with orthonormal basis of L^2 (𝒮×𝒯). It ensures that a large proportion of the four-dimensional covariance function C(s,t;u,v) can be explained by the top several terms, but is computationally demanding and not suitable for asymmetric bivariate arguments which are usually assumed in repeatedly measured functional data. To address these issues, we exploit the marginal covariance functions technique to separate the bivariate process for efficiency (). In model (<ref>),  we define W_i(s,t) = η_i(s,t) + ϵ_i,1(s,t), where ϵ_i(s,t) = ϵ_i,1(s,t) + ϵ_i,2(s,t) for model identifiability;  C_𝒲(s,u) = C_𝒮(s,u) + Γ(s,u), where C_𝒮(s,u) = ∫_𝒯C(s,t;u,t)dt is the marginal covariance with respect to 𝒮 and . Let ψ_j(s) be eigenfunction of marginal covariance function, such that C_𝒲(s,u)=∑_j=1^∞τ_jψ_j(s)ψ_j(u); {ξ_i,j(t): j ≥ 1} are random coefficient functions obtained by the projection of η_i(·,t) subject onto the direction ψ_j(s), i.e. ⟨ W(·,t), ψ_j⟩_𝒮 - ⟨ϵ_1,i(·,t), ψ_j⟩_𝒮. It is easy to see that 𝔼[ξ_i,j(t)] = 0 for t ∈𝒯 and 𝔼[⟨ξ_i,j, ξ_i,h⟩_𝒯] = τ_j1_{j=h}. We further define the covariance function Θ_j(t,v) = 𝔼[ξ_i,j(t)ξ_i,j(v)], which depends on jth eigenfunction of marginal covariance C_𝒮(s,u). Therefore, for any s ≠ u and t ≠ v, we can obtain an expression of the four-dimensional covariance function as follows, 

    C(s,t;u,v)    = 𝔼[(Y_i(s,t) - μ(s,t, 𝐗_i))(Y_i(u,v) - μ(u,v, 𝐗_i))] 
       = ∑_j=1^∞ψ_j(s)ψ_j(u)Θ_j(t,v).


Moreover, regarding two components of measurement error ϵ_i(s,t),
ϵ_i,1(s,t) with covariance function Γ(s,u)1_{t=v} depicts variation of repeated visits and white noise ϵ_i,2(s,t) is mean zero with variance φ^2. Therefore, ℰ(s,t;u,v) = {Γ(s,u) + φ^21_{s=u}}1_{t=v}, where 1_{·} is an indicator function. 


Thus, the functional response Y_i(s,t) is i.i.d. Gaussian random process with mean 𝔼[Y_i(s,t)] = μ(s,t, 𝐗_i) and variance function of Y_i(s,t) is σ^2(s,t), while for any s ≠ u and t ≠ v, covariance function is σ(s,t;u,v) = C(s,t;u,v) + ℰ(s,t;u,v).







Remark 1. The proposed bivariate model differs from existing FMEM in two key aspects. First, the expectation of the functional response in our model depends on bivariate coefficient functions, which inherently imply possible variations over the longitudinal domain. Second, our model captures complex correlation structure through the four-dimensional covariance function C(s,t;u,v), removing the constraint of the linear framework in FMEM and allowing a more flexible representation. In addition, our bivariate model addresses the pre-specification difficulty of FMEM design: whether using a random slope effect depends on hard-to-verify a priori assumptions about longitudinal correlation structure. Our data-driven approach circumvents this specification difficulty. We demonstrate this better accommodates subtle complexities like distinctions between in-school weekdays and off-school weekend patterns for our motivating examples. 

Remark 2. The four-dimensional covariance function in Equation <ref> is decomposed under a weaker assumption than strong separability, which is defined as C(s,t;u,v) = C_𝒮(s,u)C_𝒯(t,v) (). We approximate C_𝒯(t,v) within each score function ξ_i,j(t) and corresponding covariance function Θ_j(t,v). The representation relaxes the restriction to the scope of strong separability and is theoretically sound because it delivers near-optimality under the appropriate assumptions (). 

Remark 3. Compared with FMEM in longitudinal functional studies, our coefficient functions are expanded from curves to surfaces, indicating two-dimensional effects on functional responses. It is deemed necessary for repeatedly measured functional data varying in two domains. Despite this, we can still approximate the univariate effect by taking the average of the bivariate coefficient function along either domain. 

Suppose that the sampling points {s_r, r=1,…,R} and {t_l, l=1,…,L} are prefixed in this study with design densities f_𝒮(s) and f_𝒯(t) such that ∫_s_1^s_r f_𝒮(s)ds  = r/R for r ≥ 1 and similarly ∫_0^t_l f_𝒯(t)dt  = l/L for l ≥ 1. Both densities are continuous second-order differentiable with the support 𝒮 or 𝒯, which are uniformly bounded away from zero and infinity. Suppose β_p(s,t) = β_p(s)β_p(t) and only take account into 𝒯, the satisfication of identification conditions 𝔼[β_p(s)] = 1 implies β_p(t) = ∫_𝒮β_p(s,t)f_𝒮(s)ds. Hence, given an estimator of bivariate coefficient function β̂_p(s,t), univariate function estimator can be obtained β̂_p(t) = R^-1∑_r=1^Rβ̂_p(s_r,t)
with covariance estimator Cov{β̂_p(t), β̂_p(v)} =R^-2∑_r_1∑_r_2Cov{β̂_p(s_r_1,t),β̂_p(s_r_2,v)}.




 §.§ Estimation of Model Components



We propose a computationally efficient three-stage model estimation procedure. First, we estimate the fixed effects using bivariate pointwise and post-smoothing (pointwise-smoothing) estimation. We also represent the four-dimensional covariance function C(s,t;u,v) as a combination of eigenfunctions and B-spline basis functions for faster estimation. This is achieved by leveraging the decomposition of the marginal covariance. We assume densely and regularly sampled grids for both domains. Estimation with randomly sampled designs is discussed in Section <ref>. Our approach also lends itself well to parallelization, further accelerating the entire process.

Suppose design points s_r and t_l satisfy densities in Remark 3 , we let a pointwise response vector denoted by 𝐘_s_r,t_l and design matrix 𝐗_s_r,t_l = (𝐗_1, s_r,t_l^T, …, 𝐗_N,s_r,t_l^T)^T, where 𝐗_i, s_r,t_l = (1, x_i,1(s_r,t_l), …, x_i,P(s_r,t_l))^T. Our estimation procedure is as follows.

Step I: Bivariate Pointwise Estimation. We reform a linear model 
𝐘_s_r,t_l = 𝐗_s_r,t_lβ_s_r,t_l + 𝐞_s_r,t_l, where β_s_r,t_l = (β_0(s_r,t_l), …, β_P(s_r, t_l))^T and 𝐞_s_r, t_l∼ N(0, σ^2_s_r, t_l𝐈_N), under the mutual independence assumption, which is commonly adopted in other fixed-effect estimation approach (). We use ordinary least squares estimator β̃_s_r,t_l = (𝐗^T_s_r,t_l𝐗_s_r,t_l)^-1𝐗^T_s_r,t_l𝐘_s_r,t_l for initial estimation because it shares nice statistical properties and low computational cost. Considering the bivariate pointwise conditions, the covariance matrix estimates of coefficient functions can be obtained for any s_r_1, s_r_2, t_l_1, t_l_2,

    Cov{β̃_s_r_1,t_l_1, β̃_s_r_2,t_l_2} =σ_s_r_1,t_l_1;s_r_2,t_l_2𝐇_s_r_1,t_l_1𝐇^T_s_r_2,t_l_2,

where 𝐇_s_r,t_l = (𝐗^T_s_r,t_l𝐗_s_r,t_l)^-1𝐗^T_s_r,t_l and  the estimates for σ_s_r_1,t_l_1;s_r_2,t_l_2 will be provided in step III later. Note that for pairwise observation (s_r,t_l),
the standard error estimate of the linear model which is denoted by σ̃^2_s_r,t_l is the raw estimate to σ^2_s_r,t_l.

Step II: Bivariate Smoothing. As the raw estimate is obtained, bivariate smoothing is required to refine it by integrating neighboring temporal information. The post-smoothers for the bivariate process are rich, for example, bivariate P-splines (), thin plate regression splines (), tensor product smooths (), and sandwich (). 

Here we illustrate the use of sandwich smoother due to its computational efficiency with nice asymptotic properties. Denote raw pth estimated bivariate coefficient function by matrix β̃_p= (β̃_p(s_r, t_l))_R × L, the refined estimator is simply expressed in a closed-form solution such that β̂_p = 𝐒_2β̃_p𝐒_1, where 𝐒_1 and 𝐒_2 are smoother matrices for 𝒯 and 𝒮 respectively, utilizing P-splines in different prespecified number of knots K_R and K_L, respectively. 

With the help of the bivariate smoother, the variability of covariance estimator σ̃^2_s_r,t_l can be further diminished. Applying the sandwich smoother on standard error matrix 𝐑 = (σ̃^2_s_r,t_l)_R × L gives the final covariance estimator 𝐑 = (σ̂^2_s_r,t_l)_R × L. However, there is no guarantee that the resulting estimators are all non-negative. The issue can be handled by trimming the negative values at zero (). 

In this study, we also provide tensor product smooths in the practical implementation of bivariate smoothing. To quantify the wiggliness of the raw estimator, sandwich smoother relies on differencing matrices to account for the distance between adjacent coefficients, while tensor product smooths use common second-order derivatives penalties. The performance of two smoother will be shown in detail in simulation studies.

Step III: Covariance Estimation. To reduce the computing burden of the most time-consuming covariance estimation step, we employ flexible nonparametric methods incorporating FPCA on marginal covariance of C_𝒮(·, ·) and B-splines to approximate random functions ξ_i,j(·), tailored to imbalance number of grids for two domains. The efficient estimation procedure also consists of three stages. 



Firstly, we use the centered data, Y_i(s_r,t_l) = Y_i(s_r,t_l) - μ̂(s_r,t_l, 𝐗_i), to estimate the marginal covariance function C_𝒲(s,u), and obtain the estimates of eigenfunctions ψ_j(s) and score functions ξ_i,j(t). Specifically, given the refined estimator of coefficient functions from Step II, we pool {Y_i(·,t_l), i = 1,… N, l=1, … L} and have sample covariance

    C_𝒲(s_r_1,s_r_2) = (|𝒯|/NL)∑_i=1^N∑_l=1^LY_i(s_r_1,t_l)Y_i(s_r_2,t_l),

where 1 ≤ r_1≤ r_2≤ R. The storage memory and computation of C_𝒲(s_r_1,s_r_2) is light because the number of longitudinal grids is relatively small. We obtain the resulting positive semi-definite covariance function estimates C_𝒲(s_r_1,s_r_2) by kernel-based local linear smoothing in PACE algorithm, which excludes the diagonal entries due to the inflation by the effect of white noise (). Let {ψ̂_j(s): j = 1,…, J} be the estimated eigenfunctions and ξ̃_𝒲,i,j(t_l) = ∫_𝒮Y_i(s,t_l)ψ̂_j(s)ds be the estimated score function, which is used as the approximant to ξ̃_i,j(t_l). The number of components J is determined by a fraction of variance explained (FVE), the threshold of which is set 0.99. Additionally, the variance of white noise φ^2 can be estimated as the average difference between C_𝒲(s_r_1,s_r_2) and C_𝒲(s_r_1,s_r_2). 

Secondly, we estimate the marginal covariance functions Θ_j(t, v) by “observed" functional data ξ̃_i,j(t). Suppose each functional data has B-splines basis expansion ξ̃_i,j(t_l)=𝐁^T_j(t_l)𝐛_i,j, where 𝐁_j(t_l) = (B_j1(t_l), …, B_jK(t_l))^T and 𝐛_i,j = (b_i,j1, …, b_i,jK)^T,  B_jk(t_l) is kth B-splines basis function of jth principal component, K is the number of basis functions. Let K be the same for all j, the basis functions 𝐁(t_l) do not rely on j. The covariance estimator of ξ̃_i,j(t_l), denoted by Θ_j(t_l_1, t_l_2), can be obtained nonparametrically

    Θ_j(t_l_1, t_l_2) = N^-1∑_i=1^Nξ̂_i,j(t_l_1)ξ̂_i,j(t_l_2) = N^-1∑_i=1^N𝐁^T(t_l_1)(𝐛̂_i,j^T⊗𝐛̂_i,j)𝐁(t_l_2).


The choice of the value of K depends on a tradeoff between capturing variations adequately and ensuring computational efficiency. To ensure that the majority of variations are captured by large enough number of basis functions, we also consider the computational efficiency of the basis function expansion. Compared with FPCA, which has the computational complexity O(NL^2 + L^3) for each ξ̃_i,j(t)  (), the usage of B-splines requires O(NLK^2) computations. It implies that when K < (L+L^2/N)^1/2, our approach offers a lighter computational cost. Naturally, we can set K = K_L and K_L = min{L/2, 35} is recommended by <cit.>. Thus, in practical usage, we suggest K_L = min{(L+L^2/N)^1/2, L/2, 35}.

Finally, given Equation (<ref>), we obtain the estimator of four-dimensional covariance function C(s_r_1,t_l_1;s_r_2,t_l_2) as follows. 

    C(s_r_1,t_l_1;s_r_2,t_l_2) =N^-1𝐁^T(t_l_1) {∑_j=1^Jψ̂_j(s_r_1)ψ̂_j(s_r_2)∑_i=1^N𝐛̂_i,j^T⊗𝐛̂_i,j}𝐁(t_l_2),

and therefore, σ̂_s_r_1,t_l_1; s_r_2,t_l_2 = C(s_r_1,t_l_1; s_r_2,t_l_2) + σ̂^2_s_r_1,t_l_11_{s_r_1 = s_r_2, t_l_1 = t_l_2}. 



 §.§ Inference Procedure


Here we show the construction of pointwise and simultaneous confidence bands. These inference procedures require the estimator of Cov{β̃_s_r_1,t_l_1, β̃_s_r_2,t_l_2}, which can be explicitly obtained by Equation (<ref>) and (<ref>). However, this sample covariance function is wiggly because it relies on a pointwise estimator, although the B-splines smoothing technique controls the roughness to some degree. We still need a covariance estimator of refined coefficient Cov{β̂_s_r_1,t_l_1, β̂_s_r_2,t_l_2} in the form of the sandwich smoother, which is coincidentally equivalent to the covariance smoothing approach (). Let β̃_p = vec(β̃_p) indicate a matrix is stacked by column and Cov{β̃_p, β̃_p} be a RL × RL covariance matrix estimates of pth bivariate coefficient function formed by pth diagonal element of Cov{β̃_s_r_1,t_l_1,β̃_s_r_2,t_l_2} for all s_r_1, s_r_2, t_l_1, t_l_2. By tensor product properties β̂_p = (𝐒_1⊗𝐒_2)β̃_p,  the ultimate four-dimensional covariance estimator of the pth bivariate coefficient function is 

    Var{β̂_p(s_r,t_l)} = e_s_r,t_l^T(𝐒_1⊗𝐒_2)Cov{β̃_p, β̃_p}(𝐒_1⊗𝐒_2)^Te_s_r,t_l.
 
where e_s_r,t_l denotes RL-dimensional unit vector with 1 at the (l-1)L+rth entry.

The analytic inference for two-dimensional fixed effects using confidence bands is straightforward. Depending on the pointwise variability for every pair of (s_r,t_l) estimated by Equation (<ref>), the ± 2 standard error surfaces can be constructed by 

    β̂_p(s_r,t_l) ± 2Var{β̂_p(s_r,t_l)}^1/2.

Note that our estimator is biased so that the standard error surface is also called a 95% pointwise confidence bands (PCB) if neglecting the effect of bias term based on the nice approximation property of sandwich smoother (). While PCB is efficient to construct in analytical form, using it for statistical inference can be flawed because it ignores the inherent correlation of functional data and results in false positives ().

2
    


To address the issue, we refer to simultaneous confidence bands (SCB) to account for correlation, which are commonly constructed using nonparametric bootstrap approaches. One method involves multivariate normal simulations (), with high computational cost due to the dimensionality equalling the sampling density of the functional domain. We use parameter simulations based on the number of B-splines basis functions and reduce to a tractable number of dimensions (). The details of the bootstrap algorithm are provided in the Algorithm <ref>. Our bootstrap SCB algorithm makes construction computationally practical by leveraging B-splines technique. In addition, it can also serve as a data-driven inference tool for formal global tests about the coefficient functions without relying on distributional assumptions (). 




 §.§ Asymptotic Results


In this section, we derive the asymptotic distribution of our pointwise-smoothing estimator by showing the asymptotic bias and variance structure. Asymptotics of our estimator is established based on the properties of the least square estimator and sandwich smoother which is equivalent to a bivariate kernel regression estimator with a product kernel,  (RLh_Rh_L)^-1∑_r,lβ_p(s_r,t_l)H_m_R{h^-1_R(s-s_r)}H_m_L{h^-1_L(t-t_l)}, where H_m is the equivalent kernel for univariate penalized splines, h_R and h_L are the bandwidths (). The kernel function H_m is symmetric and bounded. For simplicity, our results are for the case of equally spaced design points and knots. For notation convenience, a ∼ b means a/b converges to 1. 

We first derive the asymptotic bias in the interior points. Let m_R and m_L are difference orders of differencing matrices, m_T = 4m_Rm_L + m_R + m_L for notation simplicity.

Propostion 1. Suppose conditions (a)-(d) and (g) are satisfied, further assume K_R∼ c_R(RL)^b_1, K_L∼ c_L(RL)^b_2, with b_1 > (m_R + 1)m_L/m_T, b_2 > (m_L + 1)m_R/m_T, h_R∼ d_R(RL)^-m_L/m_T, h_L∼ d_L(RL)^-m_R/m_T for some positive constants c_R, c_L, d_R, d_L. Then, for any (s,t) ∈ (0,1)^2, we have 
    
    bias{β̂_p(s,t)}    = (-1)^m_R + 1d^2m_R_R∂^2m_R/∂ s^2m_Rβ_p(s,t) + (-1)^m_L + 1d^2m_L_L∂^2m_L/∂ t^2m_Lβ_p(s,t) + o(h^2m_L_L).

The bias remians same with the sandwich smoother (), containing componests from 𝒮 and 𝒯. Noted that the bias converges at a slower rate at the boundary than in the interior (), the proof of which is ignored here. 

To derive the asymptotic variance of the estimator, we assume the covariates are identically and independently distributed as well as time-invariant. For simplicity of illustration, we also assume no missing values. 

Proposition 2. Suppose conditions (a)-(g) are satisfied and under the conditions of Proposition 1, when N →∞, we have

    var{β̂_p(s,t)} =  2(RLh_Rh_LN)^-1ω_pσ^2(s,t)κ(H_m_R)κ(H_m_L) + o(h_L^4m_L).

where ω_p is the (p,p)th entry of Ω = 𝔼(𝐗_1,s_r,t_l𝐗_1,s_r,t_l^T) and κ(H_m) = ∫ H^2_m(u)du. The proposition implies that the asymptotic variance structure of our estimator has an extra component because of the existence of pointwise least square estimator, compared to the sandwich smoother only (). Additionally, it also shows that the correlation influence of two points can be ignored, similarly with kernel regression estimator (). 

Based on the asymptotic bias and variance structures, and proposition 1 of <cit.>, the corresponding asymptotic distribution of our estimator is given by 

    (RL)^2m_Rm_L/m_T(β̂_p(s,t) - β_p(s,t)) → N(α_p(s,t), V_p(s,t)),

in distribution as R →∞ and L →∞, N →∞, where α_p(s,t) = (-1)^m_R + 1d^2m_R_R∂^2m_R/∂ s^2m_Rβ_p(s,t) + (-1)^m_L + 1d^2m_L_L∂^2m_L/∂ t^2m_Lβ_p(s,t) and V_p(s,t) = 2(RLh_Rh_LN)^-1ω_pσ^2(s,t)κ(H_m_R)κ(H_m_L).



§ SIMULATION



We conduct extensive simulation studies to evaluate the performance of the proposed estimation and inference procedure. Our method is examined not only in bivariate but also univariate perspectives, as other competing FMEM estimation methods often only consider fixed effects over functional domain.

The bivariate functional model is simulated as follows:

    Y_i(s,t) = β_0(s,t) + X_i(s)β_1(s,t) + γ_i0(t) + z_i(s)γ_i1(t) + ϵ_i(s, t),   (s, t) ∈ [0, 1]^2.

The fixed-effect covariates are generated from X_i(s) ∼ N(0, 4) and z_i(s) = 6(s-0.5)^2 + N(0,ρ^2), where ρ represents how noisy the signal of repeatedly measured visits is. For bivariate coefficient functions, we take account of two different types as shown in Figure <ref> of the supplementary material. The first scenario (S1) presents a continuous non-differentiable bivariate function with local zero regions (sparse), while the second scenario (S2) presents a smooth bivariate function. The random effects are simulated as γ_il(t) = a_i1ϕ_1l(t) + a_i2ϕ_2l(t), l=0,1. We use the scaled orthonormal functions

    [ ϕ_1(t); ϕ_2(t) ]∝[ 1.5 - sin(2π t) - cos(2π t)                   sin(4π t) ]^T   if  l=0 
    [ cos(2π t) sin(2π t) ]^T   if  l=1

to capture the subject-level fluctuations. The random coefficients are generated from a_i1∼ N(0, 2σ^2_B) and a_i2∼ N(0, σ^2_B) respectively and σ^2_B depends on the relative importance of random effect SNR_B. The measurement error ϵ_ij(s,t) ∼ N(0, σ_ϵ^2), where σ_ϵ^2 depends on signal-to-noise ratio SNR_ϵ. Here SNR_B is defined as the ratio of the standard deviation of fixed-effect and random-effect surfaces, while SNR_ϵ is the ratio of the standard deviation of all liner predictors and that of the measurement errors (). We set SNR_B = SNR_ϵ = 1.

The performance of the method is evaluated from three aspects reflecting the accuracy of estimation and inference, as well as computational efficiency. First, the estimation error is assessed by integrated squared error (ISE) to measure the difference between the estimate and the underlying truth, defined as ISE(β̂_p) = |𝒮×𝒯|^-1∫_𝒮∫_𝒯(β̂_p(s,t) - β_p(s,t))^2dt, where |𝒮×𝒯| denotes the area of the entire domain. Secondly, the proportion of pointwise surfaces wrapping the true plane in the sandwich form is computed for bivariate functional slope to evaluate inferential performance on fine grids. We use the empirical coverage probability of 95% PCB, defined by , defined by |𝒮×𝒯|^-1∫_𝒮∫_𝒯1_β̂_p(s,t) ∈PCB_p(s,t)dsdt. Additionally, to measure the width of confidence bands we also report integrated actual width (IAW), defined as |𝒮×𝒯|^-1∫_𝒮∫_𝒯{UB_p(s,t) - LB_p(s,t)}dsdt, where UB_p(·,·) and LB_p(·,·) are pointwise estimates of upper bound and lower bound respectively. When comparing with other methods for FMEM from the univariate perspective, we accommodate the above metrics with only 𝒯 direction. Computing time for the entire estimation procedure will also be counted to present the computational cost. 



 §.§ Bivariate Comparative Analysis


We compare our proposed approach (2dFMM) with the concurrent bivariate functional regression method denoted by 2dGAM, which is developed by tensor product smooths (). The following several simulation scenarios are considered. We let sample size N ∈{50, 75, 100}, the number of functional grids L ∈{100, 150, 200}, and the number of longitudinal grids R ∈{10, 15, 20}. The baseline setting is N = 50, R = 10, and L=100, where all other sample generating parameters are fixed at their baseline values when one is changed. The noise argument of the longitudinal signal is set ρ = 0.5. A total of 100 replicates are independently simulated. 



Figure <ref> shows better performance of the proposed 2dFMM compared with 2dGAM regarding estimation accuracy in most of the cases under the first scenario. It is primarily attributed to the choice of sandwich smoother depending on differencing matrices, which can tackle this non-differentiable function. When the difference of dimensions between longitudinal and functional domains reduces, 2dGAM improves and slightly outperforms our method because the technique of tensor product spline basis expansions is appropriate for symmetric grids. The coverage probability of 95% PCB for 2dFMM approaches the nominal level, while 2dGAM is far below it due to the ignorance of the four-dimensional correlation structure, leading to the too-narrow width of confidence bands. Since the 2dGAM model estimation depends on a representation of the generalized additive model, the computing time is not surprisingly much longer. The storage of intermediate large matrices is also problematic as it easily runs out of memory even when the sample generating parameters are not very large. 




Figure <ref> displays the results of the second scenario of the slope coefficient function, while we use tensor product smooths as post-smoother instead. The comparable performance of the two methods regarding estimation results indicates that the smoother is sufficient to compensate for the violation of the independence assumption underneath the pointwise technique when encountering continuous functions. The choice of tensor product smooths also solves the problem of symmetric numbers of longitudinal and functional grids that sandwich smoother has. Despite similar estimation accuracy, our method is still plausible given the nice inferential behaviors and low computational cost. 

In addition, we also perform the empirical coverage probability of 95% pointwise and simultaneous confidence bands in different scenarios. The simultaneous confidence bands reach the nominal level.  The detailed comparisons between them are shown in Table <ref> of the supplementary file.



 §.§ Univariate Comparative Analysis


In longitudinal functional data analysis, fixed effects in FMEM framework are often evaluated only over the functional domain. Several methods including functional additive mixed models (FAMM), fast univariate inference (FUI), and fixed-effect inference for longitudinal functional data (FILF) allowing between-visit correlations are considered for comparisons, while FUI and FILF are designed for simpler computation (). They are prespecified to incorporate random slope covariates using longitudinal time points. However, the performance of FAMM is not shown because it is similar to FUI while taking dramatically longer computing time and narrower confidence bands (). 

All approaches are evaluated in two cases: (i) bivariate functional slope is retained in the true model, but only the marginal effect over the functional domain is examined. (ii) bivariate functional intercept and slope are shrunk to univariate ones in the true model, i.e. β_p(t) = |𝒮|^-1∫_𝒮β_p(s,t)ds, p=0,1. Let the noise argument ρ∈{0.5, 2, 6} control the magnitude of the longitudinal correlation in the true model. A total of 100 replicates are independently simulated. 




Figure <ref>(a) presents the performance of different methods under case (i). It can be seen that the estimation accuracy of 2dFMM outperforms the other approaches at any strength of noise argument because of allowing two-dimensional effects. FMEM approaches are worse as the noise argument tends to be larger, indicating nearly no longitudinal correlation structures and therefore resulting in the misspecification of the functional random slope model. For confidence bands, our method is robust in terms of coverage level and actual width among all kinds of noise arguments, while FILF obtains high coverage level at the expense of wider confidence bands. On the other hand, Figure <ref>(b) reveals the disadvantage of 2dFMM in case (ii) where a strong longitudinal correlation structure and univariate effect appear, which precisely meet the prespecified FMEM framework. The relatively low estimation accuracy is due to the misrecognition of the strong random longitudinal correlated signals as quantities of fixed effects. Given the reason, it is not surprising that the 2dFMM is better with heavier longitudinal noise, while others decline for the misspecification issue again. Despite the disadvantage, our method still presents decent coverage close to the nominal level in any case, even though it is slightly conservative in unfavorable conditions. The patterns remain unchanged for larger sample generating parameters (see Figure <ref>-<ref> of the supplementary material). 

































Table <ref> shows the computing time of these three methods with varying sample sizes and numbers of functional grids. The advantage of 2dFMM facing high-dimensional functional data is attributed to the pointwise-smoothing technique and marginal decomposition of the covariance function, while parallel implementation can further accelerate the process. Remarkably, the computation cost of our approach is not sensitive to either parameter, which causes particular concern in the analysis of our motivating data. However, FUI slows down given a considerable number of functional grids when encountering at least a moderate sample size. It is due to the fact that the estimation procedure for x for each pair of functional time points is computationally challenging once having a larger sample size. In addition, FILF uses the generalized additive model to refine the initial estimator, resulting in the same storage and cost limitation as 2dGAM. 






§ APPLICATION



In this section, we apply our proposed method to two studies for illustration. The first study uses motivating accelerometer data to examine the intraday and interday dynamic associations between adolescents' physical activity and their demographic characteristics, family socioeconomic status, and physical and mental health assessments. The second study uses a public environmental dataset to assess the association between electricity demand and temperature.

To present the statistical inferences of bivariate coefficient functions, we define a new metric Î_p(s,t) = β̂_p(s,t)1_{LB_p(s,t) > 0  or UB_p(s,t) < 0} to quantify and interpret the dynamicity of the associations. Heatmaps of the significance evaluation Î_p(s,t) will be displayed: white regions indicate no significant effects, red indicates significantly positive effects, and blue indicates significantly negative effects, while the darkness of the colors represents the magnitude of the effects. 




 §.§ Application to Shanghai School Adolescent Physical Activity Data


In the Shanghai school adolescent study, we collected each subject's demographic information including a binary indicator of gender, grade from 7th to 12th, annual family income level, and mother's education level. Annual family income contains 7 levels from “≤ 10k" to “≥ 300k", while mother’s education level includes 8 levels from “not graduated from primary school" to “at least master degree", both of which are treated as ordinal variables. Mental health screening questionnaires were also conducted during measuring periods. Several common self-report health measurements are included: Depression Anxiety Stress Scales (DASS), Anticipatory and Consummatory Interpersonal Pleasure Scale (ACIPS), Emotion Regulation Questionnaire (ERQ) including two metrics: cognitive reappraisal and expressive suppression. All health assessment results are numeric to present the subjective psychological conditions within different aspects. 

We first regress the subjects' profiles to demographic and socioeconomic covariates in the baseline model.   To avoid the collinearity of mental health outcomes, each covariate is added to the baseline model separately. The functional response consists of a 16,191 × 1,440 dimensional matrix, where each row corresponds to “day of week” and each column corresponds to “time of day”.







Figure <ref> shows the results of the baseline model, where only demographic data and socioeconomic status are incorporated. The evaluation heatmap of the estimated intercept function is consistent with the overall student physical activity patterns. It is observed that girls usually fall asleep later but wake up earlier, while boys are more vigorous in the daytime; yet it is surprising that girls essentially have more intense activity than boys on Saturday night and Sunday. Grade is an alternative characterization for age but has more substantial indications about school workloads and pressure. From the heatmap, elder students are critically more sedentary most of the time, except at 6 a.m., 6 p.m., and midnight, suggesting longer studying time due to laborious tasks. As for two socioeconomic covariates, the effect of family income level is weaker than mother's education level. As mohter's education level increases, students get up and go to bed later, probably because of more demanding assignments from family or richer material well-being such as commuting by private cars and more entertaining options at night. In addition, the results of the univariate effects along the time of day for only weekday or weekend in the baseline model are illustrated in Figure <ref> of the supplementary material.



Figure <ref> demonstrates inferential heatmaps for self-reported health assessment results assessed by well-known questionnaires. The associations between physical activity and scores of subjective depression, anxiety, or stress (DASS) measures share similar patterns. Students tend to be more active at the weekend midnight if they have more severe symptoms of mental health. Nevertheless, being sedentary during 6 p.m. after school or on weekends is because a student who is considered symptomatic for DASS is likely reluctant to exercise. Anticipatory and Consummatory Interpersonal Pleasure Scale (AICPS) evaluates one's pleasant sensations. It is expected that students who own high ACIPS tend to be more active in the daytime, especially on weekends, revealing that happy students are typically more motivated and energetic. In addition, emotion regulation strategies (ERQ) measure respondents' propensity to adjust their emotions, comprising cognitive reappraisal (ERQ(CR)) and expressive suppression (ERQ(ES)). For example, heavily relying on expressive suppression leads to substantial physical activity reduction on weekends, while cognitive reappraisal mainly affects student behaviors during the weekend daytime. In addition, the results of the univariate effects for health assessment results along the time of day for only weekdays or weekends are illustrated in Figure <ref> of the supplementary material.



 §.§ Application to Electricity Demand Data


In Adelaide, Australia, summer electricity demand exhibits high volatility and a strong correlation with temperature. Several studies have examined this relationship under various temperature conditions (). The electricity demand and temperature data, accessible through the R  package (), span from July 6, 1997, to March 31, 2007, providing half-hourly recordings for each day of the week. Our investigation focuses on whether electricity demand is associated with temperature and the impact of weekends as a binary indicator under our proposed modeling framework. The dataset comprises 63 two-dimensional samples, each representing a weekday of a year. The response variable is electricity demand, measured in megawatts, recorded at half-hourly intervals within each week of the year, resulting in 48 functional grids over a day and 52 longitudinal grids over a year. 





Figure <ref> demonstrates the inferential heatmaps for intercept, temperature, and weekend. The intercept shows that electricity demand was generally higher during the mid-year weeks (Australian winter). During summer months, temperature had a strong positive effect on electricity demand, particularly from 10 a.m. to 8 p.m., suggesting a combination of factors including residential cooling needs and the population's daily activity patterns. Conversely, the mid-year weeks (Australian winter) show decreased demand with rising temperatures. This negative association may be attributed to Adelaide's pleasant winter temperatures, reducing the need for heating. Additionally, compared to weekdays, weekend corresponds to lower electricity demands throughout the entire year, especially in the morning around 6 a.m., which demonstrates the habit of waking up late on weekends. Figure <ref> presents the univariate effects of time of day and week of year separately. While they provide strong supporting evidence for bivariate effects, the interplay of two temporal directions is overlooked. For example, the trend of prolonged positive effect was observed around week 10 and week 40 at 3 p.m. only from a bivariate perspective, possibly because work-related and other activities offset the climate effect.  



§ CONCLUSION



Motivated by daily activity profiles obtained by wearable device technologies, our work has provided an effective method for analyzing functional data over a series of times, which is common for longitudinal studies and known as repeatedly measured functional data. The analysis of this type of functional data is increasingly concerned in large-scale medical health research and many other studies such as biomedicine and environmental sciences. The development of efficient and flexible statistical tools for analyzing such data with ultra-high dimensionality and complex longitudinal-functional structure is an urgent problem. In this study, we aim to establish the two-dimensional functional mixed-effect model and efficient fixed-effect inference. We illustrate our method with the analysis of daily adolescent physical activity profiles and hourly electricity demands data, exploring their associations with various covariates of interests. 

The proposed fast three-stage estimation procedure sufficiently reduces computing costs for big data, especially with large-scale samples and functional grids. While we focus on the situation where the longitudinal sampling design is dense and regular in this study, real data is sometimes irregularly sampled. To address this, the estimation of bivariate coefficient functions can be adjusted by taking the average inside equal-size rectangular bins (). Moreover, the covariance function can be estimated utilizing a local-linear smoothing approach and functional principal components analysis through conditional expectation (), while the remaining two-dimensional fixed-effect inference procedure is unchanged. Therefore, our method promotes the use of wearable devices in health research and has wider applications for longitudinal studies and spatial analysis. 




SUPPLEMENTARY MATERIAL






  Supplementary Material: The Supplementary Materials contain results of additional simulation studies, additional analyses of Shanghai School Adolescent Physical Activity Data study, and the proof of the propositions.

  R-package for method and simulations: All code for model implementation and simulation is available at <https://github.com/Cheng-0621/2dFMM>.





agsm


