HW 19 analysis. 4/2/08 ======================== > unlist(lapply(OrchardSprays,class)) decrease rowpos colpos treatment "numeric" "numeric" "numeric" "factor" > attach(OrchardSprays) > fit1 = lm(decrease ~ treatment) > summary(fit1)$coef Estimate Std. Error t value Pr(>|t|) (Intercept) 4.625 7.253328 0.6376383 5.263078e-01 treatmentB 3.000 10.257755 0.2924617 7.710145e-01 treatmentC 20.625 10.257755 2.0106739 4.918440e-02 treatmentD 30.375 10.257755 2.9611742 4.488728e-03 treatmentE 58.500 10.257755 5.7030022 4.600226e-07 treatmentF 64.375 10.257755 6.2757396 5.391817e-08 treatmentG 63.875 10.257755 6.2269960 6.479213e-08 treatmentH 85.625 10.257755 8.3473430 2.080424e-11 > fit2 = lm(decrease ~ treatment + factor(rowpos) + factor(colpos)) > summary(fit2)$coef Estimate Std. Error t value Pr(>|t|) (Intercept) 29.90625 11.441621 2.6138123 1.237639e-02 treatmentB 3.00000 9.757447 0.3074575 7.600146e-01 treatmentC 20.62500 9.757447 2.1137701 4.051552e-02 treatmentD 30.37500 9.757447 3.1130069 3.328111e-03 treatmentE 58.50000 9.757447 5.9954206 4.032712e-07 treatmentF 64.37500 9.757447 6.5975248 5.487335e-08 treatmentG 63.87500 9.757447 6.5462819 6.501307e-08 treatmentH 85.62500 9.757447 8.7753486 4.767062e-11 factor(rowpos)2 -10.12500 9.757447 -1.0376690 3.053616e-01 factor(rowpos)3 -11.37500 9.757447 -1.1657762 2.502827e-01 factor(rowpos)4 -24.37500 9.757447 -2.4980919 1.648724e-02 factor(rowpos)5 -25.62500 9.757447 -2.6261992 1.199721e-02 factor(rowpos)6 -24.87500 9.757447 -2.5493348 1.453363e-02 factor(rowpos)7 -23.00000 9.757447 -2.3571739 2.315095e-02 factor(rowpos)8 -18.25000 9.757447 -1.8703663 6.841262e-02 factor(colpos)2 2.25000 9.757447 0.2305931 8.187505e-01 factor(colpos)3 -13.87500 9.757447 -1.4219908 1.624131e-01 factor(colpos)4 -8.25000 9.757447 -0.8455080 4.026185e-01 factor(colpos)5 -14.50000 9.757447 -1.4860444 1.447359e-01 factor(colpos)6 -14.00000 9.757447 -1.4348015 1.587486e-01 factor(colpos)7 -14.00000 9.757447 -1.4348015 1.587486e-01 factor(colpos)8 -2.25000 9.757447 -0.2305931 8.187505e-01 > table(treatment) treatment A B C D E F G H 8 8 8 8 8 8 8 8 > fit3 = lm(decrease ~ treatment + I(rowpos > 3) -1) > summary(fit3)$coef Estimate Std. Error t value Pr(>|t|) treatmentA 14.66146 7.355590 1.993240 5.120501e-02 treatmentB 17.66146 7.355590 2.401093 1.975440e-02 treatmentC 35.28646 7.355590 4.797230 1.269226e-05 treatmentD 45.03646 7.355590 6.122753 1.018439e-07 treatmentE 73.16146 7.355590 9.946376 6.763870e-14 treatmentF 79.03646 7.355590 10.745088 3.993877e-15 treatmentG 78.53646 7.355590 10.677112 5.066516e-15 treatmentH 100.28646 7.355590 13.634047 2.779293e-19 I(rowpos > 3)TRUE -16.05833 4.886788 -3.286071 1.772441e-03 > label = c("A","B","C","D","E","F","G","H") > detach() > contrasts(OrchardSprays$treatment) B C D E F G H A 0 0 0 0 0 0 0 B 1 0 0 0 0 0 0 C 0 1 0 0 0 0 0 D 0 0 1 0 0 0 0 E 0 0 0 1 0 0 0 F 0 0 0 0 1 0 0 G 0 0 0 0 0 1 0 H 0 0 0 0 0 0 1 > contrasts(OrchardSprays$treatment) = matrix( as.numeric(outer(1:8,1:7, function(x,y) x>y)), ncol=7, dimnames=list(label,label[2:8])) > contrasts(OrchardSprays$treatment) B C D E F G H A 0 0 0 0 0 0 0 B 1 0 0 0 0 0 0 C 1 1 0 0 0 0 0 D 1 1 1 0 0 0 0 E 1 1 1 1 0 0 0 F 1 1 1 1 1 0 0 G 1 1 1 1 1 1 0 H 1 1 1 1 1 1 1 > fit1B = lm(decrease ~ treatment, data=OrchardSprays) > summary(fit1B)$coef Estimate Std. Error t value Pr(>|t|) (Intercept) 4.625 7.253328 0.63763832 0.526307813 treatmentB 3.000 10.257755 0.29246165 0.771014549 treatmentC 17.625 10.257755 1.71821220 0.091282915 treatmentD 9.750 10.257755 0.95050037 0.345942415 treatmentE 28.125 10.257755 2.74182798 0.008187303 treatmentF 5.875 10.257755 0.57273740 0.569115065 treatmentG -0.500 10.257755 -0.04874361 0.961297055 treatmentH 21.750 10.257755 2.12034697 0.038420180 ### Siginificance tests given essentially by the last column here: ### only the E vs D and H vs G contrasts look individually ### significant, although I am sure that an overall increasing ### linear trend would be significant in addition to these two. > fit3B = lm(decrease ~ treatment + I(rowpos > 3), data=OrchardSprays) > summary(fit3B)$coef Estimate Std. Error t value Pr(>|t|) (Intercept) 14.66146 7.355590 1.99324034 0.051205014 treatmentB 3.00000 9.463224 0.31701669 0.752431199 treatmentC 17.62500 9.463224 1.86247308 0.067881174 treatmentD 9.75000 9.463224 1.03030426 0.307376845 treatmentE 28.12500 9.463224 2.97203151 0.004382031 treatmentF 5.87500 9.463224 0.62082436 0.537279564 treatmentG -0.50000 9.463224 -0.05283612 0.958053929 treatmentH 21.75000 9.463224 2.29837104 0.025366494 I(rowpos > 3)TRUE -16.05833 4.886788 -3.28607129 0.001772441 ### So the indication of high row number is highly significant! ### Variance matrix for coefficient estimators are given as ### summary(fit3B)$cov.uns*summary(fit3B)$sigma^2, similarly > round(summary(fit1B)$cov.uns*summary(fit1B)$sigma^2,3) (Intercept) treatmentB treatmentC treatmentD treatmentE treatmentF (Intercept) 52.611 -52.611 0.000 0.000 0.000 0.000 treatmentB -52.611 105.222 -52.611 0.000 0.000 0.000 treatmentC 0.000 -52.611 105.222 -52.611 0.000 0.000 treatmentD 0.000 0.000 -52.611 105.222 -52.611 0.000 treatmentE 0.000 0.000 0.000 -52.611 105.222 -52.611 treatmentF 0.000 0.000 0.000 0.000 -52.611 105.222 treatmentG 0.000 0.000 0.000 0.000 0.000 -52.611 treatmentH 0.000 0.000 0.000 0.000 0.000 0.000 treatmentG treatmentH (Intercept) 0.000 0.000 treatmentB 0.000 0.000 treatmentC 0.000 0.000 treatmentD 0.000 0.000 treatmentE 0.000 0.000 treatmentF -52.611 0.000 treatmentG 105.222 -52.611 treatmentH -52.611 105.222 ### Output text file: > write.table(OrchardSprays,"OrchSprays.txt", row.names=F) ## remove quotes using word-processor >>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>> ### Now to re-do it all in SAS. libname home "."; options nocenter nodate linesize=80; data home.OrchSpray (drop=i); infile "OrchSprays.txt"; input decr rowpos colpos treat $ ; array RX{7} ; do i = 1 to 7; RX(i)=0; end; select ; when (treat="B") RX(1)=1; when (treat="C") do i=1 to 2; RX(i)=1; end; when (treat="D") do i=1 to 3; RX(i)=1; end; when (treat="E") do i=1 to 4; RX(i)=1; end; when (treat="F") do i=1 to 5; RX(i)=1; end; when (treat="G") do i=1 to 6; RX(i)=1; end; when (treat="H") do i=1 to 7; RX(i)=1; end; otherwise; end; if decr ne . ; run; PROC REG data=home.OrchSpray COVOUT outest=varmat; model decr = RX1-RX7 ; run; Parameter Estimates Parameter Standard Variable DF Estimate Error t Value Pr > |t| Intercept 1 4.62500 7.25333 0.64 0.5263 RX1 1 3.00000 10.25776 0.29 0.7710 RX2 1 17.62500 10.25776 1.72 0.0913 RX3 1 9.75000 10.25776 0.95 0.3459 RX4 1 28.12500 10.25776 2.74 0.0082 RX5 1 5.87500 10.25776 0.57 0.5691 RX6 1 -0.50000 10.25776 -0.05 0.9613 RX7 1 21.75000 10.25776 2.12 0.0384 PROC PRINT data=varmat; where _NAME_ ne ""; Var Intercept RX1-RX7; run; Obs Intercept RX1 RX2 RX3 RX4 RX5 RX6 RX7 2 52.6108 -52.611 0.000 -0.000 0.000 0.000 -0.000 0.000 3 -52.6108 105.222 -52.611 0.000 -0.000 0.000 -0.000 0.000 4 0.0000 -52.611 105.222 -52.611 0.000 -0.000 0.000 0.000 5 -0.0000 0.000 -52.611 105.222 -52.611 -0.000 0.000 0.000 6 0.0000 -0.000 0.000 -52.611 105.222 -52.611 -0.000 0.000 7 0.0000 0.000 -0.000 -0.000 -52.611 105.222 -52.611 0.000 8 -0.0000 -0.000 0.000 0.000 -0.000 -52.611 105.222 -52.611 9 0.0000 0.000 0.000 0.000 0.000 0.000 -52.611 105.222