{"id":"01d92b5f-2a4e-4b7a-a5b3-fa4a47d079ec","arxiv_id":"2412.17327","paper_version":1,"verdict":"CONDITIONAL","confidence":"MODERATE","novelty_score":6.0,"correctness_risk":"medium","formal_verification":"none","parameter_count":3,"one_line_summary":"Spatial autoregressive function-on-function regression is estimated by projecting response and predictor curves onto functional principal components, then applying multivariate spatial autoregressive least squares.","lead":"Researchers propose a regression model for functional data collected over space, where outcomes at nearby locations influence each other through a spatial autoregressive term and covariates are entire curves. The method is tested on Brazilian COVID-19 data, where it predicts 2022 daily deaths from 2021 curves slightly better than non-spatial alternatives.","discovery_kind":"extension","skeptic_critique":{"model":"deepseek-v4-flash","headline":"Theorem 3.1 is proven only under exact finite-rank Condition C3, but the paper's own simulation design violates C3; the reported strong-dependence results and sqrt(n)-consistency claim are therefore not covered by the theorem.","rationale":"The reader's primary concern (C3) is also the load-bearing one: the exactness of the reduction and the theorem both hinge on finite-rank representation, and the simulation uses non-finite-rank truth with K chosen by a variance threshold. The contradiction between Proposition 2.1 and the alpha=0.9 design is real but less decisive, because the proposition gives a sufficient, not necessary, condition; the actual contraction can still hold. A careful stress-test should therefore center on whether the estimator tracks the true functions or only their projections. This does not destroy the paper: the algebraic reduction is useful and the empirical comparisons are suggestive, but the theoretical guarantee as stated is narrower than the simulation claims. The appropriate verdict remains CONDITIONAL, which is what the reader already gave; I recommend no change.","tokens_in":25532,"tokens_out":15262,"duration_ms":160684,"concrete_test":"Re-run the existing simulation at n=250 and n=1000 for alpha=0.9 (inverse-distance W), and in each replication compute the squared L2 distance between the true beta and its projection onto the first Kx FPC basis, and between the true rho and its projection onto the first Ky SFPC basis. If the mean ISE converges to the mean projection error rather than to zero, or if ISE(n=1000) is not roughly half of ISE(n=250), the consistency claim in Theorem 3.1 is not what the simulation is measuring. A positive control with beta and rho generated exactly as finite linear combinations of fixed known orthonormal bases, with K fixed, would confirm that the estimator achieves sqrt(n) decay when C3 holds.","verdict_should_be":"UNCHANGED","load_bearing_attack":"The central reduction (3.10)-(3.12) and Theorem 3.1 rely on Condition C3 in the supplement: the true rho(u,t) and beta(s,t) must lie exactly in the span of the first Ky SFPC and Kx FPC eigenfunctions. In the simulation, beta(s,t)=2+s+t+0.5 sin(2*pi*s*t) and rho(u,t)=alpha(1+ut)/(1+|u-t|) are smooth functions that are not finite-rank combinations of the estimated empirical eigenbases, and K is selected by a 95% variance rule, not by exact representation. The proof's step 'by C1-C3, Y(t)=..., X(s)=...' is therefore inapplicable to the advertised operating regime: the least-squares estimator is consistent for the truncated/projected parameters, not necessarily for the true functions, unless K is allowed to grow and a truncation-bias bound is supplied. No such bound or growing-K condition appears. Consequently the empirical gains reported under alpha=0.9, although plausible, cannot be cited as evidence for the consistency theorem. I do not treat the violation of Proposition 2.1's sufficient condition at alpha=0.9 as the main issue: that condition is sufficient, not necessary, and a rough row-sum bound gives spectral radius <1 for the integral operator, so the generated process may still be stable.","agreement_with_reader":"partial"},"referee_report":{"model":"deepseek-v4-flash","summary":"The paper proposes a spatial function-on-function regression (SFoFR) model, given in (2.2)-(2.3), in which the functional response at each areal unit depends on a spatial lag of the response through a bivariate kernel rho(u,t) and on a functional covariate through beta(s,t). The estimation strategy projects the response onto spatial functional principal components (SFPC) and the predictor onto classical functional principal components, rewrites the infinite-dimensional model as the multivariate spatial autoregressive model (3.12), and applies the least-squares estimator of Zhu et al. (2020) to the score vectors. The authors state sqrt(n)-consistency and asymptotic normality for the estimated functions (Theorem 3.1), support the method with simulations against FPC and FPLS competitors, and illustrate it on Brazilian COVID-19 data, with an R package available.","tokens_in":25792,"tokens_out":6805,"duration_ms":63907,"significance":"The algebraic reduction from model (2.3) to the MSAR form (3.12) is clean, and, once the exact finite-rank representation is accepted, the use of the Zhu et al. (2020) estimator is natural and computationally attractive. The paper ships an R package and conducts a genuinely out-of-sample COVID-19 evaluation (2021 training, 2022 prediction), which are strengths that deserve explicit credit. The substantive contribution, adding a SAR-type spatial dependence component to function-on-function regression, fills a real gap. However, the advertised operating regime, an infinite-dimensional model, is not what Theorem 3.1 proves: the theorem holds only under exact finite-rank conditions that the paper's own simulation design does not satisfy, and no truncation-bias analysis is supplied. The contribution is therefore conditional on a nontrivial extension.","major_comments":[{"comment":"The exact equality (3.12) and Theorem 3.1 require beta(s,t) and rho(u,t) to lie exactly in the span of the first Kx and Ky eigenfunctions. In the simulation study (Section 4), beta(s,t)=2+s+t+0.5 sin(2*pi*s*t) and rho(u,t)=alpha(1+ut)/(1+|u-t|) are smooth infinite-rank functions, and Kx and Ky are chosen by a 95% cumulative variance rule rather than by exact representation. Hence the reported estimators are least-squares estimates of a truncated or projected model, and Theorem 3.1 does not cover the reported finite-sample results. The paper must either supply a growing-K truncation-bias bound that establishes consistency for the true functions or redesign the simulations so that the data-generating parameters satisfy C3.","section":"Section 3.1 and supplement, Condition C3"},{"comment":"The text states that the response is generated with a Neumann series and that rho(u,t) satisfies ||rho||_inf < 1/||W||_inf, as required by Proposition 2.1. For alpha=0.9, ||rho||_inf = max_{u,t in [0,1]^2} 0.9(1+ut)/(1+|u-t|) = 1.8, while the row-normalized weight matrix has ||W||_inf = 1, so the stated sufficient condition is violated. The generated process may still be stable because the condition is only sufficient, but the justification written in the paper is incorrect; the authors should verify stability directly, for example through the spectral radius of the induced integral operator, or choose a rho that satisfies the stated condition.","section":"Section 4, simulation design vs. Proposition 2.1"},{"comment":"The proof imports Lemmas 2-6 from Zhu et al. (2020) and Condition C9 directly, but it does not verify the hypotheses of those lemmas when the scores are obtained from estimated eigenfunctions and when Kx and Ky are either fixed or growing. In particular, the asymptotic covariance operator is expressed in terms of the estimated basis, yet no argument shows that the uniform convergence in C1 is sufficient to control the effect of basis estimation inside the least-squares objective. This gap must be closed before the main asymptotic claim can be accepted in the functional setting.","section":"Supplement, proof of Theorem 3.1"}],"minor_comments":[{"comment":"The weak-spatial-dependence cases are labeled 'alpha = 1' in the text, but the simulations use alpha=0.1; the labels should be corrected.","section":"Section 4, discussion of Tables 1 and 2"},{"comment":"The symbol Y^T denotes the n x Ky matrix of scores rather than the transpose of the vector-valued curve Y(t); a distinct symbol such as Y_scores would avoid confusion, especially because eY=vec(Y^T) is introduced immediately afterward.","section":"Section 3.1, equations (3.10)-(3.13)"},{"comment":"The sentence 'The [package anonymized for review] package in provides...' is incomplete because of anonymization; the package name and language should be restored in the published version.","section":"Abstract"},{"comment":"The phrase 'under some regulatory conditions' should read 'regularity conditions.'","section":"Section 1, Introduction"},{"comment":"The bound is attributed to Jensen's inequality, but the step is an application of the Minkowski integral inequality followed by ||Y||_1 <= ||Y||_p on [0,1]; the result is correct, but the justification should be corrected.","section":"Supplement, proof of Proposition 2.1"},{"comment":"The caption says '100 generated sample curves,' but the simulations use n_train in {100, 250, 500, 1000}; the displayed sample size should be clarified.","section":"Figure 1 caption"}],"recommendation":"major_revision","confidential_remarks":"The paper is not circular: the COVID-19 evaluation is genuinely out-of-sample, and the estimator is inherited from Zhu et al. (2020). The main risk is the mismatch between the theorem's finite-rank condition and the simulation's infinite-rank design; I view this as fixable by adding a truncation-bias theory or by adjusting the simulation protocol, hence major revision rather than rejection."},"author_rebuttal":null,"desk_editor":{"model":"deepseek-v4-flash","letter":"Two things to know about this paper. It is the first to put a functional autoregressive spatial term and a functional predictor together in a function-on-function regression, using SFPC for the response and FPC for the predictor to reduce to a finite-dimensional multivariate SAR. The reduction in (3.10)-(3.12) is clean and correct under exact eigen-representations. And the COVID-19 evaluation is genuinely out-of-sample: train on 2021, predict 2022, with R2new of 0.765 versus 0.71-0.72 for non-spatial baselines. Real, if modest.\n\nThe soft spot is exactly where the stress-test note lands. Theorem 3.1's consistency and asymptotic normality require Condition C3: rho and beta lie exactly in the span of the chosen eigenfunctions. In the simulations, the true functions are smooth with infinite-rank representations, and K is chosen by 95% variance, not by exact representation. So the theorem does not cover the advertised operating regime. The authors need a growing-K argument with a truncation-bias bound, or a statement that the estimator is consistent for the projected parameters. Without that, the headline empirical gains under strong dependence cannot be cited as evidence for the sqrt(n) theory. That is the main load-bearing gap.\n\nSecondary point: the paper says the alpha=0.9 simulation satisfies Proposition 2.1's contraction condition, but ||rho||_inf is 1.8 there and ||W||_inf is 1. So the claim is false. As the stress-test notes, the condition is sufficient not necessary, and the process may still be stable, so this is a sloppy statement rather than a fatal one. But it needs fixing.\n\nWhat the paper does well: it positions against Hoshino (2024) and Zhu et al. (2020) honestly, includes an extensive simulation with two weight matrices and n up to 1000, and ships code (anonymized). The application is interesting and the interpretation of the estimated rho function is plausible.\n\nWho this is for: researchers working on areal functional data, especially in spatial econometrics or environmental stats, who want a ready-to-run estimator. It deserves a serious referee; the gap is fixable. I would encourage the editor to send it out, with a clear request to address the C3 issue and the contraction discrepancy.","headline":"Genuinely new combination for spatial function-on-function regression, but the asymptotics only cover an exact finite-rank case the simulations don't satisfy; fixable with a truncation-bound argument.","tokens_in":26411,"tokens_out":3335,"would_cite":true,"duration_ms":31775,"reading_group":"maybe","serious_thinker":"yes","would_accept_peer_review":true},"rs_alignment":null,"lean_confirmation":null,"pith_extraction":{"msc":["62R10","62M30","62J05","62G20"],"pacs":[],"model":"deepseek-v4-flash","headline":"The paper claims that the infinite-dimensional spatial function-on-function regression model can be rewritten exactly as a finite multivariate spatial autoregressive model, and that least-squares estimates of the coefficient functions are…","keywords":["function-on-function regression","spatial autoregressive model","functional principal component analysis","spatial dependence","multivariate spatial autoregressive model","least squares estimation","functional data","COVID-19 case study"],"falsifier":"One concrete check: in the paper's strong-dependence simulation setting, $\\rho(u,t)=0.9(1+ut)/(1+|u-t|)$ gives $\\|\\rho\\|_\\infty = 1.8$ for a row-normalized $W$ with $\\|W\\|_\\infty = 1$, violating the stated contraction condition $\\|\\rho\\|_\\infty < 1/\\|W\\|_\\infty$ from Proposition 2.1, yet the method still works — so computing the spectral radius of the corresponding spatial operator would show whether the true stability condition is much weaker than the paper's proposition claims.","tokens_in":25254,"feed_emoji":"📊","tokens_out":8270,"duration_ms":69792,"temperature":0.7,"pith_summary":"This paper introduces a spatial function-on-function regression (SFoFR) model in which the functional response at each location depends on the functional responses of its neighbors through a spatial lag integral, in addition to the usual functional covariate effect. The central claim is that projecting the response onto spatial functional principal components and the predictor onto ordinary functional principal components converts the infinite-dimensional model exactly into the finite-dimensional multivariate spatial autoregressive equation $Y^{\\top} = WY^{\\top}\\rho + X^{\\top}\\beta + e^{\\top}$. From that reduction, the paper obtains least-squares estimators of the spatial autocorrelation function $\\rho(u,t)$ and the regression coefficient function $\\beta(s,t)$, and proves they are $\\sqrt{n}$-consistent and asymptotically normal. Simulations show the estimator outperforms functional partial least squares and functional PCA when spatial dependence is moderate to strong, and an application to Brazilian COVID-19 data reports improved in-sample and out-of-sample prediction.","feed_headline":"Spatial autoregression upgrades function-on-function regression","feed_subtitle":"Projection to principal components turns it into a finite spatial autoregression with √n-consistent estimates.","key_machinery":"The machinery is principal-component projection followed by algebraic collapse: spatial functional principal component (SFPC) analysis on the response and classical functional principal component (FPC) analysis on the predictor rewrite each infinite-dimensional object using a finite score vector, and orthonormality of the eigenfunctions yields the identity $Y^{\\top} = WY^{\\top}\\rho + X^{\\top}\\beta + e^{\\top}$ (equation 3.12). This identity turns the original inverse problem into a multivariate spatial autoregressive model whose coefficient matrices $\\rho$ and $\\beta$ are estimated by least squares; multiplying the estimated scores back by the eigenfunctions gives the functional estimates $\\hat{\\rho}(u,t)$ and $\\hat{\\beta}(s,t)$.","core_discovery":"The paper establishes that estimating the two functional parameters in the SFoFR model is algebraically equivalent to estimating a finite multivariate spatial autoregressive (MSAR) model once the data are expressed in principal component coordinates. Substituting the truncated Karhunen-Loève expansions of $Y(t)$, $X(s)$, $\\rho(u,t)$, and $\\beta(s,t)$ into equation (2.3) and using the orthonormality of the basis functions collapses the infinite-dimensional model to $Y^{\\top} = WY^{\\top}\\rho + X^{\\top}\\beta + e^{\\top}$, with no loss beyond the truncation itself. The paper then imports the least-squares estimator of Zhu et al. (2020) for the vectorized MSAR equation and, under regularity conditions that include the true functions lying in the span of the retained eigenfunctions, proves the recovered functional estimates are $\\sqrt{n}$-consistent and converge weakly to a Gaussian process.","pith_inferences":["The reduction is purely algebraic, so the same projection-to-scores strategy should extend to spatial error models, spatial Durbin forms, multiple functional predictors, and scalar covariates, provided the corresponding multivariate spatial estimator exists.","The sup-norm contraction condition in Proposition 2.1 is sufficient but likely not necessary; in the paper's own $\\alpha=0.9$ setting the condition is violated, yet the Neumann series simulations converge, suggesting the operative requirement is a spectral-radius condition on the spatial operator.","A direct finite-sample check of the theorem would be to simulate under condition C3, construct Gaussian-process confidence bands using the derived covariance operator, and measure coverage; any substantial undercoverage would point to a missing term in the covariance formula.","If users choose $K$ too small, the estimator will converge to the projection of the truth rather than the truth itself; comparing estimates across increasing $K$ offers a practical diagnostic for truncation bias."],"forward_implications":["Because the infinite-dimensional model collapses to a finite MSAR equation, estimation reduces to a least-squares computation on score matrices, making the method computationally feasible for large areal datasets.","The $\\sqrt{n}$-consistency and asymptotic normality results give a covariance operator from which confidence bands for $\\rho(u,t)$ and $\\beta(s,t)$ can be constructed, enabling inference about the strength of spatial dependence.","In the paper's simulations, SFoFR attains lower integrated squared error for $\\beta(s,t)$ than functional partial least squares and functional PCA when spatial dependence is moderate to strong ($\\alpha=0.5$ and $\\alpha=0.9$), while the comparators degrade substantially.","On the Brazilian COVID-19 data, the SFoFR model yields higher $R^2$ and $R^2_{\\text{new}}$ values than FPLS and FPC, and its estimated spatial autocorrelation function tracks the functional Moran's I pattern across the two years."],"supporting_citations":[{"why":"Provides the least-squares estimator and the invertibility, identification, and asymptotic lemmas for the multivariate spatial autoregressive model used after truncation.","marker":"Zhu et al. (2020)"},{"why":"Supplies the spatial functional principal component analysis used to decompose the functional response.","marker":"Khoo et al. (2023)"},{"why":"Contributes the functional spatial autoregressive model and the Neumann series expansion the paper uses to generate spatially correlated functional responses.","marker":"Hoshino (2024)"},{"why":"Gives the functional partial least squares method used as the baseline comparator in simulations.","marker":"Zhou (2021)"},{"why":"Introduced the function-on-function regression model that this paper extends with spatial dependence.","marker":"Ramsay and Dalzell (1991)"}],"fun_headline_variants":["Spatial autoregression for functional responses","Functional data regression with spatial dependence","Spatial functional regression via principal components","Spatial function-on-function regression: a PCA approach"],"cache_read_input_tokens":3200,"weakest_assumption_plain":"The theoretical claims require the true spatial autocorrelation and regression coefficient functions to lie exactly in the span of the first estimated principal component functions; if that fails, truncation bias makes the estimator consistent for an approximation rather than the truth.","fun_headline_variants_meta":{"raw":{"variants":["Spatial autoregression for functional responses","Functional data regression with spatial dependence","Spatial functional regression via principal components","Spatial function-on-function regression: a PCA approach"]},"model":"deepseek-v4-flash","effort":"low","cost_usd":0.000588,"raw_usage":{"total_tokens":2750,"prompt_tokens":926,"completion_tokens":1824,"prompt_tokens_details":{"cached_tokens":384},"prompt_cache_hit_tokens":384,"prompt_cache_miss_tokens":542,"completion_tokens_details":{"reasoning_tokens":1768}},"tokens_in":542,"tokens_out":1824,"duration_ms":15336,"temperature":1.0,"reasoning_tokens":1768,"cache_read_input_tokens":384,"cache_creation_input_tokens":0},"cache_creation_input_tokens":0},"created_at":"2026-08-11T05:35:47.582548+00:00","model_set":{"reader":"deepseek-v4-flash"},"falsifier":"One concrete check: in the paper's strong-dependence simulation setting, $\\rho(u,t)=0.9(1+ut)/(1+|u-t|)$ gives $\\|\\rho\\|_\\infty = 1.8$ for a row-normalized $W$ with $\\|W\\|_\\infty = 1$, violating the stated contraction condition $\\|\\rho\\|_\\infty < 1/\\|W\\|_\\infty$ from Proposition 2.1, yet the method still works — so computing the spectral radius of the corresponding spatial operator would show whether the true stability condition is much weaker than the paper's proposition claims.","supporting_citations":[],"review_version":1}