如果你查看WRS2::rmmcp的源代码，前几行显示了错误的原因。
WRS2::rmmcp
function (y, groups, blocks, tr = 0.2, alpha = 0.05)
{
cols1 <- deparse(substitute(y))
cols2 <- deparse(substitute(groups))
cols3 <- deparse(substitute(blocks))
dat <- data.frame(y, groups, blocks)
colnames(dat) <- c(cols1, cols2, cols3)
#...
#...
}
deparse(substitute()) 代码在与rlang::exec 一起使用时无法按预期工作。我们可以对列名进行硬编码以使其与rlang::exec 一起工作。
tmp <- function (y, groups, blocks, tr = 0.2, alpha = 0.05)
{
cols1 <- 'col1' #Change
cols2 <- 'col2' #Change
cols3 <- 'col3' #Change
dat <- data.frame(y, groups, blocks)
colnames(dat) <- c(cols1, cols2, cols3)
cl <- match.call()
x <- reshape(dat, idvar = cols3, timevar = cols2, direction = "wide")[-1]
grp <- c(1:length(x))
con = 0
dif = TRUE
flagcon = F
if (!is.matrix(x))
x <- matl(x)
if (!is.matrix(x))
stop("Data must be stored in a matrix or in list mode.")
con <- as.matrix(con)
J <- ncol(x)
xbar <- vector("numeric", J)
x <- elimna(x)
nval <- nrow(x)
h1 <- nrow(x) - 2 * floor(tr * nrow(x))
df <- h1 - 1
for (j in 1:J) xbar[j] <- mean(x[, j], tr)
if (sum(con^2 != 0))
CC <- ncol(con)
if (sum(con^2) == 0)
CC <- (J^2 - J)/2
ncon <- CC
if (alpha == 0.05) {
dvec <- c(0.05, 0.025, 0.0169, 0.0127, 0.0102, 0.00851,
0.0073, 0.00639, 0.00568, 0.00511)
if (ncon > 10) {
avec <- 0.05/c(11:ncon)
dvec <- c(dvec, avec)
}
}
if (alpha == 0.01) {
dvec <- c(0.01, 0.005, 0.00334, 0.00251, 0.00201, 0.00167,
0.00143, 0.00126, 0.00112, 0.00101)
if (ncon > 10) {
avec <- 0.01/c(11:ncon)
dvec <- c(dvec, avec)
}
}
if (alpha != 0.05 && alpha != 0.01)
dvec <- alpha/c(1:ncon)
if (sum(con^2) == 0) {
flagcon <- T
psihat <- matrix(0, CC, 5)
dimnames(psihat) <- list(NULL, c("Group", "Group", "psihat",
"ci.lower", "ci.upper"))
test <- matrix(NA, CC, 6)
dimnames(test) <- list(NULL, c("Group", "Group", "test",
"p.value", "p.crit", "se"))
temp1 <- 0
jcom <- 0
for (j in 1:J) {
for (k in 1:J) {
if (j < k) {
jcom <- jcom + 1
q1 <- (nrow(x) - 1) * winvar(x[, j], tr)
q2 <- (nrow(x) - 1) * winvar(x[, k], tr)
q3 <- (nrow(x) - 1) * wincor(x[, j], x[, k],
tr)$cov
sejk <- sqrt((q1 + q2 - 2 * q3)/(h1 * (h1 -
1)))
if (!dif) {
test[jcom, 6] <- sejk
test[jcom, 3] <- (xbar[j] - xbar[k])/sejk
temp1[jcom] <- 2 * (1 - pt(abs(test[jcom,
3]), df))
test[jcom, 4] <- temp1[jcom]
psihat[jcom, 1] <- j
psihat[jcom, 2] <- k
test[jcom, 1] <- j
test[jcom, 2] <- k
psihat[jcom, 3] <- (xbar[j] - xbar[k])
}
if (dif) {
dv <- x[, j] - x[, k]
test[jcom, 6] <- trimse(dv, tr)
temp <- trimci(dv, alpha = alpha/CC, pr = FALSE,
tr = tr)
test[jcom, 3] <- temp$test.stat
temp1[jcom] <- temp$p.value
test[jcom, 4] <- temp1[jcom]
psihat[jcom, 1] <- j
psihat[jcom, 2] <- k
test[jcom, 1] <- j
test[jcom, 2] <- k
psihat[jcom, 3] <- mean(dv, tr = tr)
psihat[jcom, 4] <- temp$ci[1]
psihat[jcom, 5] <- temp$ci[2]
}
}
}
}
temp2 <- order(0 - temp1)
zvec <- dvec[1:ncon]
sigvec <- (test[temp2] >= zvec)
if (sum(sigvec) < ncon) {
dd <- ncon - sum(sigvec)
ddd <- sum(sigvec) + 1
zvec[ddd:ncon] <- dvec[ddd]
}
test[temp2, 5] <- zvec
if (!dif) {
psihat[, 4] <- psihat[, 3] - qt(1 - alpha/(2 * CC),
df) * test[, 6]
psihat[, 5] <- psihat[, 3] + qt(1 - alpha/(2 * CC),
df) * test[, 6]
}
}
if (sum(con^2) > 0) {
if (nrow(con) != ncol(x))
warning("The number of groups does not match the number of contrast coefficients.")
ncon <- ncol(con)
psihat <- matrix(0, ncol(con), 4)
dimnames(psihat) <- list(NULL, c("con.num", "psihat",
"ci.lower", "ci.upper"))
test <- matrix(0, ncol(con), 5)
dimnames(test) <- list(NULL, c("con.num", "test", "p.value",
"p.crit", "se"))
temp1 <- NA
for (d in 1:ncol(con)) {
psihat[d, 1] <- d
if (!dif) {
psihat[d, 2] <- sum(con[, d] * xbar)
sejk <- 0
for (j in 1:J) {
for (k in 1:J) {
djk <- (nval - 1) * wincor(x[, j], x[, k],
tr)$cov/(h1 * (h1 - 1))
sejk <- sejk + con[j, d] * con[k, d] * djk
}
}
sejk <- sqrt(sejk)
test[d, 1] <- d
test[d, 2] <- sum(con[, d] * xbar)/sejk
test[d, 5] <- sejk
temp1[d] <- 2 * (1 - pt(abs(test[d, 2]), df))
}
if (dif) {
for (j in 1:J) {
if (j == 1)
dval <- con[j, d] * x[, j]
if (j > 1)
dval <- dval + con[j, d] * x[, j]
}
temp1[d] <- trimci(dval, tr = tr, pr = FALSE)$p.value
test[d, 1] <- d
test[d, 2] <- trimci(dval, tr = tr, pr = FALSE)$test.stat
test[d, 5] <- trimse(dval, tr = tr)
psihat[d, 2] <- mean(dval, tr = tr)
}
}
test[, 3] <- temp1
temp2 <- order(0 - temp1)
zvec <- dvec[1:ncon]
sigvec <- (test[temp2, 3] >= zvec)
if (sum(sigvec) < ncon) {
dd <- ncon - sum(sigvec)
ddd <- sum(sigvec) + 1
}
test[temp2, 4] <- zvec
psihat[, 3] <- psihat[, 2] - qt(1 - test[, 4]/2, df) *
test[, 5]
psihat[, 4] <- psihat[, 2] + qt(1 - test[, 4]/2, df) *
test[, 5]
}
if (flagcon)
num.sig <- sum(test[, 4] <= test[, 5])
if (!flagcon)
num.sig <- sum(test[, 3] <= test[, 4])
fnames <- as.character(unique(groups))
psihat1 <- cbind(psihat, test[, 4:5])
result <- list(comp = psihat1, fnames = fnames, call = cl)
class(result) <- "mcp2"
result
}
请注意，我复制整个代码只是为了重现性，此函数的更改仅是第 1 3 行。
运行该函数后，您可以将其用作：
tmp(
y = WineTasting$Taste,
groups = WineTasting$Wine,
blocks = WineTasting$Taster
)
#Call:
#tmp(y = WineTasting$Taste, groups = WineTasting$Wine, blocks = WineTasting$Taster)
# psihat ci.lower ci.upper p.value p.crit sig
#Wine A vs. Wine B 0.02143 -0.02164 0.06449 0.19500 0.0500 FALSE
#Wine A vs. Wine C 0.11429 0.02148 0.20710 0.00492 0.0169 TRUE
#Wine B vs. Wine C 0.08214 0.00891 0.15538 0.00878 0.0250 TRUE
还有rlang::exec
res <- rlang::exec(
.fn = tmp,
y = WineTasting$Taste,
groups = WineTasting$Wine,
blocks = WineTasting$Taster
)
res$comp
# Group Group psihat ci.lower ci.upper p.value p.crit
#[1,] 1 2 0.02142857 -0.021636832 0.06449397 0.195004531 0.0500
#[2,] 1 3 0.11428571 0.021475579 0.20709585 0.004915566 0.0169
#[3,] 2 3 0.08214286 0.008910564 0.15537515 0.008777396 0.0250
res$fnames
#[1] "Wine A" "Wine B" "Wine C"
（虽然与rlang::exec 一起使用会破坏res$call。不知道为什么！）
运行此程序时，我遇到了could not find function elimna 或could not find function matl 之类的错误，这很奇怪，因为这些函数来自我已加载的包WRS2，但它仍然给出了错误。我必须将 https://github.com/cran/WRS2/tree/master/R 中的函数复制到我的会话中，然后它才能如上所示工作。