Запустите lm с несколькими ответами и весами

Я должен приспособить линейную модель с одной и той же модельной матрицей к нескольким ответам. Это можно легко сделать в R, указав ответ в виде матрицы вместо вектора. Таким образом, вычисления очень быстрые.

Теперь я также хотел бы добавить весовые коэффициенты к модели, которые соответствуют точности ответов. Следовательно, для каждого вектора ответа мне понадобится также различный весовой вектор. Тем не мение, lm позволяет вводить веса только как вектор, а не как матрицу. Так что я не мог ввести веса в партии и должен был бы бежать lm за каждый ответ отдельно. Таким образом, вычисления станут намного медленнее.

Есть ли способ запустить модели такого типа в пакетном режиме, не вызывая lm несколько раз?

Источник

1 ответ

Теперь я также хотел бы добавить весовые коэффициенты к модели, которые соответствуют точности ответов. Следовательно, для каждого вектора ответа мне понадобится также различный весовой вектор. Тем не менее lm позволяет вводить веса только как вектор, а не как матрицу. Так что я не мог ввести веса в партии и должен был бы бежать lm за каждый ответ отдельно. Таким образом, вычисления станут намного медленнее.

Как объясняется в разделе "Подгонка линейной модели с несколькими LHS", для эффективности "mlm" требуется общая матрица моделей для всех ответов LHS. Однако взвешенная регрессия не дает повторного использования матрицы модели, так как для другого набора весов оба ответа y и модель матрицы X нужно изменить Прочитайте R: lm() результат отличается при использовании weights аргумент и при использовании вручную взвешенных данных, чтобы увидеть, как работает взвешенная регрессия.

Есть ли способ запустить модели такого типа в пакетном режиме, не вызывая lm несколько раз?

Это зависит от того, что вы хотите. Если вам нужен полный lmObject тогда вам нужно позвонить lm каждый раз. Если вам нужны только коэффициенты, вы можете использовать .lm.fit, 2-я ссылка выше продемонстрировала использование lm.fit, в то время как использование .lm.fit почти идентичен. Простой шаблон может быть следующим:

## weighted mlm, by specifying matrix directly
## `xmat`: non-weighted model matrix, manually created from `model.matrix`
## `ymat`: non-weighted response matrix
## `wmat`: matrix of weights

## all matrices must have the same number of rows (not checked)
## `ymat` and `wmat` must have the same number of columns (not checked)
## no `NA` values in any where is allowed (not checked)
## all elements of `wmat` must be strictly positive (not checked)

wmlm <- function (xmat, ymat, wmat) {
  N <- ncol(ymat)
  wmlmList <- vector("list", length = N)
  for (j in 1:N) {
    rw <- sqrt(wmat[, j])
    wmlmList[[j]] <- .lm.fit(rw * xmat, rw * ymat[, j])
    }
  return(wmlmList)
  }

Рассмотрим простой пример его использования:

## a toy dataset of 200 data with 3 numerical variables and 1 factor variable
dat <- data.frame(x1 = rnorm(200), x2 = rnorm(200), x3 = rnorm(200), f = gl(5, 40, labels = letters[1:5]))

## consider a model `~ x1 + poly(x3, 3) + x2 * f`
## we construct the non-weighted model matrix
xmat <- model.matrix (~ x1 + poly(x3, 3) + x2 * f, dat)

## now let's assume we have 100 model responses as well as 100 sets of weights
ymat <- matrix(rnorm(200 * 100), 200)
wmat <- matrix(runif(200 * 100), 200)

## Let's call `wmlm`:
fit <- wmlm (xmat, ymat, wmat)

.lm.fit возвращает критическую информацию для дальнейшего вывода модели и полного lmObject унаследует большинство этих записей.

## take the first fitted model as an example
str(fit[[1]])
 #$ qr          : num [1:200, 1:14] -10.4116 0.061 0.0828 0.0757 0.0698 ...
 # ..- attr(*, "assign")= int [1:14] 0 1 2 2 2 3 4 4 4 4 ...
 # ..- attr(*, "contrasts")=List of 1
 # .. ..$ f: chr "contr.treatment"
 # ..- attr(*, "dimnames")=List of 2
 # .. ..$ : chr [1:200] "1" "2" "3" "4" ...
 # .. ..$ : chr [1:14] "(Intercept)" "x1" "poly(x3, 3)1" "poly(x3, 3)2" ...
 #$ coefficients: num [1:14] 0.1184 -0.0506 0.3032 0.1643 0.4269 ...
 #$ residuals   : num [1:200] -0.7311 -0.0795 -0.2495 0.4097 0.0495 ...
 #$ effects     : num [1:200] -0.351 -0.36 0.145 0.182 0.291 ...
 #$ rank        : int 14
 #$ pivot       : int [1:14] 1 2 3 4 5 6 7 8 9 10 ...
 #$ qraux       : num [1:14] 1.06 1.13 1.07 1.05 1.01 ...
 #$ tol         : num 1e-07
 #$ pivoted     : logi FALSE

Результат .lm.fit не поддерживает родовые функции, такие как summary, anova, predict, plot и т. д. Но сделать вывод о линейной модели несложно, так что вычислить себя легко (при условии, что вы знаете теорию, лежащую в основе):

  1. Таблица t-статистики для коэффициентов регрессии (от $qr);
  2. F-статистика и таблица ANOVA (от $effects);
  3. Остаточная стандартная ошибка, R-квадрат и скорректированный R-квадрат (от $residulas а также $rank).

Наконец, я предлагаю тест:

library(microbenchmark)
microbenchmark(wmlm = {wmlm (xmat, ymat, wmat);},
               lm = {for (i in 1:ncol(ymat))
                       lm(ymat[, i] ~ x1 + poly(x3, 3) + x2 * f, dat, weights = wmat[, i]);} )

#Unit: milliseconds
# expr       min        lq      mean    median       uq      max neval cld
# wmlm  20.84512  23.02756  27.29539  24.49314  25.9027  79.8894   100  a 
#   lm 400.53000 405.10622 430.09787 414.42152 442.2640 535.9144   100   b

Таким образом, увеличение в 17,25 раза (по среднему времени).

Источник
Другие вопросы по тегам