CRAN Package Check Results for Package MVA

Last updated on 2024-06-13 16:49:15 CEST.

Flavor Version Tinstall Tcheck Ttotal Status Flags
r-devel-linux-x86_64-debian-clang 1.0-8 3.02 46.23 49.25 OK
r-devel-linux-x86_64-debian-gcc 1.0-8 2.16 27.47 29.63 ERROR
r-devel-linux-x86_64-fedora-clang 1.0-8 61.79 OK
r-devel-linux-x86_64-fedora-gcc 1.0-8 70.75 OK
r-devel-windows-x86_64 1.0-8 4.00 69.00 73.00 OK
r-patched-linux-x86_64 1.0-8 2.40 45.04 47.44 OK
r-release-linux-x86_64 1.0-8 2.84 43.90 46.74 OK
r-release-macos-arm64 1.0-8 29.00 OK
r-release-macos-x86_64 1.0-8 43.00 OK
r-release-windows-x86_64 1.0-8 3.00 72.00 75.00 OK
r-oldrel-macos-arm64 1.0-8 32.00 OK
r-oldrel-macos-x86_64 1.0-8 46.00 OK
r-oldrel-windows-x86_64 1.0-8 4.00 71.00 75.00 OK

Check Details

Version: 1.0-8
Check: package dependencies
Result: NOTE Packages suggested but not available for checking: 'RLRsim', 'sem' Flavor: r-devel-linux-x86_64-debian-gcc

Version: 1.0-8
Check: tests
Result: ERROR Running ‘demo-test.R’ [3s/5s] Running the tests in ‘tests/demo-test.R’ failed. Complete output: > > 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)) > ################################################### > ### 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") Error in library("RLRsim") : there is no package called 'RLRsim' Calls: sapply ... FUN -> source -> withVisible -> eval -> eval -> library Execution halted Flavor: r-devel-linux-x86_64-debian-gcc