我会非常感谢一些帮助。我没有编码背景，我被运行WRS 2包的sppb函数时得到的错误消息所迷惑。这些函数使用 Bootstrap 执行强大的混合ANOVA。
sppba（公式= score ~ my_between_variable * my_within_variable，id = participant_code，data = df_long_T2）

Error in xmat[, k] <- x[[kv]] : 
  number of items to replace is not a multiple of replacement length

三个sppb函数都出现了同样的错误。除了sppba被替换成了sppbb和sppbi之外，这些函数看起来都一样。我甚至不知道这些函数要替换的是什么。这些函数对我来说可以处理其他数据。
所有涉及的类似乎都很好：score是数字，order_supplement和time是因子，participant_code是字符，df_long_T2是数据框。我有120个参与者，一组61个，另一组59个，每个参与者有两个观察结果。在涉及的列中没有NA。
Traceback（）只给我上面的一行代码和错误消息。
Debug（）给了我这个，我不知道该怎么做：“调试位置是近似的，因为位置不可用”

function (formula, id, data, est = "mom", avg = TRUE, nboot = 500, 
  MDIS = FALSE, ...) 
{
  if (missing(data)) {
    mf <- model.frame(formula)
  }
  else {
    mf <- model.frame(formula, data)
  }
  cl <- match.call()
  est <- match.arg(est, c("mom", "onestep", "median"), several.ok = FALSE)
  mf1 <- match.call()
  m <- match(c("formula", "data", "id"), names(mf1), 0L)
  mf1 <- mf1[c(1L, m)]
  mf1$drop.unused.levels <- TRUE
  mf1[[1L]] <- quote(stats::model.frame)
  mf1 <- eval(mf1, parent.frame())
  random1 <- mf1[, "(id)"]
  depvar <- colnames(mf)[1]
  if (all(length(table(random1)) == table(mf[, 3]))) {
    ranvar <- colnames(mf)[3]
    fixvar <- colnames(mf)[2]
  }
  else {
    ranvar <- colnames(mf)[2]
    fixvar <- colnames(mf)[3]
  }
  MC <- FALSE
  K <- length(table(mf[, ranvar]))
  J <- length(table(mf[, fixvar]))
  p <- J * K
  grp <- 1:p
  est <- get(est)
  fixsplit <- split(mf[, depvar], mf[, fixvar])
  indsplit <- split(mf[, ranvar], mf[, fixvar])
  dattemp <- mapply(split, fixsplit, indsplit, SIMPLIFY = FALSE)
  data <- do.call(c, dattemp)
  x <- data
  jp <- 1 - K
  kv <- 0
  kv2 <- 0
  for (j in 1:J) {
    jp <- jp + K
    xmat <- matrix(NA, ncol = K, nrow = length(x[[jp]]))
    for (k in 1:K) {
      kv <- kv + 1
      xmat[, k] <- x[[kv]]
    }
    xmat <- elimna(xmat)
    for (k in 1:K) {
      kv2 <- kv2 + 1
      x[[kv2]] <- xmat[, k]
    }
  }
  xx <- x
  nvec <- NA
  jp <- 1 - K
  for (j in 1:J) {
    jp <- jp + K
    nvec[j] <- length(x[[jp]])
  }
  bloc <- matrix(NA, nrow = J, ncol = nboot)
  mvec <- NA
  ik <- 0
  for (j in 1:J) {
    x <- matrix(NA, nrow = nvec[j], ncol = K)
    for (k in 1:K) {
      ik <- ik + 1
      x[, k] <- xx[[ik]]
      if (!avg) 
        mvec[ik] <- est(xx[[ik]])
    }
    tempv <- apply(x, 2, est)
    data <- matrix(sample(nvec[j], size = nvec[j] * nboot, 
      replace = TRUE), nrow = nboot)
    bvec <- matrix(NA, ncol = K, nrow = nboot)
    for (k in 1:K) {
      temp <- x[, k]
      bvec[, k] <- apply(data, 1, rmanogsub, temp, est)
    }
    if (avg) {
      mvec[j] <- mean(tempv)
      bloc[j, ] <- apply(bvec, 1, mean)
    }
    if (!avg) {
      if (j == 1) 
        bloc <- bvec
      if (j > 1) 
        bloc <- cbind(bloc, bvec)
    }
  }
  if (avg) {
    d <- (J^2 - J)/2
    con <- matrix(0, J, d)
    id <- 0
    Jm <- J - 1
    for (j in 1:Jm) {
      jp <- j + 1
      for (k in jp:J) {
        id <- id + 1
        con[j, id] <- 1
        con[k, id] <- 0 - 1
      }
    }
  }
  if (!avg) {
    MJK <- K * (J^2 - J)/2
    JK <- J * K
    MJ <- (J^2 - J)/2
    cont <- matrix(0, nrow = J, ncol = MJ)
    ic <- 0
    for (j in 1:J) {
      for (jj in 1:J) {
        if (j < jj) {
          ic <- ic + 1
          cont[j, ic] <- 1
          cont[jj, ic] <- 0 - 1
        }
      }
    }
    tempv <- matrix(0, nrow = K - 1, ncol = MJ)
    con1 <- rbind(cont[1, ], tempv)
    for (j in 2:J) {
      con2 <- rbind(cont[j, ], tempv)
      con1 <- rbind(con1, con2)
    }
    con <- con1
    if (K > 1) {
      for (k in 2:K) {
        con1 <- push(con1)
        con <- cbind(con, con1)
      }
    }
  }
  if (!avg) 
    bcon <- t(con) %*% t(bloc)
  if (avg) 
    bcon <- t(con) %*% (bloc)
  tvec <- t(con) %*% mvec
  tvec <- tvec[, 1]
  tempcen <- apply(bcon, 1, mean)
  vecz <- rep(0, ncol(con))
  bcon <- t(bcon)
  temp = bcon
  for (ib in 1:nrow(temp)) temp[ib, ] = temp[ib, ] - tempcen + 
    tvec
  bcon <- rbind(bcon, vecz)
  if (!MDIS) {
    if (!MC) 
      dv = pdis(bcon, center = tvec)
  }
  if (MDIS) {
    smat <- var(temp)
    bcon <- rbind(bcon, vecz)
    chkrank <- qr(smat)$rank
    if (chkrank == ncol(smat)) 
      dv <- mahalanobis(bcon, tvec, smat)
    if (chkrank < ncol(smat)) {
      smat <- ginv(smat)
      dv <- mahalanobis(bcon, tvec, smat, inverted = T)
    }
  }
  bplus <- nboot + 1
  sig.level <- 1 - sum(dv[bplus] >= dv[1:nboot])/nboot
  tvec1 <- data.frame(Estimate = tvec)
  if (avg) {
    tnames <- apply(combn(levels(mf[, fixvar]), 2), 2, paste0, 
      collapse = "-")
    rownames(tvec1) <- tnames
  }
  else {
    fixcomb <- apply(combn(levels(mf[, fixvar]), 2), 2, 
      paste0, collapse = "-")
    rnames <- levels(mf[, ranvar])
    tnames <- as.vector(t(outer(rnames, fixcomb, paste)))
    rownames(tvec1) <- tnames
  }
  result <- list(test = tvec1, p.value = sig.level, contrasts = con, 
    call = cl)
  class(result) <- c("spp")
  result
}

我希望得到这样的输出：

## Test statistics:
##                              Estimate
## time1-time2   0.3000
## 
## Test whether the corrresponding population parameters are the same:
## p-value: 0.37