{'fraction_non_alphanumeric': 0.02874809802842436, 'fraction_numerical': 0.01243160995823756, 'mean_word_length': 4.834907442387609, 'pattern_counts': {'":': 0, '<': 3, '<?xml version=': 0, '>': 0, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 1, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Spatial Matrix Completion for Spatially-Misaligned and High-Dimensional Air Pollution Data

In health-pollution cohort studies, accurate predictions of pollutant concentrations at new locations are needed, since the locations of fixed monitoring sites and study participants are often spatially misaligned. For multi-pollution data, principal component analysis (PCA) is often incorporated to obtain low-rank (LR) structure of the data prior to spatial prediction. Recently developed predictive PCA modifies the traditional algorithm to improve the overall predictive performance by leveraging both LR and spatial structures within the data. However, predictive PCA requires complete data or an initial imputation step. Nonparametric imputation techniques without accounting for spatial information may distort the underlying structure of the data, and thus further reduce the predictive performance. We propose a convex optimization problem inspired by the LR matrix completion framework and develop a proximal algorithm to solve it. Missing data are imputed and handled concurrently within the algorithm, which eliminates the necessity of a separate imputation step. We show that our algorithm has low computational burden and leads to reliable predictive performance as the severity of missing data increases.


Introduction
In multi-pollutant studies, a dataset is often represented as an (n × p) matrix X, in which the concentrations of p pollutants are collected at n monitoring locations. When evaluating the associations between health outcomes and exposures to air pollution, including some or all of these pollutants in a statistical model can be problematic due to correlations and potential interactions among these components. Hence, dimension reduction is often necessary to obtain a lower-dimensional representation of the original data.
Principal component analysis (PCA) (Jolliffe, 1986) is an unsupervised technique for dimension reduction that has been used in multi-pollutant analysis (Dominici et al., 2003).
PCA essentially provides a mapping of the original data X to a low-rank (LR) approximation U V T , in which U ∈ R n×q and V ∈ R p×q (q < p) are usually referred to as principal component (PC) scores and loadings. The product U V T can be considered to be the best rank-q approximation to X. One can also derive this quantity by solving min W 1 2 X − W 2 F + λ W 0 , for a value of λ dependent on q. Here W 0 denotes the number of non-zero singular values of W . Another approach is to replace this L0 penalty, . 0 , with the nuclear norm, Here W * denotes the nuclear norm of W , which is equal to the sum of its singular values.
When closed-form solutions exist for both problems, the L0 penalty is not convex. When some elements of X are missing, the extension to matrix completion with an L0 norm is NP-hard, and existing algorithms would require time doubly exponential in the matrix's dimensions (Candès and Recht, 2009). With missing data, the convex relaxation of the nuclear norm allows optimization via semidefinite programming, and leads to efficient LR matrix completion (LRMC) algorithms (Cai et al., 2010;Mazumder et al., 2010). LRMC is a powerful tool to reduce dimension while being able to utilize incomplete information.
When multiple pollutants are involved, it is often more challenging as modeling many correlated pollutant surfaces can be intractable. The goal is then to leverage both the LR and spatial structures in the data at monitoring locations in order to recover unobserved data at new locations of interest. While PCA can produce a mapping from the original data to their LR structure, LRMC can take advantage of the LR structure even when the monitoring data have missing observations. While both methods cannot, by themselves, produce estimates of the pollutant concentrations at new locations, the resulted LR structure can be incorporated in a spatial prediction model. However, these methods do not account for external geographic information and spatial correlations across neighboring locations. The lack of utilizing spatial structure in these methods may lead to poor predictive performance at new locations.
A spatially predictive PCA algorithm, PredPCA (Jandarov et al., 2017), was developed to produce LR structure embedded with spatial patterns that can be subsequently predicted well at new locations. While capable of leveraging both the LR and spatial structures of the monitoring data, in its current form PredPCA cannot handle input data with missing observations. This motivated a probabilistic version of PredPCA, ProPrPCA (Vu et al., 2020), that can handle incomplete monitoring data. Using external covariates, the unobserved data at new locations can be estimated directly from the hierarchical model assumptions, or by incorporating the LR structure of monitoring data in a separate spatial prediction model. However, similar to PredPCA, the estimation procedure of ProPrPCA is computationally burdensome, due to a combination of multiple rank-1 approximations, numerical optimization, and likelihood-base imputations.
In this paper, our objective is to develop a practical and computationally feasible of version of PCA that can be applied to the setting of spatially-misaligned multivariate data. Our proposed method, Spatial Matrix Completion, aims to (i) accommodate complex spatial missing patterns and (ii) incorporate spatial information into the LR structure of the data so that estimates at new locations can be easily and accurately derived. we leverage the LRMC framework to formulate a convex optimization problem, and derive a straight-forward algorithm to solve it using proximal gradient descent. While solving similar problem to the existing method ProPrPCA, our proposed method offers a simple algorithm that is easy to implement, computationally efficient, and potentially scalable. In addition, we also discuss the connections among existing methods (LRMC, PredPCA, and ProPrPCA) in conjunction with our proposed method.
For the rest of this paper, we first review the LRMC framework in Section 2. We then propose a spatial matrix completion problem, along with its solution and connections with existing methods in Section 3. In Section 4, we compare the performance of our proposed method to LRMC in recovering missing entries at monitoring locations in simulation studies.
In Section 5, we compare the proposed method to ProPrPCA in predicting exposure data at new locations under spatial misalignment, by reproducing the simulations and real data examples in Vu et al. (2020). Finally, we discuss the implications of these results, our contributions, and potential extensions in Section 6.
2 Review of low-rank matrix completion (LRMC)

Problem formulation
The optimization problem in (1) can be referred to as a LR matrix approximation, or the "fully observed" version of LRMC. When some entries of X are missing, the LR structure of X can be recovered by minimizing the residuals over only the observed indices, where Ω denotes the set of observed indices in X, and [P Ω (X)] ij = X ij if (i, j) ∈ Ω and zero otherwise. Mazumder et al. (2010) prove that the "fully observed" problem in (1) has a closed-form solutionŴ that uses the soft-thresholding operator,

Optimization
HereŨ DṼ T is the singular value decomposition (SVD) of X, with D = diag{σ 1 , ..., σ r }, with σ i being the i−th largest singular value of X, and r being the column rank of X.
We assume that the columns of X have been properly centered. The soft-thresholding operator (Donoho et al., 1995) This result is closely related to PCA, in that, if PCA were to be applied onto X,Ṽ would be returned as the loadings. Mazumder et al. (2010) propose an algorithm that uses proximal gradient descent (Rockafellar, 1976) to solve (2), which is given Table 1. A simple re-derivation of this LRMC algorithm is given in the Supplement. The algorithm consists of two major steps: a gradient descent update, and solving the proximal problem to (2). In particular, the gradient update is simply filling in missing entries with the corresponding entries of the current estimate.
The proximal problem turns out to be exactly the LR approximation problem (1), which has a closed-form solution using soft-thresholding.
3 The spatial matrix completion problem

Proposed optimization problem
The LRMC algorithm is a powerful tool to recover the LR structure of the data even though only a sampling of the entries is observed. Once the missing entries are filled, the LR structure, i.e. the PC scores, can be easily obtained by projecting the imputed dataX onto the direction of its first q right singular vectors. However, under spatial misalignment, the ultimate goal is to produce accurate predictions at cohort locations where pollution data is unavailable. One potential approach is to employ a multi-stage procedure in these cohort studies: 1) imputation to fill in missing elements of the data matrix, 2) dimension reduction to obtain the LR structure of the data, and 3) spatial prediction to estimate these scores at locations of interest.
In the second stage with dimension reduction, ideally we would like to identify principal directions such that the resulting PC scores would retain important characteristics and spatial structure. Having these spatial patterns, the PC scores could be predicted well at new locations in the spatial prediction stage. As a result, we propose the following convex optimization problem for the "fully-observed" scenario, where Z ∈ R n×k contains geographic covariates used in the later prediction stage and thinplate spline basis functions. The basis functions are included to capture any underlying spatial patterns that may not have been explained by other covariates.
When X has missing entries, we propose the following optimization problem min M

Optimization
First, we look at the complete data scenario with the proposed problem in (3). While M is the unknown quantity of the objective function, it is important to keep in mind that we are more interested inŴ = ZM whereM is the optimal solution for (3). This quantity is the LR approximation of X. We give the closed-form solution ofŴ in the following lemma, with the detailed proof in the Supplement.
Lemma. If (Z T Z) −1 exists, then the approximationŴ = ZM of X, whereM is the optimal solution for (3), has a closed-form expression, When X has missing entries, we prove thatŴ = ZM , whereM is the optimizer of (4), can be derived via a proximal algorithm similar to the LRMC algorithm. The steps are summarized in Table 2. It turns out that our algorithm is similar to the LRMC algorithm, with the exception of the insertion of an additional projection step involving Z. Derivation and proof are given in the Supplement. We refer to this as the Spatial Matrix Completion (SMC) algorithm.

Connection with existing methods
Another direct approach to induce spatial patterns into the PC scores was the Pred-PCA algorithm proposed by Jandarov et al. (2017). The PredPCA algorithm also employs the same matrix Z of covariates and spline basis functions in its objective function. This algorithm uses a biconvex formulation of PCA where the PCs are estimated sequentially, Here the algorithm estimates one PC at a time. The quantity Zα/ Zα 2 plays the role of the PC score, while v is the loading. PredPCA directly imposes a spatial structure on the score via Z. By doing so, PredPCA essentially aims to recover the best rank-1 approximation of X that has spatial structure embedded within the left singular vectors. Thus the objective function of PredPCA can be rewritten as an L0 problem with spatial constraints. Heuristically, this is similar to what SMC aims to achieve. The fundamental difference is that SMC utilizes a nuclear-norm penalty, while PredPCA resembles a spatial optimization problem with an L0 penalty. When some data are missing, such reformulated version of PredPCA can potentially be solved using hard-thresholding (Mazumder et al., 2010). However, unlike SMC, the lack of convexity means that convergence to an optimal solution and corresponding theoretical properties are not guaranteed.
In its current form, the sequential estimation and the lack of a mechanism to handle missing data are rather unsatisfying. The probabilistic version, ProPrPCA (Vu et al., 2020), provides a better alternative to PredPCA when missing data are present, as imputation and dimension reduction are handled simultaneously. However, ProPrPCA also involves a series of rank-1 approximations, and requires longer computational time overall due to numerical optimizations of spatial parameters and likelihood-based imputations for missing observations.
In later sections, we compare SMC to both LRMC, in terms of recovering missing entries at monitoring locations, and ProPrPCA, in terms of predicting exposure data at new locations of interest.
4 Low-rank matrix recovery in simulation studies: A comparison between SMC and LRMC

Setups
We first conduct two sets of simulation studies to compare the performance of LRMC and SMC in recovering missing entries. We generate spatially-correlated multivariate data on a dense 100 × 100 grid (N = 10, 000), using the following formula Here R o and R u represent the observed and unmeasured covariates, respectively, in addition to the data of interest X ∈ R N ×p . In practice, these covariates often include, but not limited to, Geographic Information System (GIS) variables. S is a mean-zero stationary structure with exponential spatial covariance, while E include independent white noises. The loadings V represent the mapping from the LR q−dimensional space to the higher p−dimensional space of the data. The specific simulation parameters can be found in the Supplement. The first set is a toy simulation study with q = 1 and p = 4, while the second set of simulations has q = 3 and p = 12. We consider four scenarios defined by the presence of unmeasured covariates (B u = 0 or B u = 0) and spatial correlations (S = 0 or S = 0).
In each simulation, we randomly choose 400 locations, and omit data entries at random.
The level of missing data completely at random (MCAR) for each feature in X varies from 5% to 40%, with a 5% increment. We use LRMC and SMC to recover the missing entries and evalue the mean squared errors (MSE) of these entries. The design matrix Z used in SMC includes the observed covariates R o and thin-plate spline basis functions to capture other spatial variability.
While we focus on matrix recovery in this section, in each simulation, we also pick 100 new locations at random. Under spatial misalignment, the scientific interest lies in the data at these new locations, however, they are 100% missing. Thus accurate prediction of data at these locations using the "observed" locations is necessary. As SMC assumes a mean model for the LR structure of the data, this method is capable of producing an estimate for the data at new locations where no observation on any feature of X is available, while LRMC cannot. Thus, we also evaluate the MSE of SMC at these new locations. Figure 3 shows the results of the toy simulations in which the data with four features (p = 4) are generated based on a rank-1 structure (q = 1). Overall, the mean MSE values for LRMC are always higher than those of SMC in these simulations. As the amount of missing data increases, mean MSE goes up for LRMC, but remains relatively the same for SMC even at MCAR level as high as 40%. It is notable that the variability of SMC results (in red) is the highest at 5% MCAR and decreases with more missing data. This is possibly due to the algorithm's greedy attempt to use rich patterns learned within observed entries to recover a much smaller number of missing entries. The results from the higher-dimensional simulations in which the data generated with 12 features are included in the Supplement. We observe trends similar to these results from the toy simulations.

Predictive performance under spatial misalignment:
A comparison between SMC and ProPrPCA As alluded to in Section 3.1, in health-pollution cohort studies, the goal is not to recover missing data at monitoring locations. However, accurate predictions of the total mass and chemical profile at cohort locations, where pollution data are not readily available, are often required prior to health-pollution analysis. A straightforward solution is a multi-stage approach with imputation (to fill in missing monitoring data if any), dimension reduction (to identify important pollutant profile and also reduce computational burden of later step), and spatial prediction (estimate exposures at unobserved locations). As discussed in Section 3.3, existing methods PCA and PredPCA require complete data as input. LRMC can be used in the imputation step, however, it does not leverage important information from relevant covariates that are often available in environmental studies, such as GIS variables.
ProPrPCA (Vu et al., 2020) provides likelihood-based imputation simultaneously within the dimension reduction step. This method also aims to produce a lower-dimensional representation of the original data that can be easily predicted at new locations. Our proposed SMC model offers an alternative approach, with similar goals to ProPrPCA, but using convex optimization. Thus, in this Section we reproduce the simulations and data application in Vu et al. (2020), to compare SMC directly to ProPrPCA, PredPCA, and traditional PCA.

Simulation studies
The full details of the original simulations are described in Vu et al. (2020). In summary, these simulations include p = 15 pollutant surfaces spanning a dense 100 × 100 grid and generated based on three underlying PC scores. These scores have various levels of variance contribution to the data overall and spatial predictability, i.e. how well these scores can be predicted at new locations using external covariates. In the first scenario where the orders of spatial predictability and variance contribution to the data are the same, all competing methods are expected to perform equally well when there is no missing training data. The second scenario is where the most spatially predictable score did not contribute the most variability to the data. We also consider various settings in which the training data are either complete, missing completely at random (MCAR), or missing at random (MAR), where the missing patterns are associated with external covariates.
In each of the 1000 simulations, 400 training locations and 100 test locations are chosen at random. We implement the multi-stage procedure and evaluate these competing dimension reduction methods based on their performance in the spatial prediction stage with universal kriging. When training data are complete, we compare the "fully-observed" version of SMC to PCA, PredPCA, and ProPrPCA. When the training data are either MCAR or MAR, LRMC is used in the imputation step prior to PCA and PredPCA, while neither ProPrPCA nor SMC requires this extra step. Specially, we compare the predicted scores, estimated at test locations by using universal kriging and the training scores obtained from each dimension reduction method, with the "true" scores, derived by projection of test data. Table 3 provides both the prediction R 2 used in Vu et al. (2020), which reflects the correlation between the predicted score and true score, as well as the median MSE for each PC and for the overall reconstructed data matrix at test locations in the first simulation scenario. As expected, all methods perform equally well when the training data are complete.
Under MAR, data among the first 5 pollutants are more likely to be missing at locations with extreme geographic covariates values. As discussed in the original paper, due to the data generating mechanism, this setup eventually leads to the decreases in median R 2 's for PC1 but slight increases for PC2. While the overall performance decreases under MAR, SMC gives the best result for PC1 (median R 2 = 0.74), followed by ProPrPCA (median R 2 = 0.69). This is in exchange for a slightly better performance for PC2 by ProPrPCA.
The differences between SMC and ProPrPCA are further examined in Figure 4, where we examine scatter plot of prediction R 2 across all simulations. Corresponding results for the second simulation scenario are shown in the Supplement.
Overall, we find that SMC produces results that are on par with ProPrPCA in simulations with various levels of missing data.

Assessment of computational burden
In addition, we find SMC to be more computationally efficient as the sample size increases, as shown in Figure 5. Even though ProPrPCA only takes about 10 seconds for n = 1000 on average, due to complex likelihood-based algebras, ProPrPCA is not as efficient as methods that leverage convex optimization. Even with 2 separate steps of imputation and dimension reduction, LRMC+PredPCA is faster at roughly 0.1 second run time on average. SMC proves to be the most efficient at approximately 0.01 second on average with the same sample size. The increase in computational burden jumps from 0.1 (n = 100) to 10 seconds (n = 1000) on average for ProPrPCA. In contrast, the burden only increases 10 fold from 0.001 to 0.01 for SMC. Although 10 seconds may not seem to be a big deal in practice for one dataset, this has a potential implication for larger datasets in practice that span multiple years or with complicated missing patterns that can complicate the likelihood calculations used in ProPrPCA. Meanwhile SMC has the potential to scale and be further optimized thanks to using similar structure to LRMC.

Data application
Next we apply SMC to the same data application used in Vu et al. (2020), which were collected nationally by the Chemical Speciation Network (CSN), a sub-network of the  (Vu et al., 2020;Jandarov et al., 2017;Keller et al., 2017), we use the 2010 CSN data, which after data processing described in Vu et al. (2020), consists of 221 CSN sites, only 130 sites of which have complete data for all 21 components (overall 32.1% missing data). We additionally look into similar data collected in 2011, which includes 208 sites in total, only 128 sites of which have complete data (overall 30.1% missing data).
We first evaluate how the pollutant profile varies across dimension reduction methods.
We also assess how the estimated loadings and corresponding PC scores change with different methods and whether all sites or just complete sites are used. The results are given in the Supplement. Overall, we find that the estimated loadings are noticeably different when including sites with missing data. SMC gives results that shares similarities with both PredPCA and ProPrPCA.
Next we evaluate the predictive performance via leave-one-site-out cross-validation. In each iteration, we leave one site with complete data on all 21 components as test data, while performing dimension reduction and building spatial prediction model based on training data. For a fair comparison, we use the same universal kring model with exponential covariance structure for the spatial prediction step. The training data in each iteration may consist of either the remaining complete site ("complete" training data), or all remaining site ("full" training data). Similar to Vu et al. (2020), we focus on the prediction for each PC instead of reconstructing the data matrix, because in practice, these PCs would represent the pollutant profile, which could be used as an effect modifier in a model investigating the associations between total mass of particulate matter and health outcomes.
The results for 2010 CSN data are shown in Table 4. When using only complete sites for training data, PCA has acceptable performance for PC2 and PC3, however, gives poor results for PC1 (R 2 = 0.24). Meanwhile, the other three methods give the same results and better R 2 for both PC1 and PC3. When training data consist of all remaining sites, SMC is better than PredPCA overall for all PCs. SMC has the same performance as ProPrPCA for PC1, but there is a slight trade-off in performance between PC2 and PC3. SMC has a substantially higher R 2 for PC2 compared to ProPrPCA, but worse performance for PC3.
These results may be due to the switching between PC2 and PC3 observed in the pollutant profile, as shown in the Supplement. The overall performance for 2011 data, as given in the Supplement, is better than 2010 across all methods. When using all available sites, the performance of SMC is more similar to ProPrPCA in this dataset compared to 2010 data.
For a fair comparison across all methods, we use the same pre-specified spatial information encoded in Z and the same universal kriging model in the spatial prediction stage for all PCs. While simplifying the comparisons, using the same Z matrix to model all PCs may not be effective, and the covariance structure might have not been correctly specified in these data applications. Potential solutions, which are beyond the scope of this paper, include an adaptive selection of features to be included in Z (Bose et al., 2018), or more sophisticated non-stationary models for spatial prediction instead of universal kriging.

Discussion
In this paper, we focus on problems arising in health-pollution cohort studies, in which multi-pollutant data are often spatially misaligned and have a large number of missing observations. The ultimate goal is to develop a dimension reduction technique that is similar to PCA but able to leverage spatial structure and external information to obtain accurate predictions at new locations of interest. The scientific motivation includes the ability to characterize the pollutant profile across locations, and to identify effect modifiers for the health-pollution associations of interest. For example, many studies on fine particulate matter (PM 2.5 ) have shown evidence that the associations between health outcomes and PM 2.5 total mass can be significantly modified by the PM 2.5 chemical composition (Krall et al., 2013;Zanobetti et al., 2014;Kioumourtzoglou et al., 2015;Wang et al., 2017;Keller et al., 2018).
We formulate a convex optimization problem based on the existing idea of LRMC with nuclear-norm penalization. We show that a closed-form solution exists when the original data are fully observed. In addition, we also derive a proximal algorithm to solve the problem when some elements of the data are missing. In simulations with generated spatial data, we show that SMC outperforms LRMC in recovering missing data entries. We also illustrate how SMC can produce reliable estimates of the entire pollutant profile at new locations while LRMC simply cannot. In addition, we show the merits of SMC under spatial misalignment typically seen in air pollution cohort studies by reproducing and comparing to simulated results and data applications given in Vu et al. (2020).
A slight complication of SMC compared to other methods is that it requires the penalty parameter λ. In our current algorithm, the choice of λ is based on a small grid search to reach the desired rank q of the low-rank approximation. However, the grid search does not have a major impact on computation time, as shown in simulation results. Computation time can also be shortened using warm starts (Mazumder et al., 2010).
Under the setting of spatially misaligned and multivariate air pollution data, SMC can produce LR structure with spatial patterns and impute for missing data with considerations of external geographic and spatial information, as illustrated in Table 5. SMC is also able to estimate all PCs simultaneously, whereas the ProPrPCA model requires the PCs to be estimated sequentially. Under SMC, estimated PCs are guaranteed to be orthogonal, and thus considered to be uncorrelated, which is one of the desirable properties of PCA. It is important to note that for ProPrPCA, when data is missing, the parameter estimation and data imputation are separate. The loadings are estimated based only on the observed data.
The data is then imputed with consideration of Z, and projected onto the directions of the loadings to derive the PC scores. One can practically use different Z matrices for the estimation and imputation procedures. This can potentially be more beneficial and more accurate, particularly when there is reasonable evidence to believe that the missing mechanism only depends on a subset of covariates and spline terms included in Z. Meanwhile, imputation and dimension reduction are essentially intertwined in the SMC algorithm, and there is no flexibility in modifying information used for imputing data only. The results discussed in Section 5 show that ProPrPCA consistently produces better results, although SMC follows very closely. The SMC algorithm offers a faster, more compact and elegant alternative to ProPrPCA.
In its current form proposed by Jandarov et al. (2017), an imputation step is required prior to PredPCA. LRMC is a useful method to fill in the missing data, but it does not account for spatial structure while imputing the data. Using LRMC prior to PredPCA may distort the underlying structure of the data even before dimension reduction, and thus worsen the predictive performance.
In recent literature, LRMC has been employed in various problems involving spatially correlated data. For example, LRMC is used in video denoising (Ji et al., 2010), seismic data reconstruction (Yang et al., 2013), and imaging recovery (Shin et al., 2014). (2017) develops a two-phase matrix-completion-based procedure with spatial and temporal considerations to recover corrupted weather data. These methods are intriguing and work well for the purpose of handling correlated missing data. However, none of these approaches directly impose spatial patterns into the the LR structure of the data. To the best of our knowledge, no other method has utilize the LRMC framework in such a way that is relevant to spatially-misaligned multi-pollutant data.
While initially focusing on health-pollution cohort studies, our proposed framework can be applicable to other fields where spatial misalignment motivates a separate exposure model. In such cases, the design matrix Z can be modified to incorporate whatever covariates are necessary, with spline terms that represents various structures not limited to just spatial correlations. Future work includes potential extension to misaligned spatio-temporal data.

Supporting information
Additional information and supporting material for this article is available online at the journal's website. Data used in this paper are available upon request through the MESA Air team at the University of Washington.   Input X, q, λ, and t max Initialize W (0) = 0, t = 1 while not converged or t < t max dõ Input X, Z, q, λ, and t max while not converged or t < t max dȏ