Using recursive algorithms for the efficient identification of smoothing spline ANOVA models

22 Pages·2010·0.72 MB·English

AStAAdvStatAnal(2010)94:367–388
DOI10.1007/s10182-010-0148-8
ORIGINAL PAPER
Using recursive algorithms for the efficient
identification of smoothing spline ANOVA models
MarcoRatto·AndreaPagano
Received:4September2009/Accepted:26July2010
©Springer-Verlag2010
Abstract In this paper we present a unified discussion of different approaches to
theidentificationofsmoothingsplineanalysisofvariance(ANOVA)models:(i)the
“classical”approach(inthelineofWahbainSplineModelsforObservationalData,
1990; Gu in Smoothing Spline ANOVA Models, 2002; Storlie et al. in Stat. Sin.,
2011)and(ii)theState-DependentRegression(SDR)approachofYounginNonlin-
earDynamicsandStatistics(2001).Thelatterisanonparametricapproachwhichis
very similar to smoothing splines and kernel regression methods, but based on re-
cursive filtering and smoothing estimation (the Kalman filter combined with fixed
interval smoothing). We will show that SDR can be effectively combined with the
“classical”approachtoobtainamoreaccurateandefficientestimationofsmoothing
splineANOVAmodelstobeappliedforemulationpurposes.Wewillalsoshowthat
suchanapproachcancomparefavorablywithkriging.
Keywords SmoothingsplineANOVAmodels·Recursivealgorithms·Backfitting·
Sensitivityanalysis
1 Introduction
In the analysis of data from computer experiments, approximation models are built
toemulatethebehavioroflargecomputationalmodels.Smoothingsplinemodelsare
(cid:2)
M.Ratto( )·A.Pagano
JointResearchCentre,EuropeanCommission,TP361IPSC,ViaFermi2749,21027Ispra(VA),Italy
e-mail:[email protected]
A.Pagano
e-mail:[email protected]
Presentaddress:
A.Pagano
Euro-MediterraneanCentreforClimateChange,Milano,Italy
368 M.Ratto,A.Pagano
a useful tool for this kind of analysis (Levy and Steinberg 2010). In this paper we
presentaunifieddiscussionofdifferentapproachestotheidentificationofsmoothing
splineanalysisofvariance(ANOVA)models.The“classical”approachtosmoothing
splineANOVAmodelsisalongthesamelinesastheworkofWahba(1990)andGu
(2002).Recently,Storlieetal.(2011)presentedtheAdaptiveCOmponentSelection
and Selection Operator (ACOSSO), “a new regularization method for simultaneous
modelfittingandvariableselectioninnonparametricregressionmodelsintheframe-
workofsmoothingsplineANOVA.”ThismethodisanimprovementtoCOSSO(Lin
and Zhang 2006), penalizing the sum of component norms, instead of the squared
normemployedinthetraditionalsmoothingsplinemethod.InACOSSO,anadaptive
weight is used in the COSSO penalty which allows for more flexibility in estimat-
ingimportantfunctionalcomponentswhilegivingaheavierpenaltyto unimportant
functionalcomponents.
Ina“parallel”streamofresearch,usingtheState-Dependent(Parameter)Regres-
sion(SDR)approachofYoung(2001),Rattoetal.(2007)havedevelopedanonpara-
metric approach, very similar to smoothing splines and kernel regression methods,
based on recursive filtering and smoothing estimation (the Kalman filter combined
with fixed interval smoothing). Such a recursive least-squares implementation has
somekeycharacteristics:(a)itispluggedwithoptimalmaximumlikelihoodestima-
tion, thus allowing for an objective estimation of the smoothing hyper-parameters,
and(b)itprovidesgreaterflexibilityinadaptingtolocaldiscontinuities,heavynon-
linearity,andheteroscedasticerrorterms.Theuseofrecursivealgorithmsinsmooth-
ing splines is not new in statistical literature: the works of Weinert et al. (1980)
andWeckerandAnsley(1983)demonstratedtheapplicabilityofastochasticframe-
work for recursive computation of smoothing splines. However, such works were
limitedtotheunivariatecase,whilethesubsequenthistoryoftensorproductsmooth-
ing splines developed in the “standard” nonrecursive form. The recursive approach
of Young (2001) can be seen as an extension and generalization of such earlier pa-
pers,whichisapplicabletothemultivariatecase,andRattoetal.(2007)discussedits
possibleextensionfortheestimationofinteractiontermsoftheANOVA.
Thepurposesofthispaperare:
1. todevelopaformalcomparisonanddemonstrateequivalencesbetweenthe“clas-
sical” tensor product smoothing spline with reproducing kernel Hilbert space
(RKHS) algebra and the SDR, extending the results derived in Weinert et al.
(1980)andWeckerandAnsley(1983)tothemultivariatecase;
2. todiscusstheadvantagesanddisadvantagesoftheseapproaches;
3. toproposeaunifiedapproachtosmoothingsplineANOVAmodelsthatcombines
thebestofthediscussedmethods:theuseoftherecursivealgorithmsinparticular
can be very effective in detecting the important functional components, adding
valuableinformationintheACOSSOframework;atthesametime,suchacom-
bined approach improves the first extension of the SDR approach to interaction
termsproposedbyRattoetal.(2007);
4. tocompareresultsofanalysisofcomputerexperimentscarriedoutusingkriging-
based techniques, e.g., Design and Analysis of Computer Experiments (DACE)
andtheproposedunifiedmethod.
RecursivealgorithmsforsmoothingsplineANOVAmodels 369
Concerning item 1 above, for the sake of parsimony we concentrate here on the
specialcaseofcubicsmoothingsplines.Thiscaseis,infact,themostwidelyapplied
splinemethodaswellasthebasiccoreintheACOSSOframework.However,simi-
larlytothediscussioninWeckerandAnsley(1983),afullclassofequivalencescan
bederivedbetweentheSDRrecursiveapproachandpolynomialsplinesofanyorder.
In this context, we also note that Young and Pedregal (1999) already discussed the
equivalence between the recursive and en bloc formulation of the smoothing prob-
lem,whentheenbloccaseiscastinthe“discrete”form,asintheHodrick–Prescott
filter(HodrickandPrescott1981;Leser1961).Hereweextendsuchequivalenceto
the“continuous”integralformtypicaloftheRKHSalgebra.
2 State-dependentregressionsandsmoothingsplines
2.1 Additivemodels
We denote the generic mapping as z(X) and assume without loss of generality that
X ∈ [0,1]p, where p is the number of input variables. The simplest example of
smoothingsplinemappingestimationofzistheadditivemodel:
(cid:2)p
f(X)=f + f (X ). (1)
0 j j
j=1
Toestimatef wecanuseamultivariate(cubic)smoothingsplineminimizationprob-
lem,thatis,givenλ=(λ ,...,λ ),findtheminimizerf(X )of:
1 p k
(cid:5)
1 (cid:2)N (cid:3) (cid:4) (cid:2)p 1(cid:6) (cid:7)
z −f(X ) 2+ λ f(cid:3)(cid:3)(X ) 2dX , (2)
N k k j j j j
k=1 j=1 0
whereaMonteCarlo(MC)sampleofdimensionN isassumed.
This statistical estimation problem requires the estimation of the p hyper-
parameters λ (also denoted as smoothing parameters). Various ways of doing that
j
areavailableintheliterature,byapplyinggeneralizedcross-validation(GCV),gen-
eralized maximumlikelihood(GML) procedures and so on (see, e.g., Wahba1990;
Gu 2002). Here we discuss the recursive estimation approach, where the additive
model is put into the State-DependentParameter Regression(SDR) form of Young
(2001)asappliedtotheestimationofANOVAmodelsbyRattoetal.(2007).
We highlight here the key features of Young’s recursive algorithms of SDR by
considering the case of p =1 and z(X)=g(X)+e, with e∼N(0,σ2); i.e., we
rewrite the smoothing problem as z =s +e , where k =1,...,N and s is the
k k k k
estimateofg(X ).Tomaketherecursiveapproachmeaningful,theMCsampleneeds
k
to be sorted in ascending order with respect to X: i.e., k and k−1 subscripts are
adjacent elements under such ordering (see Fig. 1), implying 0≤X <X <···<
1 2
X <···<X ≤1.
k N
To recursively estimate the s in SDR, it is necessary to characterize it in some
k
stochastic manner, borrowing from nonstationary time series processes. In general
370 M.Ratto,A.Pagano
Fig.1 Exampleofthesorted
k-ordering:therecursive
algorithmspanstheoutputdata
zkfromthelefttotherightof
theplot
thisisaccomplishedbyassumingthattheevolutionofs followsonememberofthe
k
generalizedrandomwalk(GRW)classonnonstationaryrandomsequences(see,e.g.,
YoungandNg1989andNgandYoung1990).
In the present context, the integrated random walk (IRW) process provides the
samesmoothingpropertiesofacubicspline,intheoverallstate-spaceformulation:
ObservationEquation: z =s +e
k k k
StateEquations: sk =sk−1+dk−1 (3)
dk =dk−1+ηk
whered isthe“slope”ofs ,η ∼N(0,σ2)andη (“systemdisturbance”insystems
k k k η k
terminology)isassumedtobeindependentofthe“observationnoise”e ∼N(0,σ2).
k
Within the framework of the analysis of computer experiments, it seems neces-
sarytodiscusstheterme in(3).Normalityandindependencearestrictlyappropri-
k
atewhenthereisobservationalerrorinthedatabutcanbereasonableforsmoothing
observeddataevenincomputerexperimentsbecausetherecanbeapplicationswhere
the“computed”valueisproducedwithsomeerrororvariability,dueto,e.g.,conver-
gence of numerical algorithms. Moreover, e descends naturally from the standard
k
smoothingsplineformulation(2):suchapenalizedleast-squaresregressionrulesout
a perfect fit for f(X ), thus implying an “observation noise” linked to the nonzero
k
residual term (z −f(X )). This residual reflects the fact that the present emula-
k k
tion approach is based on identifying a truncated ANOVA expansion that approxi-
matesz(X).Thisisdonebyincludinga“small”subsetofq ANOVAterms(mainef-
fectsandinteractions)thatarestatisticallyidentifiablefromtheavailableMCsample.
Thus,e canbeseenasthesumofallthetermsthatarenotincludedbyashrinkage
k
procedure.ThissetofdroppedANOVAtermsusuallyincludesaverylargenumber
of elements (namely 2p −q, where q (cid:6)2p), which are orthogonal (independent)
by definition. It does not seem out of place to model the sum of a large number of
independentvariablesinstatisticalterms(CentralLimitTheorem).Wewillseethat
theinclusionofthis“error”term,ratherthanbeingadrawbackofthismethod,turns
RecursivealgorithmsforsmoothingsplineANOVAmodels 371
outtobeanadvantage(seetheexamplesinSect.3),sinceitimpliesthatemulation
(andtherefore“prediction”atuntriedX values)isperformedonlyusingstatistically
significantANOVAterms,enhancingtherobustnessofout-of-sampleperformances.
GiventheascendingorderingoftheMCsample,s canbeestimatedbyusingthe
k
recursive Kalman filter (KF) and the associated recursive fixed interval smoothing
(FIS)algorithm(see,e.g.,Kalman1960;Young1999fordetails).
First, it is necessary to optimize the hyper-parameters associated with the state-
space model (3), namely the white noise variances σ2 and σ2. In fact, by a simple
η
reformulationoftheKFandFISalgorithms,theIRWmodelcanbeentirelycharac-
terized by one noise variance ratio (NVR) hyper-parameter, where NVR=σ2/σ2.
η
ThisNVRvalueis,ofcourse,unknownapriori andneedstobeoptimized:forex-
ample, in the above references, this is accomplished by maximum likelihood (ML)
optimizationusingpredictionerrordecomposition(Schweppe1965).TheNVRplays
theinverseroleofasmoothingparameter:thesmallertheNVR,thesmootherthees-
timate of s (and at the limit NVR=0, s will be a straight line). Given the NVR,
k k
theFISalgorithmthenyieldsanestimate sˆk|N of sk ateachdatasample,anditcan
be seen that the sˆk|N from the IRW process is the equivalent of f(Xk) in the cubic
smoothing spline model. At the same time, the recursive procedures provide, in a
natural way, standard errors of the estimated sˆk|N that allow for the testing of their
relativesignificance.
We need to clarify here the meaning of the ML optimization in this recursive
context. We first observe that, to avoid a perfect fit solution, a penalty term is used
inthe“classical”smoothingsplineestimates.The“penalty”appearsintheobjective
function (GCV, GML, etc.) which is used to optimize the λ’s. In this way one may
limitthe“degreesoffreedom”ofthesplinemodel.
Forexample,inGCV,wehavetofindaλthatminimizes
(cid:8)
(z −f (X ))2
GCV =1/N· k k λ k , (4)
λ (1−df(λ)/N)2
wheredf ∈[0,N]denotesthe“degreesoffreedom”asafunctionofλ.Similarly,in
theSDRnotation,welookforanNVRminimizing
(cid:8)
GCV =1/N· k(zk−sˆk|N)2 , (5)
NVR (1−df(NVR)/N)2
where the “degrees of freedom” df depend on the NVR (by the above-mentioned
equivalence).
Remark1 Withoutthepenaltyterm,theoptimumwouldalwaysbeattainedatλ=0
(orNVR→∞),i.e.,perfectfit.
A perfect fit “optimal” solution would never happen within the SDR recursive
context. In this case penalty is directly associated with the estimate procedure, by
thefactthatthepredictionerrordecomposition(ML)estimateisbasedonthefiltered
estimatesˆk|k−1=sk−1+dk−1andnotonthesmoothedestimatesˆk|N.Namely,inML
372 M.Ratto,A.Pagano
Fig.2 Pictureoftheonestep
aheadpredictionsoftheIRW
processinthecaseof
NVR→∞.Thestars/dotted
linedenotethezkdataseries,
whilethesquaresdenotetheone
stepaheadpredictionsˆk|k−1
estimation,wefindtheNVRthatmaximizesthelog-likelihoodfunctionL,where:
(cid:2)N (cid:3) (cid:4)
ˆ
−2·log(L)=const+ log(1+Pk|k−1)+(N−2)·log σ2 ,
k=3
σˆ2= 1 (cid:2)N (zk−sˆk|k−1)2, (6)
N−2k=3 (1+Pk|k−1)
ˆ
where σ2 is the “weighted average” of the squared innovations (i.e., the prediction
erroroftheIRWmodel),andPk|k−1 istheonestepaheadforecasterrorofthestate
sˆk|k−1 provided by the KF (both Pk|k−1 and sˆk|k−1 are functions of NVR). Since
sˆk|k−1 isbasedonlyontheinformationcontainedinthesamplevalues[1,...,k−1]
(whilesmoothedestimatesusetheentireinformationset[1,...,N]),itcanbeeasily
seenthatthelimitNVR→∞(λ→0)isnolongera“perfectfit”situation,sincea
zero variance for ek implies sˆk|k−1=sk−1+dk−1=zk−1+dk−1; i.e., the one step
aheadpredictionofzk isgivenbythelinearextrapolationoftheadjacentvaluezk−1,
thusimplyinganonzeropredictionerrorinthislimitcase.
ThisisfurtherexemplifiedinFig.2:thesquaresintheplotsdenotetheonestep
ahead prediction sˆk|k−1 and the arrows show the linear extrapolation mechanism of
theIRWprocesswhenNVR→∞.Suchapredictiondepartsconsiderablynotonly
froma“perfectfit”situationbutalsofroma“decentfit,”implyingthattheMLesti-
matewillautomaticallypenalizethiskindofsituationandprovidethe“right”value
fortheNVR(seealsoWeckerandAnsley1983foradiscussionofthisMLestimator
inthesmoothingsplinecontext).
TocompletetheequivalencebetweentheSDRandcubicsplineformulations,we
needtolinktheNVRestimatedbytheMLproceduretothesmoothingparametersλ.
It can be easily verified that by setting λ=1/(NVR·N4), and with evenly spaced
Xk values,the f(Xk) estimateinthecubicsmoothingsplinemodelequalsthe sˆk|N
estimate from the IRW process. The present results are also in line with the cited
RecursivealgorithmsforsmoothingsplineANOVAmodels 373
workofWeckerandAnsley(1983),whoassume,however,adifferentstochasticform
fors ,namelyanARIMA(0,2,2)process.
k
As mentioned in the Introduction, this is not the only possible equivalence be-
tween recursive algorithms and polynomial splines. For example, one can verify
that assuming a random walk stochastic process for s with some ML optimized
k
NVR(RW) is equivalent to estimating a linear spline with smoothing parameter
λ=1/(NVR(RW)·N2)(seealsoWeckerandAnsley1983).
The most interesting aspect of the SDR approach is that it is not limited to the
univariatecase,butcanbeeffectivelyextendedtothemostrelevantmultivariateone.
In the general additive case (1), for example, the recursive procedure needs to be
applied,inturn,foreachtermfj(Xj,k)=sˆj,k|N,requiringadifferentsortingstrategy
foreachsˆj,k|N.Hence,the“backfitting”procedure,asdescribedinYoung(2000)and
Young (2001), is exploited (see Appendix A.1). This procedure provides both ML
estimates of all NVRj’s and the smoothed estimates of the additive terms sˆj,k|N. It
canbeeasilyverifiedthattheequivalencebetweenthe λ’sandNVRs presentedfor
p=1alsoholdsfortheadditivemodelwithp>1.So,theestimatedNVR ’scanbe
j
convertedintoλ valuesusingλ =1/(NVR ·N4),allowingustoputtheadditive
j j j
modelintothestandardcubicsplineform.
2.2 ANOVAmodelswithinteractionfunctions
Theadditivemodelconcept(1)canbegeneralizedtoincludetwo-way(andhigher)
interaction functions via the functional ANOVA decomposition (Wahba 1990; Gu
2002).Forexample,wecanlet
(cid:2)p (cid:2)p
f(X)=f + f (X )+ f (X ,X ). (7)
0 j j j,i j i
j=1 j<i
In the ANOVA smoothing spline context, corresponding optimization problems
withinteractionfunctionsandtheirsolutionscanbeobtainedconvenientlywiththe
reproducing kernel Hilbert space (RKHS) approach (see Wahba 1990). In the SDR
context, we propose here to formalize an interaction function as the product of two
states s ·s , each of them characterized by an IRW stochastic process. Hence, the
1 2
estimationofasingleinteractiontermz∗(X )=f(X ,X )+e isexpressedas
k 1,k 2,k k
ObservationEquation: z∗=sI ·sI +e
k 1,k 2,k k
StateEquations:(j =1,2) sjI,k =sjI,k−1+djI,k−1 (8)
dI =dI +ηI
j,k j,k−1 j,k
wherez∗isthemodeloutputafterthemaineffectsaretakenout,I =1,2isthemulti-
indexdenotingtheinteractiontermunderestimation,andηI ∼N(0,σ2 ).Thetwo
j,k ηI
j
termssI areestimatediterativelybyrunningtherecursiveprocedureinturn;i.e.,
j,k
– takeaninitialestimateofsI andsI byregressingzwiththeproductofsimple
1,k 2,k
linearorquadraticpolynomialsp (X )·p (X )andsetsI,0=p (X );
1 1 2 2 j,k j j,k
374 M.Ratto,A.Pagano
– iteratei=1,2:
– fixsI,i−1 andestimateNVRI andsI,i usingtherecursiveprocedure;
2,k 1 1,k
– fixsI,i andestimateNVRI andsI,i usingtherecursiveprocedure;
1,k 2 2,k
– theproductsI,2·sI,2obtainedaftertheseconditerationprovidestherecursiveSDR
1,k 2,k
estimateoftheinteractionfunction.
Thelatterstoppingcriterionisaconvenientchoicetolimitthecomputationtime,
andisduetotheobservationthattheestimateoftheinteractiontermneverchanged
toomuchinanysubsequentiteration.
Unfortunately, in the case of interaction functions we cannot derive an explicit
andfullequivalencebetweenSDRandcubicsplinesofthetypementionedforfirst-
orderANOVAterms.Therefore,inordertobeabletoexploittheestimationresults
inthecontextofasmoothingsplineANOVAmodel,weproposetotakeadifferent
approach,similartotheACOSSOcase.
2.3 VeryshortsummaryofACOSSO
Wemaketheusualassumptionthatf ∈F,whereF is(cid:9)anRKHS.ThespaceF can
bewrittenasanorthogonaldecompositionF ={1}⊕{ q F },whereeachF is
j=1 j j
itselfanRKHSandj =1,...,q spansANOVAtermsofvariousorders.Typically,q
includesthemaineffectsplusrelevantinteractionterms.
Wereformulate(2)forthegeneralcasewithinteractionsusingthefunctionf that
minimizes
(cid:2)N (cid:3) (cid:4) (cid:2)q (cid:10) (cid:10)
1 z −f(X ) 2+λ 1 (cid:10)Pjf(cid:10)2 , (9)
k k 0 F
N θ
k=1 j=1 j
where Pjf is the orthogonal projection of f onto F and the q-dimensional vec-
j
tor θ of smoothing parameters needs to be optimized somehow. This is typically
j
a formidable problem, and in the simplest case θ is set to one, with the single λ
j 0
estimatedbyGCVorGML.
Problem (9) also poses the issue of selection of F terms: this is tackled rather
j
effectivelywithintheCOSSO/ACOSSOframework.
The COSSO (Lin and Zhang 2006) penalizes the sum of norms, using a Least
AbsoluteShrinkageandSelectionOperator(LASSO)-typepenalty(Tibshirani1996)
for the ANOVA model, which allows us to identify the informative predictor terms
F withanestimateoff thatminimizes
j
(cid:2)N (cid:3) (cid:4) (cid:2)Q (cid:10) (cid:10)
1 z −f(X ) 2+λ (cid:10)Pjf(cid:10) (10)
k k F
N
k=1 j=1
usingasinglesmoothingparameterλ,andwhereQincludesallANOVAtermstobe
potentiallyincludedinf,e.g.,withatruncationuptosecond-orthird-orderinterac-
tions.
RecursivealgorithmsforsmoothingsplineANOVAmodels 375
ItcanbeshownthattheCOSSOestimateisalsotheminimizerof
(cid:2)N (cid:3) (cid:4) (cid:2)Q (cid:10) (cid:10)
1 z −f(X ) 2+ 1 (cid:10)Pjf(cid:10)2 (11)
k k F
N θ
k=1 j=1 j
(cid:8)
Q
subjectto 1/θ <M (wherethereisa1-1mappingbetweenM andλ).Sowe
j=1 j
can think of the COSSO penalty as the traditional smoothing spline penalty plus a
penaltyontheQsmoothingparametersusedforeachcomponent.TheLASSO-type
penalty has the effect of setting some of the functional components (F ’s) equal to
j
zero (e.g., the variable X or the interaction (X ,X ) is not in the model); thus, it
j j i
“automatically”selectstheappropriatesubsetq oftermsoutoftheQ“candidates.”
The key property of COSSO is that with one single smoothing parameter (λ or M)
itprovidesproperestimatesofallθ parameters:therefore,itconsiderablyimproves
j
problem(9)withθ =1(stillwithonesinglesmoothingparameterλ )andismuch
j 0
morecomputationallyefficientthanthefullproblem(9)withoptimizedθ ’s.
j
IntheadaptiveCOSSO(ACOSSO)ofStorlieetal.(2011),f ∈F minimizes
(cid:2)N (cid:3) (cid:4) (cid:2)q (cid:10) (cid:10)
1 z −f(X ) 2+λ w (cid:10)Pjf(cid:10) , (12)
k k j F
N
k=1 j=1
where 0<w ≤∞ are weights that depend on an initial estimate of f˜, either us-
j
ing (9) with θ = 1 or the COSSO estimate (10). The adaptive weights are ob-
j
tained as w = (cid:11)Pjf˜(cid:11)−γ, typically with γ = 2 and the L norm (cid:11)Pjf˜(cid:11) =
((cid:11) (Pjf˜(X))j2dX)1/2.ThLe2useofadaptiveweightsimprovesthe2predictivecapabLi2lity
ofANOVAmodelswithrespecttotheCOSSOcase.
2.4 CombiningSDRandACOSSOforinteractionfunctions
ThereisanobviouswayofexploitingtheSDRidentificationandestimationstepsin
theACOSSOframework:namely,theSDRestimatesofadditiveandinteractionfunc-
˜
tiontermscanbetakenastheinitialf usedtocomputetheweightsintheACOSSO.
However, this would be a minimal approach, whereas the SDR identification and
estimationprovidesmoredetailedinformationaboutf termsthatisworthexploit-
j
ing.WedefineK(cid:12)j(cid:13) tobethereproducingkernel(r.k.)ofanadditivetermFj ofthe
ANOVAdecompositionofthespaceF.Inthecubicsplinecase,thisisconstructedas
thesumoftwotermsK(cid:12)j(cid:13)=K01(cid:12)j(cid:13)⊕K1(cid:12)j(cid:13),whereK01(cid:12)j(cid:13)isther.k.oftheparametric
(linear)partandK1(cid:12)j(cid:13) isther.k.ofthepurelynonparametricpart.Thesecond-order
interactiontermsareconstructedasthetensorproductofthefirst-orderterms,fora
totaloffourelements,i.e.,
K(cid:12)i,j(cid:13)=(K01(cid:12)i(cid:13)⊕K1(cid:12)i(cid:13))⊗(K01(cid:12)j(cid:13)⊕K1(cid:12)j(cid:13))
=(K01(cid:12)i(cid:13)⊗K01(cid:12)j(cid:13))⊕(K01(cid:12)i(cid:13)⊗K1(cid:12)j(cid:13))⊕(K1(cid:12)i(cid:13)⊗K01(cid:12)j(cid:13))⊕(K1(cid:12)i(cid:13)⊗K1(cid:12)j(cid:13)).
(13)
376 M.Ratto,A.Pagano
Ingeneral,consideringproblem(9),oneshouldattributeaspecificcoefficientθ(cid:12)·(cid:13)
to each single element of the r.k. of F (see, e.g., Gu 2002, Chap. 3), i.e. two θ’s
j
foreachmaineffect,four θ’sforeachtwo-wayinteraction,andsoon.Infact,each
Fj would be optimally fitted by opportunely choosing weights in the sum of K(cid:12)·,·(cid:13)
elements.This,however,makestheestimationproblemrathercomplex,so,usually,
thetensorproduct(13)isdirectlyused,withouttuningtheweightsofeachelement
ofthesum.ThisstrategyisalsoappliedinACOSSO.
Instead,weproposetouseSDRestimatesofinteractiontosettheweights.
Inparticular,wecanseethattheSDRestimateoftheinteraction(8)isgivenbythe
productoftwounivariatecubicsplines.So,onecaneasilydecomposeeachestimated
sˆI into the sum of a linear (sˆI ) and nonparametric term (sˆI ). This provides a
j 01(cid:12)j(cid:13) 1(cid:12)j(cid:13)
decompositionoftheSDRinteractionoftheform
sˆI ·sˆI =sˆI sˆI +sˆI sˆI +sˆI sˆI +sˆI sˆI , (14)
i j 01(cid:12)i(cid:13) 01(cid:12)j(cid:13) 01(cid:12)i(cid:13) 1(cid:12)j(cid:13) 1(cid:12)i(cid:13) 01(cid:12)j(cid:13) 1(cid:12)i(cid:13) 1(cid:12)j(cid:13)
whichcanbethoughtofasaproxyofthefourelementsofther.k.ofthesecond-order
tensorproductcubicspline.
This suggests that a natural use of the SDR identification and estimation in the
ACOSSOframeworkistoapplyspecificweightstoeachelementofther.k.K(cid:12)·,·(cid:13) in
(13). In particular, the weights are the L norms of each of the four elements esti-
2
matedin(14).Wewillshowintheexamplesthatthischoiceleadstosignificantim-
provementintheaccuracyofANOVAmodelswithrespecttotheoriginalACOSSO
approach.
2.5 Krigingmethod:theDACEMatlabtoolbox
DACE(Lophavenetal.2002)isaMatlabtoolboxusedtoconstructkrigingapprox-
imationmodelsonthebasisofdatafromcomputerexperiments.Oncewehavethis
approximatemodel,wecanuseitasasurrogatemodel(emulator,meta-model).
WebrieflyhighlightthemainfeaturesofDACE.
Keepingwhereverpossiblethesamenotationthatweusedpreviously,thekriging
modelcanbeexpressedasaregression
zˆ(X)=β f (X)+···+β f (X)+ζ(X), (15)
1 1 q q
wheref ,i=1,...,qaredeterministicregressionterms,β aretherelatedregression
i i
coefficients, and ζ is a zero-mean random process whose variance depends on the
process variance ω2 and on the correlation R(v,w) between ζ(v) and ζ(w). The
toolboxprovidesasetofcorrelationfunctionsdefinedas
(cid:12)
R(θ,v,w)= R (θ ,w −v ).
j j j j
j=1:p
In particular, for the generalized exponential correlation function, used in the next
section,onehas
(cid:3) (cid:4)
Rj(θj,wj −vj)=exp −θj|wj −vj|θp+1 .