Factor Selection and Structural Identification in the Interaction ANOVA

28 Pages·2011·0.29 MB·English

Factor Selection and Structural Identification
in the Interaction ANOVA Model
Justin Post and Howard D. Bondell
Abstract
When faced with categorical predictors and a continuous response, the objective of
analysisoftenconsistsoftwotasks: Findingwhichfactorsareimportantanddetermining
whichlevelsofthefactorsdiffersignificantlyfromoneanother. Oftentimesthesetasksare
done separately using Analysis of Variance (ANOVA) followed by a post-hoc hypothesis
testingproceduresuchasTukey’sHonestlySignificantDifferencetest. Wheninteractions
between factors are included in the model the collapsing of levels of a factor becomes a
moredifficultproblem. Forinterpretability, whencollapsingtwolevelsofafactoritmay
notmakesensetocollapsethemaineffectdifferencetozerowhilestillhavinginteraction
differences between those levels nonzero. This structure between the main effects and
interactions in a model is similar to the idea of heredity used in regression models. This
paper introduces a new method for accomplishing both of the common analysis tasks si-
multaneouslyinaninteractionmodelwhilealsoadheringtotheheredity-typeconstrainton
themodel. Agroupnormpenaltyisplacedonparametersthatencourageslevelsoffactors
tocollapseandentirefactorstobesettozero. TheprocedureiscalledGASH-ANOVAfor
grouping and selection using heredity in ANOVA. It is shown that the procedure has the
oracle property implying that asymptotically it performs as well as if the exact structure
were known beforehand. We also discuss the application of a modified GASH-ANOVA
estimatorforestimatinginteractionsintheunreplicatedcase. Simulationstudiesshowthat
theGASH-ANOVAprocedureoutperformsposthochypothesistestingproceduresaswell
assimilarmethodsthatdonotincludeastructuralconstraint. Themethodisalsoillustrated
usingarealdataexample.
1
1 Introduction
Considerthecommoncaseofacontinuousresponsevariableandcategoricalpredictors
(factors). AconventionalwaytojudgetheimportanceofthesefactorsistouseAnalysisof
Variance(ANOVA).Oncefactorsaredeemedimportant,tocheckwhichlevelsofthefac-
torsdifferfromoneanotherthenextstepisoftentodoaposthocanalysissuchasTukey’s
honestly significantly different test, Fisher’s least significant difference test, or pairwise
comparisonsoflevelsusingaBonferroniadjustmentoraBenjamini-Hochberg(Benjamini
and Hochberg, 1995) type adjustment. Rather than carry out these two tasks separately,
a technique called CAS-ANOVA, for collapsing and shrinkage in ANOVA (Bondell and
Reich,2009),hasbeendevelopedtoperformtheseactionssimultaneously. Thisprocedure
is a constrained or penalized regression technique. Much of the recent variable selection
literature is of this form and examples include the least absolute shrinkage and selection
operator (LASSO) (Tibshirani, 1996), the Elastic Net (EN) (Zou and Hastie, 2005), the
SmoothlyClippedAbsoluteDeviationpenalty(SCAD)(FanandLi,2001),andtheOctag-
onal Shrinkage and Clustering Algorithm for Regression (OSCAR) (Bondell and Reich,
2008). The CAS-ANOVA procedure places an L constraint directly on the pairwise dif-
1
ferencesineachfactorallowinganentirefactortobezeroedoutwhilealsoallowinglevels
within a factor to collapse (be set equal) into groups. Not only does the CAS-ANOVA
procedure accomplish both tasks at once, the nature of the penalty requires the levels of
eachfactortobecollapsedintonon-overlappinggroups. Thismethodwasshowntohave
the oracle property, implying that its performance is asymptotically equivalent to having
the true grouping structure known beforehand and performing the standard least squares
ANOVAfittothiscollapseddesign.
AlimitationoftheCAS-ANOVAmethodisthatitwasdevelopedassumingamainef-
fectonlymodelandoftenaninteractionmodelismorerealistic. IntheinteractionANOVA
model, when determining whether two levels of a factor should collapse we must look at
more than just their respective main effects. It is not appropriate to claim the two levels
should collapse if there is an interaction term that differs between the two levels for any
given level of another factor. This leads to the idea that only levels whose interaction ef-
fects are all identical, should have their main effects also collapsed together. Thus, if we
canenforcethistypeofstructureonourmodelwewillbeabletoaccomplishbothtasksof
2
ouranalysis. Thenotionofthismodelstructureissimilartotheideaofheredityoftenused
in regression models when higher order terms are involved. For a regression model such
asE(y ) = β +β x +β x +β x x +β x2 ,wherey isacontinuousresponsewith
i 0 1 i1 2 i2 3 i1 i2 4 i1 i
x and x continuous predictors, strong heredity implies x2 should only appear in the
i1 i2 i1
modelifx appearsinthemodelandtheinteractiontermx x shouldonlyappearinthe
i1 i1 i2
modelifbothx andx appearinthemodel. Theheredityprincipleisimportantbecause
i1 i2
it aids in interpretation of the model and ensures that the model makes sense structurally.
This type of constraint on the structure of the predictors in a regression model has seen
much use. For example, Choi, N., Li,W. and Zhu, J. (2010) used heredity to extend the
LASSOvariableselectiontechniquetoincludeinteractionterms,Yuan,M.,Joseph,V.and
Zou,H.(2009)usedheredityinvariableselectionwiththenon-negativegarrote,Yuan,M.,
Joseph,V.andLin,Y.(2007)usedtheconstraintwiththeLARSalgorithm,andChipman
(1996)usedtheconstraintintheBayesianvariableselectioncontext. However,allofthese
approachesdealwithcontinuouspredictors, inthataninteractionisasingletermderived
from the product of two other predictors. In the interaction ANOVA model, interaction
terms are not single terms, they arise as products of groups of variables, and thus need to
betreateddifferently.
This paper develops a method for use in the interaction ANOVA model to simultane-
ously perform the two main goals of analysis. The method utilizes a weighted penalty
thatenforcestheheredity-typestructureonthemodel. ThenewmethodiscalledGASH-
ANOVAforGroupingAndSelectionusingHeredityinANOVA.Weshowthattheoracle
property holds for the GASH-ANOVA procedure, in that asymptotically it performs as
wellasiftheexactstructurewereknownbeforehand. Thispropertyalsoimpliesthat,even
after the possible dimension reduction, asymptotic inference can be based on standard
ANOVAtheory. WealsodiscusstheuseofanunweightedversionoftheGASH-ANOVA
estimatortoestimateinteractiontermsintheunreplicatedANOVAmodel.
The rest of the article is organized as follows: In section 2, notation is introduced
and the related CAS-ANOVA procedure is reviewed. In section 3.1, the GASH-ANOVA
procedure is introduced. The extension to the unreplicated case is discussed in section
3.2. Insection4,theasymptoticpropertiesoftheGASH-ANOVAsolutionarepresented.
Section5discussescomputationandtuningfortheGASH-ANOVAmethod. Toillustrate
the method, section 6 gives simulation studies and their results and section 7 shows the
3
method’susefulnessonarealdataexample. Finally,section8givesadiscussion.
2 Notation and CAS-ANOVA
2.1 Notation
Tosimplifynotation, considerthetwofactorANOVAmodelwithfactorAandfactor
B having a and b levels respectively. The extension to more than two factors is straight-
forward and will be discussed shortly. We denote the number of observations at each
(cid:80)
level combination by n and denote the total sample size by n = n . Note that
ij i,j ij
for a balanced design we have n = abn . We use the matrix representation of the
ij
linear model, y = Xθ + (cid:15), where y is the n × 1 vector of responses, X is the typical
n×poverparameterizedANOVAdesignmatrixconsistingofzerosandonescorrespond-
ing to the combination of levels for each observation (where p = 1 + a + b + ab), θ
is a p × 1 vector of parameters, and (cid:15) is a n × 1 vector of error terms. The parame-
ter vector consists of the mean and three stacked vectors, θ = (µ,αT,βT,(αβ)T)T,
where µ is the intercept, αT = (α ,α ,...α ), βT = (β ,β ,...β ), and (αβ)T =
1 2 a 1 2 b
((αβ) ,(αβ) ,...(αβ) ,(αβ) ,...,(αβ) ,...,(αβ) ). Here, α corresponds to the
11 12 1b 21 2b ab i
maineffectoflevelioffactorA,β correspondstothemaineffectoflevelj offactorB,
j
and(αβ) correspondstotheinteractioneffectoflevelioffactorAandlevelj offactor
ij
B.Theordinaryleastsquares(OLS)solutioncanbewrittenasfollows:
θ(cid:98)OLS = argminθ||y−Xθ||2 (1)
a a
(cid:88) (cid:88)
subjectto α = 0, (αβ) = 0forallj,
i ij
i=1 i=1
b b
(cid:88) (cid:88)
β = 0, (αβ) = 0foralli,
j ij
j=1 j=1
wheretheconstraintsontheparametersareaptlynamedthesumtozeroconstraints. Note
thatn mustbeatleasttwoinordertoestimatetheinteractiontermswhenusingOLS.
ij
2.2 CAS-ANOVA
In order to shrink the coefficients and to perform variable selection in the additive
model,the(adaptive)CAS-ANOVA(BondellandReich,2009)procedureplacesaweighted
constraintdirectlyonthepairwisedifferencesofthelevelsofeachfactor. Inthetwofactor
4
maineffectonlymodel,whichhasthesamesetupastheabovemodelignoringtheinter-
actiontermsandconstraintsontheinteractionterms, theCAS-ANOVAprocedureaddsa
constraintthatisoftheform
(cid:88) (km) (cid:88) (km)
w |α −α |+ w |β −β | ≤ t,
α,CAS k m β,CAS k m
1≤k<m≤a 1≤k<m≤b
(km) (km)
where t > 0 is a tuning constant and w and w are scaled adaptive LASSO
α,CAS β,CAS
(Zou,2006)typeweightsforthepairoflevelsk andmofeachfactor.
2.3 NoHeredityMethod
We denote the extension of the CAS-ANOVA procedure to the interaction ANOVA
(cid:91)
modelasthe‘NoHeredity’(NH)method. Letα(cid:98)i,OLS,β(cid:98)j,OLS,and(αβ)ij,OLS betheOLS
estimates of the corresponding parameters found using equation (1). The NH method’s
penaltyisgivenby
(cid:88) (km) (cid:88) (km)
w |α −α |+ w |β −β |
α,NH k m β,NH k m
1≤k<m≤a 1≤k<m≤b
(cid:88) (cid:88) (km,lm)
+ w |(αβ) −(αβ) |
αβ,NH km lm
1≤m≤b1≤k<l≤a
(cid:88) (cid:88) (km,kl)
+ w |(αβ) −(αβ) | ≤ t,
αβ,NH km kl
1≤k≤a1≤m≤l<b
wheretheweightsonthemaineffectdifferencesaregivenby
(cid:12) (cid:12)−1
wα(k,NmH) = |α(cid:98)k,OLS −α(cid:98)m,OLS|−1 andwβ(k,NmH) = (cid:12)(cid:12)β(cid:98)k,OLS −β(cid:98)m,OLS(cid:12)(cid:12)
andtheweightsplacedontheinteractiondifferencesaregivenby
w(km,kl) = (cid:12)(cid:12)((cid:91)αβ) −((cid:91)αβ) (cid:12)(cid:12)−1
αβ,NH (cid:12) km,OLS kl,OLS(cid:12)
and
w(km,lm) = (cid:12)(cid:12)((cid:91)αβ) −((cid:91)αβ) (cid:12)(cid:12)−1.
αβ,NH (cid:12) km,OLS lm,OLS(cid:12)
Asymptotically the NH method may perform well, although it is unlikely to perform
well at the two main tasks of choosing significant factors and collapsing levels in small
samples because of the lack of important heredity-type structure discussed in section
1. In problems where the factors have a large number of levels, the NH method will
have difficulty collapsing levels due to the sizable number of interaction differences that
5
need to be set to zero. Using our notation we can now describe the procedure for col-
lapsing two levels of a factor when interactions are present in more detail. In the two
factor situation considered here, our procedure implies that level k of factor A should
only be collapsed to level m of factor A if the main effect difference, α −α , and all
m k
the interaction differences between factor A and B involving level k and m of factor A,
(αβ) −(αβ) ,(αβ) −(αβ) ,...,(αβ) −(αβ) ,havebeensettozero. Thereisa
m1 k1 m2 k2 mb kb
chancefortheNHmethodtoselectfactorsandcollapselevelsbutnothingforcesthistobe
thecase,astraynonzerointeractiondifferencecanpreventthecollapsingfromoccurring.
Thus, for interpretability of the model and to accomplish our two main tasks of analysis
simultaneously,thisextensionoftheCAS-ANOVAprocedureisnotidealwithinteractions
present. Insection6and7,theNHmethodiscomparedtotheGASH-ANOVAprocedure.
3 GASH-ANOVA
3.1 Method
For computation and the statement of the theoretical results of the GASH-ANOVA
estimator it is more convenient to reparametrize to a full rank design matrix using a ref-
erencelevelasabaseline. Thislessensthenumberofparametersandconstraintsneeded.
Thus, from this point forward we will use the full rank design. We choose the first level
of each factor as the reference level, although this choice is arbitrary as the levels can
be relabeled. Define the new design matrix by X∗ and the new parameter vector by
θ∗ = (µ∗,α∗T,β∗T,(αβ)∗T)T,whereµ∗ = µ+α +β +(αβ) ,
1 1 11
α∗T = (α∗,...,α∗) = (α −α ,...,α −α ),
2 a 2 1 a 1
β∗ isa(b−1)×1vectordefinedsimilarlyforfactorB,and
(αβ)∗T = ((αβ)∗ ,(αβ)∗ ,...,(αβ)∗ ,(αβ)∗ ,...,(αβ)∗ ,...,(αβ)∗ ),
22 23 2b 32 3b ab
where(αβ∗) = (αβ) −(αβ) .
ij ij 11
To achieve the automatic factor selection and collapsing of levels in the interaction
model the GASH-ANOVA approach uses a weighted heredity-type constraint. To en-
courage the collapsing of levels, an infinity norm constraint is placed on (overlapping)
groups of pairwise differences belonging to different levels of each factor. In detail, we
form G = (cid:0)a(cid:1) + (cid:0)b(cid:1) groups where each group contains a main effect difference be-
2 2
tweentwolevelsofafactoralongwithallinteractiondifferencesthatinvolvethosesame
6
two levels. We denote each group of parameters by either φ ,1 ≤ i < j ≤ a or
α,ij
φ ,1 ≤ i < j ≤ b, where φ = (α∗,(αβ)∗ ,(αβ)∗ ,...,(αβ)∗ ) for 2 ≤ j ≤ a and
β,ij α,1j j j2 j3 jb
φ = (α∗ −α∗,(αβ)∗ −(αβ)∗ ,(αβ)∗ −(αβ)∗ ,...,(αβ)∗ −(αβ)∗) for 1 < i <
α,ij j i j2 i2 j3 i3 jb ib
j ≤ aandφ isdefinedsimilarly. Notethatthesegroupssharesomeinteractionterms.
β,ij
By judicious choice of overlapping groups, two main effects of a factor can be set equal
to one another only if all of the interactions for those two levels are also set equal and,
with probability one, an interaction difference is only present if the corresponding main
effect differences are also present. Thus, the GASH-ANOVA procedure adheres to our
heredity-type structure which encourags levels of each factor to be estimated with exact
equalityandentirefactorstobesettozero. Thisoverlappinggrouppenaltyonthediffer-
ences is related to the family of Composite Absolute Penalties (CAP) (Zhang, P., Rocha,
G.andYu, B.,2009). However, theCAPtreatscoefficientsthemselvesasgroups, notthe
differencesofcoefficientsasgroups.
TheGASH-ANOVAsolutioncanbewrittenindetailasfollows:
θ(cid:98)∗ = argminθ∗||y−X∗θ∗||2 (2)
(cid:80) (ij) (cid:80) (ij)
subjectto w max{|φ |}+ w max{|φ |} ≤ t,
1≤i<j≤a α α,ij 1≤i<j≤b β β,ij
(ij) (ij)
wherew andw areadaptiveweights,t > 0isatuningconstant,
α β
|φ | = (|α∗|,|(αβ)∗ |,|(αβ)∗ |,...,|(αβ)∗ |)T
α,ij j j2 j3 jb
for2 ≤ j ≤ aand
|φ∗ | = (|α∗−α∗|,|(αβ)∗ −(αβ)∗ |,|(αβ)∗ −(αβ)∗ |,...,|(αβ)∗ −(αβ)∗|)T
α,ij j i j2 i2 j3 i3 jb ib
for1 < i < j ≤ a,and|φ |issimilarlydefined. Usingequation(1)toobtaintheOLS
β,ij
(cid:110)(cid:12) (cid:12)(cid:111)
solution, the weight wα(ij) is given by (max (cid:12)(cid:12)φ(cid:98)α,ij,OLS(cid:12)(cid:12) )−1, where φ(cid:98)α,ij,OLS denotes
(ij)
the use of the OLS estimate for the differences and the form of the weight w is given
β
similarly. TheadaptiveweightsallowfortheasymptoticpropertiesoftheGASH-ANOVA
proceduregiveninsection4.
It may be of interest to not only be able to collapse entire levels of a factor, but to
also be able to collapse individual interaction differences. If all interaction differences
for all factors are set to zero then we are left with the additive model. To accomplish the
7
collapsing of individual interaction differences we can explicitly add these terms to the
penalty. Thenewpenaltywouldbegivenby
(cid:88) w(ij)max{|φ |}+ (cid:88) w(ij)max{|φ |}
α α,ij β β,ij
1≤i<j≤a 1≤i<j≤b
(cid:88) (cid:88) (km,lm)
+ w |(αβ) −(αβ) |
αβ,NH km lm
1≤m≤b1≤k<l≤a
(cid:88) (cid:88) (km,kl)
+ w |(αβ) −(αβ) | ≤ t,
αβ,NH km kl
1≤k≤a1≤m≤l<b
wheretheweightsareasdefinedpreviously. Theasymptotictheoryforthispenaltygiven
insection4shouldstillhold. Ifonedesiredevenmorecontroloftheinteractiondifferences
wecouldleavetheoriginalGASH-ANOVApenaltyaloneandcreateasecondpenaltywith
its own tuning parameter, say t , that involved only the interaction differences (explicitly
2
givenbythelasttwosumsofthepreviouspenalty). However,inthiscasewewouldneed
tofitalatticeofpointsoverourtuningparameterstofindasolution. Thisgreatlyincreases
thenumberofGASH-ANOVAsolutionstocompute.
When more than two factors are included in the model the method follows directly.
Withmorethantwofactorstheideaofcollapsingtwolevelsofafactorremainsthesame.
In order to collapse two levels, we need the main effect difference and any interaction
differences that involve those two levels to be set to zero. Thus, we need only augment
theφvectorswithallinteractiondifferencesnecessaryforcollapsingthetwolevelsofthe
givenfactor. Ifweassumeinteractionsoforderthreeorgreaterarenull,theφvectorsneed
only be augmented with all two-way interaction differences that involve the given levels
ofthefactor. Ifweallowforallhigherorderinteractions,theφvectorneedstoincludeall
higherorderinteractiondifferencesbetweenthelevels.
3.2 InvestigatingInteractionsintheUnreplicatedCase
Assume there is only one observation for each level combination, i.e. n = 1 for all
ij
i,j. This is common case when using a Randomized Complete Block Design (RCBD).
Using OLS for this case one usually assumes there is no interaction between the factors
(orthefactorandtheblock)astherearenotenoughdegreesoffreedomtoinvestigateand
testinteractioneffects. Anumberofsolutionshavebeenproposedtoinvestigatetheinter-
actionwhennoreplicationispresent(seeFranck,C.,Osborne,J.,andNielsen,D.(2011)
foradetailedreviewandcomparisonofmethods).
8
WecanuseamodifiedGASH-ANOVAproceduretoinvestigateinteractionsintheun-
replicatedcase. Thismodifiedversionusesanunweightedpenaltyenablingustoestimate
our model fits. The optional penalty that includes explicit penalization of the interaction
differencesdiscussedintheprevioussectioncouldalsobeappliedhere.
One issue of note for the unweighted GASH-ANOVA procedure is that of scaling. In
penalized regression procedures, having the effects on a comparable scale is important to
ensure the penalization is done equally. In most penalization methods, the design matrix
is scaled as to have unit L norm. However, since we have differences of parameters in
2
our penalty, the way to standardize the variables is not clear. The usual GASH-ANOVA
procedureremediesthisissuebyusingadaptiveweights,essentiallyplacingthepenalized
termsonthesamescale.
4 Asymptotic Properties
WheninvestigatingtheasymptoticpropertiesoftheGASH-ANOVAestimatorweas-
sume that each φ group is either truly all zero (collapsed) or all differences in that group
aretrulynonzero. Thisimpliesthatifamaineffectdifferenceinaφgroupistrulynonzero
then,providedtheotherfactorhasatleasttwodistinctlevels,thecorrespondinginteraction
differencesinthatφgroupmustallbenonzeroaswell.
Let A = {(i,j) : α (cid:54)= α } and A = {(i,j) : β (cid:54)= β } be defined as the set of in-
α i j β i j
dices for the main effect differences of each factor that are truly nonzero and let A =
α,n
(cid:110) (cid:111)
{(i,j) : αˆ (cid:54)= αˆ } and A = (i,j) : βˆ (cid:54)= βˆ be defined as the set of indices for
i j β,n i j
each factor whose main effect differences are estimated as nonzero. For the pairwise
differences indexed by A and A , let η be the vector of those pairwise differ-
α β Aα,Aβ
ences along with their corresponding interaction differences. Notice that the sets A
α
and A contain the indices for the truly significant level and factor structure. If this in-
β
formation were known a priori, the solution would be estimated by collapsing down to
this structure and then conducting the usual ANOVA analysis. Define η as this
(cid:101)Aα,Aβ
so called ‘oracle’ estimator of η . It is well known that under standard conditions
Aα,Aβ
(cid:16) (cid:17)
n−1/2 η −η → N(0,Σ). Let η denote the GASH-ANOVA estimator
(cid:101)Aα,Aβ Aα,Aβ (cid:98)Aα,Aβ
ofη . Theorem1givenbelowshowsthattheGASH-ANOVAobtainstheoracleprop-
Aα,Aβ
erty.
ThetheoremismosteasilystatedwhenwerewritetheGASH-ANOVAcriterioninits
9
correspondingLagrangianformulation:
 
∗ ||y−X∗θ∗||2+λn(cid:80)1≤i<j≤a w√α(inj)max{|φα,ij|}
θ(cid:98) = argminθ∗ +λn(cid:80)1≤i<j≤b w√β(inj)max{|φβ,ij|} .
Notethatthereisaone-to-onecorrespondencewiththetuningparametertandλ .
n
Theorem 1: Suppose that λn → ∞and √λnn → 0. The GASH-ANOVA estimator θ(cid:98)
anditscorrespondingestimatorofthepairwisedifferencesη hasthefollowingproperties:
(cid:98)
a)P(A = A ) → 1andP(A = A ) → 1
α,n α β,n β
b)n−1/2(η −η ) → N(0,Σ)
(cid:98)Aα,Aβ Aα,Aβ
TheproofofTheorem1isgivenintheappendix.
Theoraclepropertystatesthatthemethoddeterminesthecorrectstructureofthemodel
with probability tending to one. Additionally, it tells us that one can create a new design
matrixcorrespondingtothereducedmodelstructureselectedandconductinferenceusing
thestandardasymptoticvarianceobtainedfromOLSestimationonthatdesign. Notethat
thissecondlevelofinferencemaynotbenecessary,dependingonthegoalsofone’sstudy.
5 Computation and Tuning
TheGASH-ANOVAproblemcanbeexpressedasaquadraticprogrammingproblem.
Define
ζ = Mθ∗ = (µ∗,α∗T,ξT,β∗T,ξT,(αβ)∗T,ξT ,ξT )T,
α β αβ,A αβ,B
whereξ andξ arevectorscontainingthemaineffectpairwisedifferencesforeachfac-
α β
tor that do not involve the baseline level and ξ and ξ are vectors containing the
αβ,A αβ,B
interaction pairwise differences of interest for factor A and factor B, respectively, that do
notinvolvethebaselinelevels. Thematrix
 
1 0 0 0
 
 0 M 0 0 
 1 
M =  
 
0 0 M 0
 2 
 
0 0 0 M
3
needed to create this new parameter vector is block diagonal. The first block (a scalar)
corresponds to µ∗. The second block corresponds to factor A, M = (cid:2)I DT(cid:3)T , and
1 a−1 1
consists of an identity matrix of size a−1 and a matrix D of ±1 that creates ξ that is
1 α
ofdimension(cid:0)a−1(cid:1)×(a−1). Thethirdblock,M ,isdefinedlikewiseforfactorB.The
2 2
10