Optimal Estimation of Large Functional and Longitudinal Data by Using Functional Linear Mixed Model

Checking the adequacy of functional linear quantile regression model

Modern medical devices are increasingly producing complex data that could offer deeper insights into physiological mechanisms of underlying diseases. One type of complex data that arises frequently in medical imaging studies is functional data, whose sampling unit is a smooth continuous function. In this work, with the goal of establishing the scientific validity of experiments involving modern medical imaging devices, we focus on the problem of evaluating reliability and reproducibility of multiple functional data that are measured on the same subjects by different methods (i.e. different technologies or raters). Specifically, we develop a series of intraclass correlation coefficient and concordance correlation coefficient indices that can assess intra-method, inter-method, and total (intra + inter) agreement based on multivariate multilevel functional data consisting of replicated functional data measurements produced by each of the different methods. For efficient estimation, the proposed indices are expressed using variance components of a multivariate multilevel functional mixed effect model, which can be smoothly estimated by functional principal component analysis. Extensive simulation studies are performed to assess the finite-sample properties of the estimators. The proposed method is applied to evaluate the reliability and reproducibility of renogram curves produced by a high-tech radionuclide image scan used to non-invasively detect kidney obstruction.

Functional data analysis has attracted considerable interest and is facing new challenges, one of which is the increasingly available data in a streaming manner. In this article we develop an online nonparametric method to dynamically update the estimates of mean and covariance functions for functional data. The kernel-type estimates can be decomposed into two sufficient statistics depending on the data-driven bandwidths. We propose to approximate the future optimal bandwidths by a sequence of dynamically changing candidates and combine the corresponding statistics across blocks to form the updated estimation. The proposed online method is easy to compute based on the stored sufficient statistics and the current data block. We derive the asymptotic normality and, more importantly, the relative efficiency lower bounds of the online estimates of mean and covariance functions. This provides insight into the relationship between estimation accuracy and computational cost driven by the length of candidate bandwidth sequence. Simulations and real data examples are provided to support such findings. Supplementary materials for this article are available online.

The covariance function of the stochastic component of various soil properties plays an essential role in the estimation and interpolation of spatial phenomena using kriging techniques and in the solution of transport problems in heterogeneous media using stochastic continuum transport modeling. In this study the effect of the sampling network on the estimates of covariance functions of two‐dimensional, isotropic, second‐order stationary processes, characterized by different underlying correlation scales, was analyzed. Estimates of the covariance function based on a standard systematic sampling network were compared with estimates of the covariance functions based on a modified sampling network. Results of these comparisons showed that the modified sampling network was better suited than the systematic network to resolving the statistical properties of the underlying stochastic process at small separation distances. This may result in an improvement in structure identification. Analyses of the effect of the sampling network used for estimating the covariance function on the results of conditional inference procedures such as kriging showed that the kriging estimates were essentially insensitive to the geometric configuration of the sampling network. Uncertainties about these estimates were reduced when using the modified sampling network, however. In the case of unconditional inferences, such as with stochastic analyses of transport problems in heterogeneous media, predictions of the head variance provided by using estimates of the parameters of the covariance function of the log saturated hydraulic conductivity based on the modified sampling network fit the theoretical head variance better than those based on the systematic sampling network. The effect of the sampling network used for estimating the covariance function on the predicted head variance was significant, particularly when the underlying correlation scale was relatively small, when the porous medium texture became finer, and when the hydraulic conductivity variations were treated as being three‐dimensional.

Data with spatio-temporal structure can arise in many contexts, therefore a considerable interest in modelling these data has been generated, but the complexity of spatio-temporal models, together with the size of the dataset, results in a challenging task. The modelization is even more complex in presence of multivariate data. Since some modelling problems are more natural to think through in functional terms, even if only a finite number of observations is available, treating the data as functional can be useful (Berrendero et al. in Comput Stat Data Anal 55:2619–2634, 2011). Although in Ramsay and Silverman (Functional data analysis, 2nd edn. Springer, New York, 2005) the case of multivariate functional data is also contemplated, they do not cope with more than one dimension (only time is considered as covariate) in the domain of the considered functions. In estimating functional data through smoothing methods, a proper framework for incorporating space-time structures can be found in the generalized additive models (GAM), while classical dimension reduction techniques for functional data lead to functional principal component analysis (FPCA). In a previous work Ruggieri et al. (J Appl Stat 40(4):795–807, 2013) extended temporal FPCA, that is FPCA on data modelled as functions of the one-dimensional time, to multivariate (more than one variable, in our case pollutants) context. In this paper the computational aspects of FPCA are extended to more than one dimension: space (long, lat) and/or space-time; moreover, multidimensional (spatial, spatio-temporal) FPCA is extended to multivariate case. In order to provide a generalization of FPCA to multidimensional (spatio-temporal) and simultaneously multivariate data, we link GAM models together with the approach proposed by Ramsay and Silverman (2005). The paper describes all the computational details useful to implement this approach, while its effectiveness will be shown by a multivariate spatio-temporal environmental dataset, whose structure is actually very common in literature.

The functional linear model is an important extension of the classical regression model allowing for scalar responses to be modeled as functions of stochastic processes. Yet, despite the usefulness and popularity of the functional linear model in recent years, most treatments, theoretical and practical alike, suffer either from (i) lack of resistance towards the many types of anomalies one may encounter with functional data or (ii) biases resulting from the use of discretely sampled functional data instead of completely observed data. To address these deficiencies, the first class of robust functional regression estimators for partially observed functional data is introduced and studied. The proposed broad class of estimators is based on thin-plate splines with a novel, computationally efficient quadratic penalty, is easily implementable and enjoys good theoretical properties under weak assumptions. It is shown that, in the incomplete data setting, both the sample size and discretization error of the processes determine the asymptotic rate of convergence of functional regression estimators and the latter cannot be ignored. These theoretical properties remain valid even with multi-dimensional random fields acting as predictors and random smoothing parameters. The effectiveness of the proposed class of estimators in practice is demonstrated by means of a simulation study and a real-data example.

SummaryWe consider the problem of testing functional constraints in a class of functional concurrent linear models where both the predictors and the response are functional data measured at discrete time points. We propose test procedures based on the empirical likelihood with bias-corrected estimating equations to conduct both pointwise and simultaneous inferences. The asymptotic distributions of the test statistics are derived under the null and local alternative hypotheses, where sparse and dense functional data are considered in a unified framework. We find a phase transition in the asymptotic null distributions and the orders of detectable alternatives from sparse to dense functional data. Specifically, the tests proposed can detect alternatives of √n-order when the number of repeated measurements per curve is of an order larger than nη0 with n being the number of curves. The transition points η0 for pointwise and simultaneous tests are different and both are smaller than the transition point in the estimation problem. Simulation studies and real data analyses are conducted to demonstrate the methods proposed.

I will present three projects that are related to the modeling of covariance structures on the Euclidean space and on sphere. Firstly, we propose a method to model isotropic random field on sphere, where a tapered Matern covariance function is used to capture the local behavior while a nonparametric expansion controls the behavior at large distances. A model selection procedure based on residual sum of squares with penalization is used to reduce over fitting. Secondly, we address the issue of modeling axially symmetric spatial random fields on sphere with a kernel convolution approach. The observed random field is generated by convolving a latent uncorrelated random field with a class of Matern type kernel functions. By allowing the parameters in the kernel functions to vary with locations, we are able to generate a flexible class of covariance functions and capture the nonstationary properties. We use pre-computation tables to speed up the computation. For regular grid data on sphere, the block circulant property of the covariance matrix enables us to use Fast Fourier Transform (FFT) to get the determinant and the inverse efficiently. Thirdly, we proposed a semiparametric variogram estimating method through its spectral representation to model the intrinsically stationary random on R. The low frequency part of the spectral density is estimated by solving a regularized inverse problem through quadratic programming. The behavior at high frequencies, however, is modeled via a parametric tail in the form of a power decaying function. The power parameter in the tail is estimated by a log likelihood method. All proposed methodologies are supplemented with simulation studies and real data analyses.

Inference in functional linear quantile regression

Genetics plays a role in age-related macular degeneration (AMD), a common cause of blindness in the elderly. There is a need for powerful methods for carrying out region-based association tests between a dichotomous trait like AMD and genetic variants on family data. Here, we apply our new generalized functional linear mixed models (GFLMM) developed to test for gene-based association in a set of AMD families. Using common and rare variants, we observe significant association with two known AMD genes: CFH and ARMS2. Using rare variants, we find suggestive signals in four genes: ASAH1, CLEC6A, TMEM63C, and SGSM1. Intriguingly, ASAH1 is down-regulated in AMD aqueous humor, and ASAH1 deficiency leads to retinal inflammation and increased vulnerability to oxidative stress. These findings were made possible by our GFLMM which model the effect of a major gene as a fixed mean, the polygenic contributions as a random variation, and the correlation of pedigree members by kinship coefficients. Simulations indicate that the GFLMM likelihood ratio tests (LRTs) accurately control the Type I error rates. The LRTs have similar or higher power than existing retrospective kernel and burden statistics. Our GFLMM-based statistics provide a new tool for conducting family-based genetic studies of complex diseases. Supplementary materials for this article, including a standardized description of the materials available for reproducing the work, are available as an online supplement.

Longitudinal designs for biomarker traits are very appealing for a number of reasons. First, the collection of repeated biomarker measurements over time allows investigators to characterize both the overall and individual trajectories in these measurements. This helps investigators to understand overall effects as well as how individuals vary around these average effects. Distinguishing individual heterogeneity from average trend is impossible without longitudinal data. This distinction is important for understanding individual between-subject in its own right, as well as in developing predictors of disease outcome from longitudinal marker data. There are many statistical issues that are unique to longitudinal data analysis ranging from correctly specifying realistic models which characterize average changes over time and individual variation, accounting for observations that are taken at unequally spaced observation times which may differ between individuals, to accounting for missing data due to occasionally missed measurements or loss to follow-up. When the reason for missing data is not related to the underlying observation process, most longitudinal data analysis methods can be routinely applied. However, complex modelling is necessary when the reason for missing an observation or dropout is directly related to the underlying outcome process [1–2]. Statistical techniques that appropriately analyze continuous and discrete longitudinal data have been discussed by many authors. Diggle, et al. [3] provides a good summary of these techniques for the practitioner, while Fitzmaurice et al. [4] discusses recent advances from a more technical perspective. The three papers in this section present innovative approaches for solving important methodological problems in biomarker studies. Most standard longitudinal data analysis techniques assume that the measurement process is independent of the outcome process. However, when the outcome is related to the underlying observation process, failure to properly account for the observation process can result in bias [5]. Schildcrout et al.[6] proposed an approach for regression modelling for longitudinal binary data where the outcome is measured with probabilities determined by an auxiliary variable that is related, but not identical to the binary outcome. Modelling multivariate longitudinal biomarker data is an increasingly common issue with the advent of an increasing number of important biomarker techniques (e.g., multiple cytokines in a multiplex assay or multiple reproductive hormones that describe complex biological phenomenon). Although various authors have proposed techniques for analyzing multivariate longitudinal data [7–8 and others], little work has focused on estimating multivariate linear mixed models under known constraints in the relationship between the multivariate outcomes. Roy et al [9] proposes such an approach for analyzing multivariate reproductive hormones over the menstrual cycle, where established relationships have solid biological justification. Finally, the collection of longitudinal biomarker data as disease predictors is becoming an increasingly important design in diagnostic medicine. Pfeiffer et al. [10] proposes flexible non-parametric methodology for combining biomarkers that exploits the longitudinal structure in formulating a score that is used for predicting a binary outcome.

Abstract. This paper is devoted to the modeling of altitude-dependent patterns of ozone variations over time. Umkehr ozone profiles (quarter of Umkehr layer) from 1978 to 2011 are investigated at two locations: Boulder (USA) and Arosa (Switzerland). The study consists of two statistical stages. First we approximate ozone profiles employing an appropriate basis. To capture primary modes of ozone variations without losing essential information, a functional principal component analysis is performed. It penalizes roughness of the function and smooths excessive variations in the shape of the ozone profiles. As a result, data-driven basis functions (empirical basis functions) are obtained. The coefficients (principal component scores) corresponding to the empirical basis functions represent dominant temporal evolution in the shape of ozone profiles. We use those time series coefficients in the second statistical step to reveal the important sources of the patterns and variations in the profiles. We estimate the effects of covariates – month, year (trend), quasi-biennial oscillation, the solar cycle, the Arctic oscillation, the El Niño/Southern Oscillation cycle and the Eliassen–Palm flux – on the principal component scores of ozone profiles using additive mixed effects models. The effects are represented as smooth functions and the smooth functions are estimated by penalized regression splines. We also impose a heteroscedastic error structure that reflects the observed seasonality in the errors. The more complex error structure enables us to provide more accurate estimates of influences and trends, together with enhanced uncertainty quantification. Also, we are able to capture fine variations in the time evolution of the profiles, such as the semi-annual oscillation. We conclude by showing the trends by altitude over Boulder and Arosa, as well as for total column ozone. There are great variations in the trends across altitudes, which highlights the benefits of modeling ozone profiles.

Penalized PCA approaches for B-spline expansions of smooth functional data

ABSTRACTIn a plant science Root Image Study, the process of seedling roots bending in response to gravity is recorded using digital cameras, and the bending rates are modeled as functional plant phenotype data. The functional phenotypes are collected from seeds representing a large variety of genotypes and have a three-level nested hierarchical structure, with seeds nested in groups nested in genotypes. The seeds are imaged on different days of the lunar cycle, and an important scientific question is whether there are lunar effects on root bending. We allow the mean function of the bending rate to depend on the lunar day and model the phenotypic variation between genotypes, groups of seeds imaged together, and individual seeds by hierarchical functional random effects. We estimate the covariance functions of the functional random effects by a fast penalized tensor product spline approach, perform multi-level functional principal component analysis (FPCA) using the best linear unbiased predictor of the principal component scores, and improve the efficiency of mean estimation by iterative decorrelation. We choose the number of principal components using a conditional Akaike information criterion and test the lunar day effect using generalized likelihood ratio test statistics based on the marginal and conditional likelihoods. We also propose a permutation procedure to evaluate the null distribution of the test statistics. Our simulation studies show that our model selection criterion selects the correct number of principal components with remarkably high frequency, and the likelihood-based tests based on FPCA have higher power than a test based on working independence. Supplementary materials for this article are available online.

In this paper, we consider rank estimation for partial functional linear regression models based on functional principal component analysis. The proposed rank-based method is robust to outliers in the errors and highly efficient under a wide range of error distributions. The asymptotic properties of the resulting estimators are established under some regularity conditions. A simulation study conducted to investigate the finite sample performance of the proposed estimators shows that the proposed rank method performs well. Furthermore, the proposed methodology is illustrated by means of an analysis of the Berkeley growth data.