multreg

Multiple Regression by Iterative Fitting(EDTTS, Ch.7) We can use the same principle with more variables:(1a) Fit y using x1(2a) Fit residuals from step (1a) to x2(1b) Fit residuals from step (2a) to x1(2b) Fit residuals from step (1b) to x2(1c) Fit residuals from step (2b) to x1(2c) Fit residuals from step (1c) to x2and so on, until the coefficients are effectively zero.When fitting by least squares, the final residuals will be the same as the least squares residuals.Final coefficients must be calculated:Intercept: add intercept terms together from all regressions for x1:add coefficients of x1from (1a), (2a), (3a), ...for x2:add coefficients of x2from (1b), (2b), (3b), ... NB:This is a very slow way to obtain least squares fitA more direct way of finding this least squares fit that does not require iteration is to adjust each variable for the suc-cessive ones as described in EDTTS (Ch.7).We use this same "adjusted iterative fitting" when we replace the LS fit by the RR fit.Multiple Regression by Adjusted Iterative Fitting "Sweep out" one variable at a timeBasic idea:(1) Adjust y for x1→y⋅1(2) Adjust x2for x1→x2⋅1(3) Adjust y⋅1for x2⋅1Now y is adjusted for both x1and x2.Also, see how much of x1is in x2.See how this results in LS line:(1) Fit y to x1and form residuals:ˆy=c0+c1x1y⋅1=y−c1x1c1=ni=1Σ(x i1−x1)y i/ni=1Σ(x i1−x1)2(2) Fit x2to x1and form residuals:ˆx2=d0+d1x1x2⋅1=x2−d1x1d1=ni=1Σ(x i1−x1)x i2/ni=1Σ(x i1−x1)2(3) Fit y⋅1to x2⋅1and form residuals:ˆy⋅1=c0+c2x2⋅1y⋅12=y⋅1−c2x2⋅1c2=ni=1Σ(x2⋅1,i−x2⋅1)y⋅1,i/ni=1Σ(x2⋅1,i−x2⋅1)2Least squares fit isˆy=b0+b1x1+b2x2,whereb2=c2b1=c1−d1c2b0=mean(y i−b1x i1−b2x i2).Note:If you plot y⋅12against x1,the least squares slope of the line will be zero; likewise for x2.Everywhere above,replace "LS fit" with "RR line fit", but the "Note"no longer applies, so must add iterative steps: (4) Plot y⋅12against x1.If slope is not "small", adjust c1slope accordingly,forming new c1,y⋅1,c2,y⋅12. (5) Plot y⋅12against x2.If slope is not "small", adjust c1slope accordingly,return to (4).Otherwise, go to (6).(6) Obtainfinal residuals from step (5);obtain final fit from either(data−residuals)or from assembling all coefficients from all stepsor fromˆy=b0+b1+x1+b2x2whereb2=c2b1=c1−d1c2b0=median(y1−b1x i1−b2x i2).EDTTS, p.251: adjusted iterative fitting:(1) Fit y to x1=Poverty:y=87.54+1. 72x1+y⋅1(2) Fit x4=Infant mortality rate to x1:x4=89.82+1. 15x1+x4⋅1(3) Fit y⋅1to x4⋅1:y⋅1=−4. 34−0. 329x4⋅1+y⋅14 Assemble:y=87.54+1. 72x1−4. 34−0. 329(x4−89.82−1. 15x1)+y⋅14 =11 2.75+2. 10x1−0. 329x4+y⋅14Note:The correlation between x1and x4is high -- 0.57 --so unadjusted least squares fitting requires nine steps to converge to the least squares solution!(When the correla-tion is only 0.30, as for x1and x2,only six steps are needed.) We reduce the need for extensive iteration by adjusting variables as we go along.Same decomposition as before:data=fit1+residual1=fit2+residual2=fit3+residual3with final decompositiondata=(fit1+fit2+fit3)+residual3Here,fit2is obtained by fitting x2.1to residual1,and so on, as variables are used in turn.The results are not the same as fitting x2to residual1,unless x1is uncorrelated with x2 (there is no component of x1in x2).General procedure, > 2 independent variables:(1) Fit y to x1→y=c10+c1x1+y⋅1(2) Fit x2to x1→x2=d20+d21x1+x2⋅1(1) + (2)→(3):(3) Fit y⋅1to x2⋅1→y⋅1=c21+c2x2⋅1+y⋅12(4) Fit x3to x1and x2⋅1→x3=d30+d31x1+d32x2⋅1+x3⋅12(3) + (4)→(5)(5) Fit y⋅12to x3⋅12→y⋅12=c30+c3x3⋅12+y⋅123(6) Fit x4to x1,x2⋅1,x3⋅12→x4=d40+d41x1+d42x2⋅1+d43x3⋅12+x4⋅123(5) + (6)→(7):(7) Fit y⋅123to x4⋅123→y⋅123=c40+c4x4⋅123+y⋅1234(8) Fit x5to x1,x2⋅1,x3⋅12,x4⋅123→(etc.)Assemble coefficients (for intercepts,x1,x2,etc.)Strategies for removing variables:(1) Plot y versus x k,in turn; choose that x k for whichrelationship looks strongest.Note: some plots maysuggest transforming x kfirst; e.g.,1/income,log(rate).(2) Regress (rank of y i)on(rank of x ik), and choose thatx k for which the slope coefficient is largest (use ofranks reduces impact of outliers on the LS regression-- but there could be many ties for the "most signifi-cant" if there are only few data points and severalvariables have monotonic relationships with y)Some important facts to remember:(1) Thefit depends on what variables are in the equation.Omission of a variable can make drastic changes infit -- even rev e rse sign of coefficients(cf. fit of y⋅some to x k⋅some).(2) The LS coefficient b k can be obtained by regressingy⋅(all but k)on x k⋅(all but k).So a plot of y⋅(all but k)onx k⋅(all but k)is more informative than y on x k.Robust Regression: "Iteratively reweighted LS"y=Xβ+epislon=fit+roughˆβ=(X′X)−1X′ylinear function of y=(y1,... ,y n)′Weighted least squares:ˆβ=(X′WX)−1X′WyWWeight matrix W=diag(w1,...,w n)w i=BIG if residual is small(point i is "consistent" with "model")w i=little if resid-ual is large(point i is "inconsistent" with "model")Deciding on BIG and little residuals:•s_r=measure of spread of residuals (e.g., ...)•BIG residual if|resid|>6⋅s_r:w i=[1−(resid/6s2_r)]2|resid|<6s r(*)=0otherwiseIteratively re-weighted LS:1. Fit OLS→ˆβ(0)→r(0)=y−Xˆβ(0)2Calculate W(1 )[i,i]using r(0)i for resid for w i(*)3. Fitˆβvia WLS (weighted LS):→ˆβ(1 )→r(1 )=y−Xˆβ(1 )4Calculate W(2 )[i,i]using r(1 )i for resid for w i5. Fitˆβvia WLS using W(2 ):→ˆβ(2 )→r(2 )=y−Xˆβ(2 )|beta(k+1)−ˆβ(k)|<τ(tolerance)Iterate untilˆNotes:1. Starting with OLS is NOT robust, but it’s easy.2. Eqn(*) is called bisquare weight function:w i=1when residual=0w i=0when |residual|>6s_rfalls smoothly to zero3. The value "6" is a "tuning parameter"monly4≤c≤6.4. Any weight function that smoothly falls to zeroyields robust estimate ofβ.Biweight regression isespecially efficient.5. Easy to do computationally.See rlm(MASS):rlm(formula, data, weights, ..., subset, na.action,method = c("M", "MM", "model.frame"), ...)rlm(x, y,weights, w = rep(1, nrow(x)), init = "ls",-10-psi = psi.huber,method = c("M", "MM"), ... )method: currently either M-estimation or MM-estimationor (for the ’formula’ method only) find the model frame. MM-estimation is M-estimation with Tukey’s biweight initialized by a specific S-estimator.See ’Details’Psi functions are supplied for the Huber,Hampel and Tukey bisquare proposals as ’psi.huber’, ’psi.hampel’, ‘psi.bisquare’. Selecting ’method = "MM"’ selects a specific set of options which ensures that the estimator has a high breakdown point. The initial set of coefficients and the final scale are selected by an S-estimator with ’k0 = 1.548’; this gives(for n >> p) breakdown point 0.5. The final estimator isan M-estimator with Tukey’s biweight and fixed scale that will inherit this breakdown point provided ’c > k0’; thisis true for the default value of ’c’ that corresponds to95% relative efficiency at the normal.Case weights arenot supported for ’method = "MM"’.。