HW2 Solutions Log 2/26/03 > Strataux <- cbind.data.frame(factor(c(rep("Lo",24),rep("Hi",24))), factor(rep(c(rep("Lo",12),rep("Hi",12)),2)), factor(rep(c(rep("Lo",6),rep("Hi",6)),4)), factor(rep(c(rep(0,3),rep(1,3)),8)), factor(rep(c("Mob","Sng","Apt"),16))) > Strataux <-cbind.data.frame(Strataux,apply(StratMR[,,1]* StratMR[,,3],2,sum,na.rm=T),apply(StratMR[,,2]*StratMR[,,3], 2,sum,na.rm=T), apply(StratMR[,,3],2,sum,na.rm=T)) > names(Strataux) <- c("FB","FOWN","FSPOU","SNGF","HTYP","Nobs", "Nprdict","Cellct") ### Part (A) > Strataux FB FOWN FSPOU SNGF HTYP Nobs Nprdict Cellct 1 Lo Lo Lo 0 Mob 280409 2.798101e+05 480654 2 Lo Lo Lo 0 Sng 5826088 5.777823e+06 7747696 3 Lo Lo Lo 0 Apt 9327630 9.343762e+06 16565293 4 Lo Lo Lo 1 Mob 11647 1.177485e+04 24346 5 Lo Lo Lo 1 Sng 1510825 1.525286e+06 2220227 6 Lo Lo Lo 1 Apt 354218 3.499646e+05 647147 7 Lo Lo Hi 0 Mob 16171 1.616454e+04 26242 8 Lo Lo Hi 0 Sng 418904 4.260484e+05 529834 9 Lo Lo Hi 0 Apt 618742 6.070303e+05 905528 10 Lo Lo Hi 1 Mob 4335 4.100371e+03 8403 11 Lo Lo Hi 1 Sng 160829 1.624760e+05 234272 12 Lo Lo Hi 1 Apt 34419 3.270576e+04 51489 13 Lo Hi Lo 0 Mob 657803 6.614505e+05 1071545 14 Lo Hi Lo 0 Sng 1967305 1.973422e+06 2584631 15 Lo Hi Lo 0 Apt 900182 9.070258e+05 1520398 16 Lo Hi Lo 1 Mob 146947 1.473461e+05 268254 17 Lo Hi Lo 1 Sng 7670939 7.645492e+06 10124852 18 Lo Hi Lo 1 Apt 939319 9.551800e+05 1590689 19 Lo Hi Hi 0 Mob 669282 6.678720e+05 1105674 20 Lo Hi Hi 0 Sng 1334325 1.347321e+06 1828940 21 Lo Hi Hi 0 Apt 569624 5.622483e+05 936940 22 Lo Hi Hi 1 Mob 649754 6.508337e+05 1102193 23 Lo Hi Hi 1 Sng 19524109 1.952132e+07 24493384 24 Lo Hi Hi 1 Apt 949088 9.527199e+05 1493406 25 Hi Lo Lo 0 Mob 12555 1.200543e+04 30292 26 Hi Lo Lo 0 Sng 690491 7.116195e+05 1183169 27 Hi Lo Lo 0 Apt 1218300 1.215556e+06 2776401 28 Hi Lo Lo 1 Mob 2591 2.222638e+03 5912 29 Hi Lo Lo 1 Sng 430159 4.383632e+05 785382 30 Hi Lo Lo 1 Apt 88085 7.969790e+04 195271 31 Hi Lo Hi 0 Mob 45 4.272667e+01 76 32 Hi Lo Hi 0 Sng 3222 3.722507e+03 5051 33 Hi Lo Hi 0 Apt 19886 1.905023e+04 29978 34 Hi Lo Hi 1 Mob 15 1.126832e+01 33 35 Hi Lo Hi 1 Sng 365 3.668574e+02 541 36 Hi Lo Hi 1 Apt 16 2.459659e+01 43 37 Hi Hi Lo 0 Mob 20889 1.966105e+04 45111 38 Hi Hi Lo 0 Sng 95809 9.636233e+04 156109 39 Hi Hi Lo 0 Apt 33686 3.332685e+04 73556 40 Hi Hi Lo 1 Mob 19928 1.889220e+04 42565 41 Hi Hi Lo 1 Sng 1134877 1.138112e+06 1793386 42 Hi Hi Lo 1 Apt 79383 7.589318e+04 171281 43 Hi Hi Hi 0 Mob 811 6.243272e+02 1356 44 Hi Hi Hi 0 Sng 1688 1.756731e+03 2834 45 Hi Hi Hi 0 Apt 170 1.382131e+02 266 46 Hi Hi Hi 1 Mob 1173 1.060271e+03 2281 47 Hi Hi Hi 1 Sng 66302 6.699401e+04 88801 48 Hi Hi Hi 1 Apt 880 7.702011e+02 1295 > motif() attach(Strataux) par(mfrow=c(1,2), cex=1) pltchar <- rev(c(0,5,1,18)) > for (i in 1:2) { plot(Nprdict/Cellct,Nobs/Cellct,type="n", xlab="Predicted MR rates", ylab="Observed MR rates", main=paste("\n","for FOWN =",levels(FOWN)[i])) for (j in 1:2) { tmpind <- FOWN==levels(FOWN)[i] & SNGF==levels(SNGF)[j] points((Nprdict/Cellct)[tmpind],(Nobs/Cellct)[tmpind], pch=pltchar[2*i+j-2]) } legend(locator(), legend=c("SNGF=Hi","SNGF=Lo"), marks=pltchar[2*i-2+1:2]) } > par(mfg=c(1,1,1,1)) title("Observed vs Predicted MR by FOWN and SNGF") #### part (B) > Bigstrat <- (1:48)[Cellct>= 5000] > SmallFrm <- cbind.data.frame(Strataux,Nobs/Cellct)[Bigstrat[ordstrat],] > names(SmallFrm)[9] <- "MRrate" ### completes part (C) > Regions <- list(NEng=c(7,20,21,31,35,40,47), MidAtP=c(8,9,12, 22,32,39,46), NoCen=c(14,15,23,25,29,36,42,49,50), MidW=c(16,17,27,30,37), South=c(1,3,10,11,18,19,24,28,41,43,44), West=c(2,4,5,6,13,26,33,34,38,45,48,51)) > par(mfrow=c(3,2)) > unlist(lapply(Regions,length)) NEng MidAtP NoCen MidW South West 7 7 9 5 11 12 > for (i in 1:6) { inds <- (1:(length(Regions[[i]])*39))[!is.na(c(StratMR[ Regions[[i]],Bigstrat,2]))] plot(factor(rep(dimnames(StratMR)[[1]][ Regions[[i]]],39))[inds], c(StratMR[ Regions[[i]],Bigstrat,1]/StratMR[ Regions[[i]],Bigstrat,2])[inds], xlab="", ylab="Observed over Predicted", main= paste("Boxplot for",names(Regions)[i])) } ## Part (D) ** NOTE: inds vector used to remove missing values ### Resulting Boxplot for part (D) shows ratios of Nobs (1st level of ## 3rd dim of StratMR) over Nprdict (2nd level of 3rd dim of StratMR) ## Corresponding boxplot for array aggregated over covariates ## is much less interesting !! ** No regional patterns are very persuasive so far !! #### Part (E) > plot(Strataux[Bigstrat,"HTYP"],Strataux[Bigstrat,"Nobs"]/ Strataux[Bigstrat,"Nprdict"], xlab="Housing Type", ylab="Obs/Predicted MR Count", main="Boxplot of Model Discrepancies by HTYP") ## Now produce similar boxplot by different covariates, by region > ByRegion <- function(Covnam = "HTYP") { par(mfrow = c(3, 2)) for(j in 1:6) { reglgth <- length(Regions[[j]]) inds <- (1:(reglgth * 39))[!is.na(c(StratMR[ Regions[[j]], Bigstrat, 2]))] plot(c(matrix(rep(Strataux[Bigstrat, Covnam], reglgth), nrow = reglgth, byrow = T))[inds], c(StratMR[Regions[[j]], Bigstrat, 1]/StratMR[ Regions[[j]], Bigstrat, 2])[inds], xlab = "", ylab = "Observed over Predicted", main = paste( "Boxplot for", Covnam, "in", names(Regions)[j])) } } > ByRegion("HTYP") > ByRegion("FB") > ByRegion("FOWN") > ByRegion("FSPOU") > ByRegion("SNGF"). ## Of these, the one for HTYP is most interesting and shows: ratios are very tight for Apt and Sng in South, for Sng in MidAtP. In Neng there is a stratum with Sng ratio >>1, otherwise, the extremes are usually <<1 for Sng. Mob ratios always very widely dispersed. Ratios for Apt skewed up in MidAtP, down in West. ## Patterns not too clear, but that's life !!