R Under development (unstable) (2025-01-28 r87664 ucrt) -- "Unsuffered Consequences" Copyright (C) 2025 The R Foundation for Statistical Computing Platform: x86_64-w64-mingw32/x64 R is free software and comes with ABSOLUTELY NO WARRANTY. You are welcome to redistribute it under certain conditions. Type 'license()' or 'licence()' for distribution details. R is a collaborative project with many contributors. Type 'contributors()' for more information and 'citation()' on how to cite R or R packages in publications. Type 'demo()' for some demos, 'help()' for on-line help, or 'help.start()' for an HTML browser interface to help. Type 'q()' to quit R. > > library("MVA") Loading required package: HSAUR2 Loading required package: tools > > sapply(demo(package = "MVA")$results[, "Item"], demo, character.only = TRUE, + ask = FALSE) demo(Ch-CA) ---- ~~~~~ > ### R code from vignette source 'Ch-CA.Rnw' > > ################################################### > ### code chunk number 1: setup > ################################################### > library("MVA") > set.seed(280875) > library("lattice") > lattice.options(default.theme = + function() + standard.theme("pdf", color = FALSE)) > if (file.exists("deparse.R")) { + if (!file.exists("figs")) dir.create("figs") + source("deparse.R") + options(prompt = "R> ", continue = "+ ", width = 64, + digits = 4, show.signif.stars = FALSE, useFancyQuotes = FALSE) + + options(SweaveHooks = list(onefig = function() {par(mfrow = c(1,1))}, + twofig = function() {par(mfrow = c(1,2))}, + figtwo = function() {par(mfrow = c(2,1))}, + threefig = function() {par(mfrow = c(1,3))}, + figthree = function() {par(mfrow = c(3,1))}, + fourfig = function() {par(mfrow = c(2,2))}, + sixfig = function() {par(mfrow = c(3,2))}, + nomar = function() par("mai" = c(0, 0, 0, 0)))) + } > ################################################### > ### code chunk number 2: ch:CA:data > ################################################### > measure <- + structure(list(V1 = 1:20, V2 = c(34L, 37L, 38L, 36L, 38L, 43L, + 40L, 38L, 40L, 41L, 36L, 36L, 34L, 33L, 36L, 37L, 34L, 36L, 38L, + 35L), V3 = c(30L, 32L, 30L, 33L, 29L, 32L, 33L, 30L, 30L, 32L, + 24L, 25L, 24L, 22L, 26L, 26L, 25L, 26L, 28L, 23L), V4 = c(32L, + 37L, 36L, 39L, 33L, 38L, 42L, 40L, 37L, 39L, 35L, 37L, 37L, 34L, + 38L, 37L, 38L, 37L, 40L, 35L)), .Names = c("V1", "V2", "V3", + "V4"), class = "data.frame", row.names = c(NA, -20L)) > measure <- measure[,-1] > names(measure) <- c("chest", "waist", "hips") > measure$gender <- gl(2, 10) > levels(measure$gender) <- c("male", "female") > ################################################### > ### code chunk number 3: ch:CA-scplot > ################################################### > library("mvtnorm") > dat <- rbind(rmvnorm(25, mean = c(3,3)), + rmvnorm(20, mean = c(10, 8)), + rmvnorm(10, mean = c(20, 1))) > plot(abs(dat), xlab = expression(x[1]), ylab = expression(x[2])) > ################################################### > ### code chunk number 4: ch:CA-dist > ################################################### > set.seed(29) > x1 <- c(0.7, 0.8, 0.85, 0.9, 1.1, 1, 0.95) > x <- c(x1, x1 + 1.5) > y1 <- sample(x1) > y <- c(y1, y1 + 1) > plot(x, y, main = "single") > lines(c(1, 0.7 + 1.5), c(1.1, 0.7 + 1), col = "grey") > set.seed(29) > x1 <- c(0.7, 0.8, 0.85, 0.9, 1.1, 1, 0.95) > x <- c(x1, x1 + 1.5) > y1 <- sample(x1) > y <- c(y1, y1 + 1) > plot(x, y, main = "complete") > lines(c(0.7, 2.5), c(0.7, 1.1 + 1), col = "grey") > set.seed(29) > x1 <- c(0.7, 0.8, 0.85, 0.9, 1.1, 1, 0.95) > x <- c(x1, x1 + 1.5) > y1 <- sample(x1) > y <- c(y1, y1 + 1) > plot(x, y, main = "average") > for (i in 1:7) { + for (j in 8:14) lines(x[c(i, j)], y[c(i, j)], col = rgb(0.1, 0.1, 0.1, 0.1)) + } > ################################################### > ### code chunk number 5: ch:CA:measure (eval = FALSE) > ################################################### > ## (dm <- dist(measure[, c("chest", "waist", "hips")])) > > > ################################################### > ### code chunk number 6: ch:CA:measure > ################################################### > dm <- dist(measure[, c("chest", "waist", "hips")]) > round(dm, 2) 1 2 3 4 5 6 7 8 9 10 11 12 2 6.16 3 5.66 2.45 4 7.87 2.45 4.69 5 4.24 5.10 3.16 7.48 6 11.00 6.08 5.74 7.14 7.68 7 12.04 5.92 7.00 5.00 10.05 5.10 8 8.94 3.74 4.00 3.74 7.07 5.74 4.12 9 7.81 3.61 2.24 5.39 4.58 3.74 5.83 3.61 10 10.10 4.47 4.69 5.10 7.35 2.24 3.32 3.74 3.00 11 7.00 8.31 6.40 9.85 5.74 11.05 12.08 8.06 7.48 10.25 12 7.35 7.07 5.48 8.25 6.00 9.95 10.25 6.16 6.40 8.83 2.24 13 7.81 8.54 7.28 9.43 7.55 12.08 11.92 7.81 8.49 10.82 2.83 2.24 14 8.31 11.18 9.64 12.45 8.66 14.70 15.30 11.18 11.05 13.75 3.74 5.20 15 7.48 6.16 4.90 7.07 6.16 9.22 9.00 4.90 5.74 7.87 3.61 1.41 16 7.07 6.00 4.24 7.35 5.10 8.54 9.11 5.10 5.00 7.48 3.00 1.41 17 7.81 7.68 6.71 8.31 7.55 11.40 10.77 6.71 7.87 9.95 3.74 2.24 18 6.71 6.08 4.58 7.28 5.39 9.27 9.49 5.39 5.66 8.06 2.83 1.00 19 9.17 5.10 4.47 5.48 7.07 6.71 5.74 2.00 4.12 5.10 6.71 4.69 20 7.68 9.43 7.68 10.82 7.00 12.41 13.19 9.11 8.83 11.53 1.41 3.00 13 14 15 16 17 18 19 2 3 4 5 6 7 8 9 10 11 12 13 14 3.74 15 3.00 6.40 16 3.61 6.40 1.41 17 1.41 5.10 2.24 3.32 18 2.83 5.83 1.00 1.00 2.45 19 6.40 9.85 3.46 3.74 5.39 4.12 20 2.45 2.45 4.36 4.12 3.74 3.74 7.68 > ################################################### > ### code chunk number 7: ch:CA:measure:plot > ################################################### > plot(cs <- hclust(dm, method = "single")) > plot(cc <- hclust(dm, method = "complete")) > plot(ca <- hclust(dm, method = "average")) > ################################################### > ### code chunk number 8: ch:CA:measure:plotplot > ################################################### > body_pc <- princomp(dm, cor = TRUE) > layout(matrix(1:6, nr = 2), height = c(2, 1)) > plot(cs <- hclust(dm, method = "single"), main = "Single") > abline(h = 3.8, col = "lightgrey") > xlim <- range(body_pc$scores[,1]) > plot(body_pc$scores[,1:2], type = "n", xlim = xlim, ylim = xlim, + xlab = "PC1", ylab = "PC2") > lab <- cutree(cs, h = 3.8) > text(body_pc$scores[,1:2], labels = lab, cex=0.6) > plot(cc <- hclust(dm, method = "complete"), main = "Complete") > abline(h = 10, col = "lightgrey") > plot(body_pc$scores[,1:2], type = "n", xlim = xlim, ylim = xlim, + xlab = "PC1", ylab = "PC2") > lab <- cutree(cc, h = 10) > text(body_pc$scores[,1:2], labels = lab, cex=0.6) > plot(ca <- hclust(dm, method = "average"), main = "Average") > abline(h = 7.8, col = "lightgrey") > plot(body_pc$scores[,1:2], type = "n", xlim = xlim, ylim = xlim, + xlab = "PC1", ylab = "PC2") > lab <- cutree(ca, h = 7.8) > text(body_pc$scores[,1:2], labels = lab, cex=0.6) > ################################################### > ### code chunk number 9: ch:CA:measure:pca (eval = FALSE) > ################################################### > ## body_pc <- princomp(dm, cor = TRUE) > ## xlim <- range(body_pc$scores[,1]) > ## plot(body_pc$scores[,1:2], type = "n", > ## xlim = xlim, ylim = xlim) > ## lab <- cutree(cs, h = 3.8) > ## text(body_pc$scores[,1:2], labels = lab, cex = 0.6) > > > ################################################### > ### code chunk number 10: ch:CA:jet:tab > ################################################### > jet <- + structure(list(V1 = c(82L, 89L, 101L, 107L, 115L, 122L, 127L, + 137L, 147L, 166L, 174L, 175L, 177L, 184L, 187L, 189L, 194L, 197L, + 201L, 204L, 255L, 328L), V2 = c(1.468, 1.605, 2.168, 2.054, 2.467, + 1.294, 2.183, 2.426, 2.607, 4.567, 4.588, 3.618, 5.855, 2.898, + 3.88, 0.455, 8.088, 6.502, 6.081, 7.105, 8.548, 6.321), V3 = c(3.3, + 3.64, 4.87, 4.72, 4.11, 3.75, 3.97, 4.65, 3.84, 4.92, 3.82, 4.32, + 4.53, 4.48, 5.39, 4.99, 4.5, 5.2, 5.65, 5.4, 4.2, 6.45), V4 = c(0.166, + 0.154, 0.177, 0.275, 0.298, 0.15, 0, 0.117, 0.155, 0.138, 0.249, + 0.143, 0.172, 0.178, 0.101, 0.008, 0.251, 0.366, 0.106, 0.089, + 0.222, 0.187), V5 = c(0.1, 0.1, 2.9, 1.1, 1, 0.9, 2.4, 1.8, 2.3, + 3.2, 3.5, 2.8, 2.5, 3, 3, 2.64, 2.7, 2.9, 2.9, 3.2, 2.9, 2), + V6 = c(0L, 0L, 1L, 0L, 1L, 0L, 1L, 0L, 0L, 1L, 0L, 0L, 1L, + 0L, 1L, 0L, 1L, 1L, 1L, 1L, 0L, 1L)), .Names = c("V1", "V2", + "V3", "V4", "V5", "V6"), class = "data.frame", row.names = c(NA, + -22L)) > colnames(jet) <- c("FFD", "SPR", "RGF", "PLF", "SLF", "CAR") > rownames(jet) <- c("FH-1", "FJ-1", "F-86A", "F9F-2", "F-94A", "F3D-1", "F-89A", + "XF10F-1", "F9F-6", "F-100A", "F4D-1", "F11F-1", + "F-101A", "F3H-2", "F-102A", "F-8A", "F-104B", + "F-105B", "YF-107A", "F-106A", "F-4B", "F-111A") > jet$CAR <- factor(jet$CAR, labels = c("no", "yes")) > toLatex(HSAURtable(jet), pcol = 1, + caption = "Jet fighters data.", + label = "ch:CA:jet:tab") \index{jet data@\Robject{jet} data} \begin{center} \begin{longtable}{ rrrrrr } \caption{\Robject{jet} data. Jet fighters data. \label{ch:CA:jet:tab}} \\ \hline \Robject{FFD} & \Robject{SPR} & \Robject{RGF} & \Robject{PLF} & \Robject{SLF} & \Robject{CAR} \\ \hline \endfirsthead \caption[]{\Robject{jet} data (continued).} \\ \hline \Robject{FFD} & \Robject{SPR} & \Robject{RGF} & \Robject{PLF} & \Robject{SLF} & \Robject{CAR} \\ \hline \endhead 82 & 1.468 & 3.30 & 0.166 & 0.10 & no \\ 89 & 1.605 & 3.64 & 0.154 & 0.10 & no \\ 101 & 2.168 & 4.87 & 0.177 & 2.90 & yes \\ 107 & 2.054 & 4.72 & 0.275 & 1.10 & no \\ 115 & 2.467 & 4.11 & 0.298 & 1.00 & yes \\ 122 & 1.294 & 3.75 & 0.150 & 0.90 & no \\ 127 & 2.183 & 3.97 & 0.000 & 2.40 & yes \\ 137 & 2.426 & 4.65 & 0.117 & 1.80 & no \\ 147 & 2.607 & 3.84 & 0.155 & 2.30 & no \\ 166 & 4.567 & 4.92 & 0.138 & 3.20 & yes \\ 174 & 4.588 & 3.82 & 0.249 & 3.50 & no \\ 175 & 3.618 & 4.32 & 0.143 & 2.80 & no \\ 177 & 5.855 & 4.53 & 0.172 & 2.50 & yes \\ 184 & 2.898 & 4.48 & 0.178 & 3.00 & no \\ 187 & 3.880 & 5.39 & 0.101 & 3.00 & yes \\ 189 & 0.455 & 4.99 & 0.008 & 2.64 & no \\ 194 & 8.088 & 4.50 & 0.251 & 2.70 & yes \\ 197 & 6.502 & 5.20 & 0.366 & 2.90 & yes \\ 201 & 6.081 & 5.65 & 0.106 & 2.90 & yes \\ 204 & 7.105 & 5.40 & 0.089 & 3.20 & yes \\ 255 & 8.548 & 4.20 & 0.222 & 2.90 & no \\ 328 & 6.321 & 6.45 & 0.187 & 2.00 & yes \\ \hline \end{longtable} \end{center} > ################################################### > ### code chunk number 11: ch:CA:jet:hclust > ################################################### > X <- scale(jet[, c("SPR", "RGF", "PLF", "SLF")], + center = FALSE, scale = TRUE) > dj <- dist(X) > plot(cc <- hclust(dj), main = "Jets clustering") > cc Call: hclust(d = dj) Cluster method : complete Distance : euclidean Number of objects: 22 > ################################################### > ### code chunk number 12: ch:CA:jet:hclust:plot > ################################################### > X <- scale(jet[, c("SPR", "RGF", "PLF", "SLF")], + center = FALSE, scale = TRUE) > dj <- dist(X) > plot(cc <- hclust(dj), main = "Jets clustering") > cc Call: hclust(d = dj) Cluster method : complete Distance : euclidean Number of objects: 22 > ################################################### > ### code chunk number 13: ch:CA:jet:PCA > ################################################### > pr <- prcomp(dj)$x[, 1:2] > plot(pr, pch = (1:2)[cutree(cc, k = 2)], + col = c("black", "darkgrey")[jet$CAR], + xlim = range(pr) * c(1, 1.5)) > legend("topright", col = c("black", "black", + "darkgrey", "darkgrey"), + legend = c("1 / no", "2 / no", "1 / yes", "2 / yes"), + pch = c(1:2, 1:2), title = "Cluster / CAR", bty = "n") > ################################################### > ### code chunk number 14: ch:CA:crime:tab > ################################################### > `crime` <- + structure(c(2, 2.2, 2, 3.6, 3.5, 4.6, 10.7, 5.2, 5.5, 5.5, 6, + 8.9, 11.3, 3.1, 2.5, 1.8, 9.2, 1, 4, 3.1, 4.4, 4.9, 9, 31, 7.1, + 5.9, 8.1, 8.6, 11.2, 11.7, 6.7, 10.4, 10.1, 11.2, 8.1, 12.8, + 8.1, 13.5, 2.9, 3.2, 5.3, 7, 11.5, 9.3, 3.2, 12.6, 5, 6.6, 11.3, + 8.6, 4.8, 14.8, 21.5, 21.8, 29.7, 21.4, 23.8, 30.5, 33.2, 25.1, + 38.6, 25.9, 32.4, 67.4, 20.1, 31.8, 12.5, 29.2, 11.6, 17.7, 24.6, + 32.9, 56.9, 43.6, 52.4, 26.5, 18.9, 26.4, 41.3, 43.9, 52.7, 23.1, + 47, 28.4, 25.8, 28.9, 40.1, 36.4, 51.6, 17.3, 20, 21.9, 42.3, + 46.9, 43, 25.3, 64.9, 53.4, 51.1, 44.9, 72.7, 31, 28, 24, 22, + 193, 119, 192, 514, 269, 152, 142, 90, 325, 301, 73, 102, 42, + 170, 7, 16, 51, 80, 124, 304, 754, 106, 41, 88, 99, 214, 367, + 83, 208, 112, 65, 80, 224, 107, 240, 20, 21, 22, 145, 130, 169, + 59, 287, 135, 206, 343, 88, 106, 102, 92, 103, 331, 192, 205, + 431, 265, 176, 235, 186, 434, 424, 162, 148, 179, 370, 32, 87, + 184, 252, 241, 476, 668, 167, 99, 354, 525, 319, 605, 222, 274, + 408, 172, 278, 482, 285, 354, 118, 178, 243, 329, 538, 437, 180, + 354, 244, 286, 521, 401, 103, 803, 755, 949, 1071, 1294, 1198, + 1221, 1071, 735, 988, 887, 1180, 1509, 783, 1004, 956, 1136, + 385, 554, 748, 1188, 1042, 1296, 1728, 813, 625, 1225, 1340, + 1453, 2221, 824, 1325, 1159, 1076, 1030, 1461, 1787, 2049, 783, + 1003, 817, 1792, 1845, 1908, 915, 1604, 1861, 1967, 1696, 1162, + 1339, 2347, 2208, 2697, 2189, 2568, 2758, 2924, 2822, 1654, 2574, + 2333, 2938, 3378, 2802, 2785, 2801, 2500, 2049, 1939, 2677, 3008, + 3090, 2978, 4131, 2522, 1358, 2423, 2846, 2984, 4373, 1740, 2126, + 2304, 1845, 2305, 3417, 3142, 3987, 3314, 2800, 3078, 4231, 3712, + 4337, 4074, 3489, 4267, 4163, 3384, 3910, 3759, 164, 228, 181, + 906, 705, 447, 637, 776, 354, 376, 328, 628, 800, 254, 288, 158, + 439, 120, 99, 168, 258, 272, 545, 975, 219, 169, 208, 277, 430, + 598, 193, 544, 267, 150, 195, 442, 649, 714, 215, 181, 169, 486, + 343, 419, 223, 478, 315, 402, 762, 604, 328), .Dim = c(51L, 7L + ), .Dimnames = list(c("ME", "NH", "VT", "MA", "RI", "CT", "NY", + "NJ", "PA", "OH", "IN", "IL", "MI", "WI", "MN", "IA", "MO", "ND", + "SD", "NE", "KS", "DE", "MD", "DC", "VA", "WV", "NC", "SC", "GA", + "FL", "KY", "TN", "AL", "MS", "AR", "LA", "OK", "TX", "MT", "ID", + "WY", "CO", "NM", "AZ", "UT", "NV", "WA", "OR", "CA", "AK", "HI" + ), c("Murder", "Rape", "Robbery", "Assault", "Burglary", "Theft", + "Vehicle"))) > crime <- as.data.frame(crime) > toLatex(HSAURtable(crime), pcol = 1, rownames = TRUE, + caption = "Crime data.", + label = "ch:CA:crime:tab") \index{crime data@\Robject{crime} data} \begin{center} \begin{longtable}{l rrrrrrr } \caption{\Robject{crime} data. Crime data. \label{ch:CA:crime:tab}} \\ \hline & \Robject{Murder} & \Robject{Rape} & \Robject{Robbery} & \Robject{Assault} & \Robject{Burglary} & \Robject{Theft} & \Robject{Vehicle} \\ \hline \endfirsthead \caption[]{\Robject{crime} data (continued).} \\ \hline & \Robject{Murder} & \Robject{Rape} & \Robject{Robbery} & \Robject{Assault} & \Robject{Burglary} & \Robject{Theft} & \Robject{Vehicle} \\ \hline \endhead ME & 2.0 & 14.8 & 28 & 102 & 803 & 2347 & 164 \\ NH & 2.2 & 21.5 & 24 & 92 & 755 & 2208 & 228 \\ VT & 2.0 & 21.8 & 22 & 103 & 949 & 2697 & 181 \\ MA & 3.6 & 29.7 & 193 & 331 & 1071 & 2189 & 906 \\ RI & 3.5 & 21.4 & 119 & 192 & 1294 & 2568 & 705 \\ CT & 4.6 & 23.8 & 192 & 205 & 1198 & 2758 & 447 \\ NY & 10.7 & 30.5 & 514 & 431 & 1221 & 2924 & 637 \\ NJ & 5.2 & 33.2 & 269 & 265 & 1071 & 2822 & 776 \\ PA & 5.5 & 25.1 & 152 & 176 & 735 & 1654 & 354 \\ OH & 5.5 & 38.6 & 142 & 235 & 988 & 2574 & 376 \\ IN & 6.0 & 25.9 & 90 & 186 & 887 & 2333 & 328 \\ IL & 8.9 & 32.4 & 325 & 434 & 1180 & 2938 & 628 \\ MI & 11.3 & 67.4 & 301 & 424 & 1509 & 3378 & 800 \\ WI & 3.1 & 20.1 & 73 & 162 & 783 & 2802 & 254 \\ MN & 2.5 & 31.8 & 102 & 148 & 1004 & 2785 & 288 \\ IA & 1.8 & 12.5 & 42 & 179 & 956 & 2801 & 158 \\ MO & 9.2 & 29.2 & 170 & 370 & 1136 & 2500 & 439 \\ ND & 1.0 & 11.6 & 7 & 32 & 385 & 2049 & 120 \\ SD & 4.0 & 17.7 & 16 & 87 & 554 & 1939 & 99 \\ NE & 3.1 & 24.6 & 51 & 184 & 748 & 2677 & 168 \\ KS & 4.4 & 32.9 & 80 & 252 & 1188 & 3008 & 258 \\ DE & 4.9 & 56.9 & 124 & 241 & 1042 & 3090 & 272 \\ MD & 9.0 & 43.6 & 304 & 476 & 1296 & 2978 & 545 \\ DC & 31.0 & 52.4 & 754 & 668 & 1728 & 4131 & 975 \\ VA & 7.1 & 26.5 & 106 & 167 & 813 & 2522 & 219 \\ WV & 5.9 & 18.9 & 41 & 99 & 625 & 1358 & 169 \\ NC & 8.1 & 26.4 & 88 & 354 & 1225 & 2423 & 208 \\ SC & 8.6 & 41.3 & 99 & 525 & 1340 & 2846 & 277 \\ GA & 11.2 & 43.9 & 214 & 319 & 1453 & 2984 & 430 \\ FL & 11.7 & 52.7 & 367 & 605 & 2221 & 4373 & 598 \\ KY & 6.7 & 23.1 & 83 & 222 & 824 & 1740 & 193 \\ TN & 10.4 & 47.0 & 208 & 274 & 1325 & 2126 & 544 \\ AL & 10.1 & 28.4 & 112 & 408 & 1159 & 2304 & 267 \\ MS & 11.2 & 25.8 & 65 & 172 & 1076 & 1845 & 150 \\ AR & 8.1 & 28.9 & 80 & 278 & 1030 & 2305 & 195 \\ LA & 12.8 & 40.1 & 224 & 482 & 1461 & 3417 & 442 \\ OK & 8.1 & 36.4 & 107 & 285 & 1787 & 3142 & 649 \\ TX & 13.5 & 51.6 & 240 & 354 & 2049 & 3987 & 714 \\ MT & 2.9 & 17.3 & 20 & 118 & 783 & 3314 & 215 \\ ID & 3.2 & 20.0 & 21 & 178 & 1003 & 2800 & 181 \\ WY & 5.3 & 21.9 & 22 & 243 & 817 & 3078 & 169 \\ CO & 7.0 & 42.3 & 145 & 329 & 1792 & 4231 & 486 \\ NM & 11.5 & 46.9 & 130 & 538 & 1845 & 3712 & 343 \\ AZ & 9.3 & 43.0 & 169 & 437 & 1908 & 4337 & 419 \\ UT & 3.2 & 25.3 & 59 & 180 & 915 & 4074 & 223 \\ NV & 12.6 & 64.9 & 287 & 354 & 1604 & 3489 & 478 \\ WA & 5.0 & 53.4 & 135 & 244 & 1861 & 4267 & 315 \\ OR & 6.6 & 51.1 & 206 & 286 & 1967 & 4163 & 402 \\ CA & 11.3 & 44.9 & 343 & 521 & 1696 & 3384 & 762 \\ AK & 8.6 & 72.7 & 88 & 401 & 1162 & 3910 & 604 \\ HI & 4.8 & 31.0 & 106 & 103 & 1339 & 3759 & 328 \\ \hline \end{longtable} \end{center} > ################################################### > ### code chunk number 15: ch:CA:crime:plot > ################################################### > plot(crime, pch = ".", cex = 1.5) > ################################################### > ### code chunk number 16: ch:CA:crime:outlier > ################################################### > subset(crime, Murder > 15) Murder Rape Robbery Assault Burglary Theft Vehicle DC 31 52.4 754 668 1728 4131 975 > ################################################### > ### code chunk number 17: ch:CA:crime:plot2 > ################################################### > plot(crime, pch = c(".", "+")[(rownames(crime) == "DC") + 1], cex = 1.5) > ################################################### > ### code chunk number 18: ch:CA:crime:var > ################################################### > sapply(crime, var) Murder Rape Robbery Assault Burglary Theft 23.20215 212.31228 18993.37020 22004.31294 177912.83373 582812.83843 Vehicle 50007.37490 > ################################################### > ### code chunk number 19: ch:CA:crime:stand > ################################################### > rge <- sapply(crime, function(x) diff(range(x))) > crime_s <- sweep(crime, 2, rge, FUN = "/") > sapply(crime_s, var) Murder Rape Robbery Assault Burglary Theft Vehicle 0.02578017 0.05687124 0.03403775 0.05439933 0.05277909 0.06411424 0.06516672 > ################################################### > ### code chunk number 20: ch:CA:crime:wss > ################################################### > n <- nrow(crime_s) > wss <- rep(0, 6) > wss[1] <- (n - 1) * sum(sapply(crime_s, var)) > for (i in 2:6) + wss[i] <- sum(kmeans(crime_s, + centers = i)$withinss) > plot(1:6, wss, type = "b", xlab = "Number of groups", + ylab = "Within groups sum of squares") > ################################################### > ### code chunk number 21: ch:CA:crimes:k2 > ################################################### > kmeans(crime_s, centers = 2)$centers * rge Murder Rape Robbery Assault Burglary Theft Vehicle 1 10.359091 567.6133 628.0876 562.400515 52.25841 726.112 1991.5395 2 9.965621 259.7625 311.3398 8.893949 378.99966 1560.437 253.6552 > ################################################### > ### code chunk number 22: ch:CA:crime:PCA > ################################################### > crime_pca <- prcomp(crime_s) > plot(crime_pca$x[, 1:2], + pch = kmeans(crime_s, centers = 2)$cluster) > ################################################### > ### code chunk number 23: ca:CA:pottery:dist > ################################################### > pottery_dist <- dist(pots <- scale(pottery[, colnames(pottery) != "kiln"], + center = FALSE)) > library("lattice") > levelplot(as.matrix(pottery_dist), xlab = "Pot Number", + ylab = "Pot Number") > ################################################### > ### code chunk number 24: ch:CA:pottery:distplot > ################################################### > trellis.par.set(standard.theme(color = FALSE)) > plot(levelplot(as.matrix(pottery_dist), xlab = "Pot Number", ylab = "Pot Number", + scales = list(x = list(draw = FALSE), y = list(draw = FALSE)))) > ################################################### > ### code chunk number 25: ch:CA:pottery:wss > ################################################### > n <- nrow(pots) > wss <- rep(0, 6) > wss[1] <- (n - 1) * sum(sapply(pots, var)) > for (i in 2:6) + wss[i] <- sum(kmeans(pots, + centers = i)$withinss) > plot(1:6, wss, type = "b", xlab = "Number of groups", + ylab = "Within groups sum of squares") > ################################################### > ### code chunk number 26: ch:CA:pots:PCA > ################################################### > pots_pca <- prcomp(pots) > plot(pots_pca$x[, 1:2], + pch = kmeans(pots, centers = 3)$cluster) > ################################################### > ### code chunk number 27: ch:CA:pottery:cluster > ################################################### > set.seed(29) > pottery_cluster <- kmeans(pots, centers = 3)$cluster > xtabs(~ pottery_cluster + kiln, data = pottery) kiln pottery_cluster 1 2 3 4 5 1 21 0 0 0 0 2 0 0 0 5 5 3 0 12 2 0 0 > ################################################### > ### code chunk number 28: ch:CA:thompson > ################################################### > cnt <- c("Iceland", "Norway", "Sweden", "Finland", "Denmark", "UK", "Eire", + "Germany", "Netherlands", "Belgium", "Switzerland", "France", "Spain", + "Portugal", "Italy", "Greece", "Yugoslavia", "Albania", "Bulgaria", "Romania", + "Hungary", "Czechia", "Slovakia", "Poland", "CIS", "Lithuania", "Latvia", "Estonia") > thomson <- expand.grid(answer = factor(c("no", "yes")), + question = factor(paste("Q", 1:6, sep = "")), + country = factor(cnt, levels = cnt)) > thomson$Freq <- c( + 0, 5, 0, 5, 0, 4, 0, 5, 0, 5, 0, 5, + 1, 6, 1, 5, 0, 6, 0, 5, 0, 4, 1, 4, + 0, 11, 4, 7, 0, 7, 0, 11, 5, 5, 3, 6, + 0, 6, 2, 4, 0, 6, 0, 6, 1, 5, 2, 4, + 1, 12, 4, 9, 0, 12, 3, 9, 7, 4, 6, 7, + 7, 12, 2, 16, 0, 20, 1, 19, 9, 10, 0, 17, + 0, 1, 1, 2, 0, 3, 2, 0, 2, 0, 0, 3, + 0, 14, 0, 13, 0, 13, 2, 12, 11, 2, 1, 13, + 0, 8, 0, 8, 0, 8, 1, 7, 2, 5, 1, 7, + 2, 0, 0, 2, 0, 2, 1, 1, 2, 0, 0, 2, + 0, 5, 0, 5, 0, 4, 2, 2, 5, 0, 0, 4, + 7, 3, 1, 7, 3, 5, 8, 2, 10, 0, 1, 7, + 11, 1, 0, 12, 2, 8, 5, 6, 11, 0, 0, 11, + 5, 1, 0, 6, 2, 4, 3, 3, 6, 0, 0, 6, + 8, 7, 0, 15, 1, 13, 9, 6, 13, 2, 0, 15, + 7, 1, 0, 8, 3, 5, 7, 1, 8, 0, 0, 7, + 11, 4, 0, 15, 7, 8, 11, 4, 15, 0, 0, 14, + 3, 2, 2, 3, 3, 2, 3, 2, 3, 3, 2, 3, + 3, 0, 0, 3, 2, 1, 3, 0, 3, 0, 0, 3, + 7, 0, 0, 6, 6, 1, 6, 1, 6, 1, 0, 7, + 4, 1, 0, 5, 1, 4, 5, 0, 5, 0, 0, 5, + 18, 2, 0, 20, 17, 3, 20, 0, 20, 0, 0, 20, + 13, 0, 1, 14, 14, 0, 16, 0, 13, 0, 15, 0, + 18, 0, 0, 19, 13, 5, 17, 2, 17, 0, 0, 19, + 7, 0, 1, 6, 5, 2, 7, 0, 7, 0, 1, 6, + 8, 0, 0, 8, 8, 0, 8, 0, 8, 0, 0, 8, + 5, 0, 0, 5, 5, 0, 5, 0, 5, 0, 0, 5, + 2, 2, 0, 3, 0, 3, 3, 0, 3, 0, 0, 3) > ttab <- xtabs(Freq ~ country + answer + question, data = thomson) > thomsonprop <- prop.table(ttab, c(1,3))[,"yes",] > plot(1:(22 * 6), rep(-1, 22 * 6), + ylim = c(-nlevels(thomson$country), -1), type = "n", + axes = FALSE, xlab = "", ylab = "") > for (q in 1:6) { + tmp <- ttab[,,q] + xstart <- (q - 1) * 22 + 1 + y <- -rep(1:nrow(tmp), rowSums(tmp)) + x <- xstart + unlist(sapply(rowSums(tmp), function(i) 1:i)) + pch <- unlist(apply(tmp, 1, function(x) c(rep(19, x[2]), rep(1, x[1])))) + points(x, y, pch = pch) + } > axis(2, at = -(1:nlevels(thomson$country)), labels = levels(thomson$country), + las = 2, tick = FALSE, line = 0) > mtext(text = paste("Question", 1:6), 3, at = 22 * (0:5), adj = 0) > ################################################### > ### code chunk number 29: ch:CA:thompsonMC > ################################################### > library("mclust") Package 'mclust' version 6.1.1 Type 'citation("mclust")' for citing this R package in publications. Attaching package: 'mclust' The following object is masked from 'package:mvtnorm': dmvnorm > (mc <- Mclust(thomsonprop)) 'Mclust' model object: (VEV,2) Available components: [1] "call" "data" "modelName" "n" [5] "d" "G" "BIC" "loglik" [9] "df" "bic" "icl" "hypvol" [13] "parameters" "z" "classification" "uncertainty" > ################################################### > ### code chunk number 30: ch:CA:thompsonMC:plot > ################################################### > plot(mc, thomsonprop, what = "BIC", col = "black") > ################################################### > ### code chunk number 31: ch:CA:thompsonMC > ################################################### > cl <- mc$classification > nm <- unlist(sapply(1:3, function(i) names(cl[cl == i]))) > ttab <- ttab[nm,,] > plot(1:(22 * 6), rep(-1, 22 * 6), + ylim = c(-nlevels(thomson$country), -1), type = "n", + axes = FALSE, xlab = "", ylab = "") > for (q in 1:6) { + tmp <- ttab[,,q] + xstart <- (q - 1) * 22 + 1 + y <- -rep(1:nrow(tmp), rowSums(tmp)) + x <- xstart + unlist(sapply(rowSums(tmp), function(i) 1:i)) + pch <- unlist(apply(tmp, 1, function(x) c(rep(19, x[2]), rep(1, x[1])))) + points(x, y, pch = pch) + } > axis(2, at = -(1:nlevels(thomson$country)), labels = dimnames(ttab)[[1]], + las = 2, tick = FALSE, line = 0) > mtext(text = paste("Question", 1:6), 3, at = 22 * (0:5), adj = 0) > abline(h = -cumsum(table(cl))[-3] - 0.5, col = "grey") > text(-c(0.75, 0.75, 0.75), -cumsum(table(cl)) + table(cl)/2, + label = paste("Cluster", 1:3), srt = 90, pos = 1) > ################################################### > ### code chunk number 32: ch:CA:neighbor > ################################################### > library("flexclust") Loading required package: grid Loading required package: modeltools Loading required package: stats4 > library("mvtnorm") > set.seed(290875) > x <- rbind(rmvnorm(n = 20, mean = c(0, 0), + sigma = diag(2)), + rmvnorm(n = 20, mean = c(3, 3), + sigma = 0.5 * diag(2)), + rmvnorm(n = 20, mean = c(7, 6), + sigma = 0.5 * (diag(2) + 0.25))) > k <- cclust(x, k = 5, save.data = TRUE) > plot(k, hull = FALSE, col = rep("black", 5), xlab = "x", ylab = "y") > ################################################### > ### code chunk number 33: ch:CA:pottery:neighbor > ################################################### > k <- cclust(pots, k = 3, save.data = TRUE) > plot(k, project = prcomp(pots), hull = FALSE, col = rep("black", 3), + xlab = "PC1", ylab = "PC2") > ################################################### > ### code chunk number 34: ch:CA:art:stripes1 > ################################################### > set.seed(912345654) > x <- rbind(matrix(rnorm(100, sd = 0.5), ncol= 2 ), + matrix(rnorm(100, mean =4, sd = 0.5), ncol = 2), + matrix(rnorm(100, mean =7, sd = 0.5), ncol = 2), + matrix(rnorm(100, mean =-1.0, sd = 0.7), ncol = 2), + matrix(rnorm(100, mean =-4.0, sd = 1.0), ncol = 2)) > c5 <- cclust(x, 5, save.data = TRUE) > stripes(c5, type = "second", col = 1) > ################################################### > ### code chunk number 35: ch:CA:art:stripes2 > ################################################### > set.seed(912345654) > x <- rbind(matrix(rnorm(100, sd = 2.5), ncol = 2), + matrix(rnorm(100, mean = 3, sd = 0.5), ncol = 2), + matrix(rnorm(100, mean = 5, sd = 0.5), ncol = 2), + matrix(rnorm(100, mean = -1.0, sd = 1.5), ncol = 2), + matrix(rnorm(100, mean = -4.0, sd = 2.0), ncol = 2)) > c5 <- cclust(x, 5, save.data = TRUE) > stripes(c5, type = "second", col = 1) > ################################################### > ### code chunk number 36: ch:CA:pottery:stripes > ################################################### > set.seed(15) > c5 <- cclust(pots, k = 3, save.data = TRUE) > stripes(c5, type = "second", col = "black") demo(Ch-EFA) ---- ~~~~~~ > ### R code from vignette source 'Ch-EFA.Rnw' > > ################################################### > ### code chunk number 1: setup > ################################################### > library("MVA") > set.seed(280875) > library("lattice") > lattice.options(default.theme = + function() + standard.theme("pdf", color = FALSE)) > if (file.exists("deparse.R")) { + if (!file.exists("figs")) dir.create("figs") + source("deparse.R") + options(prompt = "R> ", continue = "+ ", width = 64, + digits = 4, show.signif.stars = FALSE, useFancyQuotes = FALSE) + + options(SweaveHooks = list(onefig = function() {par(mfrow = c(1,1))}, + twofig = function() {par(mfrow = c(1,2))}, + figtwo = function() {par(mfrow = c(2,1))}, + threefig = function() {par(mfrow = c(1,3))}, + figthree = function() {par(mfrow = c(3,1))}, + fourfig = function() {par(mfrow = c(2,2))}, + sixfig = function() {par(mfrow = c(3,2))}, + nomar = function() par("mai" = c(0, 0, 0, 0)))) + } > ################################################### > ### code chunk number 2: ch:EFA:data > ################################################### > d <- + c(0.447, + 0.422, 0.619, + 0.435, 0.604, 0.583, + 0.114, 0.068, 0.053, 0.115, + 0.203, 0.146, 0.139, 0.258, 0.349, + 0.091, 0.103, 0.110, 0.122, 0.209, 0.221, + 0.082, 0.063, 0.066, 0.097, 0.321, 0.355, 0.201, + 0.513, 0.445, 0.365, 0.482, 0.186, 0.315, 0.150, 0.154, + 0.304, 0.318, 0.240, 0.368, 0.303, 0.377, 0.163, 0.219, 0.534, + 0.245, 0.203, 0.183, 0.255, 0.272, 0.323, 0.310, 0.288, 0.301, 0.302, + 0.101, 0.088, 0.074, 0.139, 0.279, 0.367, 0.232, 0.320, 0.204, 0.368, 0.340, + 0.245, 0.199, 0.184, 0.293, 0.278, 0.545, 0.232, 0.314, 0.394, 0.467, 0.392, 0.511) > druguse <- diag(13) / 2 > druguse[upper.tri(druguse)] <- d > druguse <- druguse + t(druguse) > rownames(druguse) <- colnames(druguse) <- c("cigarettes", "beer", "wine", "liquor", "cocaine", + "tranquillizers", "drug store medication", "heroin", + "marijuana", "hashish", "inhalants", "hallucinogenics", "amphetamine") > ################################################### > ### code chunk number 3: ch:EFA:life:tab > ################################################### > "life" <- + structure(.Data = list(c(63., 34., 38., 59., 56., 62., 50., 65., 56., 69., 65., 64., 56., 60., 61., 49., 59., 63., 59., 65., 65., 64., + 64., 67., 61., 68., 67., 65., 59., 58., 57.) + , c(51., 29., 30., 42., 38., 44., 39., 44., 46., 47., 48., 50., 44., 44., 45., 40., 42., 44., 44., 48., 48., 63., + 43., 45., 40., 46., 45., 46., 43., 44., 46.) + , c(30., 13., 17., 20., 18., 24., 20., 22., 24., 24., 26., 28., 25., 22., 22., 22., 22., 23., 24., 28., 26., 21., + 21., 23., 21., 23., 23., 24., 23., 24., 28.) + , c(13., 5., 7., 6., 7., 7., 7., 7., 11., 8., 9., 11., 10., 6., 8., 9., 6., 8., 8., 14., 9., 7., 6., 8., 10., 8., + 8., 9., 10., 9., 9.) + , c(67., 38., 38., 64., 62., 69., 55., 72., 63., 75., 68., 66., 61., 65., 65., 51., 61., 67., 63., 68., 67., 68., + 68., 74., 67., 75., 74., 71., 66., 62., 60.) + , c(54., 32., 34., 46., 46., 50., 43., 50., 54., 53., 50., 51., 48., 45., 49., 41., 43., 48., 46., 51., 49., 47., + 47., 51., 46., 52., 51., 51., 49., 47., 49.) + , c(34., 17., 20., 25., 25., 28., 23., 27., 33., 29., 27., 29., 27., 25., 27., 23., 22., 26., 25., 29., 27., 25., + 24., 28., 25., 29., 28., 28., 27., 25., 28.) + , c(15., 6., 7., 8., 10., 14., 8., 9., 19., 10., 10., 11., 12., 9., 10., 8., 7., 9., 8., 13., 10., 9., 8., 10., 11., + 10., 10., 10., 12., 10., 11.) + ) + , class = "data.frame" + , names = c("m0", "m25", "m50", "m75", "w0", "w25", "w50", "w75") + , row.names = c("Algeria", "Cameroon", "Madagascar", "Mauritius", "Reunion", "Seychelles", "South Africa (C)", "South Africa (W)", + "Tunisia", "Canada", "Costa Rica", "Dominican Rep.", "El Salvador", "Greenland", "Grenada", "Guatemala", + "Honduras", "Jamaica", "Mexico", "Nicaragua", "Panama", "Trinidad (62)", "Trinidad (67)", + "United States (66)", "United States (NW66)", "United States (W66)", "United States (67)", "Argentina", + "Chile", "Colombia", "Ecuador") + ) > toLatex(HSAURtable(life), pcol = 1, rownames = TRUE, + caption = "Life expectancies for different countries by age and gender.", + label = "ch:EFA:life:tab") \index{life data@\Robject{life} data} \begin{center} \begin{longtable}{l rrrrrrrr } \caption{\Robject{life} data. Life expectancies for different countries by age and gender. \label{ch:EFA:life:tab}} \\ \hline & \Robject{m0} & \Robject{m25} & \Robject{m50} & \Robject{m75} & \Robject{w0} & \Robject{w25} & \Robject{w50} & \Robject{w75} \\ \hline \endfirsthead \caption[]{\Robject{life} data (continued).} \\ \hline & \Robject{m0} & \Robject{m25} & \Robject{m50} & \Robject{m75} & \Robject{w0} & \Robject{w25} & \Robject{w50} & \Robject{w75} \\ \hline \endhead Algeria & 63 & 51 & 30 & 13 & 67 & 54 & 34 & 15 \\ Cameroon & 34 & 29 & 13 & 5 & 38 & 32 & 17 & 6 \\ Madagascar & 38 & 30 & 17 & 7 & 38 & 34 & 20 & 7 \\ Mauritius & 59 & 42 & 20 & 6 & 64 & 46 & 25 & 8 \\ Reunion & 56 & 38 & 18 & 7 & 62 & 46 & 25 & 10 \\ Seychelles & 62 & 44 & 24 & 7 & 69 & 50 & 28 & 14 \\ South Africa (C) & 50 & 39 & 20 & 7 & 55 & 43 & 23 & 8 \\ South Africa (W) & 65 & 44 & 22 & 7 & 72 & 50 & 27 & 9 \\ Tunisia & 56 & 46 & 24 & 11 & 63 & 54 & 33 & 19 \\ Canada & 69 & 47 & 24 & 8 & 75 & 53 & 29 & 10 \\ Costa Rica & 65 & 48 & 26 & 9 & 68 & 50 & 27 & 10 \\ Dominican Rep. & 64 & 50 & 28 & 11 & 66 & 51 & 29 & 11 \\ El Salvador & 56 & 44 & 25 & 10 & 61 & 48 & 27 & 12 \\ Greenland & 60 & 44 & 22 & 6 & 65 & 45 & 25 & 9 \\ Grenada & 61 & 45 & 22 & 8 & 65 & 49 & 27 & 10 \\ Guatemala & 49 & 40 & 22 & 9 & 51 & 41 & 23 & 8 \\ Honduras & 59 & 42 & 22 & 6 & 61 & 43 & 22 & 7 \\ Jamaica & 63 & 44 & 23 & 8 & 67 & 48 & 26 & 9 \\ Mexico & 59 & 44 & 24 & 8 & 63 & 46 & 25 & 8 \\ Nicaragua & 65 & 48 & 28 & 14 & 68 & 51 & 29 & 13 \\ Panama & 65 & 48 & 26 & 9 & 67 & 49 & 27 & 10 \\ Trinidad (62) & 64 & 63 & 21 & 7 & 68 & 47 & 25 & 9 \\ Trinidad (67) & 64 & 43 & 21 & 6 & 68 & 47 & 24 & 8 \\ United States (66) & 67 & 45 & 23 & 8 & 74 & 51 & 28 & 10 \\ United States (NW66) & 61 & 40 & 21 & 10 & 67 & 46 & 25 & 11 \\ United States (W66) & 68 & 46 & 23 & 8 & 75 & 52 & 29 & 10 \\ United States (67) & 67 & 45 & 23 & 8 & 74 & 51 & 28 & 10 \\ Argentina & 65 & 46 & 24 & 9 & 71 & 51 & 28 & 10 \\ Chile & 59 & 43 & 23 & 10 & 66 & 49 & 27 & 12 \\ Colombia & 58 & 44 & 24 & 9 & 62 & 47 & 25 & 10 \\ Ecuador & 57 & 46 & 28 & 9 & 60 & 49 & 28 & 11 \\ \hline \end{longtable} \end{center} > ################################################### > ### code chunk number 4: ch:EFA:life > ################################################### > sapply(1:3, function(f) + factanal(life, factors = f, method ="mle")$PVAL) objective objective objective 1.879555e-24 1.911514e-05 4.578204e-01 > ################################################### > ### code chunk number 5: ch:EFA:life3 > ################################################### > factanal(life, factors = 3, method ="mle") Call: factanal(x = life, factors = 3, method = "mle") Uniquenesses: m0 m25 m50 m75 w0 w25 w50 w75 0.005 0.362 0.066 0.288 0.005 0.011 0.020 0.146 Loadings: Factor1 Factor2 Factor3 m0 0.964 0.122 0.226 m25 0.646 0.169 0.438 m50 0.430 0.354 0.790 m75 0.525 0.656 w0 0.970 0.217 w25 0.764 0.556 0.310 w50 0.536 0.729 0.401 w75 0.156 0.867 0.280 Factor1 Factor2 Factor3 SS loadings 3.375 2.082 1.640 Proportion Var 0.422 0.260 0.205 Cumulative Var 0.422 0.682 0.887 Test of the hypothesis that 3 factors are sufficient. The chi square statistic is 6.73 on 7 degrees of freedom. The p-value is 0.458 > ################################################### > ### code chunk number 6: ch:EFA:life:scores > ################################################### > (scores <- factanal(life, factors = 3, method = "mle", + scores = "regression")$scores) Factor1 Factor2 Factor3 Algeria -0.258062561 1.90095771 1.91581631 Cameroon -2.782495791 -0.72340014 -1.84772224 Madagascar -2.806428187 -0.81158820 -0.01210318 Mauritius 0.141004934 -0.29028454 -0.85862443 Reunion -0.196352142 0.47429917 -1.55046466 Seychelles 0.367371307 0.82902375 -0.55214085 South Africa (C) -1.028567629 -0.08065792 -0.65421971 South Africa (W) 0.946193522 0.06400408 -0.91995289 Tunisia -0.862493550 3.59177195 -0.36442148 Canada 1.245304248 0.29564122 -0.27342781 Costa Rica 0.508736247 -0.50500435 1.01328707 Dominican Rep. 0.106044085 0.01111171 1.83871599 El Salvador -0.608155779 0.65100820 0.48836431 Greenland 0.235114220 -0.69123901 -0.38558654 Grenada 0.132008172 0.25241049 -0.15220645 Guatemala -1.450336359 -0.67765804 0.65911906 Honduras 0.043253249 -1.85175707 0.30633182 Jamaica 0.462124701 -0.51918493 0.08032855 Mexico -0.052332675 -0.72020002 0.44417800 Nicaragua 0.268974443 0.08407227 1.70568388 Panama 0.442333434 -0.73778272 1.25218728 Trinidad (62) 0.711367053 -0.95989475 -0.21545329 Trinidad (67) 0.787286051 -1.10729029 -0.51958264 United States (66) 1.128331259 0.16389896 -0.68177046 United States (NW66) 0.400058903 -0.36230253 -0.74299137 United States (W66) 1.214345385 0.40877239 -0.69225320 United States (67) 1.128331259 0.16389896 -0.68177046 Argentina 0.731344988 0.24811968 -0.12817725 Chile 0.009751528 0.75222637 -0.49198911 Colombia -0.240602517 -0.29543613 0.42919600 Ecuador -0.723451797 0.44246371 1.59164974 > ################################################### > ### code chunk number 7: ch:EFA:life:3d > ################################################### > cex <- 0.8 > plot(scores[,1], scores[,2], type = "n", xlab = "Factor 1", ylab = "Factor 2") > text(scores[,1], scores[,2], abbreviate(rownames(life), 5), cex = cex) > plot(scores[,1], scores[,3], type = "n", xlab = "Factor 1", ylab = "Factor 3") > text(scores[,1], scores[,3], abbreviate(rownames(life), 5), cex = cex) > plot(scores[,2], scores[,3], type = "n", xlab = "Factor 2", ylab = "Factor 3") > text(scores[,2], scores[,3], abbreviate(rownames(life), 5), cex = cex) > ################################################### > ### code chunk number 8: ch:EFA:druguse:plot > ################################################### > ord <- order.dendrogram(as.dendrogram(hclust(dist(druguse)))) > panel.corrgram <- + function(x, y, z, subscripts, at, + level = 0.9, label = FALSE, ...) + { + require("ellipse", quietly = TRUE) + x <- as.numeric(x)[subscripts] + y <- as.numeric(y)[subscripts] + z <- as.numeric(z)[subscripts] + zcol <- level.colors(z, at = at, col.regions = grey.colors, ...) + for (i in seq(along = z)) { + ell <- ellipse(z[i], level = level, npoints = 50, + scale = c(.2, .2), centre = c(x[i], y[i])) + panel.polygon(ell, col = zcol[i], border = zcol[i], ...) + } + if (label) + panel.text(x = x, y = y, lab = 100 * round(z, 2), cex = 0.8, + col = ifelse(z < 0, "white", "black")) + } > print(levelplot(druguse[ord, ord], at = do.breaks(c(-1.01, 1.01), 20), + xlab = NULL, ylab = NULL, colorkey = list(space = "top"), + scales = list(x = list(rot = 90)), + panel = panel.corrgram, label = TRUE)) Attaching package: 'ellipse' The following object is masked from 'package:flexclust': pairs The following object is masked from 'package:graphics': pairs > ################################################### > ### code chunk number 9: ch:EFA:drugs > ################################################### > sapply(1:6, function(nf) + factanal(covmat = druguse, factors = nf, + method = "mle", n.obs = 1634)$PVAL) objective objective objective objective objective objective 0.000000e+00 9.786000e-70 7.363910e-28 1.794578e-11 3.891743e-06 9.752967e-02 > ################################################### > ### code chunk number 10: ch:EFA:drugs > ################################################### > (factanal(covmat = druguse, factors = 6, + method = "mle", n.obs = 1634)) Call: factanal(factors = 6, covmat = druguse, n.obs = 1634, method = "mle") Uniquenesses: cigarettes beer wine 0.563 0.368 0.374 liquor cocaine tranquillizers 0.412 0.681 0.522 drug store medication heroin marijuana 0.785 0.669 0.318 hashish inhalants hallucinogenics 0.005 0.541 0.620 amphetamine 0.005 Loadings: Factor1 Factor2 Factor3 Factor4 Factor5 Factor6 cigarettes 0.494 0.407 0.110 beer 0.776 0.112 wine 0.786 liquor 0.720 0.121 0.103 0.115 0.160 cocaine 0.519 0.132 0.158 tranquillizers 0.130 0.564 0.321 0.105 0.143 drug store medication 0.255 0.372 heroin 0.532 0.101 0.190 marijuana 0.429 0.158 0.152 0.259 0.609 0.110 hashish 0.244 0.276 0.186 0.881 0.194 0.100 inhalants 0.166 0.308 0.150 0.140 0.537 hallucinogenics 0.387 0.335 0.186 0.288 amphetamine 0.151 0.336 0.886 0.145 0.137 0.187 Factor1 Factor2 Factor3 Factor4 Factor5 Factor6 SS loadings 2.301 1.415 1.116 0.964 0.676 0.666 Proportion Var 0.177 0.109 0.086 0.074 0.052 0.051 Cumulative Var 0.177 0.286 0.372 0.446 0.498 0.549 Test of the hypothesis that 6 factors are sufficient. The chi square statistic is 22.41 on 15 degrees of freedom. The p-value is 0.0975 > ################################################### > ### code chunk number 11: ch:EFA:drugdiff > ################################################### > pfun <- function(nf) { + fa <- factanal(covmat = druguse, factors = nf, + method = "mle", n.obs = 1634) + est <- tcrossprod(fa$loadings) + diag(fa$uniquenesses) + ret <- round(druguse - est, 3) + colnames(ret) <- rownames(ret) <- + abbreviate(rownames(ret), 3) + ret + } > pfun(6) cgr ber win lqr ccn trn dsm hrn mrj hsh inh cgr 0.000 -0.001 0.014 -0.018 0.010 0.001 -0.020 -0.004 0.001 0 0.010 ber -0.001 0.000 -0.002 0.004 0.004 -0.011 -0.001 0.007 0.002 0 -0.004 win 0.014 -0.002 0.000 -0.001 -0.001 -0.005 0.008 0.008 -0.004 0 -0.007 lqr -0.018 0.004 -0.001 0.000 -0.008 0.021 -0.006 -0.018 0.003 0 0.012 ccn 0.010 0.004 -0.001 -0.008 0.000 0.000 0.008 0.004 -0.004 0 -0.003 trn 0.001 -0.011 -0.005 0.021 0.000 0.000 0.006 -0.004 -0.004 0 0.002 dsm -0.020 -0.001 0.008 -0.006 0.008 0.006 0.000 -0.015 0.008 0 0.004 hrn -0.004 0.007 0.008 -0.018 0.004 -0.004 -0.015 0.000 0.006 0 -0.002 mrj 0.001 0.002 -0.004 0.003 -0.004 -0.004 0.008 0.006 0.000 0 -0.006 hsh 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0 0.000 inh 0.010 -0.004 -0.007 0.012 -0.003 0.002 0.004 -0.002 -0.006 0 0.000 hll -0.005 0.005 -0.001 -0.005 -0.008 -0.008 -0.002 0.020 0.003 0 -0.002 amp 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0 0.000 hll amp cgr -0.005 0 ber 0.005 0 win -0.001 0 lqr -0.005 0 ccn -0.008 0 trn -0.008 0 dsm -0.002 0 hrn 0.020 0 mrj 0.003 0 hsh 0.000 0 inh -0.002 0 hll 0.000 0 amp 0.000 0 > ################################################### > ### code chunk number 12: ch:opt > ################################################### > op <- options(width = 150, prompt = " R> ") > pfun2 <- pfun > pfun <- function(...) { + x <- pfun2(...) + rownames(x) <- paste(" ", rownames(x)) + x + } > ################################################### > ### code chunk number 13: ch:EFA:drufdiff34 > ################################################### > pfun(3) cgr ber win lqr ccn trn dsm hrn mrj hsh inh hll amp cgr 0.000 -0.001 0.009 -0.013 0.011 0.009 -0.011 -0.004 0.003 -0.027 0.039 -0.017 0.002 ber -0.001 0.000 -0.002 0.002 0.002 -0.014 0.000 0.005 -0.001 0.019 -0.002 0.009 -0.007 win 0.009 -0.002 0.000 0.000 -0.002 -0.004 0.012 0.013 0.001 -0.017 -0.007 0.004 0.002 lqr -0.013 0.002 0.000 0.000 -0.008 0.024 -0.017 -0.020 -0.001 0.014 -0.002 -0.015 0.006 ccn 0.011 0.002 -0.002 -0.008 0.000 0.031 0.038 0.082 -0.002 0.041 0.023 -0.030 -0.075 trn 0.009 -0.014 -0.004 0.024 0.031 0.000 -0.021 0.026 -0.002 -0.016 -0.038 -0.058 0.044 dsm -0.011 0.000 0.012 -0.017 0.038 -0.021 0.000 0.021 0.007 -0.040 0.113 0.000 -0.038 hrn -0.004 0.005 0.013 -0.020 0.082 0.026 0.021 0.000 0.006 -0.035 0.031 -0.005 -0.049 mrj 0.003 -0.001 0.001 -0.001 -0.002 -0.002 0.007 0.006 0.000 0.001 0.003 -0.002 -0.002 hsh -0.027 0.019 -0.017 0.014 0.041 -0.016 -0.040 -0.035 0.001 0.000 -0.035 0.034 0.010 inh 0.039 -0.002 -0.007 -0.002 0.023 -0.038 0.113 0.031 0.003 -0.035 0.000 0.007 -0.015 hll -0.017 0.009 0.004 -0.015 -0.030 -0.058 0.000 -0.005 -0.002 0.034 0.007 0.000 0.041 amp 0.002 -0.007 0.002 0.006 -0.075 0.044 -0.038 -0.049 -0.002 0.010 -0.015 0.041 0.000 > pfun(4) cgr ber win lqr ccn trn dsm hrn mrj hsh inh hll amp cgr 0.000 -0.001 0.008 -0.012 0.009 0.008 -0.015 -0.007 0.001 -0.023 0.037 -0.020 0.000 ber -0.001 0.000 -0.001 0.001 0.000 -0.016 -0.002 0.003 -0.001 0.018 -0.005 0.006 0.000 win 0.008 -0.001 0.000 0.000 -0.001 -0.005 0.012 0.014 0.001 -0.020 -0.008 0.001 0.000 lqr -0.012 0.001 0.000 0.000 -0.004 0.029 -0.015 -0.015 -0.001 0.018 0.001 -0.010 -0.001 ccn 0.009 0.000 -0.001 -0.004 0.000 0.024 -0.014 0.007 -0.003 0.035 -0.022 -0.028 0.000 trn 0.008 -0.016 -0.005 0.029 0.024 0.000 -0.020 0.027 -0.001 0.001 -0.032 -0.028 0.001 dsm -0.015 -0.002 0.012 -0.015 -0.014 -0.020 0.000 -0.018 0.003 -0.042 0.090 0.008 0.000 hrn -0.007 0.003 0.014 -0.015 0.007 0.027 -0.018 0.000 0.003 -0.037 -0.001 0.005 0.000 mrj 0.001 -0.001 0.001 -0.001 -0.003 -0.001 0.003 0.003 0.000 0.000 0.001 -0.002 0.000 hsh -0.023 0.018 -0.020 0.018 0.035 0.001 -0.042 -0.037 0.000 0.000 -0.031 0.055 -0.001 inh 0.037 -0.005 -0.008 0.001 -0.022 -0.032 0.090 -0.001 0.001 -0.031 0.000 0.021 0.000 hll -0.020 0.006 0.001 -0.010 -0.028 -0.028 0.008 0.005 -0.002 0.055 0.021 0.000 0.000 amp 0.000 0.000 0.000 -0.001 0.000 0.001 0.000 0.000 0.000 -0.001 0.000 0.000 0.000 > ################################################### > ### code chunk number 14: ch:opt > ################################################### > options(op) demo(Ch-LME) ---- ~~~~~~ > ### R code from vignette source 'Ch-LME.Rnw' > > ################################################### > ### code chunk number 1: setup > ################################################### > library("MVA") > set.seed(280875) > library("lattice") > lattice.options(default.theme = + function() + standard.theme("pdf", color = FALSE)) > if (file.exists("deparse.R")) { + if (!file.exists("figs")) dir.create("figs") + source("deparse.R") + options(prompt = "R> ", continue = "+ ", width = 64, + digits = 4, show.signif.stars = FALSE, useFancyQuotes = FALSE) + + options(SweaveHooks = list(onefig = function() {par(mfrow = c(1,1))}, + twofig = function() {par(mfrow = c(1,2))}, + figtwo = function() {par(mfrow = c(2,1))}, + threefig = function() {par(mfrow = c(1,3))}, + figthree = function() {par(mfrow = c(3,1))}, + fourfig = function() {par(mfrow = c(2,2))}, + sixfig = function() {par(mfrow = c(3,2))}, + nomar = function() par("mai" = c(0, 0, 0, 0)))) + } > ################################################### > ### code chunk number 2: setup > ################################################### > library("nlme") > ################################################### > ### code chunk number 3: ch:LME:ex > ################################################### > exd <- data.frame(ID = factor(1:6), Group = gl(2, 3), + matrix(rpois(6 * 4, lambda = 10), nrow = 6)) > colnames(exd)[-(1:2)] <- paste("Day", c(1, 2, 5, 7), sep = ".") > ex_wide <- exd > ################################################### > ### code chunk number 4: ch:LME:wide > ################################################### > ex_wide ID Group Day.1 Day.2 Day.5 Day.7 1 1 1 15 15 10 7 2 2 1 10 9 11 12 3 3 1 8 7 6 9 4 4 2 11 8 13 7 5 5 2 11 12 11 11 6 6 2 12 12 6 10 > ################################################### > ### code chunk number 5: ch:LME:wide > ################################################### > reshape(ex_wide, direction = "long", idvar = "ID", + varying = colnames(ex_wide)[-(1:2)]) ID Group time Day 1.1 1 1 1 15 2.1 2 1 1 10 3.1 3 1 1 8 4.1 4 2 1 11 5.1 5 2 1 11 6.1 6 2 1 12 1.2 1 1 2 15 2.2 2 1 2 9 3.2 3 1 2 7 4.2 4 2 2 8 5.2 5 2 2 12 6.2 6 2 2 12 1.5 1 1 5 10 2.5 2 1 5 11 3.5 3 1 5 6 4.5 4 2 5 13 5.5 5 2 5 11 6.5 6 2 5 6 1.7 1 1 7 7 2.7 2 1 7 12 3.7 3 1 7 9 4.7 4 2 7 7 5.7 5 2 7 11 6.7 6 2 7 10 > ################################################### > ### code chunk number 6: ch:LME:timber:tab > ################################################### > "timber" <- + matrix(c(0., 0., 0., 0., 0., 0., 0., 0., 2.3799999999999999, 2.6899999999999999, 2.8500000000000001, 2.46, + 2.9700000000000002, 3.96, 3.1699999999999999, 3.3599999999999999, 4.3399999999999999, 4.75, + 4.8899999999999997, 4.2800000000000002, 4.6799999999999997, 6.46, 5.3300000000000001, 5.4500000000000002, + 6.6399999999999997, 7.04, 6.6100000000000003, 5.8799999999999999, 6.6600000000000001, 8.1400000000000006, + 7.1399999999999997, 7.0800000000000001, 8.0500000000000007, 9.1999999999999993, 8.0899999999999999, + 7.4299999999999997, 8.1099999999999994, 9.3499999999999996, 8.2899999999999991, 8.3200000000000003, + 9.7799999999999994, 10.94, 9.7200000000000006, 8.3200000000000003, 9.6400000000000006, 10.720000000000001, + 9.8599999999999994, 9.9100000000000001, 10.970000000000001, 12.23, 11.029999999999999, 9.9199999999999999, + 11.06, 11.84, 11.07, 11.06, 12.050000000000001, 13.19, 12.140000000000001, 11.1, 12.25, 12.85, + 12.130000000000001, 12.210000000000001, 12.98, 14.08, 13.18, 12.23, 13.35, 13.83, 13.15, 13.16, 13.94, + 14.66, 14.119999999999999, 13.24, 14.539999999999999, 14.85, 14.09, 14.050000000000001, 14.74, + 15.369999999999999, 15.09, 14.19, 15.529999999999999, 15.789999999999999, 15.109999999999999, + 14.960000000000001, 16.129999999999999, 16.890000000000001, 16.68, 16.07, 17.379999999999999, + 17.390000000000001, 16.690000000000001, 16.239999999999998, 17.98, 17.780000000000001, 17.940000000000001, + 17.43, 18.760000000000002, 18.440000000000001, 17.690000000000001, 17.34, 19.52, 18.41, 18.219999999999999, + 18.359999999999999, 19.809999999999999, 19.460000000000001, 18.710000000000001, 18.23, 19.969999999999999, + 18.969999999999999, 19.399999999999999, 18.93, 20.620000000000001, 20.050000000000001, 19.539999999999999, + 18.870000000000001) + , nrow = 8, ncol = 15) > slippage <- c((0:10)/10, seq(from = 1.2, to = 1.8, by = 0.2)) > colnames(timber) <- paste("s", slippage, sep = "") > timber <- as.data.frame(timber) > timber$specimen <- factor(paste("spec", 1:nrow(timber), sep = "")) > timber.dat <- reshape(timber, direction = "long", idvar = "specimen", + varying = matrix(colnames(timber)[1:15], nr = 1), + timevar = "slippage") > names(timber.dat)[3] <- "loads" > timber.dat$slippage <- slippage[timber.dat$slippage] > timber <- timber.dat > toLatex(HSAURtable(timber), pcol = 2, + caption = paste("Data giving loads needed for a given slippage", + "in eight specimens of timber, with data in ``long'' form."), + label = "ch:LME:timber:tab", rownames = FALSE) \index{timber data@\Robject{timber} data} \begin{center} \begin{longtable}{ rrr|rrr } \caption{\Robject{timber} data. Data giving loads needed for a given slippage in eight specimens of timber, with data in ``long'' form. \label{ch:LME:timber:tab}} \\ \hline \Robject{specimen} & \Robject{slippage} & \Robject{loads} & \Robject{specimen} & \Robject{slippage} & \Robject{loads} \\ \hline \endfirsthead \caption[]{\Robject{timber} data (continued).} \\ \hline \Robject{specimen} & \Robject{slippage} & \Robject{loads} & \Robject{specimen} & \Robject{slippage} & \Robject{loads} \\ \hline \endhead spec1 & 0.0 & 0.00 & spec5 & 0.7 & 12.25 \\ spec2 & 0.0 & 0.00 & spec6 & 0.7 & 12.85 \\ spec3 & 0.0 & 0.00 & spec7 & 0.7 & 12.13 \\ spec4 & 0.0 & 0.00 & spec8 & 0.7 & 12.21 \\ spec5 & 0.0 & 0.00 & spec1 & 0.8 & 12.98 \\ spec6 & 0.0 & 0.00 & spec2 & 0.8 & 14.08 \\ spec7 & 0.0 & 0.00 & spec3 & 0.8 & 13.18 \\ spec8 & 0.0 & 0.00 & spec4 & 0.8 & 12.23 \\ spec1 & 0.1 & 2.38 & spec5 & 0.8 & 13.35 \\ spec2 & 0.1 & 2.69 & spec6 & 0.8 & 13.83 \\ spec3 & 0.1 & 2.85 & spec7 & 0.8 & 13.15 \\ spec4 & 0.1 & 2.46 & spec8 & 0.8 & 13.16 \\ spec5 & 0.1 & 2.97 & spec1 & 0.9 & 13.94 \\ spec6 & 0.1 & 3.96 & spec2 & 0.9 & 14.66 \\ spec7 & 0.1 & 3.17 & spec3 & 0.9 & 14.12 \\ spec8 & 0.1 & 3.36 & spec4 & 0.9 & 13.24 \\ spec1 & 0.2 & 4.34 & spec5 & 0.9 & 14.54 \\ spec2 & 0.2 & 4.75 & spec6 & 0.9 & 14.85 \\ spec3 & 0.2 & 4.89 & spec7 & 0.9 & 14.09 \\ spec4 & 0.2 & 4.28 & spec8 & 0.9 & 14.05 \\ spec5 & 0.2 & 4.68 & spec1 & 1.0 & 14.74 \\ spec6 & 0.2 & 6.46 & spec2 & 1.0 & 15.37 \\ spec7 & 0.2 & 5.33 & spec3 & 1.0 & 15.09 \\ spec8 & 0.2 & 5.45 & spec4 & 1.0 & 14.19 \\ spec1 & 0.3 & 6.64 & spec5 & 1.0 & 15.53 \\ spec2 & 0.3 & 7.04 & spec6 & 1.0 & 15.79 \\ spec3 & 0.3 & 6.61 & spec7 & 1.0 & 15.11 \\ spec4 & 0.3 & 5.88 & spec8 & 1.0 & 14.96 \\ spec5 & 0.3 & 6.66 & spec1 & 1.2 & 16.13 \\ spec6 & 0.3 & 8.14 & spec2 & 1.2 & 16.89 \\ spec7 & 0.3 & 7.14 & spec3 & 1.2 & 16.68 \\ spec8 & 0.3 & 7.08 & spec4 & 1.2 & 16.07 \\ spec1 & 0.4 & 8.05 & spec5 & 1.2 & 17.38 \\ spec2 & 0.4 & 9.20 & spec6 & 1.2 & 17.39 \\ spec3 & 0.4 & 8.09 & spec7 & 1.2 & 16.69 \\ spec4 & 0.4 & 7.43 & spec8 & 1.2 & 16.24 \\ spec5 & 0.4 & 8.11 & spec1 & 1.4 & 17.98 \\ spec6 & 0.4 & 9.35 & spec2 & 1.4 & 17.78 \\ spec7 & 0.4 & 8.29 & spec3 & 1.4 & 17.94 \\ spec8 & 0.4 & 8.32 & spec4 & 1.4 & 17.43 \\ spec1 & 0.5 & 9.78 & spec5 & 1.4 & 18.76 \\ spec2 & 0.5 & 10.94 & spec6 & 1.4 & 18.44 \\ spec3 & 0.5 & 9.72 & spec7 & 1.4 & 17.69 \\ spec4 & 0.5 & 8.32 & spec8 & 1.4 & 17.34 \\ spec5 & 0.5 & 9.64 & spec1 & 1.6 & 19.52 \\ spec6 & 0.5 & 10.72 & spec2 & 1.6 & 18.41 \\ spec7 & 0.5 & 9.86 & spec3 & 1.6 & 18.22 \\ spec8 & 0.5 & 9.91 & spec4 & 1.6 & 18.36 \\ spec1 & 0.6 & 10.97 & spec5 & 1.6 & 19.81 \\ spec2 & 0.6 & 12.23 & spec6 & 1.6 & 19.46 \\ spec3 & 0.6 & 11.03 & spec7 & 1.6 & 18.71 \\ spec4 & 0.6 & 9.92 & spec8 & 1.6 & 18.23 \\ spec5 & 0.6 & 11.06 & spec1 & 1.8 & 19.97 \\ spec6 & 0.6 & 11.84 & spec2 & 1.8 & 18.97 \\ spec7 & 0.6 & 11.07 & spec3 & 1.8 & 19.40 \\ spec8 & 0.6 & 11.06 & spec4 & 1.8 & 18.93 \\ spec1 & 0.7 & 12.05 & spec5 & 1.8 & 20.62 \\ spec2 & 0.7 & 13.19 & spec6 & 1.8 & 20.05 \\ spec3 & 0.7 & 12.14 & spec7 & 1.8 & 19.54 \\ spec4 & 0.7 & 11.10 & spec8 & 1.8 & 18.87 \\ \hline \end{longtable} \end{center} > ################################################### > ### code chunk number 7: ch:LME:plasma:tab > ################################################### > "plasma" <- + matrix(c(4.2999999999999998, 3.7000000000000002, 4., 3.6000000000000001, 4.0999999999999996, 3.7999999999999998, + 3.7999999999999998, 4.4000000000000004, 5., 3.7000000000000002, 3.7000000000000002, 4.4000000000000004, + 4.7000000000000002, 4.2999999999999998, 5., 4.5999999999999996, 4.2999999999999998, 3.1000000000000001, + 4.7999999999999998, 3.7000000000000002, 5.4000000000000004, 3., 4.9000000000000004, 4.7999999999999998, + 4.4000000000000004, 4.9000000000000004, 5.0999999999999996, 4.7999999999999998, 4.2000000000000002, + 6.5999999999999996, 3.6000000000000001, 4.5, 4.5999999999999996, 3.2999999999999998, 2.6000000000000001, + 4.0999999999999996, 3., 3.7999999999999998, 2.2000000000000002, 3., 3.8999999999999999, 4., + 3.1000000000000001, 2.6000000000000001, 3.7000000000000002, 3.1000000000000001, 3.2999999999999998, + 4.9000000000000004, 4.4000000000000004, 3.8999999999999999, 3.1000000000000001, 5., 3.1000000000000001, + 4.7000000000000002, 2.5, 5., 4.2999999999999998, 4.2000000000000002, 4.2999999999999998, 4.0999999999999996, + 4.5999999999999996, 3.5, 6.0999999999999996, 3.3999999999999999, 4., 4.4000000000000004, 3., + 2.6000000000000001, 3.1000000000000001, 2.2000000000000002, 2.1000000000000001, 2., 2.3999999999999999, + 2.7999999999999998, 3.3999999999999999, 2.8999999999999999, 2.6000000000000001, 3.1000000000000001, + 3.2000000000000002, 3., 4.0999999999999996, 3.8999999999999999, 3.1000000000000001, 3.2999999999999998, + 2.8999999999999999, 3.2999999999999998, 3.8999999999999999, 2.2999999999999998, 4.0999999999999996, + 4.7000000000000002, 4.2000000000000002, 4., 4.5999999999999996, 4.5999999999999996, 3.7999999999999998, + 5.2000000000000002, 3.1000000000000001, 3.7000000000000002, 3.7999999999999998, 2.6000000000000001, + 1.8999999999999999, 2.2999999999999998, 2.7999999999999998, 3., 2.6000000000000001, 2.5, 2.1000000000000001, + 3.3999999999999999, 2.2000000000000002, 2.2999999999999998, 3.2000000000000002, 3.2999999999999998, + 2.6000000000000001, 3.7000000000000002, 3.8999999999999999, 3.1000000000000001, 2.6000000000000001, + 2.7999999999999998, 2.7999999999999998, 4.0999999999999996, 2.2000000000000002, 3.7000000000000002, + 4.5999999999999996, 3.3999999999999999, 4., 4.0999999999999996, 4.4000000000000004, 3.6000000000000001, + 4.0999999999999996, 2.7999999999999998, 3.2999999999999998, 3.7999999999999998, 2.2000000000000002, + 2.8999999999999999, 2.8999999999999999, 2.8999999999999999, 3.6000000000000001, 3.7999999999999998, + 3.1000000000000001, 3.6000000000000001, 3.2999999999999998, 1.5, 2.8999999999999999, 3.7000000000000002, + 3.2000000000000002, 2.2000000000000002, 3.7000000000000002, 3.7000000000000002, 3.1000000000000001, + 2.6000000000000001, 2.2000000000000002, 2.8999999999999999, 2.7999999999999998, 2.1000000000000001, + 3.7000000000000002, 4.7000000000000002, 3.5, 3.2999999999999998, 3.3999999999999999, 4.0999999999999996, + 3.2999999999999998, 4.2999999999999998, 2.1000000000000001, 2.3999999999999999, 3.7999999999999998, 2.5, + 3.2000000000000002, 3.1000000000000001, 3.8999999999999999, 3.3999999999999999, 3.6000000000000001, + 3.3999999999999999, 3.7999999999999998, 3.6000000000000001, 2.2999999999999998, 2.2000000000000002, + 4.2999999999999998, 4.2000000000000002, 2.5, 4.0999999999999996, 4.2000000000000002, 3.1000000000000001, + 1.8999999999999999, 3.1000000000000001, 3.6000000000000001, 3.7000000000000002, 2.6000000000000001, + 4.0999999999999996, 3.7000000000000002, 3.3999999999999999, 4.0999999999999996, 4.2000000000000002, 4., + 3.1000000000000001, 3.7999999999999998, 2.3999999999999999, 2.2999999999999998, 3.6000000000000001, + 3.3999999999999999, 3.1000000000000001, 3.8999999999999999, 3.7999999999999998, 3.6000000000000001, 3., + 3.5, 4., 4., 2.7000000000000002, 3.1000000000000001, 3.8999999999999999, 3.7000000000000002, + 2.3999999999999999, 4.7000000000000002, 4.7999999999999998, 3.6000000000000001, 2.2999999999999998, 3.5, + 4.2999999999999998, 3.5, 3.2000000000000002, 4.7000000000000002, 3.6000000000000001, 3.7999999999999998, + 4.2000000000000002, 4.4000000000000004, 3.7999999999999998, 3.5, 4.2000000000000002, 2.5, + 3.1000000000000001, 3.7999999999999998, 4.4000000000000004, 3.8999999999999999, 4., 4., 3.7000000000000002, + 3.5, 3.7000000000000002, 3.8999999999999999, 4.2999999999999998, 2.7999999999999998, 3.8999999999999999, + 4.7999999999999998, 4.2999999999999998, 3.3999999999999999, 4.9000000000000004, 5., 4., 2.7000000000000002, + 3.6000000000000001, 4.4000000000000004, 3.7000000000000002, 3.5, 4.9000000000000004, 3.8999999999999999, + 4., 4.2999999999999998, 4.9000000000000004, 3.7999999999999998, 3.8999999999999999, 4.7999999999999998, + 3.5, 3.2999999999999998, 3.7999999999999998) + , nrow = 33, ncol = 8 + , dimnames = list(character(0) + , c("T0.0", "T0.5", "T1.0", "T1.5", "T2.0", "T2.5", "T3.0", "T4.0") + ) + ) > time <- c(0, 0.5, 1, 1.5, 2, 3, 4) > plasma <- as.data.frame(plasma) > plasma$Subject <- factor(paste("id", + formatC(1:nrow(plasma), format = "g", width = 2, flag = "0"), sep = "")) > plasma$group <- factor(c(rep("control", 20), rep("obese", 13))) > plasma <- reshape(plasma, direction = "long", idvar = "Subject", + varying = matrix(colnames(plasma)[1:8], nr = 1), + timevar = "time") > colnames(plasma)[4] <- "plasma" > toLatex(HSAURtable(plasma), pcol = 2, + caption = "Plasma inorganic phosphate levels from 33 subjects, with data in ``long'' form.", + label = "ch:LME:plasma:tab", rownames = FALSE) \index{plasma data@\Robject{plasma} data} \begin{center} \begin{longtable}{ rrrr|rrrr } \caption{\Robject{plasma} data. Plasma inorganic phosphate levels from 33 subjects, with data in ``long'' form. \label{ch:LME:plasma:tab}} \\ \hline \Robject{Subject} & \Robject{group} & \Robject{time} & \Robject{plasma} & \Robject{Subject} & \Robject{group} & \Robject{time} & \Robject{plasma} \\ \hline \endfirsthead \caption[]{\Robject{plasma} data (continued).} \\ \hline \Robject{Subject} & \Robject{group} & \Robject{time} & \Robject{plasma} & \Robject{Subject} & \Robject{group} & \Robject{time} & \Robject{plasma} \\ \hline \endhead id01 & control & 1 & 4.3 & id01 & control & 5 & 2.2 \\ id02 & control & 1 & 3.7 & id02 & control & 5 & 2.9 \\ id03 & control & 1 & 4.0 & id03 & control & 5 & 2.9 \\ id04 & control & 1 & 3.6 & id04 & control & 5 & 2.9 \\ id05 & control & 1 & 4.1 & id05 & control & 5 & 3.6 \\ id06 & control & 1 & 3.8 & id06 & control & 5 & 3.8 \\ id07 & control & 1 & 3.8 & id07 & control & 5 & 3.1 \\ id08 & control & 1 & 4.4 & id08 & control & 5 & 3.6 \\ id09 & control & 1 & 5.0 & id09 & control & 5 & 3.3 \\ id10 & control & 1 & 3.7 & id10 & control & 5 & 1.5 \\ id11 & control & 1 & 3.7 & id11 & control & 5 & 2.9 \\ id12 & control & 1 & 4.4 & id12 & control & 5 & 3.7 \\ id13 & control & 1 & 4.7 & id13 & control & 5 & 3.2 \\ id14 & control & 1 & 4.3 & id14 & control & 5 & 2.2 \\ id15 & control & 1 & 5.0 & id15 & control & 5 & 3.7 \\ id16 & control & 1 & 4.6 & id16 & control & 5 & 3.7 \\ id17 & control & 1 & 4.3 & id17 & control & 5 & 3.1 \\ id18 & control & 1 & 3.1 & id18 & control & 5 & 2.6 \\ id19 & control & 1 & 4.8 & id19 & control & 5 & 2.2 \\ id20 & control & 1 & 3.7 & id20 & control & 5 & 2.9 \\ id21 & obese & 1 & 5.4 & id21 & obese & 5 & 2.8 \\ id22 & obese & 1 & 3.0 & id22 & obese & 5 & 2.1 \\ id23 & obese & 1 & 4.9 & id23 & obese & 5 & 3.7 \\ id24 & obese & 1 & 4.8 & id24 & obese & 5 & 4.7 \\ id25 & obese & 1 & 4.4 & id25 & obese & 5 & 3.5 \\ id26 & obese & 1 & 4.9 & id26 & obese & 5 & 3.3 \\ id27 & obese & 1 & 5.1 & id27 & obese & 5 & 3.4 \\ id28 & obese & 1 & 4.8 & id28 & obese & 5 & 4.1 \\ id29 & obese & 1 & 4.2 & id29 & obese & 5 & 3.3 \\ id30 & obese & 1 & 6.6 & id30 & obese & 5 & 4.3 \\ id31 & obese & 1 & 3.6 & id31 & obese & 5 & 2.1 \\ id32 & obese & 1 & 4.5 & id32 & obese & 5 & 2.4 \\ id33 & obese & 1 & 4.6 & id33 & obese & 5 & 3.8 \\ id01 & control & 2 & 3.3 & id01 & control & 6 & 2.5 \\ id02 & control & 2 & 2.6 & id02 & control & 6 & 3.2 \\ id03 & control & 2 & 4.1 & id03 & control & 6 & 3.1 \\ id04 & control & 2 & 3.0 & id04 & control & 6 & 3.9 \\ id05 & control & 2 & 3.8 & id05 & control & 6 & 3.4 \\ id06 & control & 2 & 2.2 & id06 & control & 6 & 3.6 \\ id07 & control & 2 & 3.0 & id07 & control & 6 & 3.4 \\ id08 & control & 2 & 3.9 & id08 & control & 6 & 3.8 \\ id09 & control & 2 & 4.0 & id09 & control & 6 & 3.6 \\ id10 & control & 2 & 3.1 & id10 & control & 6 & 2.3 \\ id11 & control & 2 & 2.6 & id11 & control & 6 & 2.2 \\ id12 & control & 2 & 3.7 & id12 & control & 6 & 4.3 \\ id13 & control & 2 & 3.1 & id13 & control & 6 & 4.2 \\ id14 & control & 2 & 3.3 & id14 & control & 6 & 2.5 \\ id15 & control & 2 & 4.9 & id15 & control & 6 & 4.1 \\ id16 & control & 2 & 4.4 & id16 & control & 6 & 4.2 \\ id17 & control & 2 & 3.9 & id17 & control & 6 & 3.1 \\ id18 & control & 2 & 3.1 & id18 & control & 6 & 1.9 \\ id19 & control & 2 & 5.0 & id19 & control & 6 & 3.1 \\ id20 & control & 2 & 3.1 & id20 & control & 6 & 3.6 \\ id21 & obese & 2 & 4.7 & id21 & obese & 6 & 3.7 \\ id22 & obese & 2 & 2.5 & id22 & obese & 6 & 2.6 \\ id23 & obese & 2 & 5.0 & id23 & obese & 6 & 4.1 \\ id24 & obese & 2 & 4.3 & id24 & obese & 6 & 3.7 \\ id25 & obese & 2 & 4.2 & id25 & obese & 6 & 3.4 \\ id26 & obese & 2 & 4.3 & id26 & obese & 6 & 4.1 \\ id27 & obese & 2 & 4.1 & id27 & obese & 6 & 4.2 \\ id28 & obese & 2 & 4.6 & id28 & obese & 6 & 4.0 \\ id29 & obese & 2 & 3.5 & id29 & obese & 6 & 3.1 \\ id30 & obese & 2 & 6.1 & id30 & obese & 6 & 3.8 \\ id31 & obese & 2 & 3.4 & id31 & obese & 6 & 2.4 \\ id32 & obese & 2 & 4.0 & id32 & obese & 6 & 2.3 \\ id33 & obese & 2 & 4.4 & id33 & obese & 6 & 3.6 \\ id01 & control & 3 & 3.0 & id01 & control & 7 & 3.4 \\ id02 & control & 3 & 2.6 & id02 & control & 7 & 3.1 \\ id03 & control & 3 & 3.1 & id03 & control & 7 & 3.9 \\ id04 & control & 3 & 2.2 & id04 & control & 7 & 3.8 \\ id05 & control & 3 & 2.1 & id05 & control & 7 & 3.6 \\ id06 & control & 3 & 2.0 & id06 & control & 7 & 3.0 \\ id07 & control & 3 & 2.4 & id07 & control & 7 & 3.5 \\ id08 & control & 3 & 2.8 & id08 & control & 7 & 4.0 \\ id09 & control & 3 & 3.4 & id09 & control & 7 & 4.0 \\ id10 & control & 3 & 2.9 & id10 & control & 7 & 2.7 \\ id11 & control & 3 & 2.6 & id11 & control & 7 & 3.1 \\ id12 & control & 3 & 3.1 & id12 & control & 7 & 3.9 \\ id13 & control & 3 & 3.2 & id13 & control & 7 & 3.7 \\ id14 & control & 3 & 3.0 & id14 & control & 7 & 2.4 \\ id15 & control & 3 & 4.1 & id15 & control & 7 & 4.7 \\ id16 & control & 3 & 3.9 & id16 & control & 7 & 4.8 \\ id17 & control & 3 & 3.1 & id17 & control & 7 & 3.6 \\ id18 & control & 3 & 3.3 & id18 & control & 7 & 2.3 \\ id19 & control & 3 & 2.9 & id19 & control & 7 & 3.5 \\ id20 & control & 3 & 3.3 & id20 & control & 7 & 4.3 \\ id21 & obese & 3 & 3.9 & id21 & obese & 7 & 3.5 \\ id22 & obese & 3 & 2.3 & id22 & obese & 7 & 3.2 \\ id23 & obese & 3 & 4.1 & id23 & obese & 7 & 4.7 \\ id24 & obese & 3 & 4.7 & id24 & obese & 7 & 3.6 \\ id25 & obese & 3 & 4.2 & id25 & obese & 7 & 3.8 \\ id26 & obese & 3 & 4.0 & id26 & obese & 7 & 4.2 \\ id27 & obese & 3 & 4.6 & id27 & obese & 7 & 4.4 \\ id28 & obese & 3 & 4.6 & id28 & obese & 7 & 3.8 \\ id29 & obese & 3 & 3.8 & id29 & obese & 7 & 3.5 \\ id30 & obese & 3 & 5.2 & id30 & obese & 7 & 4.2 \\ id31 & obese & 3 & 3.1 & id31 & obese & 7 & 2.5 \\ id32 & obese & 3 & 3.7 & id32 & obese & 7 & 3.1 \\ id33 & obese & 3 & 3.8 & id33 & obese & 7 & 3.8 \\ id01 & control & 4 & 2.6 & id01 & control & 8 & 4.4 \\ id02 & control & 4 & 1.9 & id02 & control & 8 & 3.9 \\ id03 & control & 4 & 2.3 & id03 & control & 8 & 4.0 \\ id04 & control & 4 & 2.8 & id04 & control & 8 & 4.0 \\ id05 & control & 4 & 3.0 & id05 & control & 8 & 3.7 \\ id06 & control & 4 & 2.6 & id06 & control & 8 & 3.5 \\ id07 & control & 4 & 2.5 & id07 & control & 8 & 3.7 \\ id08 & control & 4 & 2.1 & id08 & control & 8 & 3.9 \\ id09 & control & 4 & 3.4 & id09 & control & 8 & 4.3 \\ id10 & control & 4 & 2.2 & id10 & control & 8 & 2.8 \\ id11 & control & 4 & 2.3 & id11 & control & 8 & 3.9 \\ id12 & control & 4 & 3.2 & id12 & control & 8 & 4.8 \\ id13 & control & 4 & 3.3 & id13 & control & 8 & 4.3 \\ id14 & control & 4 & 2.6 & id14 & control & 8 & 3.4 \\ id15 & control & 4 & 3.7 & id15 & control & 8 & 4.9 \\ id16 & control & 4 & 3.9 & id16 & control & 8 & 5.0 \\ id17 & control & 4 & 3.1 & id17 & control & 8 & 4.0 \\ id18 & control & 4 & 2.6 & id18 & control & 8 & 2.7 \\ id19 & control & 4 & 2.8 & id19 & control & 8 & 3.6 \\ id20 & control & 4 & 2.8 & id20 & control & 8 & 4.4 \\ id21 & obese & 4 & 4.1 & id21 & obese & 8 & 3.7 \\ id22 & obese & 4 & 2.2 & id22 & obese & 8 & 3.5 \\ id23 & obese & 4 & 3.7 & id23 & obese & 8 & 4.9 \\ id24 & obese & 4 & 4.6 & id24 & obese & 8 & 3.9 \\ id25 & obese & 4 & 3.4 & id25 & obese & 8 & 4.0 \\ id26 & obese & 4 & 4.0 & id26 & obese & 8 & 4.3 \\ id27 & obese & 4 & 4.1 & id27 & obese & 8 & 4.9 \\ id28 & obese & 4 & 4.4 & id28 & obese & 8 & 3.8 \\ id29 & obese & 4 & 3.6 & id29 & obese & 8 & 3.9 \\ id30 & obese & 4 & 4.1 & id30 & obese & 8 & 4.8 \\ id31 & obese & 4 & 2.8 & id31 & obese & 8 & 3.5 \\ id32 & obese & 4 & 3.3 & id32 & obese & 8 & 3.3 \\ id33 & obese & 4 & 3.8 & id33 & obese & 8 & 3.8 \\ \hline \end{longtable} \end{center} > ################################################### > ### code chunk number 8: ch:LME:timber:lm > ################################################### > summary(lm(loads ~ slippage, data = timber)) Call: lm(formula = loads ~ slippage, data = timber) Residuals: Min 1Q Median 3Q Max -3.5158 -0.9810 0.4075 1.2982 2.4907 Coefficients: Estimate Std. Error t value Pr(>|t|) (Intercept) 3.5158 0.2644 13.30 <2e-16 *** slippage 10.3726 0.2835 36.59 <2e-16 *** --- Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1 Residual standard error: 1.65 on 118 degrees of freedom Multiple R-squared: 0.919, Adjusted R-squared: 0.9183 F-statistic: 1339 on 1 and 118 DF, p-value: < 2.2e-16 > ################################################### > ### code chunk number 9: ch:LME:timber:plot (eval = FALSE) > ################################################### > ## xyplot(loads ~ slippage | specimen, data = timber, > ## layout = c(4, 2)) > > > ################################################### > ### code chunk number 10: ch:LME:timber:plot > ################################################### > plot(xyplot(loads ~ slippage | specimen, data = timber, + layout = c(4, 2))) > ################################################### > ### code chunk number 11: ch:LME:timber:lme > ################################################### > timber.lme <- lme(loads ~ slippage, + random = ~1 | specimen, + data = timber, method = "ML") > timber.lme1 <- lme(loads ~ slippage, + random = ~slippage | specimen, + data = timber, method = "ML") > ################################################### > ### code chunk number 12: ch:LME:timber:LRT > ################################################### > library("RLRsim") > exactRLRT(timber.lme) Using restricted likelihood evaluated at ML estimators. Refit with method="REML" for exact results. simulated finite sample distribution of RLRT. (p-value based on 10000 simulated values) data: RLRT = 3.2185e-07, p-value = 0.4419 > ################################################### > ### code chunk number 13: ch:LME:timber:anova > ################################################### > anova(timber.lme, timber.lme1) Model df AIC BIC logLik Test L.Ratio p-value timber.lme 1 4 466.6631 477.8131 -229.3316 timber.lme1 2 6 470.6631 487.3881 -229.3316 1 vs 2 2.251318e-08 1 > ################################################### > ### code chunk number 14: ch:LME:timber:lme:summary > ################################################### > summary(timber.lme1) Linear mixed-effects model fit by maximum likelihood Data: timber AIC BIC logLik 470.6631 487.3881 -229.3316 Random effects: Formula: ~slippage | specimen Structure: General positive-definite, Log-Cholesky parametrization StdDev Corr (Intercept) 3.229387e-04 (Intr) slippage 1.089836e-05 0 Residual 1.635842e+00 Fixed effects: loads ~ slippage Value Std.Error DF t-value p-value (Intercept) 3.515758 0.2644018 111 13.29703 0 slippage 10.372598 0.2834685 111 36.59172 0 Correlation: (Intr) slippage -0.822 Standardized Within-Group Residuals: Min Q1 Med Q3 Max -2.1492039 -0.5996683 0.2491274 0.7935870 1.5225692 Number of Observations: 120 Number of Groups: 8 > ################################################### > ### code chunk number 15: ch:LME:timber:lme:predict > ################################################### > timber$pred1 <- predict(timber.lme1) > ################################################### > ### code chunk number 16: ch:LME:timber:plotpredict (eval = FALSE) > ################################################### > ## pfun <- function(x, y) { > ## panel.xyplot(x, y[1:length(x)]) > ## panel.lines(x, y[1:length(x) + length(x)], lty = 1) > ## } > ## plot(xyplot(cbind(loads, pred1) ~ slippage | specimen, > ## data = timber, panel = pfun, layout = c(4, 2), > ## ylab = "loads")) > > > ################################################### > ### code chunk number 17: ch:LME:timber:plotpredict > ################################################### > pfun <- function(x, y) { + panel.xyplot(x, y[1:length(x)]) + panel.lines(x, y[1:length(x) + length(x)], lty = 1) + } > plot(xyplot(cbind(loads, pred1) ~ slippage | specimen, data = timber, + panel = pfun, layout = c(4, 2), ylab = "loads")) > ################################################### > ### code chunk number 18: ch:LME:timber:quad > ################################################### > timber.lme2 <- lme(loads ~ slippage + I(slippage^2), + random = ~slippage | specimen, + data = timber, method = "ML") > anova(timber.lme1, timber.lme2) Model df AIC BIC logLik Test L.Ratio p-value timber.lme1 1 6 470.6631 487.3881 -229.33156 timber.lme2 2 7 213.5439 233.0564 -99.77197 1 vs 2 259.1192 <.0001 > ################################################### > ### code chunk number 19: ch:LME:timber:quad:plot > ################################################### > timber$pred2 <- predict(timber.lme2) > plot(xyplot(cbind(loads, pred2) ~ slippage | specimen, + data = timber, panel = pfun, layout = c(4, 2), + ylab = "loads")) > ################################################### > ### code chunk number 20: ch:LME:plasma:plot (eval = FALSE) > ################################################### > ## x <- reshape(plasma, direction = "wide", timevar = "time", > ## idvar = "Subject", v.names = "plasma") > ## X <- as.matrix(x[,-(1:2)]) > ## parallelplot(~ X | group, data = x, horizontal = FALSE, > ## col = "black", scales = list(x = list(labels = 1:8)), > ## ylab = "Plasma inorganic phosphate", > ## xlab = "Time (hours after oral glucose challenge)") > > > ################################################### > ### code chunk number 21: ch:LME:plasma:plot > ################################################### > x <- reshape(plasma, direction = "wide", timevar = "time", + idvar = "Subject", v.names = "plasma") > X <- as.matrix(x[,-(1:2)]) > plot(parallelplot(~ X | group, data = x, horizontal = FALSE, + col = "black", scales = list(x = list(labels = 1:8)), + ylab = "Plasma inorganic phosphate", + xlab = "Time (hours after oral glucose challenge)")) > ################################################### > ### code chunk number 22: ch:LME:plasma:splom > ################################################### > plot(splom(~ x[, grep("plasma", colnames(x))] | group, data = x, + cex = 1.5, pch = ".", pscales = NULL, varnames = 1:8)) > ################################################### > ### code chunk number 23: ch:LME:plasma:lme > ################################################### > plasma.lme1 <- lme(plasma ~ time + I(time^2) + group, + random = ~ time | Subject, + data = plasma, method = "ML") > summary(plasma.lme1) Linear mixed-effects model fit by maximum likelihood Data: plasma AIC BIC logLik 390.4724 419.08 -187.2362 Random effects: Formula: ~time | Subject Structure: General positive-definite, Log-Cholesky parametrization StdDev Corr (Intercept) 0.69771527 (Intr) time 0.09382651 -0.7 Residual 0.38480107 Fixed effects: plasma ~ time + I(time^2) + group Value Std.Error DF t-value p-value (Intercept) 4.879775 0.17090659 229 28.552295 0.0000 time -0.803283 0.05075453 229 -15.826820 0.0000 I(time^2) 0.084668 0.00520763 229 16.258479 0.0000 groupobese 0.437054 0.18588679 31 2.351181 0.0253 Correlation: (Intr) time I(t^2) time -0.641 I(time^2) 0.457 -0.923 groupobese -0.428 0.000 0.000 Standardized Within-Group Residuals: Min Q1 Med Q3 Max -2.771508466 -0.548688207 -0.002764541 0.564434549 2.889632993 Number of Observations: 264 Number of Groups: 33 > ################################################### > ### code chunk number 24: ch:LME:plasma:lm > ################################################### > summary(lm(plasma ~ time + I(time^2) + group, data = plasma)) Call: lm(formula = plasma ~ time + I(time^2) + group, data = plasma) Residuals: Min 1Q Median 3Q Max -1.63231 -0.44010 0.03475 0.47499 2.01696 Coefficients: Estimate Std. Error t value Pr(>|t|) (Intercept) 4.857611 0.166860 29.112 < 2e-16 *** time -0.803283 0.083350 -9.637 < 2e-16 *** I(time^2) 0.084668 0.009041 9.365 < 2e-16 *** groupobese 0.493317 0.084788 5.818 1.74e-08 *** --- Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1 Residual standard error: 0.6731 on 260 degrees of freedom Multiple R-squared: 0.3278, Adjusted R-squared: 0.32 F-statistic: 42.26 on 3 and 260 DF, p-value: < 2.2e-16 > ################################################### > ### code chunk number 25: ch:LME:plasma:predictplot1 > ################################################### > pfun <- function(x, y, subscripts, groups) { + panel.xyplot(x, y[1:length(x)], pch = c(1:2)[groups[subscripts]]) + panel.lines(x, y[1:length(x) + length(x)], lty = 1) + } > plasma$pred1 <- predict(plasma.lme1) > plot(xyplot(cbind(plasma, pred1) ~ time | Subject, data = plasma, groups = group, + type = "b", + pch = c(1, 2), ylab = "Plasma inorganic phosphate", + xlab = "Time (hours after oral glucose challenge)", panel = pfun)) > ################################################### > ### code chunk number 26: ch:LME:plasma:inter > ################################################### > plasma.lme2 <- lme(plasma ~ time*group +I(time^2), + random = ~time | Subject, + data = plasma, method = "ML") > ################################################### > ### code chunk number 27: ch:LME:plasma:interaov > ################################################### > anova(plasma.lme1, plasma.lme2) Model df AIC BIC logLik Test L.Ratio p-value plasma.lme1 1 8 390.4724 419.0800 -187.2362 plasma.lme2 2 9 383.3149 415.4984 -182.6575 1 vs 2 9.157478 0.0025 > ################################################### > ### code chunk number 28: ch:LME:plasma:summary > ################################################### > summary(plasma.lme2) Linear mixed-effects model fit by maximum likelihood Data: plasma AIC BIC logLik 383.3149 415.4984 -182.6574 Random effects: Formula: ~time | Subject Structure: General positive-definite, Log-Cholesky parametrization StdDev Corr (Intercept) 0.64189651 (Intr) time 0.07626155 -0.631 Residual 0.38480121 Fixed effects: plasma ~ time * group + I(time^2) Value Std.Error DF t-value p-value (Intercept) 4.659307 0.17806323 228 26.166589 0.0000 time -0.759215 0.05178067 228 -14.662140 0.0000 groupobese 0.996703 0.25482663 31 3.911300 0.0005 I(time^2) 0.084668 0.00521767 228 16.227177 0.0000 time:groupobese -0.111864 0.03476379 228 -3.217818 0.0015 Correlation: (Intr) time gropbs I(t^2) time -0.657 groupobese -0.564 0.181 I(time^2) 0.440 -0.907 0.000 time:groupobese 0.385 -0.264 -0.683 0.000 Standardized Within-Group Residuals: Min Q1 Med Q3 Max -2.72435961 -0.53605316 -0.01071299 0.58567768 2.95029238 Number of Observations: 264 Number of Groups: 33 > ################################################### > ### code chunk number 29: ch:LME:plasma:predictplot2 > ################################################### > plasma$pred2 <- predict(plasma.lme2) > plot(xyplot(cbind(plasma, pred2) ~ time | Subject, data = plasma, groups = group, + type = "b", + pch = c(1, 2), ylab = "Plasma inorganic phosphate", + xlab = "Time (hours after oral glucose challenge)", panel = pfun)) > ################################################### > ### code chunk number 30: ch:LME:plasma:random > ################################################### > res.int <- random.effects(plasma.lme2)[,1] > res.slope <- random.effects(plasma.lme2)[,2] > ################################################### > ### code chunk number 31: ch:LME:plasma:fitplot > ################################################### > qqnorm(res.int,ylab="Estimated random intercepts",main="Random intercepts") > qqnorm(res.slope,ylab="Estimated random slopes",main="Random slopes") > resids<-resid(plasma.lme2) > qqnorm(resids,ylab="Estimated residuals",main="Residuals") > ################################################### > ### code chunk number 32: ch:LME:BtheB:tab > ################################################### > data("BtheB", package = "HSAUR2") > toLatex(HSAURtable(BtheB), pcol = 1, + caption = "Data of a randomised trial evaluating the effects of Beat the Blues.", + label = "ch:LME:BtheB:tab") \index{BtheB data@\Robject{BtheB} data} \begin{center} \begin{longtable}{ rrrrrrrr } \caption{\Robject{BtheB} data. Data of a randomised trial evaluating the effects of Beat the Blues. \label{ch:LME:BtheB:tab}} \\ \hline \Robject{drug} & \Robject{length} & \Robject{treatment} & \Robject{bdi.pre} & \Robject{bdi.2m} & \Robject{bdi.3m} & \Robject{bdi.5m} & \Robject{bdi.8m} \\ \hline \endfirsthead \caption[]{\Robject{BtheB} data (continued).} \\ \hline \Robject{drug} & \Robject{length} & \Robject{treatment} & \Robject{bdi.pre} & \Robject{bdi.2m} & \Robject{bdi.3m} & \Robject{bdi.5m} & \Robject{bdi.8m} \\ \hline \endhead No & >6m & TAU & 29 & 2 & 2 & NA & NA \\ Yes & >6m & BtheB & 32 & 16 & 24 & 17 & 20 \\ Yes & <6m & TAU & 25 & 20 & NA & NA & NA \\ No & >6m & BtheB & 21 & 17 & 16 & 10 & 9 \\ Yes & >6m & BtheB & 26 & 23 & NA & NA & NA \\ Yes & <6m & BtheB & 7 & 0 & 0 & 0 & 0 \\ Yes & <6m & TAU & 17 & 7 & 7 & 3 & 7 \\ No & >6m & TAU & 20 & 20 & 21 & 19 & 13 \\ Yes & <6m & BtheB & 18 & 13 & 14 & 20 & 11 \\ Yes & >6m & BtheB & 20 & 5 & 5 & 8 & 12 \\ No & >6m & TAU & 30 & 32 & 24 & 12 & 2 \\ Yes & <6m & BtheB & 49 & 35 & NA & NA & NA \\ No & >6m & TAU & 26 & 27 & 23 & NA & NA \\ Yes & >6m & TAU & 30 & 26 & 36 & 27 & 22 \\ Yes & >6m & BtheB & 23 & 13 & 13 & 12 & 23 \\ No & <6m & TAU & 16 & 13 & 3 & 2 & 0 \\ No & >6m & BtheB & 30 & 30 & 29 & NA & NA \\ No & <6m & BtheB & 13 & 8 & 8 & 7 & 6 \\ No & >6m & TAU & 37 & 30 & 33 & 31 & 22 \\ Yes & <6m & BtheB & 35 & 12 & 10 & 8 & 10 \\ No & >6m & BtheB & 21 & 6 & NA & NA & NA \\ No & <6m & TAU & 26 & 17 & 17 & 20 & 12 \\ No & >6m & TAU & 29 & 22 & 10 & NA & NA \\ No & >6m & TAU & 20 & 21 & NA & NA & NA \\ No & >6m & TAU & 33 & 23 & NA & NA & NA \\ No & >6m & BtheB & 19 & 12 & 13 & NA & NA \\ Yes & <6m & TAU & 12 & 15 & NA & NA & NA \\ Yes & >6m & TAU & 47 & 36 & 49 & 34 & NA \\ Yes & >6m & BtheB & 36 & 6 & 0 & 0 & 2 \\ No & <6m & BtheB & 10 & 8 & 6 & 3 & 3 \\ No & <6m & TAU & 27 & 7 & 15 & 16 & 0 \\ No & <6m & BtheB & 18 & 10 & 10 & 6 & 8 \\ Yes & <6m & BtheB & 11 & 8 & 3 & 2 & 15 \\ Yes & <6m & BtheB & 6 & 7 & NA & NA & NA \\ Yes & >6m & BtheB & 44 & 24 & 20 & 29 & 14 \\ No & <6m & TAU & 38 & 38 & NA & NA & NA \\ No & <6m & TAU & 21 & 14 & 20 & 1 & 8 \\ Yes & >6m & TAU & 34 & 17 & 8 & 9 & 13 \\ Yes & <6m & BtheB & 9 & 7 & 1 & NA & NA \\ Yes & >6m & TAU & 38 & 27 & 19 & 20 & 30 \\ Yes & <6m & BtheB & 46 & 40 & NA & NA & NA \\ No & <6m & TAU & 20 & 19 & 18 & 19 & 18 \\ Yes & >6m & TAU & 17 & 29 & 2 & 0 & 0 \\ No & >6m & BtheB & 18 & 20 & NA & NA & NA \\ Yes & >6m & BtheB & 42 & 1 & 8 & 10 & 6 \\ No & <6m & BtheB & 30 & 30 & NA & NA & NA \\ Yes & <6m & BtheB & 33 & 27 & 16 & 30 & 15 \\ No & <6m & BtheB & 12 & 1 & 0 & 0 & NA \\ Yes & <6m & BtheB & 2 & 5 & NA & NA & NA \\ No & >6m & TAU & 36 & 42 & 49 & 47 & 40 \\ No & <6m & TAU & 35 & 30 & NA & NA & NA \\ No & <6m & BtheB & 23 & 20 & NA & NA & NA \\ No & >6m & TAU & 31 & 48 & 38 & 38 & 37 \\ Yes & <6m & BtheB & 8 & 5 & 7 & NA & NA \\ Yes & <6m & TAU & 23 & 21 & 26 & NA & NA \\ Yes & <6m & BtheB & 7 & 7 & 5 & 4 & 0 \\ No & <6m & TAU & 14 & 13 & 14 & NA & NA \\ No & <6m & TAU & 40 & 36 & 33 & NA & NA \\ Yes & <6m & BtheB & 23 & 30 & NA & NA & NA \\ No & >6m & BtheB & 14 & 3 & NA & NA & NA \\ No & >6m & TAU & 22 & 20 & 16 & 24 & 16 \\ No & >6m & TAU & 23 & 23 & 15 & 25 & 17 \\ No & <6m & TAU & 15 & 7 & 13 & 13 & NA \\ No & >6m & TAU & 8 & 12 & 11 & 26 & NA \\ No & >6m & BtheB & 12 & 18 & NA & NA & NA \\ No & >6m & TAU & 7 & 6 & 2 & 1 & NA \\ Yes & <6m & TAU & 17 & 9 & 3 & 1 & 0 \\ Yes & <6m & BtheB & 33 & 18 & 16 & NA & NA \\ No & <6m & TAU & 27 & 20 & NA & NA & NA \\ No & <6m & BtheB & 27 & 30 & NA & NA & NA \\ No & <6m & BtheB & 9 & 6 & 10 & 1 & 0 \\ No & >6m & BtheB & 40 & 30 & 12 & NA & NA \\ No & >6m & TAU & 11 & 8 & 7 & NA & NA \\ No & <6m & TAU & 9 & 8 & NA & NA & NA \\ No & >6m & TAU & 14 & 22 & 21 & 24 & 19 \\ Yes & >6m & BtheB & 28 & 9 & 20 & 18 & 13 \\ No & >6m & BtheB & 15 & 9 & 13 & 14 & 10 \\ Yes & >6m & BtheB & 22 & 10 & 5 & 5 & 12 \\ No & <6m & TAU & 23 & 9 & NA & NA & NA \\ No & >6m & TAU & 21 & 22 & 24 & 23 & 22 \\ No & >6m & TAU & 27 & 31 & 28 & 22 & 14 \\ Yes & >6m & BtheB & 14 & 15 & NA & NA & NA \\ No & >6m & TAU & 10 & 13 & 12 & 8 & 20 \\ Yes & <6m & TAU & 21 & 9 & 6 & 7 & 1 \\ Yes & >6m & BtheB & 46 & 36 & 53 & NA & NA \\ No & >6m & BtheB & 36 & 14 & 7 & 15 & 15 \\ Yes & >6m & BtheB & 23 & 17 & NA & NA & NA \\ Yes & >6m & TAU & 35 & 0 & 6 & 0 & 1 \\ Yes & <6m & BtheB & 33 & 13 & 13 & 10 & 8 \\ No & <6m & BtheB & 19 & 4 & 27 & 1 & 2 \\ No & <6m & TAU & 16 & NA & NA & NA & NA \\ Yes & <6m & BtheB & 30 & 26 & 28 & NA & NA \\ Yes & <6m & BtheB & 17 & 8 & 7 & 12 & NA \\ No & >6m & BtheB & 19 & 4 & 3 & 3 & 3 \\ No & >6m & BtheB & 16 & 11 & 4 & 2 & 3 \\ Yes & >6m & BtheB & 16 & 16 & 10 & 10 & 8 \\ Yes & <6m & TAU & 28 & NA & NA & NA & NA \\ No & >6m & BtheB & 11 & 22 & 9 & 11 & 11 \\ No & <6m & TAU & 13 & 5 & 5 & 0 & 6 \\ Yes & <6m & TAU & 43 & NA & NA & NA & NA \\ \hline \end{longtable} \end{center} > ################################################### > ### code chunk number 33: ch:LME:BtheB:plot > ################################################### > ylim <- range(BtheB[,grep("bdi", names(BtheB))], + na.rm = TRUE) > tau <- subset(BtheB, treatment == "TAU")[, + grep("bdi", names(BtheB))] > boxplot(tau, main = "Treated as Usual", ylab = "BDI", + xlab = "Time (in months)", names = c(0, 2, 3, 5, 8), + ylim = ylim) > btheb <- subset(BtheB, treatment == "BtheB")[, + grep("bdi", names(BtheB))] > boxplot(btheb, main = "Beat the Blues", ylab = "BDI", + xlab = "Time (in months)", names = c(0, 2, 3, 5, 8), + ylim = ylim) > ################################################### > ### code chunk number 34: ch:LME:BtheB:long > ################################################### > BtheB$subject <- factor(rownames(BtheB)) > nobs <- nrow(BtheB) > BtheB_long <- reshape(BtheB, idvar = "subject", + varying = c("bdi.2m", "bdi.3m", "bdi.5m", "bdi.8m"), + direction = "long") > BtheB_long$time <- rep(c(2, 3, 5, 8), rep(nobs, 4)) > ################################################### > ### code chunk number 35: ch:LME:BtheB:lme > ################################################### > BtheB_lme1 <- lme(bdi ~ bdi.pre + time + treatment + drug + + length, random = ~ 1 | subject, data = BtheB_long, + na.action = na.omit) > BtheB_lme2 <- lme(bdi ~ bdi.pre + time + treatment + drug + + length, random = ~ time | subject, data = BtheB_long, + na.action = na.omit) > ################################################### > ### code chunk number 36: ch:LME:BtheB:aov > ################################################### > anova(BtheB_lme1, BtheB_lme2) Model df AIC BIC logLik Test L.Ratio p-value BtheB_lme1 1 8 1882.903 1911.808 -933.4513 BtheB_lme2 2 10 1886.336 1922.467 -933.1681 1 vs 2 0.566461 0.7533 > ################################################### > ### code chunk number 37: ch:LME:BtheB:summary > ################################################### > summary(BtheB_lme1) Linear mixed-effects model fit by REML Data: BtheB_long AIC BIC logLik 1882.903 1911.808 -933.4513 Random effects: Formula: ~1 | subject (Intercept) Residual StdDev: 7.206358 5.028555 Fixed effects: bdi ~ bdi.pre + time + treatment + drug + length Value Std.Error DF t-value p-value (Intercept) 5.573746 2.2995442 182 2.423848 0.0163 bdi.pre 0.640351 0.0799164 92 8.012758 0.0000 time -0.701598 0.1469407 182 -4.774702 0.0000 treatmentBtheB -2.315105 1.7151530 92 -1.349795 0.1804 drugYes -2.815994 1.7729198 92 -1.588337 0.1156 length>6m 0.179015 1.6816346 92 0.106453 0.9155 Correlation: (Intr) bdi.pr time trtmBB drugYs bdi.pre -0.683 time -0.232 0.019 treatmentBtheB -0.390 0.121 0.017 drugYes -0.074 -0.236 -0.022 -0.323 length>6m -0.244 -0.241 -0.036 0.002 0.157 Standardized Within-Group Residuals: Min Q1 Med Q3 Max -2.6869948 -0.5084732 -0.0608478 0.4206731 3.8141401 Number of Observations: 280 Number of Groups: 97 > ################################################### > ### code chunk number 38: ch:LME:BtheB-missing > ################################################### > bdi <- BtheB[, grep("bdi", names(BtheB))] > plot(1:4, rep(-0.5, 4), type = "n", axes = FALSE, + ylim = c(0, 50), xlab = "Months", ylab = "BDI") > axis(1, at = 1:4, labels = c(0, 2, 3, 5)) > axis(2) > for (i in 1:4) { + dropout <- is.na(bdi[,i + 1]) + points(rep(i, nrow(bdi)) + ifelse(dropout, 0.05, -0.05), + jitter(bdi[,i]), pch = ifelse(dropout, 20, 1)) + } demo(Ch-MDS) ---- ~~~~~~ > ### R code from vignette source 'Ch-MDS.Rnw' > > ################################################### > ### code chunk number 1: setup > ################################################### > library("MVA") > set.seed(280875) > library("lattice") > lattice.options(default.theme = + function() + standard.theme("pdf", color = FALSE)) > if (file.exists("deparse.R")) { + if (!file.exists("figs")) dir.create("figs") + source("deparse.R") + options(prompt = "R> ", continue = "+ ", width = 64, + digits = 4, show.signif.stars = FALSE, useFancyQuotes = FALSE) + + options(SweaveHooks = list(onefig = function() {par(mfrow = c(1,1))}, + twofig = function() {par(mfrow = c(1,2))}, + figtwo = function() {par(mfrow = c(2,1))}, + threefig = function() {par(mfrow = c(1,3))}, + figthree = function() {par(mfrow = c(3,1))}, + fourfig = function() {par(mfrow = c(2,2))}, + sixfig = function() {par(mfrow = c(3,2))}, + nomar = function() par("mai" = c(0, 0, 0, 0)))) + } > ################################################### > ### code chunk number 2: ch:MDS:data > ################################################### > teensex <- matrix(c(21, 8, 2, 21, 9, 3, 14, 6, 4, 13, 8, 10, 8, 2, 10), nrow = 3) > rownames(teensex) <- c("No boyfriend", "Boyfriend no sex", "Boyfriend sex") > colnames(teensex) <- c("<16", "16-17", "17-18", "18-19", "19-20") > teensex <- as.table(teensex, dnn = c("Boyfriend", "Age group")) > teensex <- as.data.frame(teensex) > names(teensex) <- c("Boyfriend", "Age", "Freq") > ################################################### > ### code chunk number 3: ch:MDS:X > ################################################### > X <- matrix(c( + 3, 5, 6, 1, 4, 2, 0, 0, 7, 2, + 4, 1, 2, 1, 7, 2, 4, 6, 6, 1, + 4, 1, 0, 1, 3, 5, 1, 4, 5, 4, + 6, 7, 2, 0, 6, 1, 1, 3, 1, 3, + 1, 3, 6, 3, 2, 0, 1, 5, 4, 1), nrow = 10) > X [,1] [,2] [,3] [,4] [,5] [1,] 3 4 4 6 1 [2,] 5 1 1 7 3 [3,] 6 2 0 2 6 [4,] 1 1 1 0 3 [5,] 4 7 3 6 2 [6,] 2 2 5 1 0 [7,] 0 4 1 1 1 [8,] 0 6 4 3 5 [9,] 7 6 5 1 4 [10,] 2 1 4 3 1 > ################################################### > ### code chunk number 4: ch:MDS:dist > ################################################### > (D <- dist(X)) 1 2 3 4 5 6 7 8 2 5.196152 3 8.366600 6.082763 4 7.874008 8.062258 6.324555 5 3.464102 6.557439 8.366600 9.273618 6 5.656854 8.426150 8.831761 5.291503 7.874008 7 6.557439 8.602325 8.185353 3.872983 7.416198 5.000000 8 6.164414 8.888194 8.366600 6.928203 6.000000 7.071068 5.744563 9 7.416198 9.055385 6.855655 8.888194 6.557439 7.549834 8.831761 7.416198 10 4.358899 6.164414 7.681146 4.795832 7.141428 2.645751 5.099020 6.708204 9 2 3 4 5 6 7 8 9 10 8.000000 > ################################################### > ### code chunk number 5: ch:MDA:cmdscale > ################################################### > cmdscale(D, k = 7, eig = TRUE) $points [,1] [,2] [,3] [,4] [,5] [,6] [1,] -1.6038325 2.38060903 -2.2301092 -0.3656856 0.11536476 0.000000e+00 [2,] -2.8246377 -2.30937202 -3.9523782 0.3419185 0.33169405 2.606139e-08 [3,] -1.6908272 -5.13970089 1.2880306 0.6503227 -0.05133897 -1.951486e-08 [4,] 3.9527719 -2.43233961 0.3833746 0.6863995 -0.03460933 7.420089e-08 [5,] -3.5984894 2.75538195 -0.2551393 1.0783741 -1.26125237 4.482845e-08 [6,] 2.9520356 1.35475175 -0.1899027 -2.8211220 0.12385813 1.377060e-09 [7,] 3.4689928 0.76411068 0.3016531 1.6369166 -1.94209512 -1.234005e-08 [8,] 0.3545235 2.31408566 2.2161772 2.9240116 2.00450379 8.602551e-09 [9,] -2.9362323 -0.01279597 4.3117385 -2.5122743 -0.18911558 2.878454e-08 [10,] 1.9256952 0.32526941 -1.8734445 -1.6188611 0.90299062 1.497845e-08 [,7] [1,] 4.391725e-08 [2,] -1.326889e-08 [3,] 3.308066e-09 [4,] 1.240466e-08 [5,] -1.767287e-08 [6,] -2.610658e-08 [7,] 7.034482e-09 [8,] -5.028996e-09 [9,] 1.163985e-08 [10,] 7.553758e-09 $eig [1] 7.518716e+01 5.880560e+01 4.960516e+01 3.042789e+01 1.037419e+01 [6] 9.856468e-15 3.530811e-15 2.864197e-15 -7.294914e-15 -1.050516e-14 $x NULL $ac [1] 0 $GOF [1] 1 1 > ################################################### > ### code chunk number 6: ch:MDA:distcmd > ################################################### > max(abs(dist(X) - dist(cmdscale(D, k = 5)))) [1] 1.24345e-14 > ################################################### > ### code chunk number 7: ch:MDA:cmdscalePCA > ################################################### > max(abs(prcomp(X)$x) - abs(cmdscale(D, k = 5))) [1] 2.257569e-14 > ################################################### > ### code chunk number 8: ch:MDA:man > ################################################### > X_m <- cmdscale(dist(X, method = "manhattan"), + k = 6, eig = TRUE) > ################################################### > ### code chunk number 9: ch:MDA:man:eigen > ################################################### > (X_eigen <- X_m$eig) [1] 2.806843e+02 2.494246e+02 2.288549e+02 9.250710e+01 4.250504e+01 [6] 2.196808e+01 -2.131628e-14 -1.507023e+01 -2.804630e+01 -5.682752e+01 > ################################################### > ### code chunk number 10: ch:MDA:man:eigen2 > ################################################### > cumsum(abs(X_eigen)) / sum(abs(X_eigen)) [1] 0.2762945 0.5218182 0.7470939 0.8381542 0.8799945 0.9016190 0.9016190 [8] 0.9164536 0.9440612 1.0000000 > cumsum(X_eigen^2) / sum(X_eigen^2) [1] 0.3779304 0.6763685 0.9276127 0.9686639 0.9773307 0.9796457 0.9796457 [8] 0.9807352 0.9845085 1.0000000 > ################################################### > ### code chunk number 11: ch:MDS:airdist:tab > ################################################### > "airline.dist" <- + structure(.Data = list(c(0, 587, 1212, 701, 1936, 604, 748, 2139, 2181, 543) + , c(587, 0, 920, 940, 1745, 1188, 713, 1858, 1737, 597) + , c(1212, 920, 0, 879, 831, 1726, 1631, 949, 1021, 1494) + , c(701, 940, 879, 0, 1374, 968, 1420, 1645, 1891, 1220) + , c(1936, 1745, 831, 1374, 0, 2339, 2451, 347, 959, 2300) + , c(604, 1188, 1726, 968, 2339, 0, 1092, 2594, 2734, 923) + , c(748, 713, 1631, 1420, 2451, 1092, 0, 2571, 2408, 205) + , c(2139, 1858, 949, 1645, 347, 2594, 2571, 0, 678, 2442) + , c(218, 1737, 1021, 1891, 959, 2734, 2408, 678, 0, 2329) + , c(543, 597, 1494, 1220, 2300, 923, 205, 2442, 2329, 0) + ) + , names = c("ATL", "ORD", "DEN", "HOU", "LAX", "MIA", + "JFK", "SFO", "SEA", "IAD") + , row.names = c("ATL", "ORD", "DEN", "HOU", "LAX", "MIA", + "JFK", "SFO", "SEA", "IAD") + , class = "data.frame" + ) > airdist <- as.dist(as.matrix(airline.dist)) > tmp <- airdist > airdist <- as.data.frame(as.matrix(airdist)) > htab <- HSAURtable(airdist) > htab$data[upper.tri(htab$data)] <- " " > toLatex(htab, pcol = 1, xname = "airdist", + caption = "Airline distances between ten US cities.", + label = "ch:MDS:airdist:tab", rownames = TRUE) \index{airdist data@\Robject{airdist} data} \begin{center} \begin{longtable}{l rrrrrrrrrr } \caption{\Robject{airdist} data. Airline distances between ten US cities. \label{ch:MDS:airdist:tab}} \\ \hline & \Robject{ATL} & \Robject{ORD} & \Robject{DEN} & \Robject{HOU} & \Robject{LAX} & \Robject{MIA} & \Robject{JFK} & \Robject{SFO} & \Robject{SEA} & \Robject{IAD} \\ \hline \endfirsthead \caption[]{\Robject{airdist} data (continued).} \\ \hline & \Robject{ATL} & \Robject{ORD} & \Robject{DEN} & \Robject{HOU} & \Robject{LAX} & \Robject{MIA} & \Robject{JFK} & \Robject{SFO} & \Robject{SEA} & \Robject{IAD} \\ \hline \endhead ATL & 0 & & & & & & & & & \\ ORD & 587 & 0 & & & & & & & & \\ DEN & 1212 & 920 & 0 & & & & & & & \\ HOU & 701 & 940 & 879 & 0 & & & & & & \\ LAX & 1936 & 1745 & 831 & 1374 & 0 & & & & & \\ MIA & 604 & 1188 & 1726 & 968 & 2339 & 0 & & & & \\ JFK & 748 & 713 & 1631 & 1420 & 2451 & 1092 & 0 & & & \\ SFO & 2139 & 1858 & 949 & 1645 & 347 & 2594 & 2571 & 0 & & \\ SEA & 2181 & 1737 & 1021 & 1891 & 959 & 2734 & 2408 & 678 & 0 & \\ IAD & 543 & 597 & 1494 & 1220 & 2300 & 923 & 205 & 2442 & 2329 & 0 \\ \hline \end{longtable} \end{center} > airdist <- tmp > ################################################### > ### code chunk number 12: ch:MDS:airdist > ################################################### > airline_mds <- cmdscale(airdist, k = 6, eig = TRUE) > airline_mds$points [,1] [,2] [,3] [,4] [,5] [,6] ATL -718.4820 142.38132 -31.268404 -5.2514775 12.9215244 3.8940163 ORD -382.0846 -340.76745 -29.097077 8.5382391 9.0604048 -8.9375235 DEN 481.5988 -25.27903 -53.415844 4.2667700 -17.3575193 1.0892656 HOU -161.4740 572.87820 -2.589168 0.8043809 4.2508530 5.8408682 LAX 1203.7581 390.06234 18.236750 -15.6280052 -0.4969341 -4.3899960 MIA -1133.5570 582.11397 30.781556 4.9344072 -6.3160793 -4.5290408 JFK -1072.2837 -518.89299 34.399677 14.1980523 -3.4180108 3.0991313 SFO 1420.6236 112.51210 9.006876 17.9998946 3.2582207 1.5811509 SEA 1341.5654 -579.66890 22.493733 -7.5958487 2.1852580 1.4359808 IAD -979.6646 -335.33957 1.451902 -22.2664128 -4.0877174 0.9161472 > ################################################### > ### code chunk number 13: ch:MDS:airdistprint > ################################################### > airline_mds$points[, -9] [,1] [,2] [,3] [,4] [,5] [,6] ATL -718.4820 142.38132 -31.268404 -5.2514775 12.9215244 3.8940163 ORD -382.0846 -340.76745 -29.097077 8.5382391 9.0604048 -8.9375235 DEN 481.5988 -25.27903 -53.415844 4.2667700 -17.3575193 1.0892656 HOU -161.4740 572.87820 -2.589168 0.8043809 4.2508530 5.8408682 LAX 1203.7581 390.06234 18.236750 -15.6280052 -0.4969341 -4.3899960 MIA -1133.5570 582.11397 30.781556 4.9344072 -6.3160793 -4.5290408 JFK -1072.2837 -518.89299 34.399677 14.1980523 -3.4180108 3.0991313 SFO 1420.6236 112.51210 9.006876 17.9998946 3.2582207 1.5811509 SEA 1341.5654 -579.66890 22.493733 -7.5958487 2.1852580 1.4359808 IAD -979.6646 -335.33957 1.451902 -22.2664128 -4.0877174 0.9161472 > ################################################### > ### code chunk number 14: ch:MDS:airdist:eigen > ################################################### > (lam <- airline_mds$eig) [1] 9.581705e+06 1.686606e+06 7.736930e+03 1.466986e+03 6.523333e+02 [6] 1.851352e+02 -6.621121e-10 -6.669100e+02 -5.630189e+03 -3.524780e+04 > ################################################### > ### code chunk number 15: ch:MDS:airdist:crit > ################################################### > cumsum(abs(lam)) / sum(abs(lam)) [1] 0.8464480 0.9954429 0.9961263 0.9962559 0.9963136 0.9963299 0.9963299 [8] 0.9963888 0.9968862 1.0000000 > cumsum(lam^2) / sum(lam^2) [1] 0.9699332 0.9999859 0.9999865 0.9999865 0.9999865 0.9999865 0.9999865 [8] 0.9999865 0.9999869 1.0000000 > ################################################### > ### code chunk number 16: ch:MDS:airdist:plot > ################################################### > lim <- range(airline_mds$points[,1] * (-1)) * 1.2 > plot(airline_mds$points[,1] * (-1), airline_mds$points[,2] * (-1), + type = "n", xlab = "Coordinate 1", ylab = "Coordinate 2", + xlim = lim, ylim = lim) > text(airline_mds$points[,1] *(-1), airline_mds$points[,2] * (-1), + labels(airdist), cex = 0.7) > ################################################### > ### code chunk number 17: ch:MDA:skulls:tab > ################################################### > data("skulls", package = "HSAUR2") > toLatex(HSAURtable(skulls), pcol = 3, + caption = "Measurements of four variables taken from Egyptian skulls of five periods.", + label = "ch:MDA:skulls:tab") \index{skulls data@\Robject{skulls} data} \begin{center} \begin{longtable}{ rrrrr|rrrrr|rrrrr } \caption{\Robject{skulls} data. Measurements of four variables taken from Egyptian skulls of five periods. \label{ch:MDA:skulls:tab}} \\ \hline \Robject{epoch} & \Robject{mb} & \Robject{bh} & \Robject{bl} & \Robject{nh} & \Robject{epoch} & \Robject{mb} & \Robject{bh} & \Robject{bl} & \Robject{nh} & \Robject{epoch} & \Robject{mb} & \Robject{bh} & \Robject{bl} & \Robject{nh} \\ \hline \endfirsthead \caption[]{\Robject{skulls} data (continued).} \\ \hline \Robject{epoch} & \Robject{mb} & \Robject{bh} & \Robject{bl} & \Robject{nh} & \Robject{epoch} & \Robject{mb} & \Robject{bh} & \Robject{bl} & \Robject{nh} & \Robject{epoch} & \Robject{mb} & \Robject{bh} & \Robject{bl} & \Robject{nh} \\ \hline \endhead c4000BC & 131 & 138 & 89 & 49 & c3300BC & 137 & 136 & 106 & 49 & c200BC & 132 & 133 & 90 & 53 \\ c4000BC & 125 & 131 & 92 & 48 & c3300BC & 126 & 131 & 100 & 48 & c200BC & 134 & 134 & 97 & 54 \\ c4000BC & 131 & 132 & 99 & 50 & c3300BC & 135 & 136 & 97 & 52 & c200BC & 135 & 135 & 99 & 50 \\ c4000BC & 119 & 132 & 96 & 44 & c3300BC & 129 & 126 & 91 & 50 & c200BC & 133 & 136 & 95 & 52 \\ c4000BC & 136 & 143 & 100 & 54 & c3300BC & 134 & 139 & 101 & 49 & c200BC & 136 & 130 & 99 & 55 \\ c4000BC & 138 & 137 & 89 & 56 & c3300BC & 131 & 134 & 90 & 53 & c200BC & 134 & 137 & 93 & 52 \\ c4000BC & 139 & 130 & 108 & 48 & c3300BC & 132 & 130 & 104 & 50 & c200BC & 131 & 141 & 99 & 55 \\ c4000BC & 125 & 136 & 93 & 48 & c3300BC & 130 & 132 & 93 & 52 & c200BC & 129 & 135 & 95 & 47 \\ c4000BC & 131 & 134 & 102 & 51 & c3300BC & 135 & 132 & 98 & 54 & c200BC & 136 & 128 & 93 & 54 \\ c4000BC & 134 & 134 & 99 & 51 & c3300BC & 130 & 128 & 101 & 51 & c200BC & 131 & 125 & 88 & 48 \\ c4000BC & 129 & 138 & 95 & 50 & c1850BC & 137 & 141 & 96 & 52 & c200BC & 139 & 130 & 94 & 53 \\ c4000BC & 134 & 121 & 95 & 53 & c1850BC & 129 & 133 & 93 & 47 & c200BC & 144 & 124 & 86 & 50 \\ c4000BC & 126 & 129 & 109 & 51 & c1850BC & 132 & 138 & 87 & 48 & c200BC & 141 & 131 & 97 & 53 \\ c4000BC & 132 & 136 & 100 & 50 & c1850BC & 130 & 134 & 106 & 50 & c200BC & 130 & 131 & 98 & 53 \\ c4000BC & 141 & 140 & 100 & 51 & c1850BC & 134 & 134 & 96 & 45 & c200BC & 133 & 128 & 92 & 51 \\ c4000BC & 131 & 134 & 97 & 54 & c1850BC & 140 & 133 & 98 & 50 & c200BC & 138 & 126 & 97 & 54 \\ c4000BC & 135 & 137 & 103 & 50 & c1850BC & 138 & 138 & 95 & 47 & c200BC & 131 & 142 & 95 & 53 \\ c4000BC & 132 & 133 & 93 & 53 & c1850BC & 136 & 145 & 99 & 55 & c200BC & 136 & 138 & 94 & 55 \\ c4000BC & 139 & 136 & 96 & 50 & c1850BC & 136 & 131 & 92 & 46 & c200BC & 132 & 136 & 92 & 52 \\ c4000BC & 132 & 131 & 101 & 49 & c1850BC & 126 & 136 & 95 & 56 & c200BC & 135 & 130 & 100 & 51 \\ c4000BC & 126 & 133 & 102 & 51 & c1850BC & 137 & 129 & 100 & 53 & cAD150 & 137 & 123 & 91 & 50 \\ c4000BC & 135 & 135 & 103 & 47 & c1850BC & 137 & 139 & 97 & 50 & cAD150 & 136 & 131 & 95 & 49 \\ c4000BC & 134 & 124 & 93 & 53 & c1850BC & 136 & 126 & 101 & 50 & cAD150 & 128 & 126 & 91 & 57 \\ c4000BC & 128 & 134 & 103 & 50 & c1850BC & 137 & 133 & 90 & 49 & cAD150 & 130 & 134 & 92 & 52 \\ c4000BC & 130 & 130 & 104 & 49 & c1850BC & 129 & 142 & 104 & 47 & cAD150 & 138 & 127 & 86 & 47 \\ c4000BC & 138 & 135 & 100 & 55 & c1850BC & 135 & 138 & 102 & 55 & cAD150 & 126 & 138 & 101 & 52 \\ c4000BC & 128 & 132 & 93 & 53 & c1850BC & 129 & 135 & 92 & 50 & cAD150 & 136 & 138 & 97 & 58 \\ c4000BC & 127 & 129 & 106 & 48 & c1850BC & 134 & 125 & 90 & 60 & cAD150 & 126 & 126 & 92 & 45 \\ c4000BC & 131 & 136 & 114 & 54 & c1850BC & 138 & 134 & 96 & 51 & cAD150 & 132 & 132 & 99 & 55 \\ c4000BC & 124 & 138 & 101 & 46 & c1850BC & 136 & 135 & 94 & 53 & cAD150 & 139 & 135 & 92 & 54 \\ c3300BC & 124 & 138 & 101 & 48 & c1850BC & 132 & 130 & 91 & 52 & cAD150 & 143 & 120 & 95 & 51 \\ c3300BC & 133 & 134 & 97 & 48 & c1850BC & 133 & 131 & 100 & 50 & cAD150 & 141 & 136 & 101 & 54 \\ c3300BC & 138 & 134 & 98 & 45 & c1850BC & 138 & 137 & 94 & 51 & cAD150 & 135 & 135 & 95 & 56 \\ c3300BC & 148 & 129 & 104 & 51 & c1850BC & 130 & 127 & 99 & 45 & cAD150 & 137 & 134 & 93 & 53 \\ c3300BC & 126 & 124 & 95 & 45 & c1850BC & 136 & 133 & 91 & 49 & cAD150 & 142 & 135 & 96 & 52 \\ c3300BC & 135 & 136 & 98 & 52 & c1850BC & 134 & 123 & 95 & 52 & cAD150 & 139 & 134 & 95 & 47 \\ c3300BC & 132 & 145 & 100 & 54 & c1850BC & 136 & 137 & 101 & 54 & cAD150 & 138 & 125 & 99 & 51 \\ c3300BC & 133 & 130 & 102 & 48 & c1850BC & 133 & 131 & 96 & 49 & cAD150 & 137 & 135 & 96 & 54 \\ c3300BC & 131 & 134 & 96 & 50 & c1850BC & 138 & 133 & 100 & 55 & cAD150 & 133 & 125 & 92 & 50 \\ c3300BC & 133 & 125 & 94 & 46 & c1850BC & 138 & 133 & 91 & 46 & cAD150 & 145 & 129 & 89 & 47 \\ c3300BC & 133 & 136 & 103 & 53 & c200BC & 137 & 134 & 107 & 54 & cAD150 & 138 & 136 & 92 & 46 \\ c3300BC & 131 & 139 & 98 & 51 & c200BC & 141 & 128 & 95 & 53 & cAD150 & 131 & 129 & 97 & 44 \\ c3300BC & 131 & 136 & 99 & 56 & c200BC & 141 & 130 & 87 & 49 & cAD150 & 143 & 126 & 88 & 54 \\ c3300BC & 138 & 134 & 98 & 49 & c200BC & 135 & 131 & 99 & 51 & cAD150 & 134 & 124 & 91 & 55 \\ c3300BC & 130 & 136 & 104 & 53 & c200BC & 133 & 120 & 91 & 46 & cAD150 & 132 & 127 & 97 & 52 \\ c3300BC & 131 & 128 & 98 & 45 & c200BC & 131 & 135 & 90 & 50 & cAD150 & 137 & 125 & 85 & 57 \\ c3300BC & 138 & 129 & 107 & 53 & c200BC & 140 & 137 & 94 & 60 & cAD150 & 129 & 128 & 81 & 52 \\ c3300BC & 123 & 131 & 101 & 51 & c200BC & 139 & 130 & 90 & 48 & cAD150 & 140 & 135 & 103 & 48 \\ c3300BC & 130 & 129 & 105 & 47 & c200BC & 140 & 134 & 90 & 51 & cAD150 & 147 & 129 & 87 & 48 \\ c3300BC & 134 & 130 & 93 & 54 & c200BC & 138 & 140 & 100 & 52 & cAD150 & 136 & 133 & 97 & 51 \\ \hline \end{longtable} \end{center} > ################################################### > ### code chunk number 18: ch:MDS:skulls:maha > ################################################### > skulls_var <- tapply(1:nrow(skulls), skulls$epoch, + function(i) var(skulls[i,-1])) > S <- 0 > for (v in skulls_var) S <- S + 29 * v > (S <- S / 149) mb bh bl nh mb 20.54407159 0.03579418 0.07695749 1.955034 bh 0.03579418 22.85413870 5.06040268 2.768680 bl 0.07695749 5.06040268 23.52997763 1.102908 nh 1.95503356 2.76868009 1.10290828 9.880089 > skulls_cen <- tapply(1:nrow(skulls), skulls$epoch, + function(i) apply(skulls[i,-1], 2, mean)) > skulls_cen <- matrix(unlist(skulls_cen), + nrow = length(skulls_cen), byrow = TRUE) > skulls_mah <- apply(skulls_cen, 1, + function(cen) mahalanobis(skulls_cen, cen, S)) > skulls_mah [,1] [,2] [,3] [,4] [,5] [1,] 0.00000000 0.09354553 0.9279862 1.9330193 2.7712116 [2,] 0.09354553 0.00000000 0.7490463 1.6379863 2.2357084 [3,] 0.92798621 0.74904633 0.0000000 0.4553368 0.9359991 [4,] 1.93301928 1.63798631 0.4553368 0.0000000 0.2253343 [5,] 2.77121156 2.23570836 0.9359991 0.2253343 0.0000000 > cmdscale(skulls_mah, k = 2, eig = TRUE)$eig [1] 5.112951e+00 9.816355e-02 -5.551115e-16 -1.008636e-01 -7.777015e-01 > skulls_mds <- cmdscale(skulls_mah) > ################################################### > ### code chunk number 19: ch:MDS:skulls:plot > ################################################### > lim <- range(skulls_mds) * 1.2 > plot(skulls_mds, xlab = "Coordinate 1", ylab = "Coordinate 2", + xlim = lim, ylim = lim, type = "n") > text(skulls_mds, labels = levels(skulls$epoch), cex = 0.7) > ################################################### > ### code chunk number 20: MDS-watervoles-tab > ################################################### > data("watervoles", package = "HSAUR2") > tmp <- watervoles > colnames(watervoles) <- abbreviate(colnames(watervoles)) > watervoles <- as.data.frame(watervoles) > htab <- HSAURtable(watervoles) > htab$data[upper.tri(htab$data)] <- " " > toLatex(htab, pcol = 1, xname = "watervoles", + caption = "Water voles data-dissimilarity matrix.", + label = "MDS-watervoles-tab", + rownames = TRUE) \index{watervoles data@\Robject{watervoles} data} \begin{center} \begin{longtable}{l rrrrrrrrrrrrrr } \caption{\Robject{watervoles} data. Water voles data-dissimilarity matrix. \label{MDS-watervoles-tab}} \\ \hline & \Robject{Srry} & \Robject{Shrp} & \Robject{Yrks} & \Robject{Prth} & \Robject{Abrd} & \Robject{ElnG} & \Robject{Alps} & \Robject{Ygsl} & \Robject{Grmn} & \Robject{Nrwy} & \Robject{PyrI} & \Robject{PyII} & \Robject{NrtS} & \Robject{SthS} \\ \hline \endfirsthead \caption[]{\Robject{watervoles} data (continued).} \\ \hline & \Robject{Srry} & \Robject{Shrp} & \Robject{Yrks} & \Robject{Prth} & \Robject{Abrd} & \Robject{ElnG} & \Robject{Alps} & \Robject{Ygsl} & \Robject{Grmn} & \Robject{Nrwy} & \Robject{PyrI} & \Robject{PyII} & \Robject{NrtS} & \Robject{SthS} \\ \hline \endhead Surrey & 0.000 & & & & & & & & & & & & & \\ Shropshire & 0.099 & 0.000 & & & & & & & & & & & & \\ Yorkshire & 0.033 & 0.022 & 0.000 & & & & & & & & & & & \\ Perthshire & 0.183 & 0.114 & 0.042 & 0.000 & & & & & & & & & & \\ Aberdeen & 0.148 & 0.224 & 0.059 & 0.068 & 0.000 & & & & & & & & & \\ Elean Gamhna & 0.198 & 0.039 & 0.053 & 0.085 & 0.051 & 0.000 & & & & & & & & \\ Alps & 0.462 & 0.266 & 0.322 & 0.435 & 0.268 & 0.025 & 0.000 & & & & & & & \\ Yugoslavia & 0.628 & 0.442 & 0.444 & 0.406 & 0.240 & 0.129 & 0.014 & 0.000 & & & & & & \\ Germany & 0.113 & 0.070 & 0.046 & 0.047 & 0.034 & 0.002 & 0.106 & 0.129 & 0.000 & & & & & \\ Norway & 0.173 & 0.119 & 0.162 & 0.331 & 0.177 & 0.039 & 0.089 & 0.237 & 0.071 & 0.000 & & & & \\ Pyrenees I & 0.434 & 0.419 & 0.339 & 0.505 & 0.469 & 0.390 & 0.315 & 0.349 & 0.151 & 0.430 & 0.000 & & & \\ Pyrenees II & 0.762 & 0.633 & 0.781 & 0.700 & 0.758 & 0.625 & 0.469 & 0.618 & 0.440 & 0.538 & 0.607 & 0.000 & & \\ North Spain & 0.530 & 0.389 & 0.482 & 0.579 & 0.597 & 0.498 & 0.374 & 0.562 & 0.247 & 0.383 & 0.387 & 0.084 & 0.000 & \\ South Spain & 0.586 & 0.435 & 0.550 & 0.530 & 0.552 & 0.509 & 0.369 & 0.471 & 0.234 & 0.346 & 0.456 & 0.090 & 0.038 & 0.000 \\ \hline \end{longtable} \end{center} > watervoles <- tmp > ################################################### > ### code chunk number 21: MDS-voles-cmdscale > ################################################### > data("watervoles", package = "HSAUR2") > voles_mds <- cmdscale(watervoles, k = 7, eig = TRUE) > voles_mds$eig [1] 7.359910e-01 2.626003e-01 1.492622e-01 6.990457e-02 2.956972e-02 [6] 1.931184e-02 1.387779e-16 -1.139451e-02 -1.279569e-02 -2.849924e-02 [11] -4.251502e-02 -5.255450e-02 -7.406143e-02 -1.097833e-01 > ################################################### > ### code chunk number 22: MDS-voles-criterion1 > ################################################### > cumsum(abs(voles_mds$eig))/sum(abs(voles_mds$eig)) [1] 0.4605000 0.6248056 0.7181970 0.7619353 0.7804367 0.7925199 0.7925199 [8] 0.7996493 0.8076554 0.8254870 0.8520881 0.8849708 0.9313100 1.0000000 > ################################################### > ### code chunk number 23: MDS-voles-criterion2 > ################################################### > cumsum((voles_mds$eig)^2)/sum((voles_mds$eig)^2) [1] 0.8179213 0.9220468 0.9556875 0.9630662 0.9643865 0.9649496 0.9649496 [8] 0.9651457 0.9653929 0.9666193 0.9693486 0.9735191 0.9818014 1.0000000 > ################################################### > ### code chunk number 24: MDS-watervoles-plot > ################################################### > x <- voles_mds$points[,1] > y <- voles_mds$points[,2] > plot(x, y, xlab = "Coordinate 1", ylab = "Coordinate 2", + xlim = range(x)*1.2, type = "n") > text(x, y, labels = colnames(watervoles), cex = 0.7) > ################################################### > ### code chunk number 25: MDS-watervoles-mst > ################################################### > library("ape") > st <- mst(watervoles) > plot(x, y, xlab = "Coordinate 1", ylab = "Coordinate 2", + xlim = range(x)*1.2, type = "n") > for (i in 1:nrow(watervoles)) { + w1 <- which(st[i, ] == 1) + segments(x[i], y[i], x[w1], y[w1]) + } > text(x, y, labels = colnames(watervoles), cex = 0.7) > ################################################### > ### code chunk number 26: MDS-voting > ################################################### > library("MASS") > data("voting", package = "HSAUR2") > voting_mds <- isoMDS(voting) initial value 15.268246 iter 5 value 10.264075 final value 9.879047 converged > ################################################### > ### code chunk number 27: MDS-voting-tab > ################################################### > data("voting", package = "HSAUR2") > tmp <- voting > tmpname <- gsub("\\(.*", "", colnames(voting)) > colnames(voting) <- abbreviate(tmpname, 3) > voting <- as.data.frame(voting) > htab <- HSAURtable(voting) > htab$data[upper.tri(htab$data)] <- " " > toLatex(htab, pcol = 1, xname = "voting", + caption = paste("House of Representatives voting data;", + "(R) is short for Republican, (D) for Democrat."), + label = "MDS-voting-tab", + rownames = TRUE) \index{voting data@\Robject{voting} data} \begin{center} \begin{longtable}{l rrrrrrrrrrrrrrr } \caption{\Robject{voting} data. House of Representatives voting data; (R) is short for Republican, (D) for Democrat. \label{MDS-voting-tab}} \\ \hline & \Robject{Hnt} & \Robject{Snd} & \Robject{Hwr} & \Robject{Thm} & \Robject{Fry} & \Robject{Frs} & \Robject{Wdn} & \Robject{Roe} & \Robject{Hlt} & \Robject{Rdn} & \Robject{Mns} & \Robject{Rnl} & \Robject{Mrz} & \Robject{Dnl} & \Robject{Ptt} \\ \hline \endfirsthead \caption[]{\Robject{voting} data (continued).} \\ \hline & \Robject{Hnt} & \Robject{Snd} & \Robject{Hwr} & \Robject{Thm} & \Robject{Fry} & \Robject{Frs} & \Robject{Wdn} & \Robject{Roe} & \Robject{Hlt} & \Robject{Rdn} & \Robject{Mns} & \Robject{Rnl} & \Robject{Mrz} & \Robject{Dnl} & \Robject{Ptt} \\ \hline \endhead Hunt(R) & 0 & & & & & & & & & & & & & & \\ Sandman(R) & 8 & 0 & & & & & & & & & & & & & \\ Howard(D) & 15 & 17 & 0 & & & & & & & & & & & & \\ Thompson(D) & 15 & 12 & 9 & 0 & & & & & & & & & & & \\ Freylinghuysen(R) & 10 & 13 & 16 & 14 & 0 & & & & & & & & & & \\ Forsythe(R) & 9 & 13 & 12 & 12 & 8 & 0 & & & & & & & & & \\ Widnall(R) & 7 & 12 & 15 & 13 & 9 & 7 & 0 & & & & & & & & \\ Roe(D) & 15 & 16 & 5 & 10 & 13 & 12 & 17 & 0 & & & & & & & \\ Heltoski(D) & 16 & 17 & 5 & 8 & 14 & 11 & 16 & 4 & 0 & & & & & & \\ Rodino(D) & 14 & 15 & 6 & 8 & 12 & 10 & 15 & 5 & 3 & 0 & & & & & \\ Minish(D) & 15 & 16 & 5 & 8 & 12 & 9 & 14 & 5 & 2 & 1 & 0 & & & & \\ Rinaldo(R) & 16 & 17 & 4 & 6 & 12 & 10 & 15 & 3 & 1 & 2 & 1 & 0 & & & \\ Maraziti(R) & 7 & 13 & 11 & 15 & 10 & 6 & 10 & 12 & 13 & 11 & 12 & 12 & 0 & & \\ Daniels(D) & 11 & 12 & 10 & 10 & 11 & 6 & 11 & 7 & 7 & 4 & 5 & 6 & 9 & 0 & \\ Patten(D) & 13 & 16 & 7 & 7 & 11 & 10 & 13 & 6 & 5 & 6 & 5 & 4 & 13 & 9 & 0 \\ \hline \end{longtable} \end{center} > voting <- tmp > ################################################### > ### code chunk number 28: MDS-voting-plot > ################################################### > x <- voting_mds$points[,1] > y <- voting_mds$points[,2] > plot(x, y, xlab = "Coordinate 1", ylab = "Coordinate 2", + xlim = range(voting_mds$points[,1])*1.2, type = "n") > text(x, y, labels = colnames(voting), cex = 0.6) > voting_sh <- Shepard(voting[lower.tri(voting)], + voting_mds$points) > ################################################### > ### code chunk number 29: MDS-voting-Shepard > ################################################### > plot(voting_sh, pch = ".", xlab = "Dissimilarity", + ylab = "Distance", xlim = range(voting_sh$x), + ylim = range(voting_sh$x)) > lines(voting_sh$x, voting_sh$yf, type = "S") > ################################################### > ### code chunk number 30: ch:MDS:WWIIleaders:tab > ################################################### > WWIIleaders <- c( + 3, + 4, 6, + 7, 8, 4, + 3, 5, 6, 8, + 8, 9, 3, 9, 8, + 3, 2, 5, 7, 6, 7, + 4, 4, 3, 5, 6, 5, 4, + 8, 9, 8, 9, 6, 9, 8, 7, + 9, 9, 5, 4, 7, 8, 8, 4, 4, + 4, 5, 5, 4, 7, 2, 2, 5, 9, 5, + 7, 8, 2, 4, 7, 8, 3, 2, 4, 5, 7) > tmp <- matrix(0, ncol = 12, nrow = 12) > tmp[upper.tri(tmp)] <- WWIIleaders > tmp <- tmp + t(tmp) > rownames(tmp) <- colnames(tmp) <- c("Hitler", "Mussolini", "Churchill", + "Eisenhower", "Stalin", "Attlee", "Franco", "De Gaulle", "Mao Tse-Tung", + "Truman", "Chamberlin", "Tito") > WWIIleaders <- as.dist(tmp) > tmp <- WWIIleaders > WWIIleaders <- as.data.frame(as.matrix(WWIIleaders)) > colnames(WWIIleaders) <- abbreviate(colnames(WWIIleaders), 3) > htab <- HSAURtable(WWIIleaders) > htab$data[upper.tri(htab$data)] <- " " > toLatex(htab, pcol = 1, xname = "WWIIleaders", + caption = "Subjective distances between WWII leaders.", + label = "ch:MDS:WWIIleaders:tab", rownames = TRUE) \index{WWIIleaders data@\Robject{WWIIleaders} data} \begin{center} \begin{longtable}{l rrrrrrrrrrrr } \caption{\Robject{WWIIleaders} data. Subjective distances between WWII leaders. \label{ch:MDS:WWIIleaders:tab}} \\ \hline & \Robject{Htl} & \Robject{Mss} & \Robject{Chr} & \Robject{Esn} & \Robject{Stl} & \Robject{Att} & \Robject{Frn} & \Robject{DGl} & \Robject{MT-} & \Robject{Trm} & \Robject{Chm} & \Robject{Tit} \\ \hline \endfirsthead \caption[]{\Robject{WWIIleaders} data (continued).} \\ \hline & \Robject{Htl} & \Robject{Mss} & \Robject{Chr} & \Robject{Esn} & \Robject{Stl} & \Robject{Att} & \Robject{Frn} & \Robject{DGl} & \Robject{MT-} & \Robject{Trm} & \Robject{Chm} & \Robject{Tit} \\ \hline \endhead Hitler & 0 & & & & & & & & & & & \\ Mussolini & 3 & 0 & & & & & & & & & & \\ Churchill & 4 & 6 & 0 & & & & & & & & & \\ Eisenhower & 7 & 8 & 4 & 0 & & & & & & & & \\ Stalin & 3 & 5 & 6 & 8 & 0 & & & & & & & \\ Attlee & 8 & 9 & 3 & 9 & 8 & 0 & & & & & & \\ Franco & 3 & 2 & 5 & 7 & 6 & 7 & 0 & & & & & \\ De Gaulle & 4 & 4 & 3 & 5 & 6 & 5 & 4 & 0 & & & & \\ Mao Tse-Tung & 8 & 9 & 8 & 9 & 6 & 9 & 8 & 7 & 0 & & & \\ Truman & 9 & 9 & 5 & 4 & 7 & 8 & 8 & 4 & 4 & 0 & & \\ Chamberlin & 4 & 5 & 5 & 4 & 7 & 2 & 2 & 5 & 9 & 5 & 0 & \\ Tito & 7 & 8 & 2 & 4 & 7 & 8 & 3 & 2 & 4 & 5 & 7 & 0 \\ \hline \end{longtable} \end{center} > WWIIleaders <- tmp > ################################################### > ### code chunk number 31: ch:MDA:WWIIleadersMDS > ################################################### > (WWII_mds <- isoMDS(WWIIleaders)) initial value 20.504211 iter 5 value 15.216103 iter 5 value 15.207237 iter 5 value 15.207237 final value 15.207237 converged $points [,1] [,2] Hitler -2.5820321 -1.75960714 Mussolini -3.8806920 -1.24755866 Churchill 0.3109807 1.46671155 Eisenhower 2.9852347 2.87821624 Stalin -1.4273957 -3.75699052 Attlee -2.1067050 5.07317056 Franco -2.8589747 0.07877559 De Gaulle 0.6590874 -0.20655989 Mao Tse-Tung 4.1604539 -4.57583381 Truman 4.4961515 0.29294331 Chamberlin -2.1419835 2.75876703 Tito 2.3858746 -1.00203426 $stress [1] 15.20724 > ################################################### > ### code chunk number 32: ch:MDS:WWIIleaders:plot > ################################################### > x <- WWII_mds$points[,1] > y <- WWII_mds$points[,2] > lim <- range(c(x, y)) * 1.2 > plot(x, y, xlab = "Coordinate 1", ylab = "Coordinate 2", + xlim = lim, ylim = lim, type = "n") > text(x, y, labels = labels(WWIIleaders), cex = 0.7) > ################################################### > ### code chunk number 33: ch:MDS:teensex:tab > ################################################### > teensex <- xtabs(Freq ~ Boyfriend + Age, data = teensex) > toLatex(HSAURtable(teensex), pcol = 1, + caption = "The influence of age on relationships with boyfriends.", + label = "ch:MDS:teensex:tab", rownames = FALSE) \index{teensex data@\Robject{teensex} data} \begin{center} \begin{longtable} { rrrrrrr } \caption{\Robject{teensex} data. The influence of age on relationships with boyfriends. \label{ch:MDS:teensex:tab}} \\ & & \multicolumn{ 5 }{c}{\Robject{ Age }} \\ \Robject{ Boyfriend } & & <16 & 16-17 & 17-18 & 18-19 & 19-20 \\ & No boyfriend & 21 & 21 & 14 & 13 & 8 \\ & Boyfriend no sex & 8 & 9 & 6 & 8 & 2 \\ & Boyfriend sex & 2 & 3 & 4 & 10 & 10 \\ \end{longtable} \end{center} > ################################################### > ### code chunk number 34: ch:MDS:teensex:CA > ################################################### > D <- function(x) { + a <- t(t(x) / colSums(x)) + ret <- sqrt(colSums((a[,rep(1:ncol(x), ncol(x))] - + a[, rep(1:ncol(x), rep(ncol(x), ncol(x)))])^2 * + sum(x) / rowSums(x))) + matrix(ret, ncol = ncol(x)) + } > (dcols <- D(teensex)) [,1] [,2] [,3] [,4] [,5] [1,] 0.00000000 0.08536508 0.2574269 0.6628943 1.0738536 [2,] 0.08536508 0.00000000 0.1864428 0.5858083 1.0142343 [3,] 0.25742689 0.18644276 0.0000000 0.4066080 0.8294663 [4,] 0.66289427 0.58580829 0.4066080 0.0000000 0.5067436 [5,] 1.07385364 1.01423431 0.8294663 0.5067436 0.0000000 > (drows <- D(t(teensex))) [,1] [,2] [,3] [1,] 0.0000000 0.2035345 0.9275454 [2,] 0.2035345 0.0000000 0.9355316 [3,] 0.9275454 0.9355316 0.0000000 > ################################################### > ### code chunk number 35: ch:MDS:teensex:plot > ################################################### > r1 <- cmdscale(dcols, eig = TRUE) > c1 <- cmdscale(drows, eig = TRUE) > plot(r1$points, xlim = range(r1$points[,1], c1$points[,1]) * 1.5, + ylim = range(r1$points[,1], c1$points[,1]) * 1.5, type = "n", + xlab = "Coordinate 1", ylab = "Coordinate 2", lwd = 2) > text(r1$points, labels = colnames(teensex), cex = 0.7) > text(c1$points, labels = rownames(teensex), cex = 0.7) > abline(h = 0, lty = 2) > abline(v = 0, lty = 2) > ################################################### > ### code chunk number 36: MDS-gardenflowers-tab > ################################################### > data("gardenflowers", package = "HSAUR2") > gfnames <- attr(gardenflowers, "Labels") > attr(gardenflowers, "Labels") <- gsub(" \\(.*", "", gfnames) > tmp <- as.matrix(gardenflowers) > colnames(tmp) <- abbreviate(colnames(tmp), 3) > tmp <- as.data.frame(tmp) > tmp <- HSAURtable(tmp, xname = "gardenflowers") > tmp$data[upper.tri(tmp$data)] <- " " > toLatex(tmp, pcol = 1, + caption = "Dissimilarity matrix of $18$ species of gardenflowers.", + label = "MDS-gardenflowers-tab", + rownames = TRUE) \index{gardenflowers data@\Robject{gardenflowers} data} \begin{center} \begin{longtable}{l rrrrrrrrrrrrrrrrrr } \caption{\Robject{gardenflowers} data. Dissimilarity matrix of $18$ species of gardenflowers. \label{MDS-gardenflowers-tab}} \\ \hline & \Robject{Bgn} & \Robject{Brm} & \Robject{Cml} & \Robject{Dhl} & \Robject{F--} & \Robject{Fch} & \Robject{Grn} & \Robject{Gld} & \Robject{Hth} & \Robject{Hyd} & \Robject{Irs} & \Robject{Lly} & \Robject{L--} & \Robject{Pny} & \Robject{Pnc} & \Robject{Rdr} & \Robject{Scr} & \Robject{Tlp} \\ \hline \endfirsthead \caption[]{\Robject{gardenflowers} data (continued).} \\ \hline & \Robject{Bgn} & \Robject{Brm} & \Robject{Cml} & \Robject{Dhl} & \Robject{F--} & \Robject{Fch} & \Robject{Grn} & \Robject{Gld} & \Robject{Hth} & \Robject{Hyd} & \Robject{Irs} & \Robject{Lly} & \Robject{L--} & \Robject{Pny} & \Robject{Pnc} & \Robject{Rdr} & \Robject{Scr} & \Robject{Tlp} \\ \hline \endhead Begonia & 0.00 & & & & & & & & & & & & & & & & & \\ Broom & 0.91 & 0.00 & & & & & & & & & & & & & & & & \\ Camellia & 0.49 & 0.67 & 0.00 & & & & & & & & & & & & & & & \\ Dahlia & 0.47 & 0.59 & 0.59 & 0.00 & & & & & & & & & & & & & & \\ Forget-me-not & 0.43 & 0.90 & 0.57 & 0.61 & 0.00 & & & & & & & & & & & & & \\ Fuchsia & 0.23 & 0.79 & 0.29 & 0.52 & 0.44 & 0.00 & & & & & & & & & & & & \\ Geranium & 0.31 & 0.70 & 0.54 & 0.44 & 0.54 & 0.24 & 0.00 & & & & & & & & & & & \\ Gladiolus & 0.49 & 0.57 & 0.71 & 0.26 & 0.49 & 0.68 & 0.49 & 0.00 & & & & & & & & & & \\ Heather & 0.57 & 0.57 & 0.57 & 0.89 & 0.50 & 0.61 & 0.70 & 0.77 & 0.00 & & & & & & & & & \\ Hydrangea & 0.76 & 0.58 & 0.58 & 0.62 & 0.39 & 0.61 & 0.86 & 0.70 & 0.55 & 0.00 & & & & & & & & \\ Iris & 0.32 & 0.77 & 0.63 & 0.75 & 0.46 & 0.52 & 0.60 & 0.63 & 0.46 & 0.47 & 0.00 & & & & & & & \\ Lily & 0.51 & 0.69 & 0.69 & 0.53 & 0.51 & 0.65 & 0.77 & 0.47 & 0.51 & 0.39 & 0.36 & 0.00 & & & & & & \\ Lily-of-the-valley & 0.59 & 0.75 & 0.75 & 0.77 & 0.35 & 0.63 & 0.72 & 0.65 & 0.35 & 0.41 & 0.45 & 0.24 & 0.00 & & & & & \\ Peony & 0.37 & 0.68 & 0.68 & 0.38 & 0.52 & 0.48 & 0.63 & 0.49 & 0.52 & 0.39 & 0.37 & 0.17 & 0.39 & 0.00 & & & & \\ Pink carnation & 0.74 & 0.54 & 0.70 & 0.58 & 0.54 & 0.74 & 0.50 & 0.49 & 0.36 & 0.52 & 0.60 & 0.48 & 0.39 & 0.49 & 0.00 & & & \\ Red rose & 0.84 & 0.41 & 0.75 & 0.37 & 0.82 & 0.71 & 0.61 & 0.64 & 0.81 & 0.43 & 0.84 & 0.62 & 0.67 & 0.47 & 0.45 & 0.00 & & \\ Scotch rose & 0.94 & 0.20 & 0.70 & 0.48 & 0.77 & 0.83 & 0.74 & 0.45 & 0.77 & 0.38 & 0.80 & 0.58 & 0.62 & 0.57 & 0.40 & 0.21 & 0.00 & \\ Tulip & 0.44 & 0.50 & 0.79 & 0.48 & 0.59 & 0.68 & 0.47 & 0.22 & 0.59 & 0.92 & 0.59 & 0.67 & 0.72 & 0.67 & 0.61 & 0.85 & 0.67 & 0.00 \\ \hline \end{longtable} \end{center} demo(Ch-MVA) ---- ~~~~~~ > ### R code from vignette source 'Ch-MVA.Rnw' > > ################################################### > ### code chunk number 1: setup > ################################################### > library("MVA") > set.seed(280875) > library("lattice") > lattice.options(default.theme = + function() + standard.theme("pdf", color = FALSE)) > if (file.exists("deparse.R")) { + if (!file.exists("figs")) dir.create("figs") + source("deparse.R") + options(prompt = "R> ", continue = "+ ", width = 64, + digits = 4, show.signif.stars = FALSE, useFancyQuotes = FALSE) + + options(SweaveHooks = list(onefig = function() {par(mfrow = c(1,1))}, + twofig = function() {par(mfrow = c(1,2))}, + figtwo = function() {par(mfrow = c(2,1))}, + threefig = function() {par(mfrow = c(1,3))}, + figthree = function() {par(mfrow = c(3,1))}, + fourfig = function() {par(mfrow = c(2,2))}, + sixfig = function() {par(mfrow = c(3,2))}, + nomar = function() par("mai" = c(0, 0, 0, 0)))) + } > ################################################### > ### code chunk number 2: ch:MVA:tab:hypo > ################################################### > hypo <- + structure(list(individual = 1:10, sex = structure(c(2L, 2L, 2L, + 2L, 2L, 1L, 1L, 1L, 1L, 1L), .Label = c("Female", "Male"), class = "factor"), + age = c(21L, 43L, 22L, 86L, 60L, 16L, NA, 43L, 22L, 80L), + IQ = c(120L, NA, 135L, 150L, 92L, 130L, 150L, NA, 84L, 70L + ), depression = structure(c(2L, 1L, 1L, 1L, 2L, 2L, 2L, 2L, + 1L, 1L), .Label = c("No", "Yes"), class = "factor"), health = structure(c(3L, + 3L, 1L, 4L, 2L, 2L, 3L, 1L, 1L, 2L), .Label = c("Average", + "Good", "Very good", "Very poor"), class = "factor"), weight = c(150L, + 160L, 135L, 140L, 110L, 110L, 120L, 120L, 105L, 100L)), .Names = c("individual", + "sex", "age", "IQ", "depression", "health", "weight"), class = "data.frame", row.names = c(NA, -10L)) > toLatex(HSAURtable(hypo), pcol = 1, + caption = "Hypothetical Set of Multivariate Data.", + label = "ch:MVA:tab:hypo", rownames = FALSE) \index{hypo data@\Robject{hypo} data} \begin{center} \begin{longtable}{ rrrrrrr } \caption{\Robject{hypo} data. Hypothetical Set of Multivariate Data. \label{ch:MVA:tab:hypo}} \\ \hline \Robject{individual} & \Robject{sex} & \Robject{age} & \Robject{IQ} & \Robject{depression} & \Robject{health} & \Robject{weight} \\ \hline \endfirsthead \caption[]{\Robject{hypo} data (continued).} \\ \hline \Robject{individual} & \Robject{sex} & \Robject{age} & \Robject{IQ} & \Robject{depression} & \Robject{health} & \Robject{weight} \\ \hline \endhead 1 & Male & 21 & 120 & Yes & Very good & 150 \\ 2 & Male & 43 & NA & No & Very good & 160 \\ 3 & Male & 22 & 135 & No & Average & 135 \\ 4 & Male & 86 & 150 & No & Very poor & 140 \\ 5 & Male & 60 & 92 & Yes & Good & 110 \\ 6 & Female & 16 & 130 & Yes & Good & 110 \\ 7 & Female & NA & 150 & Yes & Very good & 120 \\ 8 & Female & 43 & NA & Yes & Average & 120 \\ 9 & Female & 22 & 84 & No & Average & 105 \\ 10 & Female & 80 & 70 & No & Good & 100 \\ \hline \end{longtable} \end{center} > ################################################### > ### code chunk number 3: ch:MVA:hypo:subset > ################################################### > hypo[1:2, c("health", "weight")] health weight 1 Very good 150 2 Very good 160 > ################################################### > ### code chunk number 4: ch:MVA:tab:measure > ################################################### > measure <- + structure(list(V1 = 1:20, V2 = c(34L, 37L, 38L, 36L, 38L, 43L, + 40L, 38L, 40L, 41L, 36L, 36L, 34L, 33L, 36L, 37L, 34L, 36L, 38L, + 35L), V3 = c(30L, 32L, 30L, 33L, 29L, 32L, 33L, 30L, 30L, 32L, + 24L, 25L, 24L, 22L, 26L, 26L, 25L, 26L, 28L, 23L), V4 = c(32L, + 37L, 36L, 39L, 33L, 38L, 42L, 40L, 37L, 39L, 35L, 37L, 37L, 34L, + 38L, 37L, 38L, 37L, 40L, 35L)), .Names = c("V1", "V2", "V3", + "V4"), class = "data.frame", row.names = c(NA, -20L)) > measure <- measure[,-1] > names(measure) <- c("chest", "waist", "hips") > measure$gender <- gl(2, 10) > levels(measure$gender) <- c("male", "female") > toLatex(HSAURtable(measure), pcol = 2, + caption = "Chest, waist, and hip measurements on 20 individuals (in inches).", + label = "ch:MVA:tab:measure", rownames = FALSE) \index{measure data@\Robject{measure} data} \begin{center} \begin{longtable}{ rrrr|rrrr } \caption{\Robject{measure} data. Chest, waist, and hip measurements on 20 individuals (in inches). \label{ch:MVA:tab:measure}} \\ \hline \Robject{chest} & \Robject{waist} & \Robject{hips} & \Robject{gender} & \Robject{chest} & \Robject{waist} & \Robject{hips} & \Robject{gender} \\ \hline \endfirsthead \caption[]{\Robject{measure} data (continued).} \\ \hline \Robject{chest} & \Robject{waist} & \Robject{hips} & \Robject{gender} & \Robject{chest} & \Robject{waist} & \Robject{hips} & \Robject{gender} \\ \hline \endhead 34 & 30 & 32 & male & 36 & 24 & 35 & female \\ 37 & 32 & 37 & male & 36 & 25 & 37 & female \\ 38 & 30 & 36 & male & 34 & 24 & 37 & female \\ 36 & 33 & 39 & male & 33 & 22 & 34 & female \\ 38 & 29 & 33 & male & 36 & 26 & 38 & female \\ 43 & 32 & 38 & male & 37 & 26 & 37 & female \\ 40 & 33 & 42 & male & 34 & 25 & 38 & female \\ 38 & 30 & 40 & male & 36 & 26 & 37 & female \\ 40 & 30 & 37 & male & 38 & 28 & 40 & female \\ 41 & 32 & 39 & male & 35 & 23 & 35 & female \\ \hline \end{longtable} \end{center} > ################################################### > ### code chunk number 5: ch:MVA:tab:pottery > ################################################### > data("pottery", package = "HSAUR2") > toLatex(HSAURtable(pottery), pcol = 1, + caption = "Romano-British pottery data.", + label = "ch:MVA:tab:pottery", rownames = FALSE) \index{pottery data@\Robject{pottery} data} \begin{center} \begin{longtable}{ rrrrrrrrrr } \caption{\Robject{pottery} data. Romano-British pottery data. \label{ch:MVA:tab:pottery}} \\ \hline \Robject{Al2O3} & \Robject{Fe2O3} & \Robject{MgO} & \Robject{CaO} & \Robject{Na2O} & \Robject{K2O} & \Robject{TiO2} & \Robject{MnO} & \Robject{BaO} & \Robject{kiln} \\ \hline \endfirsthead \caption[]{\Robject{pottery} data (continued).} \\ \hline \Robject{Al2O3} & \Robject{Fe2O3} & \Robject{MgO} & \Robject{CaO} & \Robject{Na2O} & \Robject{K2O} & \Robject{TiO2} & \Robject{MnO} & \Robject{BaO} & \Robject{kiln} \\ \hline \endhead 18.8 & 9.52 & 2.00 & 0.79 & 0.40 & 3.20 & 1.01 & 0.077 & 0.015 & 1 \\ 16.9 & 7.33 & 1.65 & 0.84 & 0.40 & 3.05 & 0.99 & 0.067 & 0.018 & 1 \\ 18.2 & 7.64 & 1.82 & 0.77 & 0.40 & 3.07 & 0.98 & 0.087 & 0.014 & 1 \\ 16.9 & 7.29 & 1.56 & 0.76 & 0.40 & 3.05 & 1.00 & 0.063 & 0.019 & 1 \\ 17.8 & 7.24 & 1.83 & 0.92 & 0.43 & 3.12 & 0.93 & 0.061 & 0.019 & 1 \\ 18.8 & 7.45 & 2.06 & 0.87 & 0.25 & 3.26 & 0.98 & 0.072 & 0.017 & 1 \\ 16.5 & 7.05 & 1.81 & 1.73 & 0.33 & 3.20 & 0.95 & 0.066 & 0.019 & 1 \\ 18.0 & 7.42 & 2.06 & 1.00 & 0.28 & 3.37 & 0.96 & 0.072 & 0.017 & 1 \\ 15.8 & 7.15 & 1.62 & 0.71 & 0.38 & 3.25 & 0.93 & 0.062 & 0.017 & 1 \\ 14.6 & 6.87 & 1.67 & 0.76 & 0.33 & 3.06 & 0.91 & 0.055 & 0.012 & 1 \\ 13.7 & 5.83 & 1.50 & 0.66 & 0.13 & 2.25 & 0.75 & 0.034 & 0.012 & 1 \\ 14.6 & 6.76 & 1.63 & 1.48 & 0.20 & 3.02 & 0.87 & 0.055 & 0.016 & 1 \\ 14.8 & 7.07 & 1.62 & 1.44 & 0.24 & 3.03 & 0.86 & 0.080 & 0.016 & 1 \\ 17.1 & 7.79 & 1.99 & 0.83 & 0.46 & 3.13 & 0.93 & 0.090 & 0.020 & 1 \\ 16.8 & 7.86 & 1.86 & 0.84 & 0.46 & 2.93 & 0.94 & 0.094 & 0.020 & 1 \\ 15.8 & 7.65 & 1.94 & 0.81 & 0.83 & 3.33 & 0.96 & 0.112 & 0.019 & 1 \\ 18.6 & 7.85 & 2.33 & 0.87 & 0.38 & 3.17 & 0.98 & 0.081 & 0.018 & 1 \\ 16.9 & 7.87 & 1.83 & 1.31 & 0.53 & 3.09 & 0.95 & 0.092 & 0.023 & 1 \\ 18.9 & 7.58 & 2.05 & 0.83 & 0.13 & 3.29 & 0.98 & 0.072 & 0.015 & 1 \\ 18.0 & 7.50 & 1.94 & 0.69 & 0.12 & 3.14 & 0.93 & 0.035 & 0.017 & 1 \\ 17.8 & 7.28 & 1.92 & 0.81 & 0.18 & 3.15 & 0.90 & 0.067 & 0.017 & 1 \\ 14.4 & 7.00 & 4.30 & 0.15 & 0.51 & 4.25 & 0.79 & 0.160 & 0.019 & 2 \\ 13.8 & 7.08 & 3.43 & 0.12 & 0.17 & 4.14 & 0.77 & 0.144 & 0.020 & 2 \\ 14.6 & 7.09 & 3.88 & 0.13 & 0.20 & 4.36 & 0.81 & 0.124 & 0.019 & 2 \\ 11.5 & 6.37 & 5.64 & 0.16 & 0.14 & 3.89 & 0.69 & 0.087 & 0.009 & 2 \\ 13.8 & 7.06 & 5.34 & 0.20 & 0.20 & 4.31 & 0.71 & 0.101 & 0.021 & 2 \\ 10.9 & 6.26 & 3.47 & 0.17 & 0.22 & 3.40 & 0.66 & 0.109 & 0.010 & 2 \\ 10.1 & 4.26 & 4.26 & 0.20 & 0.18 & 3.32 & 0.59 & 0.149 & 0.017 & 2 \\ 11.6 & 5.78 & 5.91 & 0.18 & 0.16 & 3.70 & 0.65 & 0.082 & 0.015 & 2 \\ 11.1 & 5.49 & 4.52 & 0.29 & 0.30 & 4.03 & 0.63 & 0.080 & 0.016 & 2 \\ 13.4 & 6.92 & 7.23 & 0.28 & 0.20 & 4.54 & 0.69 & 0.163 & 0.017 & 2 \\ 12.4 & 6.13 & 5.69 & 0.22 & 0.54 & 4.65 & 0.70 & 0.159 & 0.015 & 2 \\ 13.1 & 6.64 & 5.51 & 0.31 & 0.24 & 4.89 & 0.72 & 0.094 & 0.017 & 2 \\ 11.6 & 5.39 & 3.77 & 0.29 & 0.06 & 4.51 & 0.56 & 0.110 & 0.015 & 3 \\ 11.8 & 5.44 & 3.94 & 0.30 & 0.04 & 4.64 & 0.59 & 0.085 & 0.013 & 3 \\ 18.3 & 1.28 & 0.67 & 0.03 & 0.03 & 1.96 & 0.65 & 0.001 & 0.014 & 4 \\ 15.8 & 2.39 & 0.63 & 0.01 & 0.04 & 1.94 & 1.29 & 0.001 & 0.014 & 4 \\ 18.0 & 1.50 & 0.67 & 0.01 & 0.06 & 2.11 & 0.92 & 0.001 & 0.016 & 4 \\ 18.0 & 1.88 & 0.68 & 0.01 & 0.04 & 2.00 & 1.11 & 0.006 & 0.022 & 4 \\ 20.8 & 1.51 & 0.72 & 0.07 & 0.10 & 2.37 & 1.26 & 0.002 & 0.016 & 4 \\ 17.7 & 1.12 & 0.56 & 0.06 & 0.06 & 2.06 & 0.79 & 0.001 & 0.013 & 5 \\ 18.3 & 1.14 & 0.67 & 0.06 & 0.05 & 2.11 & 0.89 & 0.006 & 0.019 & 5 \\ 16.7 & 0.92 & 0.53 & 0.01 & 0.05 & 1.76 & 0.91 & 0.004 & 0.013 & 5 \\ 14.8 & 2.74 & 0.67 & 0.03 & 0.05 & 2.15 & 1.34 & 0.003 & 0.015 & 5 \\ 19.1 & 1.64 & 0.60 & 0.10 & 0.03 & 1.75 & 1.04 & 0.007 & 0.018 & 5 \\ \hline \end{longtable} \end{center} > ################################################### > ### code chunk number 6: ch:MVA:tab:exam > ################################################### > exam <- + structure(list(subject = 1:5, math = c(60L, 80L, 53L, 85L, 45L + ), english = c(70L, 65L, 60L, 79L, 80L), history = c(75L, 66L, + 50L, 71L, 80L), geography = c(58L, 75L, 48L, 77L, 84L), chemistry = c(53L, + 70L, 45L, 68L, 44L), physics = c(42L, 76L, 43L, 79L, 46L)), .Names = c("subject", + "maths", "english", "history", "geography", "chemistry", "physics" + ), class = "data.frame", row.names = c(NA, -5L)) > toLatex(HSAURtable(exam), pcol = 1, + caption = "Exam scores for five psychology students.", + label = "ch:MVA:tab:exam", rownames = FALSE) \index{exam data@\Robject{exam} data} \begin{center} \begin{longtable}{ rrrrrrr } \caption{\Robject{exam} data. Exam scores for five psychology students. \label{ch:MVA:tab:exam}} \\ \hline \Robject{subject} & \Robject{maths} & \Robject{english} & \Robject{history} & \Robject{geography} & \Robject{chemistry} & \Robject{physics} \\ \hline \endfirsthead \caption[]{\Robject{exam} data (continued).} \\ \hline \Robject{subject} & \Robject{maths} & \Robject{english} & \Robject{history} & \Robject{geography} & \Robject{chemistry} & \Robject{physics} \\ \hline \endhead 1 & 60 & 70 & 75 & 58 & 53 & 42 \\ 2 & 80 & 65 & 66 & 75 & 70 & 76 \\ 3 & 53 & 60 & 50 & 48 & 45 & 43 \\ 4 & 85 & 79 & 71 & 77 & 68 & 79 \\ 5 & 45 & 80 & 80 & 84 & 44 & 46 \\ \hline \end{longtable} \end{center} > ################################################### > ### code chunk number 7: ch:MVA:USairpollution:tab > ################################################### > data("USairpollution", package = "HSAUR2") > toLatex(HSAURtable(USairpollution), pcol = 1, + caption = paste("Air pollution in $41$ US cities."), + label = "ch:MVA:USairpollution:tab", rownames = TRUE) \index{USairpollution data@\Robject{USairpollution} data} \begin{center} \begin{longtable}{l rrrrrrr } \caption{\Robject{USairpollution} data. Air pollution in $41$ US cities. \label{ch:MVA:USairpollution:tab}} \\ \hline & \Robject{SO2} & \Robject{temp} & \Robject{manu} & \Robject{popul} & \Robject{wind} & \Robject{precip} & \Robject{predays} \\ \hline \endfirsthead \caption[]{\Robject{USairpollution} data (continued).} \\ \hline & \Robject{SO2} & \Robject{temp} & \Robject{manu} & \Robject{popul} & \Robject{wind} & \Robject{precip} & \Robject{predays} \\ \hline \endhead Albany & 46 & 47.6 & 44 & 116 & 8.8 & 33.36 & 135 \\ Albuquerque & 11 & 56.8 & 46 & 244 & 8.9 & 7.77 & 58 \\ Atlanta & 24 & 61.5 & 368 & 497 & 9.1 & 48.34 & 115 \\ Baltimore & 47 & 55.0 & 625 & 905 & 9.6 & 41.31 & 111 \\ Buffalo & 11 & 47.1 & 391 & 463 & 12.4 & 36.11 & 166 \\ Charleston & 31 & 55.2 & 35 & 71 & 6.5 & 40.75 & 148 \\ Chicago & 110 & 50.6 & 3344 & 3369 & 10.4 & 34.44 & 122 \\ Cincinnati & 23 & 54.0 & 462 & 453 & 7.1 & 39.04 & 132 \\ Cleveland & 65 & 49.7 & 1007 & 751 & 10.9 & 34.99 & 155 \\ Columbus & 26 & 51.5 & 266 & 540 & 8.6 & 37.01 & 134 \\ Dallas & 9 & 66.2 & 641 & 844 & 10.9 & 35.94 & 78 \\ Denver & 17 & 51.9 & 454 & 515 & 9.0 & 12.95 & 86 \\ Des Moines & 17 & 49.0 & 104 & 201 & 11.2 & 30.85 & 103 \\ Detroit & 35 & 49.9 & 1064 & 1513 & 10.1 & 30.96 & 129 \\ Hartford & 56 & 49.1 & 412 & 158 & 9.0 & 43.37 & 127 \\ Houston & 10 & 68.9 & 721 & 1233 & 10.8 & 48.19 & 103 \\ Indianapolis & 28 & 52.3 & 361 & 746 & 9.7 & 38.74 & 121 \\ Jacksonville & 14 & 68.4 & 136 & 529 & 8.8 & 54.47 & 116 \\ Kansas City & 14 & 54.5 & 381 & 507 & 10.0 & 37.00 & 99 \\ Little Rock & 13 & 61.0 & 91 & 132 & 8.2 & 48.52 & 100 \\ Louisville & 30 & 55.6 & 291 & 593 & 8.3 & 43.11 & 123 \\ Memphis & 10 & 61.6 & 337 & 624 & 9.2 & 49.10 & 105 \\ Miami & 10 & 75.5 & 207 & 335 & 9.0 & 59.80 & 128 \\ Milwaukee & 16 & 45.7 & 569 & 717 & 11.8 & 29.07 & 123 \\ Minneapolis & 29 & 43.5 & 699 & 744 & 10.6 & 25.94 & 137 \\ Nashville & 18 & 59.4 & 275 & 448 & 7.9 & 46.00 & 119 \\ New Orleans & 9 & 68.3 & 204 & 361 & 8.4 & 56.77 & 113 \\ Norfolk & 31 & 59.3 & 96 & 308 & 10.6 & 44.68 & 116 \\ Omaha & 14 & 51.5 & 181 & 347 & 10.9 & 30.18 & 98 \\ Philadelphia & 69 & 54.6 & 1692 & 1950 & 9.6 & 39.93 & 115 \\ Phoenix & 10 & 70.3 & 213 & 582 & 6.0 & 7.05 & 36 \\ Pittsburgh & 61 & 50.4 & 347 & 520 & 9.4 & 36.22 & 147 \\ Providence & 94 & 50.0 & 343 & 179 & 10.6 & 42.75 & 125 \\ Richmond & 26 & 57.8 & 197 & 299 & 7.6 & 42.59 & 115 \\ Salt Lake City & 28 & 51.0 & 137 & 176 & 8.7 & 15.17 & 89 \\ San Francisco & 12 & 56.7 & 453 & 716 & 8.7 & 20.66 & 67 \\ Seattle & 29 & 51.1 & 379 & 531 & 9.4 & 38.79 & 164 \\ St. Louis & 56 & 55.9 & 775 & 622 & 9.5 & 35.89 & 105 \\ Washington & 29 & 57.3 & 434 & 757 & 9.3 & 38.89 & 111 \\ Wichita & 8 & 56.6 & 125 & 277 & 12.7 & 30.58 & 82 \\ Wilmington & 36 & 54.0 & 80 & 80 & 9.0 & 40.25 & 114 \\ \hline \end{longtable} \end{center} > ################################################### > ### code chunk number 8: ch:MVA:measure:cov > ################################################### > cov(measure[, c("chest", "waist", "hips")]) chest waist hips chest 6.631579 6.368421 3.000000 waist 6.368421 12.526316 3.578947 hips 3.000000 3.578947 5.944737 > ################################################### > ### code chunk number 9: ch:MVA:measure:cov > ################################################### > cov(subset(measure, gender == "female")[, + c("chest", "waist", "hips")]) chest waist hips chest 2.277778 2.166667 1.555556 waist 2.166667 2.988889 2.755556 hips 1.555556 2.755556 3.066667 > cov(subset(measure, gender == "male")[, + c("chest", "waist", "hips")]) chest waist hips chest 6.7222222 0.9444444 3.944444 waist 0.9444444 2.1000000 3.077778 hips 3.9444444 3.0777778 9.344444 > ################################################### > ### code chunk number 10: ch:MVA:hypo:cor > ################################################### > cor(measure[, c("chest", "waist", "hips")]) chest waist hips chest 1.0000000 0.6987336 0.4778004 waist 0.6987336 1.0000000 0.4147413 hips 0.4778004 0.4147413 1.0000000 > ################################################### > ### code chunk number 11: ch:MVA:measure:dist (eval = FALSE) > ################################################### > ## dist(scale(measure[, c("chest", "waist", "hips")], > ## center = FALSE)) > > > ################################################### > ### code chunk number 12: ch:MVA:measure:dist > ################################################### > x <- dist(scale(measure[, c("chest", "waist", "hips")], + center = FALSE)) > as.dist(round(as.matrix(x), 2)[1:12, 1:12]) 1 2 3 4 5 6 7 8 9 10 11 2 0.17 3 0.15 0.08 4 0.22 0.07 0.14 5 0.11 0.15 0.09 0.22 6 0.29 0.16 0.16 0.19 0.21 7 0.32 0.16 0.20 0.13 0.28 0.14 8 0.23 0.11 0.11 0.12 0.19 0.16 0.13 9 0.21 0.10 0.06 0.16 0.12 0.11 0.17 0.09 10 0.27 0.12 0.13 0.14 0.20 0.06 0.09 0.11 0.09 11 0.23 0.28 0.22 0.33 0.19 0.34 0.38 0.25 0.24 0.32 12 0.22 0.24 0.18 0.28 0.18 0.30 0.32 0.20 0.20 0.28 0.06 > cat("...") ... > ################################################### > ### code chunk number 13: ch:MVA:fig:dmvnorm > ################################################### > library("mvtnorm") > x <- y <- seq(from = -3, to = 3, length = 50) > dat <- as.matrix(expand.grid(x, y)) > d <- dmvnorm(dat, mean = c(0, 0), + sigma = matrix(c(1, 0.5, 0.5, 1), ncol = 2)) > d <- matrix(d, ncol = length(x)) > persp(x = x, y = y, z = d, xlab = "x1", ylab = "x2", + zlab = "f(x)") > ################################################### > ### code chunk number 14: ch:MVA:fig:cdf > ################################################### > x <- seq(from = -3, to = 3, length = 1000) > Fx <- pnorm(x) > Fy <- pnorm(x, mean = -1) > plot(x, Fx, type = "l", axes = FALSE, xlab = "", + ylab = "Cumulative distribution function") > lines(x, Fy, type = "l") > x0 <- which.min((x - 1.2)^2) > x05 <- which.min((x + 0.5)^2) > x08 <- which.min((x + 0.9)^2) > xx <- which.min(abs(Fy - Fx[x0])) > arrows(x0 = c(min(x), x[x0], x[xx], x[x08], x[x08], x[x08]), + y0 = c(Fx[x0], Fx[x0], Fy[xx], 0, Fx[x08], Fy[x08]), + x1 = c(x[x0], x[x0], x[xx], x[x08], -3, -3), + y1 = c(Fx[x0], 0, 0, Fy[x08], Fx[x08], Fy[x08])) > mtext(at = c(x[x08], x[xx], x[x0]), side = 1, line = 1, text = + c(expression(q), expression(q[2](p)), expression(q[1](p)))) > mtext(at = c(0, Fx[x08], Fy[x08], Fx[x0], 1), line = 1, side = 2, text = + c(0, expression(p[1](q)), expression(p[2](q)), expression(p), 1)) > box() > ################################################### > ### code chunk number 15: ch:MVA:fig:measure:chisq:setup1 > ################################################### > x <- measure[, c("chest", "waist", "hips")] > ################################################### > ### code chunk number 16: ch:MVA:fig:measure:chisq:setup2 > ################################################### > cm <- colMeans(x) > S <- cov(x) > ################################################### > ### code chunk number 17: ch:MVA:fig:measure:chisq:setup3 > ################################################### > d <- apply(x, MARGIN = 1, function(x) + t(x - cm) %*% solve(S) %*% (x - cm)) > ################################################### > ### code chunk number 18: ch:MVA:fig:measure:qq > ################################################### > qqnorm(measure[,"chest"], main = "chest"); qqline(measure[,"chest"]) > qqnorm(measure[,"waist"], main = "waist"); qqline(measure[,"waist"]) > qqnorm(measure[,"hips"], main = "hips"); qqline(measure[,"hips"]) > ################################################### > ### code chunk number 19: ch:MVA:fig:measure:chisq > ################################################### > plot(qchisq((1:nrow(x) - 1/2) / nrow(x), df = 3), sort(d), + xlab = expression(paste(chi[3]^2, " Quantile")), + ylab = "Ordered distances") > abline(a = 0, b = 1) > ################################################### > ### code chunk number 20: ch:MVA:fig:USairpollution:qq:setup (eval = FALSE) > ################################################### > ## layout(matrix(1:8, nc = 2)) > ## sapply(colnames(USairpollution), function(x) { > ## qqnorm(USairpollution[[x]], main = x) > ## qqline(USairpollution[[x]]) > ## }) > > > ################################################### > ### code chunk number 21: ch:MVA:fig:USairpollution:qq > ################################################### > layout(matrix(1:8, nc = 2)) > sapply(colnames(USairpollution), function(x) { + qqnorm(USairpollution[[x]], main = x) + qqline(USairpollution[[x]]) + }) $SO2 NULL $temp NULL $manu NULL $popul NULL $wind NULL $precip NULL $predays NULL > ################################################### > ### code chunk number 22: ch:MVA:fig:USairpollution:chisq > ################################################### > x <- USairpollution > cm <- colMeans(x) > S <- cov(x) > d <- apply(x, 1, function(x) t(x - cm) %*% solve(S) %*% (x - cm)) > plot(qc <- qchisq((1:nrow(x) - 1/2) / nrow(x), df = 7), + sd <- sort(d), + xlab = expression(paste(chi[7]^2, " Quantile")), + ylab = "Ordered distances", xlim = range(qc) * c(1, 1.1)) > oups <- which(rank(abs(qc - sd), ties = "random") > nrow(x) - 3) > text(qc[oups], sd[oups] - 1.5, names(oups)) > abline(a = 0, b = 1) > ################################################### > ### code chunk number 23: ex > ################################################### > s <- c(3.8778, + 2.8110, 2.1210, + 3.1480, 2.2669, 2.6550, + 3.5062, 2.5690, 2.8341, 3.2352) > S <- diag(4) > S[!lower.tri(S)] <- s > S <- S + t(S) > diag(S) <- diag(S) / 2 > writeLines(apply(S, 1, function(x) paste(paste(formatC(x, format = "f"), collapse = "&"), "\\\\"))) 3.8778&2.8110&3.1480&3.5062 \\ 2.8110&2.1210&2.2669&2.5690 \\ 3.1480&2.2669&2.6550&2.8341 \\ 3.5062&2.5690&2.8341&3.2352 \\ > ################################################### > ### code chunk number 24: ex > ################################################### > X <- matrix( + c(3, 4, 4, 6, 1, + 5, 1, 1, 7, 3, + 6, 2, 0, 2, 6, + 1, 1, 1, 0, 3, + 4, 7, 3, 6, 2, + 2, 2, 5, 1, 0, + 0, 4, 1, 1, 1, + 0, 6, 4, 3, 5, + 7, 6, 5, 1, 4, + 2, 1, 4, 3, 1), ncol = 5) > writeLines(apply(X, 1, function(x) paste(paste(x, collapse = "&"), "\\\\"))) 3&6&4&0&7 \\ 4&2&7&4&6 \\ 4&0&3&1&5 \\ 6&2&6&1&1 \\ 1&6&2&1&4 \\ 5&1&2&0&2 \\ 1&1&2&6&1 \\ 1&1&5&4&4 \\ 7&0&1&3&3 \\ 3&3&0&5&1 \\ demo(Ch-PCA) ---- ~~~~~~ > ### R code from vignette source 'Ch-PCA.Rnw' > > ################################################### > ### code chunk number 1: setup > ################################################### > library("MVA") > set.seed(280875) > library("lattice") > lattice.options(default.theme = + function() + standard.theme("pdf", color = FALSE)) > if (file.exists("deparse.R")) { + if (!file.exists("figs")) dir.create("figs") + source("deparse.R") + options(prompt = "R> ", continue = "+ ", width = 64, + digits = 4, show.signif.stars = FALSE, useFancyQuotes = FALSE) + + options(SweaveHooks = list(onefig = function() {par(mfrow = c(1,1))}, + twofig = function() {par(mfrow = c(1,2))}, + figtwo = function() {par(mfrow = c(2,1))}, + threefig = function() {par(mfrow = c(1,3))}, + figthree = function() {par(mfrow = c(3,1))}, + fourfig = function() {par(mfrow = c(2,2))}, + sixfig = function() {par(mfrow = c(3,2))}, + nomar = function() par("mai" = c(0, 0, 0, 0)))) + } > ################################################### > ### code chunk number 2: ch:PCA:data > ################################################### > bc <- c( + 0.290, + 0.202, 0.415, + -0.055, 0.285, 0.419, + -0.105, -0.376, -0.521, -0.877, + -0.252, -0.349, -0.441, -0.076, 0.206, + -0.229, -0.164, -0.145, 0.023, 0.034, 0.192, + 0.058, -0.129, -0.076, -0.131, 0.151, 0.077, 0.423) > blood_sd <- c(rblood = 0.371, plate = 41.253, wblood = 1.935, + neut = 0.077, lymph = 0.071, bilir = 4.037, + sodium = 2.732, potass = 0.297) > blood_corr <- diag(length(blood_sd)) / 2 > blood_corr[upper.tri(blood_corr)] <- bc > blood_corr <- blood_corr + t(blood_corr) > blood_cov <- blood_corr * outer(blood_sd, blood_sd, "*") > ################################################### > ### code chunk number 3: ch:PCA:blood_corr > ################################################### > blood_corr [,1] [,2] [,3] [,4] [,5] [,6] [,7] [,8] [1,] 1.000 0.290 0.202 -0.055 -0.105 -0.252 -0.229 0.058 [2,] 0.290 1.000 0.415 0.285 -0.376 -0.349 -0.164 -0.129 [3,] 0.202 0.415 1.000 0.419 -0.521 -0.441 -0.145 -0.076 [4,] -0.055 0.285 0.419 1.000 -0.877 -0.076 0.023 -0.131 [5,] -0.105 -0.376 -0.521 -0.877 1.000 0.206 0.034 0.151 [6,] -0.252 -0.349 -0.441 -0.076 0.206 1.000 0.192 0.077 [7,] -0.229 -0.164 -0.145 0.023 0.034 0.192 1.000 0.423 [8,] 0.058 -0.129 -0.076 -0.131 0.151 0.077 0.423 1.000 > ################################################### > ### code chunk number 4: ch:PCA:blood_sd > ################################################### > blood_sd rblood plate wblood neut lymph bilir sodium potass 0.371 41.253 1.935 0.077 0.071 4.037 2.732 0.297 > ################################################### > ### code chunk number 5: ch:PCA:blood:PCA > ################################################### > blood_pcacov <- princomp(covmat = blood_cov) > summary(blood_pcacov, loadings = TRUE) Importance of components: Comp.1 Comp.2 Comp.3 Comp.4 Standard deviation 41.2877486 3.880212624 2.64197339 1.624583979 Proportion of Variance 0.9856182 0.008705172 0.00403574 0.001525986 Cumulative Proportion 0.9856182 0.994323381 0.99835912 0.999885108 Comp.5 Comp.6 Comp.7 Comp.8 Standard deviation 0.353951757 2.561722e-01 8.510631e-02 2.372715e-02 Proportion of Variance 0.000072436 3.794288e-05 4.187837e-06 3.255049e-07 Cumulative Proportion 0.999957544 9.999955e-01 9.999997e-01 1.000000e+00 Loadings: Comp.1 Comp.2 Comp.3 Comp.4 Comp.5 Comp.6 Comp.7 Comp.8 rblood 0.943 0.329 plate 0.999 wblood 0.192 0.981 neut 0.758 0.650 lymph -0.649 0.760 bilir -0.961 0.195 0.191 sodium -0.193 -0.979 potass 0.329 -0.942 > blood_pcacor <- princomp(covmat = blood_corr) > summary(blood_pcacor, loadings = TRUE) Importance of components: Comp.1 Comp.2 Comp.3 Comp.4 Comp.5 Standard deviation 1.6710100 1.2375848 1.1177138 0.88227419 0.78839505 Proportion of Variance 0.3490343 0.1914520 0.1561605 0.09730097 0.07769584 Cumulative Proportion 0.3490343 0.5404863 0.6966468 0.79394778 0.87164363 Comp.6 Comp.7 Comp.8 Standard deviation 0.69917350 0.66002394 0.31996216 Proportion of Variance 0.06110545 0.05445395 0.01279697 Cumulative Proportion 0.93274908 0.98720303 1.00000000 Loadings: Comp.1 Comp.2 Comp.3 Comp.4 Comp.5 Comp.6 Comp.7 Comp.8 [1,] 0.194 0.417 0.400 0.652 0.175 0.363 0.176 0.102 [2,] 0.400 0.154 0.168 -0.848 -0.230 -0.110 [3,] 0.459 0.168 -0.274 0.251 -0.403 0.677 [4,] 0.430 -0.472 -0.171 0.169 0.118 -0.237 0.678 [5,] -0.494 0.360 -0.180 -0.139 -0.136 0.157 0.724 [6,] -0.319 -0.320 -0.277 0.633 -0.162 -0.384 0.377 [7,] -0.177 -0.535 0.410 -0.163 -0.299 0.513 0.367 [8,] -0.171 -0.245 0.709 0.198 -0.469 -0.376 > ################################################### > ### code chunk number 6: ch:PCA:blood:plot1 > ################################################### > plot(blood_pcacor$sdev^2, xlab = "Component number", + ylab = "Component variance", type = "l", main = "Scree diagram") > plot(log(blood_pcacor$sdev^2), xlab = "Component number", + ylab = "log(Component variance)", type="l", + main = "Log(eigenvalue) diagram") > ################################################### > ### code chunk number 7: ch:PCA:headsize:tab > ################################################### > "headsize" <- + matrix(c(191, 195, 181, 183, 176, 208, 189, 197, 188, 192, 179, 183, 174, 190, 188, 163, 195, 186, 181, 175, 192, 174, + 176, 197, 190, 155, 149, 148, 153, 144, 157, 150, 159, 152, 150, 158, 147, 150, 159, 151, 137, 155, 153, + 145, 140, 154, 143, 139, 167, 163, 179, 201, 185, 188, 171, 192, 190, 189, 197, 187, 186, 174, 185, 195, + 187, 161, 183, 173, 182, 165, 185, 178, 176, 200, 187, 145, 152, 149, 149, 142, 152, 149, 152, 159, 151, + 148, 147, 152, 157, 158, 130, 158, 148, 146, 137, 152, 147, 143, 158, 150) + , nrow = 25, ncol = 4 + , dimnames = list(character(0) + , c("head1", "breadth1", "head2", "breadth2"))) > x <- headsize > headsize <- as.data.frame(headsize) > toLatex(HSAURtable(headsize), pcol = 2, + caption = "Head Size Data.", + label = "ch:PCA:headsize:tab", rownames = FALSE) \index{headsize data@\Robject{headsize} data} \begin{center} \begin{longtable}{ rrrr|rrrr } \caption{\Robject{headsize} data. Head Size Data. \label{ch:PCA:headsize:tab}} \\ \hline \Robject{head1} & \Robject{breadth1} & \Robject{head2} & \Robject{breadth2} & \Robject{head1} & \Robject{breadth1} & \Robject{head2} & \Robject{breadth2} \\ \hline \endfirsthead \caption[]{\Robject{headsize} data (continued).} \\ \hline \Robject{head1} & \Robject{breadth1} & \Robject{head2} & \Robject{breadth2} & \Robject{head1} & \Robject{breadth1} & \Robject{head2} & \Robject{breadth2} \\ \hline \endhead 191 & 155 & 179 & 145 & 190 & 159 & 195 & 157 \\ 195 & 149 & 201 & 152 & 188 & 151 & 187 & 158 \\ 181 & 148 & 185 & 149 & 163 & 137 & 161 & 130 \\ 183 & 153 & 188 & 149 & 195 & 155 & 183 & 158 \\ 176 & 144 & 171 & 142 & 186 & 153 & 173 & 148 \\ 208 & 157 & 192 & 152 & 181 & 145 & 182 & 146 \\ 189 & 150 & 190 & 149 & 175 & 140 & 165 & 137 \\ 197 & 159 & 189 & 152 & 192 & 154 & 185 & 152 \\ 188 & 152 & 197 & 159 & 174 & 143 & 178 & 147 \\ 192 & 150 & 187 & 151 & 176 & 139 & 176 & 143 \\ 179 & 158 & 186 & 148 & 197 & 167 & 200 & 158 \\ 183 & 147 & 174 & 147 & 190 & 163 & 187 & 150 \\ 174 & 150 & 185 & 152 & & & & \\ \hline \end{longtable} \end{center} > headsize <- x > ################################################### > ### code chunk number 8: ch:PCA:head > ################################################### > head_dat <- headsize[, c("head1", "head2")] > colMeans(head_dat) head1 head2 185.72 183.84 > cov(head_dat) head1 head2 head1 95.29333 69.66167 head2 69.66167 100.80667 > ################################################### > ### code chunk number 9: ch:PCA:head:PCA > ################################################### > head_pca <- princomp(x = head_dat) > head_pca Call: princomp(x = head_dat) Standard deviations: Comp.1 Comp.2 12.690766 5.215406 2 variables and 25 observations. > print(summary(head_pca), loadings = TRUE) Importance of components: Comp.1 Comp.2 Standard deviation 12.6907660 5.2154059 Proportion of Variance 0.8555135 0.1444865 Cumulative Proportion 0.8555135 1.0000000 Loadings: Comp.1 Comp.2 head1 0.693 0.721 head2 0.721 -0.693 > ################################################### > ### code chunk number 10: ch:PCA:head:PCAint > ################################################### > s1 <- round(diag(cov(head_pca$scores))[1], 3) > s2 <- round(diag(cov(head_pca$scores))[2], 3) > s <- summary(head_pca) > l1 <- round(s$loadings[,1], 2) > l2 <- round(s$loadings[,2], 2) > ################################################### > ### code chunk number 11: ch:PCA:head:PCAvar > ################################################### > diag(cov(head_pca$scores)) Comp.1 Comp.2 167.76619 28.33381 > ################################################### > ### code chunk number 12: ch:PCA:head:plot1 > ################################################### > a1<-183.84-0.721*185.72/0.693 > b1<-0.721/0.693 > a2<-183.84-(-0.693*185.72/0.721) > b2<--0.693/0.721 > plot(head_dat, xlab = "First son's head length (mm)", + ylab = "Second son's head length") > abline(a1, b1) > abline(a2, b2, lty = 2) > ################################################### > ### code chunk number 13: ch:PCA:head:plot2 > ################################################### > xlim <- range(head_pca$scores[,1]) > plot(head_pca$scores, xlim = xlim, ylim = xlim) > ################################################### > ### code chunk number 14: ch:PCA:heptathlon:tab > ################################################### > data("heptathlon",package="HSAUR2") > toLatex(HSAURtable(heptathlon), pcol = 1, + caption = "Results of Olympic heptathlon, Seoul, 1988.", + label = "ch:PCA:heptathlon:tab", + rownames = TRUE) \index{heptathlon data@\Robject{heptathlon} data} \begin{center} \begin{longtable}{l rrrrrrrr } \caption{\Robject{heptathlon} data. Results of Olympic heptathlon, Seoul, 1988. \label{ch:PCA:heptathlon:tab}} \\ \hline & \Robject{hurdles} & \Robject{highjump} & \Robject{shot} & \Robject{run200m} & \Robject{longjump} & \Robject{javelin} & \Robject{run800m} & \Robject{score} \\ \hline \endfirsthead \caption[]{\Robject{heptathlon} data (continued).} \\ \hline & \Robject{hurdles} & \Robject{highjump} & \Robject{shot} & \Robject{run200m} & \Robject{longjump} & \Robject{javelin} & \Robject{run800m} & \Robject{score} \\ \hline \endhead Joyner-Kersee (USA) & 12.69 & 1.86 & 15.80 & 22.56 & 7.27 & 45.66 & 128.51 & 7291 \\ John (GDR) & 12.85 & 1.80 & 16.23 & 23.65 & 6.71 & 42.56 & 126.12 & 6897 \\ Behmer (GDR) & 13.20 & 1.83 & 14.20 & 23.10 & 6.68 & 44.54 & 124.20 & 6858 \\ Sablovskaite (URS) & 13.61 & 1.80 & 15.23 & 23.92 & 6.25 & 42.78 & 132.24 & 6540 \\ Choubenkova (URS) & 13.51 & 1.74 & 14.76 & 23.93 & 6.32 & 47.46 & 127.90 & 6540 \\ Schulz (GDR) & 13.75 & 1.83 & 13.50 & 24.65 & 6.33 & 42.82 & 125.79 & 6411 \\ Fleming (AUS) & 13.38 & 1.80 & 12.88 & 23.59 & 6.37 & 40.28 & 132.54 & 6351 \\ Greiner (USA) & 13.55 & 1.80 & 14.13 & 24.48 & 6.47 & 38.00 & 133.65 & 6297 \\ Lajbnerova (CZE) & 13.63 & 1.83 & 14.28 & 24.86 & 6.11 & 42.20 & 136.05 & 6252 \\ Bouraga (URS) & 13.25 & 1.77 & 12.62 & 23.59 & 6.28 & 39.06 & 134.74 & 6252 \\ Wijnsma (HOL) & 13.75 & 1.86 & 13.01 & 25.03 & 6.34 & 37.86 & 131.49 & 6205 \\ Dimitrova (BUL) & 13.24 & 1.80 & 12.88 & 23.59 & 6.37 & 40.28 & 132.54 & 6171 \\ Scheider (SWI) & 13.85 & 1.86 & 11.58 & 24.87 & 6.05 & 47.50 & 134.93 & 6137 \\ Braun (FRG) & 13.71 & 1.83 & 13.16 & 24.78 & 6.12 & 44.58 & 142.82 & 6109 \\ Ruotsalainen (FIN) & 13.79 & 1.80 & 12.32 & 24.61 & 6.08 & 45.44 & 137.06 & 6101 \\ Yuping (CHN) & 13.93 & 1.86 & 14.21 & 25.00 & 6.40 & 38.60 & 146.67 & 6087 \\ Hagger (GB) & 13.47 & 1.80 & 12.75 & 25.47 & 6.34 & 35.76 & 138.48 & 5975 \\ Brown (USA) & 14.07 & 1.83 & 12.69 & 24.83 & 6.13 & 44.34 & 146.43 & 5972 \\ Mulliner (GB) & 14.39 & 1.71 & 12.68 & 24.92 & 6.10 & 37.76 & 138.02 & 5746 \\ Hautenauve (BEL) & 14.04 & 1.77 & 11.81 & 25.61 & 5.99 & 35.68 & 133.90 & 5734 \\ Kytola (FIN) & 14.31 & 1.77 & 11.66 & 25.69 & 5.75 & 39.48 & 133.35 & 5686 \\ Geremias (BRA) & 14.23 & 1.71 & 12.95 & 25.50 & 5.50 & 39.64 & 144.02 & 5508 \\ Hui-Ing (TAI) & 14.85 & 1.68 & 10.00 & 25.23 & 5.47 & 39.14 & 137.30 & 5290 \\ Jeong-Mi (KOR) & 14.53 & 1.71 & 10.83 & 26.61 & 5.50 & 39.26 & 139.17 & 5289 \\ Launa (PNG) & 16.42 & 1.50 & 11.78 & 26.16 & 4.88 & 46.38 & 163.43 & 4566 \\ \hline \end{longtable} \end{center} > ################################################### > ### code chunk number 15: ch:PCA:heptathlon:recode > ################################################### > heptathlon$hurdles <- with(heptathlon, max(hurdles)-hurdles) > heptathlon$run200m <- with(heptathlon, max(run200m)-run200m) > heptathlon$run800m <- with(heptathlon, max(run800m)-run800m) > score <- which(colnames(heptathlon) == "score") > round(cor(heptathlon[,-score]), 2) hurdles highjump shot run200m longjump javelin run800m hurdles 1.00 0.81 0.65 0.77 0.91 0.01 0.78 highjump 0.81 1.00 0.44 0.49 0.78 0.00 0.59 shot 0.65 0.44 1.00 0.68 0.74 0.27 0.42 run200m 0.77 0.49 0.68 1.00 0.82 0.33 0.62 longjump 0.91 0.78 0.74 0.82 1.00 0.07 0.70 javelin 0.01 0.00 0.27 0.33 0.07 1.00 -0.02 run800m 0.78 0.59 0.42 0.62 0.70 -0.02 1.00 > plot(heptathlon[,-score]) > ################################################### > ### code chunk number 16: ch:PCA:heptathlon:scatter > ################################################### > plot(heptathlon[,-score], pch = ".", cex = 1.5) > ################################################### > ### code chunk number 17: ch:PCA:heptathlon:PNG > ################################################### > heptathlon <- heptathlon[-grep("PNG", rownames(heptathlon)),] > score <- which(colnames(heptathlon) == "score") > round(cor(heptathlon[,-score]), 2) hurdles highjump shot run200m longjump javelin run800m hurdles 1.00 0.58 0.77 0.83 0.89 0.33 0.56 highjump 0.58 1.00 0.46 0.39 0.66 0.35 0.15 shot 0.77 0.46 1.00 0.67 0.78 0.34 0.41 run200m 0.83 0.39 0.67 1.00 0.81 0.47 0.57 longjump 0.89 0.66 0.78 0.81 1.00 0.29 0.52 javelin 0.33 0.35 0.34 0.47 0.29 1.00 0.26 run800m 0.56 0.15 0.41 0.57 0.52 0.26 1.00 > ################################################### > ### code chunk number 18: ch:PCA:heptathlon:scatter2 > ################################################### > plot(heptathlon[,-score], pch = ".", cex = 1.5) > ################################################### > ### code chunk number 19: PCA-opt > ################################################### > op <- options(digits = 2) > ################################################### > ### code chunk number 20: PCA-heptathlon-pca > ################################################### > heptathlon_pca <- prcomp(heptathlon[, -score], scale = TRUE) > print(heptathlon_pca) Standard deviations (1, .., p=7): [1] 2.08 0.95 0.91 0.68 0.55 0.34 0.26 Rotation (n x k) = (7 x 7): PC1 PC2 PC3 PC4 PC5 PC6 PC7 hurdles -0.45 0.058 0.17 -0.048 -0.199 0.847 -0.070 highjump -0.31 -0.651 0.21 0.557 0.071 -0.090 0.332 shot -0.40 -0.022 0.15 -0.548 0.672 -0.099 0.229 run200m -0.43 0.185 -0.13 -0.231 -0.618 -0.333 0.470 longjump -0.45 -0.025 0.27 0.015 -0.122 -0.383 -0.749 javelin -0.24 -0.326 -0.88 -0.060 0.079 0.072 -0.211 run800m -0.30 0.657 -0.19 0.574 0.319 -0.052 0.077 > ################################################### > ### code chunk number 21: PCA-heptathlon-summary > ################################################### > summary(heptathlon_pca) Importance of components: PC1 PC2 PC3 PC4 PC5 PC6 PC7 Standard deviation 2.079 0.948 0.911 0.6832 0.5462 0.3375 0.26204 Proportion of Variance 0.618 0.128 0.119 0.0667 0.0426 0.0163 0.00981 Cumulative Proportion 0.618 0.746 0.865 0.9313 0.9739 0.9902 1.00000 > ################################################### > ### code chunk number 22: PCA-heptathlon-a1 > ################################################### > a1 <- heptathlon_pca$rotation[,1] > a1 hurdles highjump shot run200m longjump javelin run800m -0.45 -0.31 -0.40 -0.43 -0.45 -0.24 -0.30 > ################################################### > ### code chunk number 23: PCA-heptathlon-scaling > ################################################### > center <- heptathlon_pca$center > scale <- heptathlon_pca$scale > ################################################### > ### code chunk number 24: PCA-heptathlon-s1 > ################################################### > hm <- as.matrix(heptathlon[,-score]) > drop(scale(hm, center = center, scale = scale) %*% + heptathlon_pca$rotation[,1]) Joyner-Kersee (USA) John (GDR) Behmer (GDR) Sablovskaite (URS) -4.758 -3.148 -2.926 -1.288 Choubenkova (URS) Schulz (GDR) Fleming (AUS) Greiner (USA) -1.503 -0.958 -0.953 -0.633 Lajbnerova (CZE) Bouraga (URS) Wijnsma (HOL) Dimitrova (BUL) -0.382 -0.522 -0.218 -1.076 Scheider (SWI) Braun (FRG) Ruotsalainen (FIN) Yuping (CHN) 0.003 0.109 0.209 0.233 Hagger (GB) Brown (USA) Mulliner (GB) Hautenauve (BEL) 0.660 0.757 1.881 1.828 Kytola (FIN) Geremias (BRA) Hui-Ing (TAI) Jeong-Mi (KOR) 2.118 2.771 3.901 3.897 > ################################################### > ### code chunk number 25: PCA-heptathlon-s1 > ################################################### > predict(heptathlon_pca)[,1] Joyner-Kersee (USA) John (GDR) Behmer (GDR) Sablovskaite (URS) -4.758 -3.148 -2.926 -1.288 Choubenkova (URS) Schulz (GDR) Fleming (AUS) Greiner (USA) -1.503 -0.958 -0.953 -0.633 Lajbnerova (CZE) Bouraga (URS) Wijnsma (HOL) Dimitrova (BUL) -0.382 -0.522 -0.218 -1.076 Scheider (SWI) Braun (FRG) Ruotsalainen (FIN) Yuping (CHN) 0.003 0.109 0.209 0.233 Hagger (GB) Brown (USA) Mulliner (GB) Hautenauve (BEL) 0.660 0.757 1.881 1.828 Kytola (FIN) Geremias (BRA) Hui-Ing (TAI) Jeong-Mi (KOR) 2.118 2.771 3.901 3.897 > ################################################### > ### code chunk number 26: PCA-heptathlon-sdev > ################################################### > sdev <- heptathlon_pca$sdev > prop12 <- round(sum(sdev[1:2]^2)/sum(sdev^2)*100, 0) > ################################################### > ### code chunk number 27: PCA-heptathlon-pca-plot (eval = FALSE) > ################################################### > ## plot(heptathlon_pca) > > > ################################################### > ### code chunk number 28: PCA-heptathlon-pca-plot > ################################################### > plot(heptathlon_pca, main = "") > ################################################### > ### code chunk number 29: PCA-scorecor > ################################################### > cor(heptathlon$score, heptathlon_pca$x[,1]) [1] -0.99 > ################################################### > ### code chunk number 30: PCA-heptathlonscore > ################################################### > plot(heptathlon$score, heptathlon_pca$x[,1]) > ################################################### > ### code chunk number 31: ch:PCA:USairpollution:scatter > ################################################### > data("USairpollution", package = "HSAUR2") > panel.hist <- function(x, ...) { + usr <- par("usr"); on.exit(par(usr)) + par(usr = c(usr[1:2], 0, 1.5) ) + h <- hist(x, plot = FALSE) + breaks <- h$breaks; nB <- length(breaks) + y <- h$counts; y <- y/max(y) + rect(breaks[-nB], 0, breaks[-1], y, col="grey", ...) + } > USairpollution$negtemp <- USairpollution$temp * (-1) > USairpollution$temp <- NULL > pairs(USairpollution[,-1], diag.panel = panel.hist, + pch = ".", cex = 1.5) > ################################################### > ### code chunk number 32: ch:PCA:USairpollution:pca > ################################################### > cor(USairpollution[,-1]) manu popul wind precip predays negtemp manu 1.000 0.955 0.238 -0.032 0.132 0.190 popul 0.955 1.000 0.213 -0.026 0.042 0.063 wind 0.238 0.213 1.000 -0.013 0.164 0.350 precip -0.032 -0.026 -0.013 1.000 0.496 -0.386 predays 0.132 0.042 0.164 0.496 1.000 0.430 negtemp 0.190 0.063 0.350 -0.386 0.430 1.000 > usair_pca <- princomp(USairpollution[,-1], cor = TRUE) > ################################################### > ### code chunk number 33: ch:PCA:USairpollution:pcasummary > ################################################### > summary(usair_pca, loadings = TRUE) Importance of components: Comp.1 Comp.2 Comp.3 Comp.4 Comp.5 Comp.6 Standard deviation 1.48 1.22 1.18 0.87 0.338 0.1856 Proportion of Variance 0.37 0.25 0.23 0.13 0.019 0.0057 Cumulative Proportion 0.37 0.62 0.85 0.98 0.994 1.0000 Loadings: Comp.1 Comp.2 Comp.3 Comp.4 Comp.5 Comp.6 manu 0.612 0.168 0.273 0.137 0.102 0.703 popul 0.578 0.222 0.350 -0.695 wind 0.354 -0.131 -0.297 -0.869 -0.113 precip -0.623 0.505 -0.171 0.568 predays 0.238 -0.708 0.311 -0.580 negtemp 0.330 -0.128 -0.672 0.306 0.558 -0.136 > ################################################### > ### code chunk number 34: ch:PCA:USairpollution:pcaplot > ################################################### > pairs(usair_pca$scores[,1:3], ylim = c(-6, 4), xlim = c(-6, 4), + panel = function(x,y, ...) { + text(x, y, abbreviate(row.names(USairpollution)), + cex = 0.6) + bvbox(cbind(x,y), add = TRUE) + }) > ################################################### > ### code chunk number 35: ch:PCA:USairpollution:lmplot > ################################################### > out <- sapply(1:6, function(i) { + plot(USairpollution$SO2,usair_pca$scores[,i], + xlab = paste("PC", i, sep = ""), + ylab = "Sulphur dioxide concentration") + }) > ################################################### > ### code chunk number 36: ch:PCA:USairpollution:lm > ################################################### > usair_reg <- lm(SO2 ~ usair_pca$scores, + data = USairpollution) > summary(usair_reg) Call: lm(formula = SO2 ~ usair_pca$scores, data = USairpollution) Residuals: Min 1Q Median 3Q Max -23.00 -8.54 -0.99 5.76 48.76 Coefficients: Estimate Std. Error t value Pr(>|t|) (Intercept) 30.049 2.286 13.15 6.9e-15 *** usair_pca$scoresComp.1 9.942 1.542 6.45 2.3e-07 *** usair_pca$scoresComp.2 -2.240 1.866 -1.20 0.2384 usair_pca$scoresComp.3 0.375 1.936 0.19 0.8475 usair_pca$scoresComp.4 8.549 2.622 3.26 0.0025 ** usair_pca$scoresComp.5 15.176 6.753 2.25 0.0312 * usair_pca$scoresComp.6 39.271 12.316 3.19 0.0031 ** --- Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1 Residual standard error: 15 on 34 degrees of freedom Multiple R-squared: 0.67, Adjusted R-squared: 0.611 F-statistic: 11.5 on 6 and 34 DF, p-value: 5.42e-07 > ################################################### > ### code chunk number 37: PCA-heptathlon-biplot (eval = FALSE) > ################################################### > ## biplot(heptathlon_pca, col = c("gray", "black")) > > > ################################################### > ### code chunk number 38: PCA-heptathlon-biplot > ################################################### > tmp <- heptathlon[, -score] > rownames(tmp) <- abbreviate(gsub(" \\(.*", "", rownames(tmp))) > biplot(prcomp(tmp, scale = TRUE), col = c("black", "darkgray"), xlim = + c(-0.5, 0.7), cex = 0.7) > ################################################### > ### code chunk number 39: ch:PCA:headsize > ################################################### > headsize.std <- sweep(headsize, 2, + apply(headsize, 2, sd), FUN = "/") > R <- cor(headsize.std) > r11 <- R[1:2, 1:2] > r22 <- R[-(1:2), -(1:2)] > r12 <- R[1:2, -(1:2)] > r21 <- R[-(1:2), 1:2] > (E1 <- solve(r11) %*% r12 %*% solve(r22) %*%r21) head1 breadth1 head1 0.32 0.32 breadth1 0.30 0.30 > (E2 <- solve(r22) %*% r21 %*% solve(r11) %*%r12) head2 breadth2 head2 0.30 0.30 breadth2 0.32 0.32 > (e1 <- eigen(E1)) eigen() decomposition $values [1] 0.6217 0.0029 $vectors [,1] [,2] [1,] 0.73 -0.70 [2,] 0.69 0.71 > (e2 <- eigen(E2)) eigen() decomposition $values [1] 0.6217 0.0029 $vectors [,1] [,2] [1,] -0.68 -0.71 [2,] -0.73 0.71 > ################################################### > ### code chunk number 40: ch:PCA:dummy > ################################################### > p <- function(x) formatC(x, format = "f", digits = 2) > f <- function(x, add = 0) paste(ifelse(x < 0, "-", "+"), p(abs(x)), "x_", 1:length(x) + add, + collapse = "") > ff <- function(x, xname) paste(ifelse(x < 0, "-", "+"), p(abs(x)), "\\\\text{", xname, "}", + collapse = "") > ################################################### > ### code chunk number 41: ch:PCA:headsize-cor > ################################################### > girth1 <- headsize.std[,1:2] %*% e1$vectors[,1] > girth2 <- headsize.std[,3:4] %*% e2$vectors[,1] > shape1 <- headsize.std[,1:2] %*% e1$vectors[,2] > shape2 <- headsize.std[,3:4] %*% e2$vectors[,2] > (g <- cor(girth1, girth2)) [,1] [1,] -0.79 > (s <- cor(shape1, shape2)) [,1] [1,] 0.054 > ################################################### > ### code chunk number 42: ch:PCA:headsize:plot > ################################################### > plot(girth1, girth2) > plot(shape1, shape2) > ################################################### > ### code chunk number 43: ch:PCA:LAdepr:tab > ################################################### > depr <- c( + 0.212, + 0.124, 0.098, + -0.164, 0.308, 0.044, + -0.101, -0.207, -0.106, -0.208, + -0.158, -0.183, -0.180, -0.192, 0.492) > LAdepr <- diag(6) / 2 > LAdepr[upper.tri(LAdepr)] <- depr > LAdepr <- LAdepr + t(LAdepr) > rownames(LAdepr) <- colnames(LAdepr) <- c("CESD", "Health", "Gender", "Age", "Edu", "Income") > x <- LAdepr > LAdepr <- as.data.frame(LAdepr) > toLatex(HSAURtable(LAdepr), + caption = "Los Angeles Depression Data.", + label = "ch:PCA:LAdepr:tab", rownames = FALSE) \index{LAdepr data@\Robject{LAdepr} data} \begin{center} \begin{longtable}{ rrrrrr } \caption{\Robject{LAdepr} data. Los Angeles Depression Data. \label{ch:PCA:LAdepr:tab}} \\ \hline \Robject{CESD} & \Robject{Health} & \Robject{Gender} & \Robject{Age} & \Robject{Edu} & \Robject{Income} \\ \hline \endfirsthead \caption[]{\Robject{LAdepr} data (continued).} \\ \hline \Robject{CESD} & \Robject{Health} & \Robject{Gender} & \Robject{Age} & \Robject{Edu} & \Robject{Income} \\ \hline \endhead 1.000 & 0.212 & 0.124 & -0.164 & -0.101 & -0.158 \\ 0.212 & 1.000 & 0.098 & 0.308 & -0.207 & -0.183 \\ 0.124 & 0.098 & 1.000 & 0.044 & -0.106 & -0.180 \\ -0.164 & 0.308 & 0.044 & 1.000 & -0.208 & -0.192 \\ -0.101 & -0.207 & -0.106 & -0.208 & 1.000 & 0.492 \\ -0.158 & -0.183 & -0.180 & -0.192 & 0.492 & 1.000 \\ \hline \end{longtable} \end{center} > LAdepr <- x > ################################################### > ### code chunk number 44: ch:PCA:LAdepr:CCA > ################################################### > r11 <- LAdepr[1:2, 1:2] > r22 <- LAdepr[-(1:2), -(1:2)] > r12 <- LAdepr[1:2, -(1:2)] > r21 <- LAdepr[-(1:2), 1:2] > (E1 <- solve(r11) %*% r12 %*% solve(r22) %*%r21) CESD Health CESD 0.084 -0.043 Health -0.033 0.133 > (E2 <- solve(r22) %*% r21 %*% solve(r11) %*%r12) Gender Age Edu Income Gender 0.0155 -0.0015 -0.018 -0.022 Age -0.0024 0.1471 -0.040 -0.016 Edu -0.0149 -0.0260 0.025 0.025 Income -0.0212 0.0130 0.022 0.029 > (e1 <- eigen(E1)) eigen() decomposition $values [1] 0.153 0.063 $vectors [,1] [,2] [1,] 0.53 -0.91 [2,] -0.85 -0.42 > (e2 <- eigen(E2)) eigen() decomposition $values [1] 1.5e-01 6.3e-02 1.4e-17 6.9e-18 $vectors [,1] [,2] [,3] [,4] [1,] 0.0026 -0.49 -0.87 0.579 [2,] 0.9801 0.32 -0.13 -0.017 [3,] -0.1858 0.43 -0.26 -0.386 [4,] 0.0699 0.69 -0.39 0.718 demo(Ch-SEM) ---- ~~~~~~ > ### R code from vignette source 'Ch-SEM.Rnw' > > ################################################### > ### code chunk number 1: setup > ################################################### > library("MVA") > set.seed(280875) > library("lattice") > lattice.options(default.theme = + function() + standard.theme("pdf", color = FALSE)) > if (file.exists("deparse.R")) { + if (!file.exists("figs")) dir.create("figs") + source("deparse.R") + options(prompt = "R> ", continue = "+ ", width = 64, + digits = 4, show.signif.stars = FALSE, useFancyQuotes = FALSE) + + options(SweaveHooks = list(onefig = function() {par(mfrow = c(1,1))}, + twofig = function() {par(mfrow = c(1,2))}, + figtwo = function() {par(mfrow = c(2,1))}, + threefig = function() {par(mfrow = c(1,3))}, + figthree = function() {par(mfrow = c(3,1))}, + fourfig = function() {par(mfrow = c(2,2))}, + sixfig = function() {par(mfrow = c(3,2))}, + nomar = function() par("mai" = c(0, 0, 0, 0)))) + } > ################################################### > ### code chunk number 2: setup > ################################################### > library("sem") > ab <- c(0.73, + 0.70, 0.68, + 0.58, 0.61, 0.57, + 0.46, 0.43, 0.40, 0.37, + 0.56, 0.52, 0.48, 0.41, 0.72) > ability <- diag(6) / 2 > ability[upper.tri(ability)] <- ab > ability <- ability + t(ability) > rownames(ability) <- colnames(ability) <- + c("SCA","PPE","PTE","PFE","EA","CP") > alienation <- matrix(c(11.839, 6.947, 6.819, 4.783, -3.834, -21.899, + 6.947, 9.364, 5.091, 5.028, -3.889, -18.831, + 6.819, 5.091,12.532, 7.495, -3.841, -21.748, + 4.783, 5.028, 7.495, 9.986, -3.625, -18.775, + -3.834,-3.889,-3.841,-3.625, 9.60, 35.522, + -21.899,-18.831,-21.748,-18.775,35.522,450.283), + ncol = 6, byrow = TRUE) > rownames(alienation) <- colnames(alienation) <- c("Anomia67", + "Powles67", "Anomia71", "Powles71", "Educ", "SEI") > d <- + c(0.447, + 0.422, 0.619, + 0.435, 0.604, 0.583, + 0.114, 0.068, 0.053, 0.115, + 0.203, 0.146, 0.139, 0.258, 0.349, + 0.091, 0.103, 0.110, 0.122, 0.209, 0.221, + 0.082, 0.063, 0.066, 0.097, 0.321, 0.355, 0.201, + 0.513, 0.445, 0.365, 0.482, 0.186, 0.315, 0.150, 0.154, + 0.304, 0.318, 0.240, 0.368, 0.303, 0.377, 0.163, 0.219, 0.534, + 0.245, 0.203, 0.183, 0.255, 0.272, 0.323, 0.310, 0.288, 0.301, 0.302, + 0.101, 0.088, 0.074, 0.139, 0.279, 0.367, 0.232, 0.320, 0.204, 0.368, 0.340, + 0.245, 0.199, 0.184, 0.293, 0.278, 0.545, 0.232, 0.314, 0.394, 0.467, 0.392, 0.511) > druguse <- diag(13) / 2 > druguse[upper.tri(druguse)] <- d > druguse <- druguse + t(druguse) > rownames(druguse) <- colnames(druguse) <- c("cigarettes", "beer", "wine", "liquor", "cocaine", + "tranquillizers", "drug store medication", "heroin", + "marijuana", "hashish", "inhalants", "haluucinogenics", "amphetamine") > ################################################### > ### code chunk number 3: ch:SEM:ability:plot > ################################################### > ord <- order.dendrogram(as.dendrogram(hclust(dist(ability)))) > panel.corrgram <- + function(x, y, z, subscripts, at, + level = 0.9, label = FALSE, ...) + { + require("ellipse", quietly = TRUE) + x <- as.numeric(x)[subscripts] + y <- as.numeric(y)[subscripts] + z <- as.numeric(z)[subscripts] + zcol <- level.colors(z, at = at, col.regions = grey.colors, ...) + for (i in seq(along = z)) { + ell <- ellipse(z[i], level = level, npoints = 50, + scale = c(.2, .2), centre = c(x[i], y[i])) + panel.polygon(ell, col = zcol[i], border = zcol[i], ...) + } + if (label) + panel.text(x = x, y = y, lab = 100 * round(z, 2), cex = 0.8, + col = ifelse(z < 0, "white", "black")) + } > print(levelplot(ability[ord, ord], at = do.breaks(c(-1.01, 1.01), 20), + xlab = NULL, ylab = NULL, colorkey = list(space = "top"), + scales = list(x = list(rot = 90)), + panel = panel.corrgram, label = TRUE)) > ################################################### > ### code chunk number 4: ch:SEM:ability-model > ################################################### > mod <- c("Ability -> SCA, lambda1, NA", + "Ability -> PPE, lambda2, NA", + "Ability -> PTE, lambda3, NA", + "Ability -> PFE, lambda4, NA", + "Aspiration -> EA, lambda5, NA", + "Aspiration -> CP, lambda6, NA", + "Ability <-> Aspiration, rho, NA", + "SCA <-> SCA, theta1, NA", + "PPE <-> PPE, theta2, NA", + "PTE <-> PTE, theta3, NA", + "PFE <-> PFE, theta4, NA", + "EA <-> EA, theta5, NA", + "CP <-> CP, theta6, NA", + "Ability <-> Ability, NA, 1", + "Aspiration <-> Aspiration, NA, 1") > ################################################### > ### code chunk number 5: ch:SEM:ability-sem > ################################################### > ability_model <- specifyModel(text = mod) Read 15 records NOTE: it is generally simpler to use specifyEquations() or cfa() see ?specifyEquations > ability_sem <- sem(ability_model, ability, 556) > ################################################### > ### code chunk number 6: ch:SEM:ability-sem > ################################################### > mod [1] "Ability -> SCA, lambda1, NA" "Ability -> PPE, lambda2, NA" [3] "Ability -> PTE, lambda3, NA" "Ability -> PFE, lambda4, NA" [5] "Aspiration -> EA, lambda5, NA" "Aspiration -> CP, lambda6, NA" [7] "Ability <-> Aspiration, rho, NA" "SCA <-> SCA, theta1, NA" [9] "PPE <-> PPE, theta2, NA" "PTE <-> PTE, theta3, NA" [11] "PFE <-> PFE, theta4, NA" "EA <-> EA, theta5, NA" [13] "CP <-> CP, theta6, NA" "Ability <-> Ability, NA, 1" [15] "Aspiration <-> Aspiration, NA, 1" > ################################################### > ### code chunk number 7: ch:SEM:ability-summary > ################################################### > summary(ability_sem) Model Chisquare = 9.3 Df = 8 Pr(>Chisq) = 0.32 AIC = 35 BIC = -41 Normalized Residuals Min. 1st Qu. Median Mean 3rd Qu. Max. -4.41e-01 -1.87e-01 -1.80e-06 -1.31e-02 2.11e-01 5.33e-01 R-square for Endogenous Variables SCA PPE PTE PFE EA CP 0.75 0.72 0.65 0.48 0.60 0.86 Parameter Estimates Estimate Std Error z value Pr(>|z|) lambda1 0.86 0.035 24.6 3.3e-133 SCA <--- Ability lambda2 0.85 0.035 24.0 7.6e-127 PPE <--- Ability lambda3 0.81 0.036 22.1 2.3e-108 PTE <--- Ability lambda4 0.70 0.039 18.0 2.1e-72 PFE <--- Ability lambda5 0.78 0.040 19.2 3.3e-82 EA <--- Aspiration lambda6 0.93 0.039 23.6 7.6e-123 CP <--- Aspiration rho 0.67 0.031 21.5 8.6e-103 Aspiration <--> Ability theta1 0.25 0.023 10.9 1.1e-27 SCA <--> SCA theta2 0.28 0.024 11.5 7.5e-31 PPE <--> PPE theta3 0.35 0.027 13.1 4.9e-39 PTE <--> PTE theta4 0.52 0.035 14.9 4.7e-50 PFE <--> PFE theta5 0.40 0.038 10.5 1.4e-25 EA <--> EA theta6 0.14 0.044 3.2 1.6e-03 CP <--> CP Iterations = 29 > ################################################### > ### code chunk number 8: ch:SEM:ability-summary > ################################################### > su <- summary(ability_sem) > ################################################### > ### code chunk number 9: ch:SEM:ability-path > ################################################### > pathDiagram(ability_sem, file = "ability_sem", + ignore.double = FALSE, edge.labels = "both", output.type = "graphics") Running dot -Tpdf -o ability_sem.pdf -Gcharset=latin1 ability_sem.dot > ################################################### > ### code chunk number 10: ch:SEM:ability-files > ################################################### > file.remove("ability_sem.dot") [1] TRUE > ################################################### > ### code chunk number 11: ch:SEM:druguse-model > ################################################### > mod <- c("Alcohol -> Cigs, lambda1, NA", + "Alcohol -> Beer, lambda3, NA", + "Alcohol -> Wine, lambda4, NA", + "Alcohol -> Liqr, lambda6, NA", + "Cannabis -> Cigs, lambda2, NA", + "Cannabis -> Wine, lambda5, NA", + "Cannabis -> Marj, lambda12, NA", + "Cannabis -> Hash, lambda13, NA", + "Hard -> Liqr, lambda7, NA", + "Hard -> Cocn, lambda8, NA", + "Hard -> Tran, lambda9, NA", + "Hard -> Drug, lambda10, NA", + "Hard -> Hern, lambda11, NA", + "Hard -> Hash, lambda14, NA", + "Hard -> Inhl, lambda15, NA", + "Hard -> Hall, lambda16, NA", + "Hard -> Amph, lambda17, NA", + "Cigs <-> Cigs, theta1, NA", + "Beer <-> Beer, theta2, NA", + "Wine <-> Wine, theta3, NA", + "Liqr <-> Liqr, theta4, NA", + "Cocn <-> Cocn, theta5, NA", + "Tran <-> Tran, theta6, NA", + "Drug <-> Drug, theta7, NA", + "Hern <-> Hern, theta8, NA", + "Marj <-> Marj, theta9, NA", + "Hash <-> Hash, theta10, NA", + "Inhl <-> Inhl, theta11, NA", + "Hall <-> Hall, theta12, NA", + "Amph <-> Amph, theta13, NA", + "Alcohol <-> Alcohol, NA, 1", + "Cannabis <-> Cannabis, NA, 1", + "Hard <-> Hard, NA, 1", + "Alcohol <-> Cannabis, rho1, NA", + "Alcohol <-> Hard, rho2, NA", + "Cannabis <-> Hard, rho3, NA") > ################################################### > ### code chunk number 12: ch:SEM:druguse-names > ################################################### > rownames(druguse) <- colnames(druguse) <- c("Cigs", + "Beer", "Wine", "Liqr", "Cocn", "Tran", "Drug", + "Hern", "Marj", "Hash", "Inhl", "Hall", "Amph") > ################################################### > ### code chunk number 13: ch:SEM:druguse-sem > ################################################### > druguse_model <- specifyModel(text = mod) Read 36 records NOTE: it is generally simpler to use specifyEquations() or cfa() see ?specifyEquations > druguse_sem <- sem(druguse_model, druguse, 1634) > ################################################### > ### code chunk number 14: ch:SEM:druguse-sem > ################################################### > mod [1] "Alcohol -> Cigs, lambda1, NA" "Alcohol -> Beer, lambda3, NA" [3] "Alcohol -> Wine, lambda4, NA" "Alcohol -> Liqr, lambda6, NA" [5] "Cannabis -> Cigs, lambda2, NA" "Cannabis -> Wine, lambda5, NA" [7] "Cannabis -> Marj, lambda12, NA" "Cannabis -> Hash, lambda13, NA" [9] "Hard -> Liqr, lambda7, NA" "Hard -> Cocn, lambda8, NA" [11] "Hard -> Tran, lambda9, NA" "Hard -> Drug, lambda10, NA" [13] "Hard -> Hern, lambda11, NA" "Hard -> Hash, lambda14, NA" [15] "Hard -> Inhl, lambda15, NA" "Hard -> Hall, lambda16, NA" [17] "Hard -> Amph, lambda17, NA" "Cigs <-> Cigs, theta1, NA" [19] "Beer <-> Beer, theta2, NA" "Wine <-> Wine, theta3, NA" [21] "Liqr <-> Liqr, theta4, NA" "Cocn <-> Cocn, theta5, NA" [23] "Tran <-> Tran, theta6, NA" "Drug <-> Drug, theta7, NA" [25] "Hern <-> Hern, theta8, NA" "Marj <-> Marj, theta9, NA" [27] "Hash <-> Hash, theta10, NA" "Inhl <-> Inhl, theta11, NA" [29] "Hall <-> Hall, theta12, NA" "Amph <-> Amph, theta13, NA" [31] "Alcohol <-> Alcohol, NA, 1" "Cannabis <-> Cannabis, NA, 1" [33] "Hard <-> Hard, NA, 1" "Alcohol <-> Cannabis, rho1, NA" [35] "Alcohol <-> Hard, rho2, NA" "Cannabis <-> Hard, rho3, NA" > ################################################### > ### code chunk number 15: ch:SEM:druguse-path > ################################################### > pathDiagram(druguse_sem, file = "druguse_sem", + ignore.double = FALSE, edge.labels = "both", output.type = "graphics") Running dot -Tpdf -o druguse_sem.pdf -Gcharset=latin1 druguse_sem.dot > ################################################### > ### code chunk number 16: ch:SEM:druguse-files > ################################################### > file.remove("druguse_sem.dot") [1] TRUE > ################################################### > ### code chunk number 17: ch:SEM:druguse-summary > ################################################### > summary(druguse_sem) Model Chisquare = 324 Df = 58 Pr(>Chisq) = 1.2e-38 AIC = 390 BIC = -105 Normalized Residuals Min. 1st Qu. Median Mean 3rd Qu. Max. -3.045264 -0.880180 -0.000004 -0.021661 0.999052 4.577128 R-square for Endogenous Variables Cigs Beer Wine Liqr Marj Hash Cocn Tran Drug Hern Inhl Hall Amph 0.39 0.63 0.62 0.59 0.83 0.45 0.22 0.46 0.13 0.23 0.29 0.38 0.58 Parameter Estimates Estimate Std Error z value Pr(>|z|) lambda1 0.36 0.034 10.4 3.5e-25 Cigs <--- Alcohol lambda3 0.79 0.023 35.0 1.0e-268 Beer <--- Alcohol lambda4 0.88 0.038 23.3 3.1e-120 Wine <--- Alcohol lambda6 0.72 0.024 30.7 2.2e-206 Liqr <--- Alcohol lambda2 0.33 0.035 9.4 4.8e-21 Cigs <--- Cannabis lambda5 -0.15 0.037 -4.2 3.2e-05 Wine <--- Cannabis lambda12 0.91 0.030 29.9 5.7e-197 Marj <--- Cannabis lambda13 0.40 0.030 13.4 9.8e-41 Hash <--- Cannabis lambda7 0.12 0.023 5.5 5.0e-08 Liqr <--- Hard lambda8 0.46 0.026 18.1 6.4e-73 Cocn <--- Hard lambda9 0.68 0.024 28.2 9.3e-175 Tran <--- Hard lambda10 0.36 0.026 13.6 4.5e-42 Drug <--- Hard lambda11 0.48 0.026 18.6 7.4e-77 Hern <--- Hard lambda14 0.38 0.029 13.1 4.6e-39 Hash <--- Hard lambda15 0.54 0.025 21.6 2.0e-103 Inhl <--- Hard lambda16 0.62 0.024 25.2 1.6e-140 Hall <--- Hard lambda17 0.76 0.023 33.0 7.6e-239 Amph <--- Hard theta1 0.61 0.024 25.8 4.7e-147 Cigs <--> Cigs theta2 0.37 0.020 18.7 2.6e-78 Beer <--> Beer theta3 0.38 0.024 16.0 5.9e-58 Wine <--> Wine theta4 0.41 0.019 21.3 4.3e-101 Liqr <--> Liqr theta5 0.78 0.029 26.8 8.5e-159 Cocn <--> Cocn theta6 0.54 0.023 23.2 2.6e-119 Tran <--> Tran theta7 0.87 0.032 27.7 2.4e-168 Drug <--> Drug theta8 0.77 0.029 26.7 1.6e-157 Hern <--> Hern theta9 0.17 0.044 3.8 1.4e-04 Marj <--> Marj theta10 0.55 0.022 24.6 1.3e-133 Hash <--> Hash theta11 0.71 0.027 25.9 2.1e-148 Inhl <--> Inhl theta12 0.62 0.025 24.7 3.1e-134 Hall <--> Hall theta13 0.42 0.021 19.7 2.3e-86 Amph <--> Amph rho1 0.63 0.027 23.3 2.0e-120 Cannabis <--> Alcohol rho2 0.31 0.029 10.7 1.3e-26 Hard <--> Alcohol rho3 0.50 0.027 18.4 1.3e-75 Hard <--> Cannabis Iterations = 34 > ################################################### > ### code chunk number 18: ch:SEM:druguse-summary > ################################################### > su <- summary(druguse_sem) > ################################################### > ### code chunk number 19: ch:SEM:druguse-cov > ################################################### > round(druguse_sem$S - druguse_sem$C, 3) Cigs Beer Wine Liqr Cocn Tran Drug Hern Marj Hash Cigs 0.000 -0.002 0.010 -0.009 -0.015 0.015 -0.009 -0.050 0.003 -0.023 Beer -0.002 0.000 0.002 0.002 -0.047 -0.021 0.014 -0.055 -0.012 0.025 Wine 0.010 0.002 0.000 -0.004 -0.039 0.005 0.039 -0.028 -0.002 0.005 Liqr -0.009 0.002 -0.004 0.000 -0.047 0.022 -0.003 -0.069 0.009 0.029 Cocn -0.015 -0.047 -0.039 -0.047 0.000 0.035 0.042 0.100 -0.026 0.034 Tran 0.015 -0.021 0.005 0.022 0.035 0.000 -0.021 0.034 0.007 -0.014 Drug -0.009 0.014 0.039 -0.003 0.042 -0.021 0.000 0.030 -0.013 -0.045 Hern -0.050 -0.055 -0.028 -0.069 0.100 0.034 0.030 0.000 -0.063 -0.057 Marj 0.003 -0.012 -0.002 0.009 -0.026 0.007 -0.013 -0.063 0.000 -0.001 Hash -0.023 0.025 0.005 0.029 0.034 -0.014 -0.045 -0.057 -0.001 0.000 Inhl 0.094 0.068 0.075 0.065 0.020 -0.044 0.115 0.030 0.054 -0.013 Hall -0.071 -0.065 -0.049 -0.077 -0.008 -0.051 0.010 0.026 -0.077 0.010 Amph 0.033 0.010 0.032 0.026 -0.077 0.029 -0.042 -0.049 0.047 0.025 Inhl Hall Amph Cigs 0.094 -0.071 0.033 Beer 0.068 -0.065 0.010 Wine 0.075 -0.049 0.032 Liqr 0.065 -0.077 0.026 Cocn 0.020 -0.008 -0.077 Tran -0.044 -0.051 0.029 Drug 0.115 0.010 -0.042 Hern 0.030 0.026 -0.049 Marj 0.054 -0.077 0.047 Hash -0.013 0.010 0.025 Inhl 0.000 0.004 -0.022 Hall 0.004 0.000 0.039 Amph -0.022 0.039 0.000 > ################################################### > ### code chunk number 20: ch:SEM:alienation:plot > ################################################### > a <- cov2cor(alienation) > ord <- order.dendrogram(as.dendrogram(hclust(dist(a)))) > panel.corrgram <- + function(x, y, z, subscripts, at, + level = 0.9, label = FALSE, ...) + { + require("ellipse", quietly = TRUE) + x <- as.numeric(x)[subscripts] + y <- as.numeric(y)[subscripts] + z <- as.numeric(z)[subscripts] + zcol <- level.colors(z, at = at, col.regions = grey.colors, ...) + for (i in seq(along = z)) { + ell <- ellipse(z[i], level = level, npoints = 50, + scale = c(.2, .2), centre = c(x[i], y[i])) + panel.polygon(ell, col = zcol[i], border = zcol[i], ...) + } + if (label) + panel.text(x = x, y = y, lab = 100 * round(z, 2), cex = 0.8, + col = ifelse(z < 0, "white", "black")) + } > print(levelplot(a[ord, ord], at = do.breaks(c(-1.01, 1.01), 20), + xlab = NULL, ylab = NULL, colorkey = list(space = "top"), + scales = list(x = list(rot = 90)), + panel = panel.corrgram, label = TRUE)) > ################################################### > ### code chunk number 21: ch:SEM:alienation-model > ################################################### > mod <- c("SES -> Educ, NA, 1", + "SES -> SEI, lambda1, NA", + "Alienation67 -> Anomia67, NA, 1", + "Alienation67 -> Powles67, lambda2, NA", + "Alienation71 -> Anomia71, NA, 1", + "Alienation71 -> Powles71, lambda3, NA", + "SES -> Alienation67, beta1, NA", + "SES -> Alienation71, beta2, NA", + "Alienation67 -> Alienation71, beta3, NA", + "Educ <-> Educ, theta1, NA", + "SEI <-> SEI, theta2, NA", + "SES <-> SES, delta0, NA", + "Anomia67 <-> Anomia67, theta3, NA", + "Powles67 <-> Powles67, theta4, NA", + "Anomia71 <-> Anomia71, theta5, NA", + "Powles71 <-> Powles71, theta6, NA", + "Alienation67 <-> Alienation67, delta1, NA", + "Alienation71 <-> Alienation71, delta2, NA") > mod2 <- c(mod, "Anomia67 <-> Anomia71,psi,NA") > ################################################### > ### code chunk number 22: ch:SEM:alienation-sem > ################################################### > alienation_model <- specifyModel(text = mod) Read 18 records NOTE: it is generally simpler to use specifyEquations() or cfa() see ?specifyEquations > alienation_sem <- sem(alienation_model, alienation, 932) > ################################################### > ### code chunk number 23: ch:SEM:alienation-sem > ################################################### > mod [1] "SES -> Educ, NA, 1" [2] "SES -> SEI, lambda1, NA" [3] "Alienation67 -> Anomia67, NA, 1" [4] "Alienation67 -> Powles67, lambda2, NA" [5] "Alienation71 -> Anomia71, NA, 1" [6] "Alienation71 -> Powles71, lambda3, NA" [7] "SES -> Alienation67, beta1, NA" [8] "SES -> Alienation71, beta2, NA" [9] "Alienation67 -> Alienation71, beta3, NA" [10] "Educ <-> Educ, theta1, NA" [11] "SEI <-> SEI, theta2, NA" [12] "SES <-> SES, delta0, NA" [13] "Anomia67 <-> Anomia67, theta3, NA" [14] "Powles67 <-> Powles67, theta4, NA" [15] "Anomia71 <-> Anomia71, theta5, NA" [16] "Powles71 <-> Powles71, theta6, NA" [17] "Alienation67 <-> Alienation67, delta1, NA" [18] "Alienation71 <-> Alienation71, delta2, NA" > ################################################### > ### code chunk number 24: ch:SEM:alienation-path > ################################################### > pathDiagram(alienation_sem, file = "alienation_sem", + ignore.double = FALSE, edge.labels = "both", output.type = "graphics") Running dot -Tpdf -o alienation_sem.pdf -Gcharset=latin1 alienation_sem.dot > ################################################### > ### code chunk number 25: ch:SEM:alienation-files > ################################################### > file.remove("alienation_sem.dot") [1] TRUE > ################################################### > ### code chunk number 26: ch:SEM:alienation-summary > ################################################### > summary(alienation_sem) Model Chisquare = 72 Df = 6 Pr(>Chisq) = 2e-13 AIC = 102 BIC = 31 Normalized Residuals Min. 1st Qu. Median Mean 3rd Qu. Max. -1.26e+00 -2.09e-01 -6.90e-06 -1.51e-02 2.44e-01 1.33e+00 R-square for Endogenous Variables Educ SEI Alienation67 Anomia67 Powles67 Alienation71 0.69 0.42 0.32 0.66 0.66 0.58 Anomia71 Powles71 0.70 0.64 Parameter Estimates Estimate Std Error z value Pr(>|z|) lambda1 5.33 0.430 12.4 2.5e-35 SEI <--- SES lambda2 0.89 0.042 21.4 1.3e-101 Powles67 <--- Alienation67 lambda3 0.85 0.040 21.2 3.9e-100 Powles71 <--- Alienation71 beta1 -0.61 0.056 -10.9 1.6e-27 Alienation67 <--- SES beta2 -0.17 0.054 -3.2 1.2e-03 Alienation71 <--- SES beta3 0.70 0.054 13.2 1.4e-39 Alienation71 <--- Alienation67 theta1 2.94 0.499 5.9 4.1e-09 Educ <--> Educ theta2 260.93 18.239 14.3 2.0e-46 SEI <--> SEI delta0 6.66 0.641 10.4 2.4e-25 SES <--> SES theta3 4.02 0.343 11.7 1.1e-31 Anomia67 <--> Anomia67 theta4 3.19 0.272 11.7 7.5e-32 Powles67 <--> Powles67 theta5 3.70 0.373 9.9 3.5e-23 Anomia71 <--> Anomia71 theta6 3.62 0.292 12.4 2.4e-35 Powles71 <--> Powles71 delta1 5.31 0.473 11.2 3.0e-29 Alienation67 <--> Alienation67 delta2 3.74 0.387 9.7 4.8e-22 Alienation71 <--> Alienation71 Iterations = 86 > ################################################### > ### code chunk number 27: ch:SEM:alienation-summary > ################################################### > su <- summary(alienation_sem) > ################################################### > ### code chunk number 28: ch:SEM:alienation-sem2 > ################################################### > alienation_model2 <- specifyModel(text = mod) Read 18 records NOTE: it is generally simpler to use specifyEquations() or cfa() see ?specifyEquations > alienation_sem2 <- sem(alienation_model2, alienation, 932) > su <- summary(alienation_sem2) demo(Ch-Viz) ---- ~~~~~~ > ### R code from vignette source 'Ch-Viz.Rnw' > > ################################################### > ### code chunk number 1: setup > ################################################### > library("MVA") > set.seed(280875) > library("lattice") > lattice.options(default.theme = + function() + standard.theme("pdf", color = FALSE)) > if (file.exists("deparse.R")) { + if (!file.exists("figs")) dir.create("figs") + source("deparse.R") + options(prompt = "R> ", continue = "+ ", width = 64, + digits = 4, show.signif.stars = FALSE, useFancyQuotes = FALSE) + + options(SweaveHooks = list(onefig = function() {par(mfrow = c(1,1))}, + twofig = function() {par(mfrow = c(1,2))}, + figtwo = function() {par(mfrow = c(2,1))}, + threefig = function() {par(mfrow = c(1,3))}, + figthree = function() {par(mfrow = c(3,1))}, + fourfig = function() {par(mfrow = c(2,2))}, + sixfig = function() {par(mfrow = c(3,2))}, + nomar = function() par("mai" = c(0, 0, 0, 0)))) + } > ################################################### > ### code chunk number 2: ch:Viz:data > ################################################### > measure <- + structure(list(V1 = 1:20, V2 = c(34L, 37L, 38L, 36L, 38L, 43L, + 40L, 38L, 40L, 41L, 36L, 36L, 34L, 33L, 36L, 37L, 34L, 36L, 38L, + 35L), V3 = c(30L, 32L, 30L, 33L, 29L, 32L, 33L, 30L, 30L, 32L, + 24L, 25L, 24L, 22L, 26L, 26L, 25L, 26L, 28L, 23L), V4 = c(32L, + 37L, 36L, 39L, 33L, 38L, 42L, 40L, 37L, 39L, 35L, 37L, 37L, 34L, + 38L, 37L, 38L, 37L, 40L, 35L)), .Names = c("V1", "V2", "V3", + "V4"), class = "data.frame", row.names = c(NA, -20L)) > measure <- measure[,-1] > names(measure) <- c("chest", "waist", "hips") > measure$gender <- gl(2, 10) > levels(measure$gender) <- c("male", "female") > data("USairpollution", package = "HSAUR2") > ################################################### > ### code chunk number 3: ch:Viz:USairpollution:plot1mlab > ################################################### > mlab <- "Manufacturing enterprises with 20 or more workers" > plab <- "Population size (1970 census) in thousands" > ################################################### > ### code chunk number 4: ch:Viz:USairpollution:plot1setup (eval = FALSE) > ################################################### > ## plot(popul ~ manu, data = USairpollution, > ## xlab = mlab, ylab = plab) > > > ################################################### > ### code chunk number 5: ch:Viz:USairpollution:plot1 > ################################################### > plot(popul ~ manu, data = USairpollution, + xlab = mlab, ylab = plab) > ################################################### > ### code chunk number 6: ch:Viz:USairpollution:plot3setup > ################################################### > layout(matrix(c(2, 0, 1, 3), nrow = 2, byrow = TRUE), + widths = c(2, 1), heights = c(1, 2), respect = TRUE) > xlim <- with(USairpollution, range(manu)) * 1.1 > plot(popul ~ manu, data = USairpollution, cex.lab = 0.9, + xlab = mlab, ylab = plab, type = "n", xlim = xlim) > with(USairpollution, text(manu, popul, cex = 0.6, + labels = abbreviate(row.names(USairpollution)))) > with(USairpollution, hist(manu, main = "", xlim = xlim)) > with(USairpollution, boxplot(popul)) > ################################################### > ### code chunk number 7: ch:Viz:USairpollution:plot2 > ################################################### > plot(popul ~ manu, data = USairpollution, + xlab = mlab, ylab = plab) > rug(USairpollution$manu, side = 1) > rug(USairpollution$popul, side = 2) > ################################################### > ### code chunk number 8: ch:Viz:USairpollution:plot3 > ################################################### > layout(matrix(c(2, 0, 1, 3), nrow = 2, byrow = TRUE), + widths = c(2, 1), heights = c(1, 2), respect = TRUE) > xlim <- with(USairpollution, range(manu)) * 1.1 > plot(popul ~ manu, data = USairpollution, cex.lab = 0.9, + xlab = mlab, ylab = plab, type = "n", xlim = xlim) > with(USairpollution, text(manu, popul, cex = 0.6, + labels = abbreviate(row.names(USairpollution)))) > with(USairpollution, hist(manu, main = "", xlim = xlim)) > with(USairpollution, boxplot(popul)) > ################################################### > ### code chunk number 9: ch:Viz:USairpollution:plot4 > ################################################### > outcity <- match(lab <- c("Chicago", "Detroit", + "Cleveland", "Philadelphia"), rownames(USairpollution)) > x <- USairpollution[, c("manu", "popul")] > bvbox(x, mtitle = "", xlab = mlab, ylab = plab) > text(x$manu[outcity], x$popul[outcity], labels = lab, + cex = 0.7, pos = c(2, 2, 4, 2, 2)) > ################################################### > ### code chunk number 10: ch:Viz:USairpollution:cor > ################################################### > with(USairpollution, cor(manu, popul)) [1] 0.96 > outcity <- match(c("Chicago", "Detroit", + "Cleveland", "Philadelphia"), + rownames(USairpollution)) > with(USairpollution, cor(manu[-outcity], popul[-outcity])) [1] 0.8 > ################################################### > ### code chunk number 11: ch:Viz:USairpollution:chull > ################################################### > (hull <- with(USairpollution, chull(manu, popul))) [1] 9 15 41 6 2 18 16 14 7 > ################################################### > ### code chunk number 12: ch:Viz:USairpollution:chullplot > ################################################### > with(USairpollution, + plot(manu, popul, pch = 1, xlab = mlab, ylab = plab)) > with(USairpollution, + polygon(manu[hull], popul[hull], density = 15, angle = 30)) > ################################################### > ### code chunk number 13: ch:Viz:USairpollution:chullcor > ################################################### > with(USairpollution, cor(manu[-hull],popul[-hull])) [1] 0.92 > ################################################### > ### code chunk number 14: ch:Viz:USairpollution:chiplot:setup (eval = FALSE) > ################################################### > ## with(USairpollution, plot(manu, popul, > ## xlab = mlab, ylab = plab, > ## cex.lab = 0.9)) > ## with(USairpollution, chiplot(manu, popul)) > > > ################################################### > ### code chunk number 15: ch:Viz:USairpollution:chiplot > ################################################### > with(USairpollution, plot(manu, popul, + xlab = mlab, ylab = plab, + cex.lab = 0.9)) > with(USairpollution, chiplot(manu, popul)) > ################################################### > ### code chunk number 16: ch:Viz:USairpollution:plot5 > ################################################### > ylim <- with(USairpollution, range(wind)) * c(0.95, 1) > plot(wind ~ temp, data = USairpollution, + xlab = "Average annual temperature (Fahrenheit)", + ylab = "Average annual wind speed (m.p.h.)", pch = 10, + ylim = ylim) > with(USairpollution, symbols(temp, wind, circles = SO2, + inches = 0.5, add = TRUE)) > ################################################### > ### code chunk number 17: ch:Viz:USairpollution:plot6 > ################################################### > plot(wind ~ temp, data = USairpollution, + xlab = "Average annual temperature (Fahrenheit)", + ylab = "Average annual wind speed (m.p.h.)", pch = 10, + ylim = ylim) > with(USairpollution, + stars(USairpollution[,-c(2,5)], locations = cbind(temp, wind), + labels = NULL, add = TRUE, cex = 0.5)) > ################################################### > ### code chunk number 18: ch:Viz:USairpollution:plot7 > ################################################### > stars(USairpollution, cex = 0.55) > ################################################### > ### code chunk number 19: ch:Viz:USairpollution:plot8 > ################################################### > pairs(USairpollution, pch = ".", cex = 1.5) > ################################################### > ### code chunk number 20: ch:Viz:USairpollution:plot9 > ################################################### > pairs(USairpollution, + panel = function (x, y, ...) { + points(x, y, ...) + abline(lm(y ~ x), col = "grey") + }, pch = ".", cex = 1.5) > ################################################### > ### code chunk number 21: ch:Viz:USairpollution:cor > ################################################### > round(cor(USairpollution), 4) SO2 temp manu popul wind precip predays SO2 1.000 -0.434 0.645 0.494 0.095 0.054 0.370 temp -0.434 1.000 -0.190 -0.063 -0.350 0.386 -0.430 manu 0.645 -0.190 1.000 0.955 0.238 -0.032 0.132 popul 0.494 -0.063 0.955 1.000 0.213 -0.026 0.042 wind 0.095 -0.350 0.238 0.213 1.000 -0.013 0.164 precip 0.054 0.386 -0.032 -0.026 -0.013 1.000 0.496 predays 0.370 -0.430 0.132 0.042 0.164 0.496 1.000 > ################################################### > ### code chunk number 22: ch:Viz-kernel-figs > ################################################### > rec <- function(x) (abs(x) < 1) * 0.5 > tri <- function(x) (abs(x) < 1) * (1 - abs(x)) > gauss <- function(x) 1/sqrt(2*pi) * exp(-(x^2)/2) > x <- seq(from = -3, to = 3, by = 0.001) > plot(x, rec(x), type = "l", ylim = c(0,1), lty = 1, + ylab = expression(K(x))) > lines(x, tri(x), lty = 2) > lines(x, gauss(x), lty = 3) > legend("topleft", legend = c("Rectangular", "Triangular", + "Gaussian"), lty = 1:3, title = "kernel functions", + bty = "n") > ################################################### > ### code chunk number 23: ch:Viz-x-bumps-data > ################################################### > x <- c(0, 1, 1.1, 1.5, 1.9, 2.8, 2.9, 3.5) > n <- length(x) > ################################################### > ### code chunk number 24: ch:Viz-x-bumps-gaussian > ################################################### > xgrid <- seq(from = min(x) - 1, to = max(x) + 1, by = 0.01) > ################################################### > ### code chunk number 25: ch:Viz-x-bumps-bumps > ################################################### > h <- 0.4 > bumps <- sapply(x, function(a) gauss((xgrid - a)/h)/(n * h)) > ################################################### > ### code chunk number 26: ch:Viz-x-bumps-setup (eval = FALSE) > ################################################### > ## plot(xgrid, rowSums(bumps), ylab = expression(hat(f)(x)), > ## type = "l", xlab = "x", lwd = 2) > ## rug(x, lwd = 2) > ## out <- apply(bumps, 2, function(b) lines(xgrid, b)) > > > ################################################### > ### code chunk number 27: ch:Viz-x-bumps > ################################################### > plot(xgrid, rowSums(bumps), ylab = expression(hat(f)(x)), + type = "l", xlab = "x", lwd = 2) > rug(x, lwd = 2) > out <- apply(bumps, 2, function(b) lines(xgrid, b)) > ################################################### > ### code chunk number 28: ch:Viz-epakernel-fig > ################################################### > epa <- function(x, y) + ((x^2 + y^2) < 1) * 2/pi * (1 - x^2 - y^2) > x <- seq(from = -1.1, to = 1.1, by = 0.05) > epavals <- sapply(x, function(a) epa(a, x)) > persp(x = x, y = x, z = epavals, xlab = "x", ylab = "y", + zlab = expression(K(x, y)), theta = -35, axes = TRUE, + box = TRUE) > ################################################### > ### code chunk number 29: ch:Viz:CYGOB1:tab > ################################################### > toLatex(HSAURtable(CYGOB1), pcol = 3, + caption = "Energy output and surface temperature of star cluster CYG OB1.", + label = "ch:Viz:CYGOB1:tab") \index{CYGOB1 data@\Robject{CYGOB1} data} \begin{center} \begin{longtable}{ rr|rr|rr } \caption{\Robject{CYGOB1} data. Energy output and surface temperature of star cluster CYG OB1. \label{ch:Viz:CYGOB1:tab}} \\ \hline \Robject{logst} & \Robject{logli} & \Robject{logst} & \Robject{logli} & \Robject{logst} & \Robject{logli} \\ \hline \endfirsthead \caption[]{\Robject{CYGOB1} data (continued).} \\ \hline \Robject{logst} & \Robject{logli} & \Robject{logst} & \Robject{logli} & \Robject{logst} & \Robject{logli} \\ \hline \endhead 4.37 & 5.23 & 4.23 & 3.94 & 4.45 & 5.22 \\ 4.56 & 5.74 & 4.42 & 4.18 & 3.49 & 6.29 \\ 4.26 & 4.93 & 4.23 & 4.18 & 4.23 & 4.34 \\ 4.56 & 5.74 & 3.49 & 5.89 & 4.62 & 5.62 \\ 4.30 & 5.19 & 4.29 & 4.38 & 4.53 & 5.10 \\ 4.46 & 5.46 & 4.29 & 4.22 & 4.45 & 5.22 \\ 3.84 & 4.65 & 4.42 & 4.42 & 4.53 & 5.18 \\ 4.57 & 5.27 & 4.49 & 4.85 & 4.43 & 5.57 \\ 4.26 & 5.57 & 4.38 & 5.02 & 4.38 & 4.62 \\ 4.37 & 5.12 & 4.42 & 4.66 & 4.45 & 5.06 \\ 3.49 & 5.73 & 4.29 & 4.66 & 4.50 & 5.34 \\ 4.43 & 5.45 & 4.38 & 4.90 & 4.45 & 5.34 \\ 4.48 & 5.42 & 4.22 & 4.39 & 4.55 & 5.54 \\ 4.01 & 4.05 & 3.48 & 6.05 & 4.45 & 4.98 \\ 4.29 & 4.26 & 4.38 & 4.42 & 4.42 & 4.50 \\ 4.42 & 4.58 & 4.56 & 5.10 & & \\ \hline \end{longtable} \end{center} > ################################################### > ### code chunk number 30: ch:Viz:CYGOB1:plot1 > ################################################### > library("KernSmooth") KernSmooth 2.23 loaded Copyright M. P. Wand 1997-2009 > CYGOB1d <- bkde2D(CYGOB1, bandwidth = sapply(CYGOB1, dpik)) > plot(CYGOB1, xlab = "log surface temperature", + ylab = "log light intensity") > contour(x = CYGOB1d$x1, y = CYGOB1d$x2, + z = CYGOB1d$fhat, add = TRUE) > ################################################### > ### code chunk number 31: ch:Viz:CYGOB1:plot2 > ################################################### > persp(x = CYGOB1d$x1, y = CYGOB1d$x2, z = CYGOB1d$fhat, + xlab = "log surface temperature", + ylab = "log light intensity", + zlab = "density") > ################################################### > ### code chunk number 32: ch:Viz:measure:plot1:setup (eval = FALSE) > ################################################### > ## panel.hist <- function(x, ...) > ## { > ## usr <- par("usr"); on.exit(par(usr)) > ## par(usr = c(usr[1:2], 0, 1.5) ) > ## h <- hist(x, plot = FALSE) > ## breaks <- h$breaks; nB <- length(breaks) > ## y <- h$counts; y <- y/max(y) > ## rect(breaks[-nB], 0, breaks[-1], y, col = "grey", ...) > ## } > ## pairs(measure[, c("chest", "waist", "hips")], > ## diag.panel = panel.hist, > ## panel = function (x,y) { > ## data <- data.frame(cbind(x,y)) > ## par(new = TRUE) > ## den <- bkde2D(data, bandwidth = sapply(data, dpik)) > ## contour(x = den$x1, y = den$x2, > ## z = den$fhat, axes = FALSE) > ## }) > > > ################################################### > ### code chunk number 33: ch:Viz:measure:plot1 > ################################################### > panel.hist <- function(x, ...) + { + usr <- par("usr"); on.exit(par(usr)) + par(usr = c(usr[1:2], 0, 1.5) ) + h <- hist(x, plot = FALSE) + breaks <- h$breaks; nB <- length(breaks) + y <- h$counts; y <- y/max(y) + rect(breaks[-nB], 0, breaks[-1], y, col = "grey", ...) + } > pairs(measure[, c("chest", "waist", "hips")], + diag.panel = panel.hist, + panel = function (x,y) { + data <- data.frame(cbind(x,y)) + par(new = TRUE) + den <- bkde2D(data, bandwidth = sapply(data, dpik)) + contour(x = den$x1, y = den$x2, + z = den$fhat, axes = FALSE) + }) > ################################################### > ### code chunk number 34: ch:Viz:measure:plot2 > ################################################### > library("scatterplot3d") > with(measure, scatterplot3d(chest, waist, hips, + pch = (1:2)[gender], type = "h", angle = 55)) > ################################################### > ### code chunk number 35: ch:Viz:USairpollution:plot10 > ################################################### > with(USairpollution, + scatterplot3d(temp, wind, SO2, type = "h", + angle = 55)) > ################################################### > ### code chunk number 36: ch:Viz:USairpollution:plot11 > ################################################### > plot(xyplot(SO2 ~ temp| cut(wind, 2), data = USairpollution)) > ################################################### > ### code chunk number 37: ch:Viz:USairpollution:plot12 > ################################################### > pollution <- with(USairpollution, equal.count(SO2,4)) > plot(cloud(precip ~ temp * wind | pollution, panel.aspect = 0.9, + data = USairpollution)) > ################################################### > ### code chunk number 38: ch:Viz:quakes:plot1 > ################################################### > plot(xyplot(lat ~ long| cut(depth, 3), data = quakes, + layout = c(3, 1), xlab = "Longitude", + ylab = "Latitude")) > ################################################### > ### code chunk number 39: ch:Viz:quakes:plot2 > ################################################### > Magnitude <- with(quakes, equal.count(mag, 4)) > depth.ord <- with(quakes, rev(order(depth))) > quakes.ordered <- quakes[depth.ord,] > depth.breaks <- with(quakes.ordered, do.breaks(range(depth),50)) > quakes.ordered$color<-level.colors(quakes.ordered$depth,at=depth.breaks, + col.regions=grey.colors) > plot(xyplot(lat ~ long | Magnitude, data = quakes.ordered, + aspect = "iso", groups = color, cex = 2, col = "black", + panel = function(x, y, groups, ..., subscripts) { + fill <- groups[subscripts] + panel.grid(h = -1, v = -1) + panel.xyplot(x, y, pch = 21, fill = fill, ...) + }, + legend = + list(right = + list(fun = draw.colorkey, + args = list(key = list(col = gray.colors, + at = depth.breaks), + draw = FALSE))), + xlab = "Longitude", ylab = "Latitude")) > ################################################### > ### code chunk number 40: ch:Viz:quakes:plot3 > ################################################### > plot(cloud(depth ~ lat * long | Magnitude, data = quakes, + zlim = rev(range(quakes$depth)), + screen = list(z = 105, x = -70), panel.aspect = 0.9, + xlab = "Longitude", ylab = "Latitude", zlab = "Depth")) > ################################################### > ### code chunk number 41: ch:Viz:USairpollution:stalac > ################################################### > stalac(USairpollution) Ch-CA Ch-EFA Ch-LME Ch-MDS Ch-MVA Ch-PCA Ch-SEM value NULL list,2 NULL Latex,33 NULL eigen,2 summary.objectiveML,34 visible FALSE FALSE FALSE TRUE FALSE TRUE FALSE Ch-Viz value NULL visible FALSE There were 12 warnings (use warnings() to see them) > > warnings() Warning messages: 1: In par(usr) : argument 1 does not name a graphical parameter 2: In par(usr) : argument 1 does not name a graphical parameter 3: In par(usr) : argument 1 does not name a graphical parameter 4: In par(usr) : argument 1 does not name a graphical parameter 5: In par(usr) : argument 1 does not name a graphical parameter 6: In par(usr) : argument 1 does not name a graphical parameter 7: In system(cmd) : 'dot' not found 8: In system(cmd) : 'dot' not found 9: In system(cmd) : 'dot' not found 10: In par(usr) : argument 1 does not name a graphical parameter 11: In par(usr) : argument 1 does not name a graphical parameter 12: In par(usr) : argument 1 does not name a graphical parameter > > proc.time() user system elapsed 5.17 0.51 5.67