Last updated on 2025-12-06 00:48:50 CET.
| Flavor | Version | Tinstall | Tcheck | Ttotal | Status | Flags |
|---|---|---|---|---|---|---|
| r-devel-linux-x86_64-debian-clang | 1.0-9 | 2.11 | 41.37 | 43.48 | OK | |
| r-devel-linux-x86_64-debian-gcc | 1.0-9 | 1.72 | 24.72 | 26.44 | ERROR | |
| r-devel-linux-x86_64-fedora-clang | 1.0-9 | 64.30 | OK | |||
| r-devel-linux-x86_64-fedora-gcc | 1.0-9 | 12.00 | 71.02 | 83.02 | OK | |
| r-devel-windows-x86_64 | 1.0-9 | 5.00 | 74.00 | 79.00 | OK | |
| r-patched-linux-x86_64 | 1.0-9 | 2.15 | 36.57 | 38.72 | OK | |
| r-release-linux-x86_64 | 1.0-9 | 2.11 | 39.14 | 41.25 | OK | |
| r-release-macos-arm64 | 1.0-9 | OK | ||||
| r-release-macos-x86_64 | 1.0-9 | 2.00 | 65.00 | 67.00 | OK | |
| r-release-windows-x86_64 | 1.0-9 | 4.00 | 74.00 | 78.00 | OK | |
| r-oldrel-macos-arm64 | 1.0-9 | OK | ||||
| r-oldrel-macos-x86_64 | 1.0-9 | 3.00 | 69.00 | 72.00 | OK | |
| r-oldrel-windows-x86_64 | 1.0-9 | 6.00 | 82.00 | 88.00 | OK |
Version: 1.0-9
Check: tests
Result: ERROR
Running ‘demo-test.R’ [2s/2s]
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.2
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")
> 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")
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