INTRODUCTION

The propos ed me thods are motiv ate d by the inv es ti gation of diffe re nces in brain funct ional organizat ion be tween peop le with auti sm spectrum di sorde r (ASD) a nd their typ ically develop ing pe ers . The autism brain imaging data exchange (ABIDE) neuroimaging study of res ting-s tate function al m agne tic res on anc e imagin g (rfMRI ) aggr ega ted and publicly shared neuroimaging data on p articip ants w ith A S D and neurotypical c ontrols from 16 in te rn ation al im aging sites . Rel ation ships be tw e en activated brain regions of in te res t (ROIs) during rest can be characte rized b y function al c onne ctivity betw e en ROIs using rfMRI data. Function al c onne ctivity betw e en 2 ROIs is typically estim ate d from the s ubje ct-spe c i fic cor rel ation con structed from the rfMRI time series (Shehzad et al., 2014 ). To avoid the computa tional bur de n of s toring a nd modeling a la rge c onne ctivity matr ix, the rfMRI time ser ies is av erage d across v oxels in each ROI before computing the cros s-ROI correl ation and modeling its r ela tionship with covaria tes. Othe r a pproaches for qua n tifyin g co v ari ate effe cts on c orr ela tion ar e pr imar ily univar iate (see, eg, Shehzad et al., 2014 ). These s ta nda rd a pproaches a re subs ta n tially unde rpowe red to detect small effect sizes and o bs cure importa n t va ria tion within ROI s. While ABIDE da ta h av e le d to some evidence that ASD can be broadly cha racte rized as a brain dysc onne ction syndrome, findings h av e varie d across s tudies (Al aerts e t al ., 2014 ;Di Mart ino et al ., 2014 ). Computat ionally and stat ist ically efficient est imat ion of functional connectivity maps is es s enti al in identifying AS D dysc onne ctions . In this pa pe r, we sho w ho w t o estimat e the effe ct of c ov ari ates, including AS D, on ov er 6.5 mi l lion within-and betw e en-ROI c orr ela tions in appr oxima tely 4 hour s using 4 GB of RAM . In con tras t, a conve n tional a nalysis t ake s ove r 190 hours a nd 175 GB of RAM.

The firs t compone n t of our efficie n t a pproach is to leve rage the spatial depe nde nce within a nd betw e en all v o xels in R OIs using Gaus si an proces s models (Ba ne rj ee e t al., 2014 ; Cres sie, 1993 ; Cressie and Wikle, 2015 ). Denote the rfMRI outcomes { y i ( s ) } N i =1 for N indepe nde n t sa mples at one of S voxels s ∈ S. The join t dis tribution of y i (S ) is ass ume d m ultiva riate Gaussia n a nd known up to a v e ctor of pa ra mete rs of in te res t θ. Without further modeling assumptions or dimension reduction te chniques, m aximum l ikel ihood est imat ion has memory and c omputation al c omplexity of O (NS 2 ) and O (NS 3 ) , respe ctiv ely, due to the S -dimension al c ov ari anc e m atrix. For infe re nce on θ when S is large, the crux of the pro b lem is to adequate ly mode l the spatial depe nde nce without s toring or inverting a large cov ari anc e m atrix.

This pro b lem h as re c eiv e d c on siderab le a t tention (Bradley e t al., 2016 ;Heaton e t al., 2019 ;Liu e t al., 2020 ;Sun e t al., 2011 ). So lution s include, for example, the composite l ikel ihood (CL) (Lindsay, 1988 ), the nea res t-nei ghbor Gaus si an proces s ( NNGP) ( Da t ta et al., 2016 ;Finley et al., 2019 ), spectral methods (Fue n tes, 2007 ), ta pe re d c ov ari ance function s (Furrer e t al., 2006 ;Ka ufman et al., 2008 ;S tein, 2013 ), lo w -ra nk a pproximation s (Cres sie and Jo hannes s on, 2008 ; Nychka e t al., 2015 ; Zimme rma n, 1989 ), and combinations thereof (Sang and Huang, 2012 ). Stein ( 2013 ) details shortcomings of the afore me ntioned methods, which a re prima rily r ela t ed t o loss of s tatis tical efficiency due to simplification of the cov ari anc e m atrix or its in verse. M oreover, thes e me thods rem ain c omputation ally burden s ome when S is very large, a pro b lem that is further exacerbated when the cov ari ance structure is nonstationary.

We propose a new distributed and int egrat ed framework for big spati al d ata an alysis . We c onsider a rec ur sive part it ion of the spatial domain S into M nested res o lution s with dis j oint s e ts of spatial observations at each res o lution, and build local, fully spec i fied di stributed model s in each s e t at the hi ghes t spa tial r eso lution . The main technical d iffic ulty arises from integrating infe re nce from these distributed models over 2 levels of dependenc e: betw e en s e ts in each res o lution and be tw e en res o lution s. The main contribution of this paper is a rec ur sive estimator that int egrat es inference over all s e ts and res o lution s by alternating betw e en lev e ls of de pe nde nce at each in tegration s tep. We also propose a seque n tial in t egrat ed estimat or that is asymptotically equivale n t to the rec ur sive int egrat ed estimat or but re duc es the c omputation al burden of rec ur sively integratin g o ve r m ultiple res o lution s . The res ult ing mult ires o lution rec ur sive integration (MR R I) framework is flexible, stat ist ically efficient, and computationally s cal ab le through its formul ation in the MapRe duc e pa radi gm.

The rest of this paper is organized as follows. Section 2 e st ablishes the formal pro b lem s e tup, and des cr ibes the distr ibuted mode l-building ste p and the rec ur sive in tegration sche me for the proposed MR R I framework. S ection 3 e valua tes the pr oposed frameworks with simul ation s . Se ction 4 prese n ts the analysis of da ta fr om ABI DE . Proofs, addit ion al res ults, A BIDE infor mation, a nd a n R package a r e pr ovide d as s upplemen ta ry mate rials .

RECUR SIVE MODEL INTEGR ATION

FR A ME WOR K 2.1 Proble m s et-up Suppos e we o bs erve y i ( s ) = α i ( s ) + ϵ i ( s ) the i th o bs erv ation a t loca tion s ∈ R d , i = 1 , . . . , N indepe nde n tly a nd s ∈ S a s e t of S loca tions, wher e α i (•) cha racte rizes the spatial variations and ϵ i (•) is an independent normally distributed measure me n t error with mean 0 and v ari ance σ 2 independent of α i (•) . The nugge t v ari ance σ 2 is assumed fixed across i and s . We justify this assumption for our neuroimaging application below, and d isc uss extensions in Section 5 . Suppose for each r eplica te i we o bs erve q exp l anatory v ari ab les X i ( s ) ∈ R q . We de note b y S = { s j } S j=1 the s e t of location s and define

j=1 . In Section 4 , S is the s e t of voxels in the left and right prece n tr al gyri, visualiz ed in the supplemen t , a nd d = 3 . Outcomes { y i ( s ) } N i =1 consist of the thinned rfMRI time series at voxel s for ABI DE part icipa n ts passing quality con trol, whe re the thinnin g remo v es ev ery se c ond time poin t to re mov e autoc orr ela tion; autocorr ela tion plots ar e pr ovide d as s upple me n ta ry m aterials . Thinning av o ids the b ias oc casion ally introduc e d by othe r a pproaches (Mon t i , 2011 ) a nd ali gns with usual rfMRI autocorr ela t ion assumpt ions (Arbabshirani et al., 2014 ). Spec i fically, for p articip ant n ∈ { 1 , . . . , 774 } and voxel s ∈ S with a time series of length 2 T n , the outcome y i ( s ) consists of the (ce n te red a nd s ta nda rdized) rfMRI outcome at time point t ∈ { 2 , 4 , . . . , 2 T n } , where i = t + n ′ <n T n ′ indexes the pa rticipa n ts a nd time poin ts . This thinning res ults in N = 75888 obs erv ation s of y i ( s ) . While we focus on 2 ROIs for simplicity of exposition, our approach gener aliz es to multip le ROIs. Vari ab les X i ( s ) = X i consist of an interc ept, AS D status, age, sex, and the age by ASD status interaction.

We assume that α i (•) is a Gaus si an proces s with mean function µ{•; X i (•) , β} and positiv e-definite c ova riance function C α {•, •; X i (•) , γ} . We further assume that µ{S ; X i (S ) , β} is known up to a pa ra mete r β ∈ R q 1 , a nd C α {S , S ; X i (S ) , γ} is known up to a pa ra mete r γ ∈ R q 2 :

The distribution of y i (S ) is S -v ari ate Gaus si an with mean µ{S ; X i (S ) , β} a nd cova ria nce C {S , S ;

r,t=1 + σ 2 I S . We define θ = ( β, γ, σ 2 ) ∈ R p the pa ra mete r of in te res t, p = q 1 + q 2 + 1 . The s ta nda r diza t ion of the rfMRI t ime se ries jus tifies our assumption of a homogeneous σ 2 across i and s . We allow the cov ari ance to depend on X i (S ) and to exhibit some nonst ation arity. In Se ction 4 , C α borr ows fr om Pacior ek a nd Sche rvish ( 2003 ) and t ake s the form

where r i j = exp { X ⊤ i ρ( s j ) } models the effect of X i on the corr ela tion betw e en v oxe ls throu gh ρ( s j ) ∈ R q . Both ρ( s j ) and τ ( s j , s j ′ ) are spec i fied in Section 4 .

A Shift from Whole to Loc a l Data Perspectives

We propose a rec ur sive part it ion of S in to m ultip le res o lution s, with multip le s e ts a t each r es o lution . Le t ⊗2 denote the outer product of a v e ctor with itse lf, name ly a ⊗2 = a a ⊤ . Using the notation of Katzfuss ( 2017 ), denote A 0 = S and part it ion S into K 1 (dis j oint) region s {A 1 , . . . , A K 1 } that are again part it ioned into K 2 (dis j oint) s ubre gion s {A k 1 1 , . . . , M, where K 0 = 1 . For comp le tenes s, le t k 0 = 0 . An examp le of a rec ur sive part it ion of o bs erv ation s on a 2-dimensional spatial domain for M = 2 res o lutions is provided as supplemen ta ry mate ri als . For res o lution m ∈ { 1 , . . . , M} , denote S k 1 ... k m the size of A k 1 ... k m , with S 0 = S . Further define k max = max k m with the maximum taken over k m = 1 , . . . , K m and m = 1 , . . . , M. Due to that fact that

We further re c ommend that values of M and K m be chosen so that S k 1 ... k M and K m , m = 1 , . . . , M ar e r ela tiv ely sm all c ompare d to S 0 , and w e ge ne rally re c ommend S k 1 ... k M ≥ 25 .

This ide a appe ars similar to MRA models (Katzfuss, 2017 ;Katzfuss and Gong, 2020 ;Nychka et al., 2015 ). Where they int egrat e from global to local lev els, how ev er, w e build local models at the hi ghes t res o lution and rec ur sively int egrat e infe re nce from the hi ghes t to the lowe st re s o lution . Unlike MRA model s, we al s o incorporate un s tructured depe nde nce within all spatial res o lution s using the ge ne r aliz ed method of mome n ts (GMM) (Ha n s e n, 1982 ) fra me work, de scribed in Section 2.3 .

A t res o lution M for k M ∈ { 1 , . . . , K M } , the l ikel ihood of the data in set A k 1 ... k M is given by y

We model the mean of the spatial process through

The function h is ass ume d known and wi l l most often be the linear function

In Section 3 , w e inv es ti gate the pe rforma nce of our method when h is the linea r a nd lo g function s . In Se ction 4 , µ{ s j ; X i ( s j ) , β} = β for s j ∈ A k 1 ... k M due the ce n te rin g of ea ch time se ries, a nd C is the nons tationa ry cova ria nce function in ( 1 ). Due to the rec ur sive part it ioning, S k 1 ... k M is small and the l ikel ihood tractable. Denote the log-l ikel ihood in

The litera tur e is r ep le te with me thods for choosing part it ions; s ee Heaton e t al. ( 2019) for a n excelle n t re view. T he cov ari ance model ne e d only be spe c i fied within each block at the hi ghes t res o lution . We thus re c ommend th at pract it ioner s rec ur sively part it ion the domain such that the cov ari ance model within each block at the hi ghes t res o lution can be easily spec i fied, consiste n tly a nd tract ably e stim ate d, a nd homoge neous acros s b locks in the s en s e th at the outc ome within e ach block is g e ne rated from the same model. In Section 4 , we rec ur sively part it ion voxels in the left and right pre c entr al gyri separ at ely int o K 1 = 2 , K 2 = K 3 = 4 dis j oint s e ts bas ed on nea res t nei ghbors, M = 3 ; the s e ts A k 1 ... k m con si st of the union of the di s j oin t pa rt it ion s e ts in the 2 gyri at each res o lution m ∈ { 1 , . . . , M} . This strategy ens ures th at each A k 1 ... k M contains locations from both gyri, so th at gyri-spe c i fic pa ra mete rs a re ide n tifiab le in all s e ts A k 1 ... k M , k M = 1 . . . , K M . As a nothe r exa mple, a pract it ioner may choose to part it ion their dom ain s uch th at the outc ome in each block at the hi ghes t res o lution has a simp le local cov ari ance structure, such as a Gaus si an or exponential cova ria nc e function (se e Se ction 3 for examples); if only the mean is of in te res t , with the cov ari ance tr ea ted as a n uisa nce, the n a part it ion with independent o bs erv ation s within each block at the hi ghes t res o lution would be practical as the model is easy to spec i fy a nd es timate, a nd the depe nde nc e betw e en blocks is ac c ounte d for in the c ombin ation steps . The ass umption is then th at the local c ov ari ance structure within all blocks at the hi ghes t res o lution is the same (exten sion s are d isc ussed in S ection 5 ). T he unstructured cova rianc e betw e en thes e b locks is ca ptured b y the GMM, d isc ussed next.

Recu rs ive Model Integration

GMM Appr oac h

We now wish to rec ur sively int egrat e the local models at each res o lution . For m = M, . . . , 0 , let k m ∈ { 1 , . . . , K m } denote the index for the s e ts at res o lution m . We describe a rec ur sive GMM appr oach tha t r e-estima tes θ a t each r es o lution M -1 bas ed on an optimally wei gh ted form of the GMM equation from the hi ghe r res o lution M. A t res o lution M -1 , we define K M-1 est imat ing funct ion s bas ed on the s core function s fr om r es o lution M:

, ther e ar e several est imat ing funct ions for each pa ra mete r. Fo llowing Han s en's ( 1982 ) GMM, we propose to minimize the quadratic form

, with a positive se mi -definite wei gh t matrix W . Ha n s en ( 1982 ) show e d th at the optim al choic e of W is the inv erse of the c ov ari ance of k 1 ... k M-1 ( θ) , which can be consistently estim ate d by

(2) the v ari ability matrix, and

e stimate s de pendenc e betw e en sc or e functions fr om all s e ts in A k 1 ... k M-1 and thus captures depe nde nc e betw e en sets

. Thi s choice i s "Han s en" optimal in the s en s e that, for all pos sib le choices of W ,

has minim um va ria nce a mong es timators a rg min

) , subs ta n tial ly smal le r tha n inve rting the S k 1 ... k M-1 × S k 1 ... k M-1 dimension al m atrix in the dir ect evalua tion of the l ikel ihood on s e t A k 1 ... k M-1 . Thi s yield s a fas te r c omputation th an a ful l li keli hood approach on A k 1 ... k M-1 . The itera tive minimiza tion in ( 3 ), ho weve r, ca n be time-consuming because it r equir e s comput ation of the score functions k 1 ... k M at each iteration.

In the spirit of Hector and Song ( 2021), a closed-form met aestimator as ymptotically equivale n t to θ GMM in ( 3 ) tha t is mor e c omputation ally attractiv e is given by

where

) can be efficiently c ompute d by Ba rtlett's ide n tity.

When updating to the next res o lution M -2 , it is tempting to proc e e d again through the GMM approach and to st ack e st imat ing funct ions

. This would, how ev er, res ult in a v e ctor of ov eridentifie d estim ating functions on θ of dimension pK M-1 K M , and ther efor e inversion of a (pK M-1 K M ) × (pK M-1 K M ) dimension al c ova ria nc e m atrix. Rec ur sively proceed ing this way would lead to the inversion of a large (pK 1 . . . K M ) × (pK 1 . . . K M ) dimensional ma trix a t r esolution m = 0 , negating the gain in computation afforded by the part it ion. To avoid this d iffic ulty, we propose wei gh ts in the spirit of opt imal est imat ing funct ion theory ( Heyde , 1997 ). These optim ally w ei gh t GMM es timating functions and re duc e their dimension for c omputation ally tractable re c ur siv e inte gration of depe nde n t models.

Weigh ted OverIden t ified Est imat ing Funct ions

We define new est imat ing funct ion s at res o lution m = M -2 th at optim ally w ei gh t the es t imat ing funct ions fr om r es o lution M -1 . Let

) a r e r e c omputed so as to evaluate the s en sitivity and v ari ability at the estima tor fr om r es o lution M -1 . This form ulation ca n also be view e d as the optim al proje ct ion of the est imat ing function k 1 ... k M-1 ( θ) onto the p arameter sp ace of θ ( Heyde , 1997 )

, a pK M-1 dimension al v e ctor of ov e ride n tifying es t imat ing functions on θ.

, one can again define the clos ed-form me t a -e stimator θ k 1 ... k M-2 in the fashion of ( 4 ). This r equir es inversion of a (pK M-1 ) × (pK M-1 ) -dimension al m atrix, subs ta n tial ly smal le r tha n inve rting the

Recurs ive In tegrati o n P roced ure

The rec ur sive integra tion pr ocedur e upda tes thr ough Sections 2.3.1 and 2.3.2 : for each k m = 1 . . . , K m , m = M -2 , . . . , 1 , do steps 1 and 2, described next. In step 1, rec ur sively compute

In step 2, compute

where

A t res o lution m = 0 , we obtain ψ i, 0 (

where

can be c ompute d following the rec ur sive pr ocedur e. A to y exa mple i l lustra ting the pr oc e dure i s g iv en in the s upple me n t. One evalua tion a t r esolution M has computation and memory complexities of

Fin ally, the c omputa tion of θ k 1 ... k m r equir es no itera tive minimization of an o bj e ctiv e function, s ubs ta n tially reducing computation al c osts . The proc e dure can be fully run on a distributed sys te m, mea ning tha t a t no point do the data ne e d to be loade d on a ce n tral se rve r or the full S × S cov ari anc e m a trix stor ed.

While the n umbe r of locations S k 1 ... k M in blocks at the highe st re s o lution driv es the c omputation al c omplexity at the distribut ed st ep, the n umbe r of blocks K m at each res o lution drives the c omputation al c omplexity at the c ombin a tion step thr ough k max . S e veral values of K m and M m ay achiev e similar values of S k 1 ... k M , and so a pr actical consider ation concerns the choice of thes e v alues. First, we recommend choosing K j close to d ivisor s of S so as to achieve a similar number of locations in each block. Next, w e re c ommend choosing values s uch th at neither K j nor M is exc essiv ely large or small. For exa mple, the to y exa mple in the supplement part it ions a 10 × 10 grid of locations using

Finally, we show in Section 3 that differe n t choices of M and K m lead to similar perform anc es .

Mult iresolut ion Est imat ing Fu n ct ion Th eory

Let

the pa ra mete r space of θ, and denote by θ 0 the true value of θ, defined formally by as sumption s as supple me n ta ry material s . In thi s se ction, w e s tudy the asymptotic prope rties of θ r by formalizing a multires o lution estimating function theory. To do this, define popul ation version s of the est imat ing function s, their v ari ability, and their s en sitivity: for

are the v ari ability, s en sitivity, a nd Goda mbe information (Goda mbe, 1991 ) m atric es, respe ctiv ely, in

, and de-

. We assume throughout that v k 1 ... k m ( θ 0 ) and v 0 ( θ 0 ) are positive definite.

Under mild as sumption s on the score functions k 1 ... k M ( θ) , the K M estimators θ k 1 ... k M from the Mth res o lution are consis te n t for θ 0 and s emi-parame trically efficient within each

More ov e r, unde r a ppr opria te condi-tion s, θ k 1 ... k M-1 is con sis te n t a nd Ha n s e n optimal b y Hector a nd Song ( 2021 ). Finally, at each res o lution, a n es timator is derived from the optimal GMM function (Han s en, 1982 ) and from opt imal est imat ing funct ion theory ( Heyde , 1997 ). This results in a hi ghly efficie n t in t egrat ed estimat or θ r in ( 6 ) that fully uses the depe nde nce within a nd betw e en each res o lution, as shown in Theorem 1.

T heore m 1 Under assum pt io ns given in the supplement as N → ∞ , θ r in ( 6) is co nsis t ent and

The pr oof pr oc e e ds b y induction afte r es tab lishing con siste ncy a nd asymptotic normality of the inte grate d estim a tors a t res o lution M -1 . It uses a rec ur sive Taylor expansion at each res o lution to estab lish the appr opria te c onv ergenc e rate of the est imat ing funct ion . Large s amp le c onfidenc e intervals for θ 0 can be c onstructe d by c ombining The orem 1 and the following Coro ll ary, whos e proof follows from the proof of Theorem 1 and is omitted.

Coroll ary 1 Under assum pt io ns given in th e supplem ent as N → ∞ , J -1 0 ( θ r ) is a co nsis t ent es t ima t o r of the asympt ot ic cova ri a nce of θ r in ( 6 ).

Seque ntial Mode l Integratio n Fra mewo rk

The most time-c ons uming step in the rec ur sive integration proc e dur e r equir es a r ec ur sive u pdate of the wei gh ts each time a new estimator of θ is computed a t r es o lution m = M -2 , . . . , 0 , as i l lustrated on the left of Figure 1 . This is because the wei gh t matrices V k 1 ... k m ( θ k 1 ... k m +1 ) and S k 1 ... k m ( θ k 1 ... k m +1 ) are evaluated at θ k 1 ... k m +1 . We propose an alt ernat e seque n tial in tegration sche me tha t evalua tes the wei gh t at the es tima tor fr om the Mth r es o lution, θ k 1 ... k M , a consis te n t es timator of θ, i l lus trated on the ri gh t of Figure 1 . The s equenti al integra tion pr ocedur e r ep l aces the rec ur sive evaluation of the wei gh ts with a seque n ti al upd ate by computing V k 1 ... k m ( θ k 1 ... k M ) and S k 1 ... k m ( θ k 1 ... k M ) . We denote by θ s the seque n tial in t egrat ed estimat or obtained by evaluating wei gh ts using

Arguing as in Section 2.3.3 , the computation and memory complexities of computing

, respe ctiv ely. A c on s e quenc e of the proof of Theorem 1 is that θ s and θ r are asymptotically equivale n t. While θ s may be computationally more adva n t age ous, it re quires th at the MLEs from res o lution M be close to θ 0 , and may not perform as well as θ r in finite s amp les.

SIMUL ATIONS

We inves ti gat e the finit e s amp le perform anc e of the proposed rec ur siv e and se quential MR R I estimators θ r and θ s . Unless otherwise spec i fied , all simulat ions are on a Linux cluster using R linked t o Int e l's MKL librarie s with analyse s a t r es o lution M performed in parallel across K = M m =1 K m CPUs with 1 GB of FIGUR E 1 E xa mple pa rt it ion of A 0 with M = 3 res o lution s . Schem a tic pr opaga tion ste ps for comput a tion of θ k 1 ... k m for r ec ur sive (left) and seque n tial (ri gh t) in tegration sche mes.

RAM each . Stand ar d err ors and c onfidenc e intervals are calculated using the results in Theorem 1 and Coro ll ary 1. In the first s e t of simul ation s, w e c onsider a S = 400dimen sional s quar e spa tial domain S = [1 , 20] 2 = { s j } 400 j=1 , s j ∈ R 2 , with N = 10000 . The Gaus si an outcomes { y i (S ) } N i =1 a re indepe nde n tly sim ulated with mea ns { X

and Gaus si an spati al cov ari ance function

a nd n ugge t v ari anc e σ 2 spe c i fie d below. The c ov ari ates X i ( s j ) = X i consist of an intercept and 2 non spati ally v arying con tin uous va riable s inde pe nde n tly ge ne ra ted fr om a Gaus si a n dis tribution with mea n 0 a nd va ria nce 4. The true value of the pa ra mete rs is β = (0 . 3 , 0 . 6 , 0 . 8) , σ 2 = 1 . 6 , τ 2 = 3 and ρ 2 = 0 . 5 . To facilit ate e st imat ion, we est imate θ = { β ⊤ , log (τ 2 ) , log (ρ 2 ) , log (σ 2 ) } using θ r and θ s . We consider 5 rec ur sive part it ions of S . In Setting I, M = 3 with

Compa rison of infe re nce a nd run times across Settings I, II, and III shows the impact of the ordering of the rec ur sive integration on the proposed proc e dure, while Settings IV and V inves ti gate the impact of M and k max on infe re nce a nd run t ime. Division of the spat ial domain is based on nea res t nei ghbors, as i l lus trated in the exa mple of the supplement . We also estimate θ using 2 c omparativ e approaches . Spe c i fically, w e c ompare to the part it ionin g approa ch (Part.) in Heaton et al. ( 2019 ) tha t evalua tes the sum of the K log-l ikel ihoods at a grid of values of τ 2 , ρ 2 , estimates τ 2 , ρ 2 using the values tha t r eturn the hi ghes t log-l ikel ihood, then estimates β and σ 2 using the least squares es timator a nd a s amp le v ari anc e, respe ctiv ely; the imple me n tation is modified directly from the code in Heaton et al. ( 2019 ). Asymptotic s ta nda rd errors (ASEs) are estim ate d as the diagonal square root of the invers e v ari ance of the score function. We also compare to the NNGP (Finley et al., 2019 ) using the 25 nea res t nei ghbors; we imple me n t this ourse lve s using sparse m atric es in the R package "Rcpp" and par alleliz e over K CPUs to make a fair comparis on . Figure 2 plots the root mean squared error (RMSE), empirical s ta nda rd e rror (ESE), a nd AS E av erage d acros s 500 simul ation s for each pa ra mete r.

The near equality of RMS E and ES E in Figure 2 for our MR R I es timators θ r a nd θ s i l lus trates the nea r unbi as ednes s of our approach in large s amp les. Further, the near equality of ES E and AS E justifies the use of the asymptotic v ari ance form ula in Theore m 1 a nd Coro ll ary 1 in l arge s amp les . Fin ally, negligib le v ari ation is o bs erv e d acros s Se ttings I-V, i l lustrating the ro bustnes s of our approach to the chos en rec ur sive or seque n tial in tegration sche me. The s tatis tical infe re nce prope rties of NNGP are simil arly favorab le, with the NNGP a ppea ring sli gh tly more efficie n t. The pa rt it ionin g approa ch yields estimators of β that are nearly unbi as ed in l arge s amp les, whereas the estimators of the cov ari ance parame ters show bi as. Stand ard errors of the estimators are ov erestim ate d for β and underestim ate d for the cov ari ance parame ters, which is pro b lematic for s tatis tical infe re nce. Thi s i s further i l lustrated in Table 1 , which reports the 95% c onfidenc e interval c ov erag e (C P) for all e stimators av erage d across 500 simul ation s. Where as C P re aches nomin al lev els for MR R I and NNGP, r egr es sion and cov ari ance para mete rs a r e, r espe ctiv ely, ov erc ov ere d and underc ov ere d using Part.

Comput ing t imes ar e r eport ed in Table 2 . Using θ s is fast er th an re c ur sively u pdating wei gh ts with θ r . In compa rison, the part it ionin g and NNGP approa ches are 212 and 467 times slow er, respe ctiv ely, th a n our es timator θ s with K 1 = 4 , K 2 = K 3 = 2 . We also observe the increased computational burden in Settings IV and V. This phe nome non is unde rpinned b y the c omputation c omplexity orde rs of θ r a nd θ s de riv e d in Se ctions 2.3.3 and 2.5 , which show that computation is slower when M or k max are larger.

A se c ond set of simulations with the same model and S = 25600 is prese n te d in the s upplemen t . The re, w e show th at our MR R I estimators scale well when S gr ows. A thir d s e t of simul ations in the supplement inves ti g ates the log-line ar model with h the log function, S = 900 , and C α ( s j , s j ′ ; γ ) the exponential cov ari ance function . Conclusion s drawn from this simulation mirror those from the first 2 s e ts of simul ation s.

The fourth s e t of simul ation s mimics the d ata an alysis of Se ction 4 with S = 800 . We consider 2 ROIs,

FIGUR E 2 R MS E, ES E, and AS E for estim ates of each pa ra mete r using NNGP, Part., θ r , and θ s in Settings I-V (r -224 , r -242 , r -422 , r -28 , r -2222 , s -224 , s -242 , s -422 , s -28 , and s -2222 , respe ctiv ely) for the first s e t of simul ation s.

TABLE 1 95% CP in pe rce n tage for the first s e t of simul ation s:

2 , 2 , 2 θ r 0.80 (0.092) 0.66 (0.060) 0.66 (0.060) 1.2 (0.86) 0.97 (0.75) θ s 0.72 (0.084) 0.60 (0.064) 0.60 (0.064) 1.2 (0.82) 1.1 (1.2) NNGP 276 (135) Part. 127 (11)

Across 500 simulations for the first set of simulations: S = 400 , θ = { 0 . 3 , 0 . 6 , 0 . 8 , log (3) , log (0 . 5) , log (1 . 6) } .

TA BLE 3 Simulat ion metrics of θ s across 500 simulations for the fourth set of simulations: S = 800 , θ = { 0 , log (3) , 0 . 5 , 0 . 5 , 0 . 5 , 0 . 6 , 0 . 6 , 0 . 6 , log (1 . 6) } .

Pa ra me ter

1.9 1.9 1.9 -0 .57 94 ρ 12 1.9 1.9 1.8 -0 .88 94 ρ 13 1.9 1.8 1.8 -1 .80 93 ρ 21 1.9 1.9 1.9 -0 .35 94 ρ 22 1.9 1.9 1.9 -1 .20 95 ρ 23 1.9 1.9 1.9 -0 .49 94 log (σ 2 ) 1.2 1.2 1.2 -0 .51 96 v e ctor of one's, and the spati al cov ari ance function in ( 1 ) with d = 2 . The spati al v ari ance is mode led throu gh τ ( s j , s j ′ ) = τ 2 , the spa tial corr ela tion is mode led throu gh ρ(

He re, X i ∈ R q consis ts of a n in te rcept a nd 2 non spati ally v arying con tin uous va riable s inde pe nde n tly ge ne ra ted fr om a Gaussia n dis tribution with mea n 0 a nd va ria nce 1. The true values of the depe nde nce pa ra mete rs a re s e t to σ 2 = 1 . 6 , τ 2 = 3 , ρ 1 = (0 . 5 , 0 . 5 , 0 . 5) , and ρ 2 = (0 . 6 , 0 . 6 , 0 . 6) . We estimate θ = { β, log (τ 2 ) , ρ ⊤ 1 , ρ ⊤ 2 , log (σ 2 ) } using θ s using the rec ur sive part it ion of S with M = 3 : K 1 = K 2 = 2 , and K 3 = 4 . Each s e t A k 1 ... k m at res o lution m is a union of the nea res t nei ghbors in S 1 and S 2 se parate ly, so that each set A k 1 k 2 k 3 consists of 25 loca tions fr om S 1 and 25 loca tions fr om S 2 , A k 1 k 2 consists of 100 loca tions fr om S 1 and 100 locations from S 2 , and s o on . Analyses a re pa r alleliz ed across 16 CPUs with 2 GB of RAM each. Table 3 reports the RMSE, ESE, ASE, mean bias (BIAS), and CP of our MR R I estim ator θ s av erage d acros s 500 simul ation s, where

Parallelizin g o ver K = 16 C PUs, me an el aps ed time (Monte Ca rlo s ta nda r d devia tion) is 29 minutes (22 minutes). Simul ation me trics in Tab le 3 are consistent with the results from the previous simul ation s: the RMSE, ESE, and ASE a re a pproxim ately e qual, the BIAS is ne gli gible, a nd the C P re aches the nominal 95% level. Further, to i l lustrate the hi gh s tatis tical powe r of our approach, we perform a test of the hypothe se s H 0 : ρ 12 = ρ 22 v ers us H A : ρ 12 ̸ = ρ 22 using the test stat ist ic Z = ( ρ 12ρ 22ρ 0 ) { V ( ρ 12 ) + V ( ρ 22 ) -2 Cov ( ρ 12 , ρ 22 ) } -1 / 2 , which follows an appr oxima te standar d normal distribution under H 0 when ρ 0 = 0 . In the context of the analysis of Section 4 , this t est evaluat es whether ASD is as s oci a ted with differ ent spa tial corr ela tion in the left and right pre c e n tral gyri. Across the 500 simulations, the test rejects the null 100% of the time at level 0.05, i l lus trating the hi gh s tatis tical power of our approa ch. Usin g ρ 0 = ρ 12ρ 22 , the type-I error rate is 3.6% across the 500 simulations.

ESTIMATION OF BR AIN FUNCTIONAL CONNECTIVIT Y

We return to the mot ivat in g neuroimagin g app lication des cribed in Section 1 . Out of 1112 A BIDE partic ipa n ts, 774 pas s ed quality control: 379 with ASD, 647 males (335 with ASD) and 127 females (44 w ith A SD), w ith me an ag e 15 ye a rs (s ta nda rd deviation 6 years). The left and right pre c entral gyri form the prim ary motor c ortex and ar e r espon sib le for executing volun ta ry move me n ts. They a re two of the la rges t ROIs in the Ha rva rd -Oxford atl as (FMRIB Softw a re Libra ry, 2018 ) with 1786 and 1888 voxe ls, re spective ly. Many indiv iduals w ith A SD h av e motor deficits ( Jansiewicz et al., 2006 ). Atypical c onne ctiv ity w ithin the left and right pre c entral gyri may indicate that these motor deficits ar e r ela t ed t o how these 2 brain regions c oordin ate mov e me n t. The pre proce ssing pipe line of rfMRI da ta has alr eady be en describe d by Craddock et al. ( 2013 ). Pa rticipa n t-spec i fic data have been reg i s te red in to a common temp l ate space such that voxel locations are comparable betw e en p articip ants.

Define S 1 = { s j } 1888 j=1 and S 2 = { s j } 3674 j=1889 the s e t of voxels in the ri gh t a nd left pre c e n tral gyri , respect ively, and S = S 1 ∪ S 2 . Voxels in the atlas outside the brain are assigned mis sing. Fo llo wing the thinning described in Section 2.1 , we obtain indepe nde n t r eplica tes y i ( s j ) at S = 3674 voxel locations, i = 1 , . . . , 75888 . We refer to this d atas e t as the "primary" d atas e t; w e c ompare estima tes fr om this primary da t aset to e stimate s from a se c ondary d atas e t, con sisting of the excluded time poin ts, with ide n tical dimen sion s N and S . There is a pr ior i no reason to believ e th at the data distribution differs across primary and se c ondary datasets,

TA BLE 4 Est imated cov ari ate effects and s.d. in the primary and se c ondary datasets . ρ 1 ρ 2 Covariate Estimate s.d. ×10 4 Estimate s.d. ×10 4 Primary d atas e t In te rcept 0 .569 1 .44 0 .561 1 .42 ASD status -0 .00676 1 .32 -0 .0050 1 .27 Age -0 .0222 1 .07 -0 .0120 0 .932 Sex 0 .0362 1 .79 0 .0655 1 .72 Age by ASD status interaction -0 .00310 1 .38 -0 .00386 1 .25 Se c ond ary d atas e t In te rcept 0 .568 1 .43 0 .562 1 .41 ASD status -0 .00511 1 .32 -0 .00472 1 .27 Age -0 .0211 1 .07 -0 .0118 0 .931 Sex 0 .0341 1 .78 0 .0643 1 .72 Age by ASD status interaction -0 .00503 1 .38 -0 .00591 1 .25

and so comparing results across both datasets wi l l al low us to qua n tify robus tness of our analysis, a notoriously d iffic ult task in analyses of rfMRI data (Uddin et al., 2017 ). Let X i be corresponding observations of q = 5 covariates for outcome i : an interc ept, AS D status (1 for ASD, 0 for neurotypical), age (ce n te red a nd s ta nda rdized), sex (0 for m ale, 1 for fem ale), and the age (c entere d and standardize d) b y ASD s tatus in te raction. We model µ( s j ; X i , β) = β, where µ( s j ; X i , β) is the mean rfMRI time series at voxel s j . The cov ari ance C α ( s j ′ , s j ; X i , γ ) is that given in Equation ( 1 ) of Section 2.1 with d = 3 and τ ( s j , s j ′ ) = { ( τ 2 1 ) 1 ( s j ∈S 1 )+ 1 ( s j ′ ∈S 1 ) ( τ 2 2 ) 1 ( s j ∈S 2 )+ 1 ( s j ′ ∈S 2 ) } 1 / 2 . As in Section 3 , the corr ela tion is modeled thr ough ρ( s j ) = ρ 1 1 ( s j ∈ S 1 ) + ρ 2 1 ( s j ∈ S 2 ) . Thus, γ = { log ( τ 2 1 ) , log ( τ 2 2 ) , ρ ⊤ 1 , ρ ⊤ 2 } ∈ R 2+2 q and θ = { β, γ, log (σ 2 ) } , with γ the parameter of intere st de scribing the effect of cov ari ates on function al c onne ctivity. We in te rpret the effect of ASD on the corr ela tion structur e in detail in the supplement . An illustration of the correlation betw e en ROIs for various values of ρ j2 , ρ j ′ 2 is provided in the supplement .

The size of each (primary and se c ondary) outc ome d atas e t is 15 GB. To ov erc ome the c omputation al burden of a whole d atas e t analysis, we estimate θ using the seque n tial es timator θ s . To part it ion the 3-dimensional spat ial dom ain, w e re c ur sively part it ion S 1 and S 2 se parate ly into K 1 = 2 , K 2 = K 3 = 4 disj oint s e ts bas ed on nea res t nei ghbors, M = 3 . The s e ts A k 1 ... k m consist of the union of the dis j oin t pa rt it ion s e ts of S 1 and S 2 at each res o lution m ∈ { 1 , . . . , M} . This s trategy e ns ures th at each A k 1 ... k M contains loca tions fr om both S 1 and S 2 , so that gyri-spec i fic parameters ρ 1 , ρ 2 , τ 1 , τ 2 a re ide n tifiab le in all s e ts A k 1 ... k M , k M = 1 . . . , K M . A plot of the part it ionin g is pro vided in the supplement .

The analysis of the prima ry a nd s econd ary d atas e ts takes 4.28 and 3.62 hours, respe ctiv ely, whe n pa r alleliz ed across 32 CPUs with 4 GB of RAM each; a comparison of these run times is detailed in the supplement . We a t t empt ed t o obtain the MLE using either the whole primary or se c onda ry dataset b y supplying the value of θ s as a s ta rting value. The est imat ion re quire d more than 170 GB of RAM and t erminat ed without comp le ting after 190 hours. The s ta nda r d devia t ion (s.d .) of block-spec i fic es ti -mates of log (σ 2 ) across the K 1 K 2 K 3 = 32 blocks is 0.401 and 0.399 in the prima ry a nd se c ondary da tasets, r espe ctiv ely, s upporting our assumption that σ 2 does not depend on s ; comp le te s umm a ries a re av ail ab le in the supplemen t . The es tim ate d effe ct and s.d. of the interc ept, AS D status, age, sex, and age by ASD s tatus in te raction a re reported in Table 4 . In the primary d atas e t, the e stimate s (s.d.) of β, log (σ 2 ) , log (τ 2 1 ) and log (τ 2 2 ) are, respe ctiv ely, -0 . 00537 ( 2 . 03 × 10 -4 ), -4 . 05 ( 3 . 90 × 10 -4 ), -0 . 108 ( 2 . 53 × 10 -4 ), and -0 . 0699 ( 2 . 48 × 10 -4 ). In the se c ondary dataset, the e stimate s (s.d.) of β, log (σ 2 ) , log (τ 2 1 ) and log (τ 2

2 ) ar e, r espe ctiv ely, -0 . 00609 ( 2 . 03 × 10 -4 ), -4 . 05 ( 3 . 88 × 10 -4 ), -0 . 111 ( 2 . 51 × 10 -4 ), and -0 . 0709 ( 2 . 46 × 10 -4 ). As expe cte d , est imates of β are close to 0 a nd es timates of σ 2 + τ 2 j are close to 1 due to the ce n te ring a nd s ta nda rdizing. The cosine simil arity be twee n s ta nda rdized es timates of θ in the prima ry a nd s econd ary d atas e ts, that is, the cosine of the angle betw e en the 2 standardize d v e ctors θ q / { V ( θ q ) } 1 / 2 , q = 1 , . . . , 14 , is 0.999994, indicating a high de gre e of agre e me n t betw e en the standardized vectors.

From a practical perspe ctiv e, estim ates and their standard e rrors a re virtually ide n tical across prima ry a nd se c ondary d atas e ts . Tw o-s amp le Z-tests, how ev er, mostly reje ct the null that ele me n ts of θ are equal in both d atas e ts at the typical 0.05 level: thi s i s a fea tur e of the s amp le size and the high power of the tes t, rathe r tha n of true unde rlying diffe re nc es betw e en the 2 d atas e ts, and i l lustrate s we l l the chal le nges of robus tness in a nalyses of rfMRI data. We calibrate the α-level of hypothesis testing proc e dures in our a nalysis b y borrowing ideas from knock-offs (Ba rbe r a nd Ca ndès, 2015 ). We e stima te a da t a -de pe nde n t type-I err or ra te thr eshold as the 5% qua n tile of the o bs erv e d p-values of the 2-s amp le tests of equality betwee n pa ra mete rs in prima ry and se c ondary datasets . The Gaussian critical value corresponding to this 5% qua n tile is z crit = 9 . 87 .

Armed with this robust critical value, we evaluate whether the AS D m ain a nd in te raction effects a re si gnifica n tly diffe re n t across the 2 brain re gions . We perform a test of the hypothe se s H m 0 : ρ 12 = ρ 22 v ers us H m A : ρ 12 ̸ = ρ 22 and H i 0 : ρ 15 = ρ 25 versus H i A : ρ 15 ̸ = ρ 25 using the tes t s tatis tic Z in Section 3 , where ρ j2 and ρ j5 are the AS D m ain a nd in te raction effe cts, respe ctively, in ROI S j , j = 1 , 2 . The tes t s tatis tic for H m 0 v ers us H m A t ake s a value of 13.99 and 3.04 in the primary and se c ondary d atas e ts, respe ctiv ely; sinc e w e reje ct H m 0 in the prim ary but not the se c ondary dataset, w e c onclude th at the m ain AS D effe ct is not si gnifica n tly diffe re n t in the corr ela tion structur e of both brain re gions . The tes t s tatis tic for H i 0 v ers us H i A t ake s a value of 6.14 and 7.13 in the primary and se c ondary da tasets, r e spectiv ely, and w e c onclude th at the AS D b y age in te raction effect is not si gnifica n tly diffe re n t in the corr ela tion structur e of both brain re gions . An an alysis th at excludes the age by AS D status in te raction is included in the supplement and agrees with this an alysis .

Summa ries of dis ta nces d j j ′ = ( s j ′ -s j ) ⊤ ( s j ′ -s j ) ar e pr ovided in the supplement . For a male p articip ant of mean age (15.13 years), the estim ate d c orre lation structure s betw e en the ri gh t a nd left pre c e n tral gyr i, within the r i gh t pre c e n tral gyrus, and within the left pre c entral gyrus, for a p articip ant with and w ithout A SD ar e, r espe ctiv ely, 0 . 915 exp (-1 . 83 d j j ′ ) and 0 . 915 exp (-1 . 86 d j j ′ ) , 0 . 898 exp (-1 . 85 d j j ′ ) and 0 . 898 exp (-1 . 89 d j j ′ ) , and 0 . 932 exp (-1 . 82 d j j ′ ) and 0 . 932 exp (-1 . 84 d j j ′ ) . AS D m a nifes ts as hype rc onne ctivity within and betw e en the right and left pre c entral gyri. These find ings conc ur with those of Nebel e t al. ( 2014 ). Thes e results also concur with a less pow ere d an alysis th at av erag es e ach p articip ant's rfMRI time series at each voxel, then r egr es s es the p articip ants ' Pearson corre lation betw e e n the ri gh t a nd left pre c e n tral gyri onto an interc ept, AS D status, age, sex, and the age by ASD status interact ion. Est im ates (s .d.) of c ov ari ate effects from this analysis are 0.726 (0.00759), -0 . 0457 (0.0101), 0.0159 (0.00745), -0 . 0111 (0.0137), and -0 . 00293 (0.0101). Only the in te rc ept, AS D s tatus, a nd age effects are si gnifica n t at level 0.05, hi ghli gh ting the supe rior powe r of our spatial approach.

DIS CUSS ION

The proposed MR R I e stimators de pend on the choice of the rec ur sive part it ion of S. We h av e s u gge ste d th at this part it ioning be performe d s uch th at S k 1 ... k m and K m , m = 1 , . . . , M ar e r e lative ly s mall compared to S 0 and shown through simul ation s that thi s lead s to de sirable st at ist ical and c omputation al perform anc e. Nonethe le ss, the GMM is known to unde res timate the va ria nce of estimators when K m is moderately large r ela tive to N (Han s en et al., 1996 ). Thus, special care should be taken to ensure K m is r ela tively lo w -dimensional.

In this pa pe r, w e h av e allow e d θ to v ary spati ally under the cons train t that it can be consis te n tly es timated using subs e ts of the spatial domain S. While we have assumed that σ 2 does not depend on s , a model that allows σ 2 to depend on s can easily be fitted in our framework under this cons train t. Othe r settings m ay c on sider a s e ttin g in which ea ch subs e t is mode led throu gh its own, subs e t-spec i fic pa ra mete r. Futur e r esear ch should focus on the developme n t of rec ur sive and seque n tial in tegration rules for this s e tt ing with part i ally he tero geneous parame ters θ fo llo wing, for example, the work of Man s chot a nd Hector ( 2022 ) a nd Hector and Reich ( 2024 ) in spatially varying coefficient models. Such exten sion s would allow pa ra mete rs such as the n ugget v ari ance σ 2 to vary with s . A pr ior i , thes e exten sion s should fo l-low from zero-pa ddin g the wei gh t m atric es in Equations ( 4) and ( 6 ).

While our approach was motivated by a comparison of function al c onne ctivity betw e e n pa rticipa n ts with a nd w ithout A SD, the proposed methods are generally app licab le to Gaus si an proce ss mode l ing of high-d imension al im ages . Our ABIDE an alysis h as prim arily focuse d on modeling c orr ela tion betw e en the left a nd ri gh t pre c e n tral gyri as a function of the within-ROI corre lations . This an alysis is s uitable when within-ROI corr ela tion is of prima ry in te res t. Exte n sion s that model the between-ROI correla tion thr ough a mor e flexib le cov ari ance structur e ar e of inter est but beyond the scope of the prese n t work.