R: как выполнить более сложные вычисления из гребенки набора данных?

Question

Прямо сейчас, у меня есть гребенка из встроенного набора диафрагмы. До сих пор я руководствовался тем, что смог найти коэффициент lm() пары значений.

myPairs <- combn(names(iris[1:4]), 2) 

formula <- apply(myPairs, MARGIN=2, FUN=paste, collapse="~") 

model <- lapply(formula, function(x) lm(formula=x, data=iris)$coefficients[2]) 

model 

Однако, я хотел бы пойти на несколько шагов вперед и использовать коэффициент от ого(), которые будут использоваться в дальнейших расчетах. Я хотел бы сделать что-то вроде этого:

Coefficient <- lm(formula=x, data=iris)$coefficients[2] 
Spread <- myPairs[1] - coefficient*myPairs[2] 
library(tseries) 
adf.test(Spread) 

Сама процедура достаточно проста, но я не смог найти способ сделать это для каждого combn в наборе данных. (В качестве оповещения adf.test не будет применяться к таким данным, но я просто использую набор данных диафрагмы для демонстрации). Мне интересно, было бы лучше написать цикл для такой процедуры?

Answer 1

Вы можете сделать все это в пределах combn.

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

fun <- function(x) coef(lm(paste(x, collapse="~"), data=iris))[2] 
combn(names(iris[1:4]), 2, fun) 

Вы можете расширить функцию, чтобы рассчитать спрэд

fun <- function(x) { 
     est <- coef(lm(paste(x, collapse="~"), data=iris))[2] 
     spread <- iris[,x[1]] - est*iris[,x[2]] 
     adf.test(spread) 
     } 

out <- combn(names(iris[1:4]), 2, fun, simplify=FALSE) 
out[[1]] 

# Augmented Dickey-Fuller Test 

#data: spread 
#Dickey-Fuller = -3.879, Lag order = 5, p-value = 0.01707 
#alternative hypothesis: stationary 

Сравнить результаты для запуска первого вручную

est <- coef(lm(Sepal.Length ~ Sepal.Width, data=iris))[2] 
spread <- iris[,"Sepal.Length"] - est*iris[,"Sepal.Width"] 
adf.test(spread) 

# Augmented Dickey-Fuller Test 

# data: spread 
# Dickey-Fuller = -3.879, Lag order = 5, p-value = 0.01707 
# alternative hypothesis: stationary 

Answer 2

Похоже, вы хотели бы написать свою собственную функцию и вызвать его в цикле myPairs (применить):

yourfun <- function(pair){ 
    fm <- paste(pair, collapse='~') 
    coef <- lm(formula=fm, data=iris)$coefficients[2] 
    Spread <- iris[,pair[1]] - coef*iris[,pair[2]] 
    return(Spread) 
} 

Тогда вы можете вызвать эту функцию:

model <- apply(myPairs, 2, yourfun) 

Я думаю, что это самый чистый путь. Но я не знаю, что именно вы хотите сделать, поэтому я составил пример для Spread. Обратите внимание, что в моем примере вы получаете предупреждающие сообщения, поскольку столбец Species является фактором.

Answer 3

Несколько советов: я бы не назвал вещи, которые вы с тем же именем, что и встроенные функции (model, formula приходят на ум в вашей первоначальной версии).

Кроме того, вы можете упростить paste, что вы делаете - см. Ниже.

И, наконец, более общее утверждение: не чувствуйте, что все должно быть сделано в *apply. Иногда краткость и короткий код на самом деле сложнее понять, и помните, что функции *apply предлагают в лучшем случае предельную скорость, получающуюся за простой цикл for. (Это не всегда было с R, но именно в этот момент).

# Get pairs 
myPairs <- combn(x = names(x = iris[1:4]),m = 2) 

# Just directly use paste() here 
myFormulas <- paste(myPairs[1,],myPairs[2,],sep = "~") 

# Store the models themselves into a list 
# This lets you go back to the models later if you need something else 
myModels <- lapply(X = myFormulas,FUN = lm,data = iris) 

# If you use sapply() and this simple function, you get back a named vector 
# This seems like it could be useful to what you want to do 
myCoeffs <- sapply(X = myModels,FUN = function (x) {return(x$coefficients[2])}) 

# Now, you can do this using vectorized operations 
iris[myPairs[1,]] - iris[myPairs[2,]] * myCoeffs[myPairs[2,]] 

Если я правильно понимаю, я считаю, что вышеизложенное будет работать. Обратите внимание, что имена на выходе в настоящее время будут бессмысленными, вам нужно будет заменить их чем-то своим собственным дизайном (возможно, значениями myFormulas).