我mult.chart在 Rmarkdown 文件中使用 SPC 来证明概念。我只想打印图表并省略所有分解、xmv、协方差和 t2。当我使用
t <- mult.chart(na.omit(test.data), type = "t2", Xmv = Xmv, S = S, colm = colm)
该对象具有除图表之外的所有内容。
> str(t)
List of 5
$ : chr "Hotelling Control Chart"
$ ucl : num 13.8
$ t2 : num [1:154, 1] 6.1 1.11 3.13 0.66 2.26 2.13 2.02 3.45 4.17 2.41 ...
$ Xmv : num [1:4] 130.9 94.8 957.4 490.1
$ covariance: num [1:4, 1:4] 320 11 130 1000 11 0.87 4.9 32 130 4.9 ...
如何从中提取图表？
我更新了函数的代码以将 ggplot 图表添加到输出中。为了大家的利益，我发布下面的代码。
mult.chart2 <- function (type = c("chi", "t2", "mewma", "mcusum",
"mcusum2"), x, Xmv, S, colm, alpha = 0.01, lambda = 0.1,
k = 0.5, h = 5.5, phase = 1, method = "sw", ...)
{
type <- match.arg(type)
p <- ncol(x)
m <- nrow(x)
if (class(x) == "matrix" || class(x) == "data.frame")
(x <- array(data.matrix(x), c(m, p, 1)))
n <- dim(x)[3]
if (!missing(Xmv))
(phase <- 2)
x.jk <- matrix(0, m, p)
t2 <- matrix(0, m, 1)
x.jk <- apply(x, 1:2, mean)
if (missing(Xmv))
(Xmv <- colMeans(x.jk))
if (missing(S))
(S <- covariance(x, method = method))
if (missing(colm))
(colm <- nrow(x))
if (type == "chi") {
name <- paste("Chi-squared Control Chart")
for (ii in 1:m) {
t2[ii, 1] <- n * t(x.jk[ii, ] - Xmv) %*% solve(S) %*%
(x.jk[ii, ] - Xmv)
}
ucl <- qchisq(1 - alpha, p)
if (any(t2 > ucl)) {
cat("The following(s) point(s) fall outside the control limits")
t3 <- which(t2 > ucl)
print(t3)
}
}
if (type == "t2") {
name <- paste("Hotelling Control Chart")
for (ii in 1:m) {
t2[ii, 1] <- n * t(x.jk[ii, ] - Xmv) %*% solve(S) %*%
(x.jk[ii, ] - Xmv)
}
ifelse(n == 1, ifelse(phase == 1, ucl <- ((colm - 1)^2)/colm *
qbeta(1 - alpha, p/2, ((colm - p - 1)/2)), ucl <- ((p *
(colm + 1) * (colm - 1))/((colm^2) - colm * p)) *
qf(1 - alpha, p, colm - p)), ifelse(phase == 1, ucl <- (p *
(colm - 1) * (n - 1))/(colm * n - colm - p + 1) *
qf(1 - alpha, p, colm * n - colm - p + 1), ucl <- (p *
(colm + 1) * (n - 1))/(colm * n - colm - p + 1) *
qf(1 - alpha, p, colm * n - colm - p + 1)))
if (any(t2 > ucl)) {
cat("The following(s) point(s) fall outside of the control limits")
t3 <- which(t2 > ucl)
print(t3)
for (ii in 1:length(t3)) {
v = 1
k = 0
for (i in 1:p) {
k <- k + factorial(p)/(factorial(i) * factorial(p -
i))
}
q <- matrix(0, k, p + 3)
for (i in 1:p) {
a <- t(combn(p, i))
for (l in 1:nrow(a)) {
for (j in 1:ncol(a)) {
q[v, j + 3] <- a[l, j]
}
v = v + 1
}
}
for (i in 1:nrow(q)) {
b <- subset(q[i, 4:ncol(q)], q[i, 4:ncol(q)] >
0)
di <- length(b)
if (length(b) > 1) {
q[i, 1] <- n * t(Xmv[b] - x.jk[t3[ii], ][b]) %*%
solve(S[b, b]) %*% (Xmv[b] - x.jk[t3[ii],
][b])
}
else (q[i, 1] <- n * (x.jk[t3[ii], ][b] - Xmv[b])^2/S[b,
b])
ifelse(n == 1, ifelse(phase == 1, q[i, 2] <- ((colm -
1)^2)/colm * qbeta(1 - alpha, di/2, (((2 *
(colm - 1)^2)/(3 * colm - 4) - di - 1)/2)),
q[i, 2] <- ((di * (colm + 1) * (colm - 1))/((colm^2) -
colm * di)) * qf(1 - alpha, di, colm -
di)), ifelse(phase == 1, q[i, 2] <- (di *
(colm - 1) * (n - 1))/(colm * n - colm -
di + 1) * qf(1 - alpha, di, colm * n - colm -
di + 1), q[i, 2] <- (di * (colm + 1) * (n -
1))/(colm * n - colm - di + 1) * qf(1 - alpha,
di, colm * n - colm - di + 1)))
q[i, 3] <- 1 - pf(q[i, 1], di, colm - 1)
}
colnames(q) <- c("t2 decomp", "ucl",
"p-value", 1:p)
print(list(`Decomposition of` = t3[ii]))
print(round(q, 4))
}
}
}
if (type == "mewma") {
h4 <- matrix(c(8.6336, 9.6476, 10.083, 10.3114, 10.4405,
10.5152, 10.5581, 10.5816, 10.5932, 10.814, 11.8961,
12.3505, 12.5845, 12.7143, 12.788, 12.8297, 12.8524,
12.8635, 12.7231, 13.8641, 14.3359, 14.576, 14.7077,
14.7818, 14.8234, 14.846, 14.857, 14.5363, 15.7293,
16.217, 16.4629, 16.5965, 16.6711, 16.7127, 16.7352,
16.7463, 16.2634, 17.5038, 18.0063, 18.2578, 18.3935,
18.4687, 18.5105, 18.5331, 18.5442, 17.9269, 19.2113,
19.7276, 19.9845, 20.1223, 20.1982, 20.2403, 20.2631,
20.2743, 19.541, 20.8665, 21.396, 21.6581, 21.798,
21.8747, 21.9171, 21.9401, 21.9515, 21.1152, 22.4796,
23.0217, 23.2887, 23.4307, 23.5082, 23.551, 23.5742,
23.5858, 22.6565, 24.0579, 24.6119, 24.8838, 25.0278,
25.1062, 25.1493, 25.1728, 25.1846), nrow = 9)
rownames(h4) <- c(seq(0.1, 0.9, by = 0.1))
colnames(h4) <- c(1:9)
z <- matrix(0, m, p)
m1 <- rownames(h4)
m2 <- colnames(h4)
l <- lambda * 10
ucl <- h4[m1[l], m2[p - 1]]
name <- paste("MEWMA Control Chart")
for (i in 1:m) {
if (i == 1) {
z[i, ] <- lambda * (x.jk[i, ] - Xmv)
}
else {
z[i, ] <- lambda * (x.jk[i, ] - Xmv) + (1 - lambda) *
z[i - 1, ]
}
weig <- S * (lambda * (1 - ((1 - lambda)^(2 * i)))/(2 -
lambda))
t2[i, 1] <- t(z[i, ]) %*% solve(weig) %*% z[i, ]
}
}
if (type == "mcusum") {
name <- paste("MCUSUM Control Chart by Crosier (1988)")
ucl <- h
dif <- sweep(x.jk, 2, Xmv)
s <- matrix(0, m, p)
ci <- matrix(0, m, 1)
ci[1] <- sqrt(dif[1, ] %*% solve((S/n)) %*% dif[1, ])
if (ci[1] > k) {
s[1, ] <- (s[1, ] + dif[1, ]) * (1 - k/ci[1])
}
else (s[1, ] = matrix(0, ncol = p))
for (i in 2:m) {
ci[i, ] = sqrt((s[i - 1, ] + dif[i, ]) %*% solve(S/n) %*%
(s[i - 1, ] + dif[i, ]))
if (ci[i] > k) {
s[i, ] = (s[i - 1, ] + dif[i, ]) * (1 - k/ci[i])
}
else {
s[i, ] = matrix(0, ncol = p)
}
}
for (i in 1:m) {
t2[i] = sqrt(s[i, ] %*% solve((S/n)) %*% (s[i, ]))
}
}
if (type == "mcusum2") {
name <- paste("MCUSUM Control Chart by Pignatiello (1990)")
ucl <- h
dif <- sweep(x.jk, 2, Xmv)
s <- matrix(0, m, p)
l <- matrix(0, m, 1)
for (i in 1:m) {
if (i == 1) {
l[i, 1] <- 1
}
if (i > 1) {
if (t2[i - 1, 1] > 0) {
l[i, 1] <- l[i - 1, 1] + 1
}
else {
l[i, 1] <- 1
}
}
if (i == ((i - l[i, 1] + 1))) {
s[i, ] <- dif[i, ]
}
else {
s[i, ] <- colSums(dif[(i - l[i, 1] + 1):i, ])
}
t2[i, 1] <- max(0, (t(s[i, ]) %*% solve(S/n) %*%
s[i, ])^0.5 - k * l[i, 1])
}
}
t3 <- which(t2 > ucl)
# par(mar = c(4, 5, 3, 5))
# plot(t2, ylim = c(0, 1.1 * max(max(t2), ucl)), main = name,
# xlab = "Sample", ylab = expression(T^2), type = "o",
# las = 1)
# points(t3, t2[t3], col = 2)
# segments(0, ucl, m, ucl, col = 2)
# mtext(paste(" UCL=", round(ucl, 2)), side = 4, at = ucl,
# las = 2)
t2df <- data.frame(t2)
t2df$oob <- ifelse(t2df$t2 > ucl, "bad", "good")
t2df$sample <- seq(1:nrow(t2df))
p2 <- ggplot(data = t2df) +
geom_point(aes(x = sample, y = t2, color = oob)) +
scale_color_manual(values = c( "Red", "Grey20")) +
geom_path(aes(x = sample, y = t2))+
geom_hline(yintercept = ucl, color = "red") +
theme(legend.position = "none")
outList = list(
name,
ucl = round(ucl, 2),
t2 = round(t2, 2),
Xmv = round(Xmv, 2),
covariance = signif(S, 2),
plot2 = p2
)
return(outList)
}