HSB data some useful advanced plots using xyplot from lattice graphics #read in prepared Bryk data > brykD = read.table("http://statweb.stanford.edu/~rag/stat209/brykdataPrepared", header = T) ## note it is safer to use statweb.stanford.edu in read.table even if the link URL is www-stat.stanford.edu > attach(brykD) > table(school) school 1224 1288 1296 1308 1317 1358 1374 1433 1436 1461 1462 1477 1499 1637 1906 1909 1942 1946 2030 2208 2277 2305 2336 2458 2467 2526 2626 2629 2639 2651 2655 2658 47 25 48 20 48 30 28 35 44 33 57 62 53 27 53 28 29 39 47 60 61 67 47 57 52 57 38 57 42 38 52 45 2755 2768 2771 2818 2917 2990 2995 3013 3020 3039 3088 3152 3332 3351 3377 3427 3498 3499 3533 3610 3657 3688 3705 3716 3838 3881 3967 3992 3999 4042 4173 4223 47 25 55 42 43 48 46 53 59 21 39 52 38 39 45 49 53 38 48 64 51 43 45 41 54 41 52 53 46 64 44 45 4253 4292 4325 4350 4383 4410 4420 4458 4511 4523 4530 4642 4868 4931 5192 5404 5619 5640 5650 5667 5720 5761 5762 5783 5815 5819 5838 5937 6074 6089 6144 6170 58 65 53 33 25 41 32 48 58 47 63 61 34 58 28 57 66 57 45 61 53 52 37 29 25 50 31 29 56 33 43 21 6291 6366 6397 6415 6443 6464 6469 6484 6578 6600 6808 6816 6897 6990 7011 7101 7172 7232 7276 7332 7341 7342 7345 7364 7635 7688 7697 7734 7890 7919 8009 8150 35 58 60 54 30 29 57 35 56 56 44 55 49 53 33 28 44 52 53 48 51 58 56 44 51 54 32 22 51 37 47 44 8165 8175 8188 8193 8202 8357 8367 8477 8531 8627 8628 8707 8775 8800 8854 8857 8874 8946 8983 9021 9104 9158 9198 9225 9292 9340 9347 9359 9397 9508 9550 9586 49 33 30 43 35 27 14 37 41 53 61 48 48 32 32 64 36 58 51 56 55 53 31 36 19 29 57 53 47 35 29 59 > table(school[sector == "Public"]) 1224 1288 1296 1358 1374 1461 1499 1637 1909 1942 1946 2030 2336 2467 2626 2639 2651 2655 2768 2771 2818 2917 2995 3013 3088 3152 3332 3351 3377 3657 3716 3881 47 25 48 30 28 33 53 27 28 29 39 47 47 52 38 42 38 52 25 55 42 43 46 53 39 52 38 39 45 51 41 41 3967 3999 4325 4350 4383 4410 4420 4458 4642 5640 5762 5783 5815 5819 5838 5937 6089 6144 6170 6291 6397 6415 6443 6464 6484 6600 6808 6897 6990 7101 7232 7276 52 46 53 33 25 41 32 48 61 57 37 29 25 50 31 29 33 43 21 35 60 54 30 29 35 56 44 49 53 28 52 53 7341 7345 7697 7734 7890 7919 8175 8188 8202 8357 8367 8477 8531 8627 8707 8775 8854 8874 8946 8983 9158 9225 9292 9340 9397 9550 51 56 32 22 51 37 33 30 35 27 14 37 41 53 48 48 32 36 58 51 53 36 19 29 47 29 > table(school[sector == "Catholic"]) 1308 1317 1433 1436 1462 1477 1906 2208 2277 2305 2458 2526 2629 2658 2755 2990 3020 3039 3427 3498 3499 3533 3610 3688 3705 3838 3992 4042 4173 4223 4253 4292 20 48 35 44 57 62 53 60 61 67 57 57 57 45 47 48 59 21 49 53 38 48 64 43 45 54 53 64 44 45 58 65 4511 4523 4530 4868 4931 5192 5404 5619 5650 5667 5720 5761 6074 6366 6469 6578 6816 7011 7172 7332 7342 7364 7635 7688 8009 8150 8165 8193 8628 8800 8857 9021 58 47 63 34 58 28 57 66 45 61 53 52 56 58 57 56 55 33 44 48 58 44 51 54 47 44 49 43 61 32 64 56 9104 9198 9347 9359 9508 9586 55 31 57 53 35 59 > library(lattice) > library(lme4) #A. Ignore schools, and plot all Cath and Public student with an "overall" regression smoother xyplot(mathach ~ cses | sector, brykD, type = c("g", "p", "smooth"), aspect = .5) # another version all in one frame xyplot(mathach ~ cses, groups = sector, brykD, type = c("g", "p", "smooth"), aspect = .5) # now overlay straight-line regression (ignoring schools) instead of smoother xyplot(mathach ~ cses, groups = sector, brykD, type = c("g", "p", "r"), aspect = .5, col = c("black", "red")) # you can see, even ignoring schools, Pub (red) is steeper and lower in mean ach than Cath # B. Do school by school plots xyplot(mathach ~ cses | school, brykD, type = c("g","p","r"), subset = sector == "Catholic", xlab = "C_centered ses", ylab = "C_Math", aspect = .5) # this gives frame-by-frame for 70 Cath schools with straight-line fit xyplot(mathach ~ cses | school, brykD, type = c("g","p","r"), subset = sector == "Public", xlab = "centered ses", ylab = "Math", aspect = .5) # same for 90 public schools # These plots look fine if you can stretch them out across a couple big monitors. # As an alternative try doing these for the subset of the first 26 schools or sample from group of schools (e.g. Bates in mlmREV Exam data ex) > xyplot(mathach ~ cses | school, brykD[1:1154,], type = c("g","p","r"), subset = sector == "Public", xlab = "centered ses", ylab = "Math", aspect = .5) > xyplot(mathach ~ cses | school, brykD[1:1154,], type = c("g","p","r"), subset = sector == "Catholic", xlab = "centered ses", ylab = "Math", aspect = .5) # C. Collection of school fits # plots of school fits for each sector xyplot(mathach ~ cses, groups = school, brykD, type = c("r"), subset = sector == "Catholic", xlab = "C_centered ses", ylab = "C_Math", aspect = .5, col = c("black")) xyplot(mathach ~ cses, groups = school, brykD, type = c("r"), subset = sector == "Public", xlab = "P_centered ses", ylab = "P_Math", aspect = .5, col = c("black")) # or do 2 frames on a page, with or without data points xyplot(mathach ~ cses | sector, groups = school, data = brykD, type = c("g","p","r"), index = function(x,y) coef(lm(y ~ x))[1], xlab = "centered ses", ylab = "Math", aspect = .5, col = c("black")) xyplot(mathach ~ cses | sector, groups = school, data = brykD, type = c("r"), index = function(x,y) coef(lm(y ~ x))[1], xlab = "centered ses", ylab = "Math", aspect = .5, col = c("black")) # D. Confint plots for lmList fits > cathRegC = lmList(mathach ~ cses| school, subset = sector == "Catholic", data = brykD) > pubRegC = lmList(mathach ~ cses| school, subset = sector == "Public", data = brykD) > length(cathRegC); length(pubRegC) [1] 70 [1] 90 > pubcoef= coef(pubRegC) > cathcoef= coef(cathRegC) # plots that are often shown in expositions, for each sector > plot(confint(cathRegC, pool = TRUE), order = 1) > > plot(confint(pubRegC, pool = TRUE), order = 1) ----------------------- END extra plots section