The table you saw, is only printed when your type is t2 and is not stored in output. You can see in the code the table is variable q, and not returned.
I basically took the calculation part out (nasty bit of code) and wrote it into a function:
printOutlier = function(x,output,alpha){
p <- ncol(x)
m <- nrow(x)
x <- array(data.matrix(x), c(m, p, 1))
n <- dim(x)[3]
phase <- 2
x.jk <- matrix(0, m, p)
t2=output$t2
x.jk <- apply(x, 1:2, mean)
Xmv <- output$Xmv
S <- output$covariance
colm <- nrow(x)
ucl = output$ucl
t3 <- which(t2 > ucl)
res = vector("list",length(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)
names(res)[ii] <- paste(`Decomposition of` = t3[ii])
res[[ii]] <- round(q, 4)
}
return(res)
}
Now if you run your mult.chart again, and use this function, it should give you the table:
library(MSQC)
set.seed(111)
a<- runif(400,0,1)
b<- matrix(a, nrow=100,ncol=4)
output <- mult.chart(type="t2", alpha=0.07,b)
printOutlier(b,output,0.07)