Быстрые пост-специальные вычисления с использованием R

Question

У меня есть большой набор данных, который я хотел бы выполнить после специального расчета:

dat = as.data.frame(matrix(runif(10000*300), ncol = 10000, nrow = 300))

dat$group = rep(letters[1:3], 100)

Вот мой код:

start <- Sys.time()

vars <- names(dat)[-ncol(dat)] 

aov.out <- lapply(vars, function(x) {
        lm(substitute(i ~ group, list(i = as.name(x))), data = dat)})

TukeyHSD.out <- lapply(aov.out, function(x) TukeyHSD(aov(x)))

Sys.time() - start

Разница во времени 4.033335 минут

Это займет около 4 минут, есть ли более эффективные и элегантные способы выполнения post hoc с использованием R?

большое спасибо

Answer 1

Ваш пример слишком велик. Для иллюстрации идеи я использую маленький.

set.seed(0)
dat = as.data.frame(matrix(runif(2*300), ncol = 2, nrow = 300))
dat$group = rep(letters[1:3], 100)

Зачем звонишь aov на встроенной модели "лм"? Это в основном соответствует той же модели.

Сначала прочитайте " Подгонка линейной модели с несколькими LHS". lm рабочая лошадка aovТаким образом, вы можете передать несколько формул LHS aov, Модель имеет класс c("maov", "aov", "mlm", "lm"),

response_names <- names(dat)[-ncol(dat)]
form <- as.formula(sprintf("cbind(%s) ~ group", toString(response_names)))
fit <- do.call("aov", list(formula = form, data = quote(dat)))

Теперь проблема в том, что нет метода "Маов" для TuckyHSD, Так что нам нужен взлом.

TuckyHSD опирается на остатки подогнанной модели. В c("aov", "lm") В случае, если остатки - это вектор, но в c("maov", "aov", "mlm", "lm") случай это матрица. Следующее демонстрирует взлом.

aov_hack <- fit
aov_hack[c("coefficients", "fitted.values")] <- NULL  ## don't need them
aov_hack[c("contrasts", "xlevels")] <- NULL  ## don't need them either
attr(aov_hack$model, "terms") <- NULL  ## don't need it
class(aov_hack) <- c("aov", "lm")  ## drop "maov" and "mlm"
## the following elements are mandatory for `TukeyHSD`
## names(aov_hack)
#[1] "residuals"   "effects"     "rank"        "assign"      "qr"         
#[6] "df.residual" "call"        "terms"       "model" 

N <- length(response_names)  ## number of response variables
result <- vector("list", N)
for (i in 1:N) {
  ## change response variable in the formula
  aov_hack$call[[2]][[2]] <- as.name(response_names[i])
  ## change residuals
  aov_hack$residuals <- fit$residuals[, i]
  ## change effects
  aov_hack$effects <- fit$effects[, i]
  ## change "terms" object and attribute
  old_tm <- terms(fit)  ## old "terms" object in the model
  old_tm[[2]] <- as.name(response_names[i])  ## change response name in terms
  new_tm <- terms.formula(formula(old_tm))  ## new "terms" object
  aov_hack$terms <- new_tm  ## replace `aov_hack$terms`
  ## replace data in the model frame
  aov_hack$model[1] <- data.frame(fit$model[[1]][, i])
  names(aov_hack$model)[1] <- response_names[i]
  ## run `TukeyHSD` on `aov_hack`
  result[[i]] <- TukeyHSD(aov_hack)
  }

result[[1]]  ## for example
#  Tukey multiple comparisons of means
#    95% family-wise confidence level
#
#Fit: aov(formula = V1 ~ group, data = dat)
#
#$group
#            diff        lwr        upr     p adj
#b-a -0.012743870 -0.1043869 0.07889915 0.9425847
#c-a -0.022470482 -0.1141135 0.06917254 0.8322109
#c-b -0.009726611 -0.1013696 0.08191641 0.9661356

Я использовал цикл "для". Замените его lapply если ты хочешь.