options ps=54 ls=80; title 'RPD - Example 8.5 - Response Surface'; proc iml; xy = { 1 -1 -1 -1 53, 1 -1 -1 1 54, 1 -1 1 -1 40, 1 -1 1 1 37, 1 1 -1 -1 84, 1 1 -1 1 76, 1 1 1 -1 40, 1 1 1 1 50, 1 0 0 0 50, 1 1.215 0 0 61, 1 -1.215 0 0 54, 1 0 1.215 0 39, 1 0 -1.215 0 67, 1 0 0 1.215 44, 1 0 0 -1.215 61, 2 -1 -1 -1 50, 2 -1 -1 1 42, 2 -1 1 -1 31, 2 -1 1 1 28, 2 1 -1 -1 57, 2 1 -1 1 78, 2 1 1 -1 49, 2 1 1 1 54, 2 0 0 0 50, 2 1.215 0 0 76, 2 -1.215 0 0 45, 2 0 1.215 0 33, 2 0 -1.215 0 54, 2 0 0 1.215 45, 2 0 0 -1.215 38, 3 -1.2500 -1.8867 -0.6350 46, 3 0.8600 -2.2200 -0.4250 66, 3 1.0000 -2.2400 -0.3100 68, 3 2.1165 -2.4167 -0.1450 75, 3 2.5825 -2.4900 -0.0800 75, 3 3.2475 -2.6667 0.0800 68, 3 1.1760 -1.3333 0 78, 3 1.4700 -1.6667 0 93, 3 1.7640 -2.0000 0 96, 3 2.0580 -2.3333 0 66}; y = xy(|,5|); x01 = j(15,1,1)//j(25,1,0); x02 = j(15,1,0)//j(15,1,1)//j(10,1,0); x03 = j(30,1,0)//j(10,1,1); x1 = xy(|,2|); x2 = xy(|,3|); x3 = xy(|,4|); x11 = x1#x1; x22 = x2#x2; x33 = x3#x3; x12 = x1#x2; x13 = x1#x3; x23 = x2#x3; xf = x01||x02||x03||x1||x2||x3||x11||x22||x33||x12||x13||x23; betaf = inv(xf` * xf) * xf` * y; yhatf = xf * betaf; ssef = (y-yhatf)` * (y-yhatf); dfef = nrow(xf)-ncol(xf); s2f = ssef/dfef; varbetaf = s2f*inv(xf` * xf); sebetaf = sqrt(diag(varbetaf)*j(ncol(xf),1,1)); tbetaf = betaf/sebetaf; print betaf sebetaf tbetaf; print s2f dfef; xr = x01||x02||x03||x1||x2||x22; betar = inv(xr` * xr) * xr` * y; yhatr = xr * betar; sser = (y-yhatr)` * (y-yhatr); dfer = nrow(xr)-ncol(xr); s2r = sser/dfer; varbetar = s2r*inv(xr` * xr); sebetar = sqrt(diag(varbetar)*j(ncol(xr),1,1)); tbetar = betar/sebetar; print betar sebetar tbetar; print s2r dfer; ffullred = ((sser-ssef)/(dfer-dfef))/(ssef/dfef); ffrcrit=finv(.95,dfer-dfef,dfef); probffr=1-probf(ffullred,dfer-dfef,dfef); print ffullred ffrcrit probffr; run; quit;