Автоматическая подгонка кривой в R

Есть ли какой-нибудь пакет, который автоматически строит кривую, используя множество простых моделей?
Под простыми моделями я подразумеваю:

  • ax+b
  • ax^2+bx+c
  • а*лог(х) + б
  • a*x^n+b
  • ax/(1+bx)
  • ax^n/(1+bx^n)
  • ...

Лучше всего иметь функцию, которая принимает два векторных параметра X и Y и возвращает список подогнанных простых моделей с их SSE.


person Tomek Tarczynski    schedule 05.07.2012    source источник
comment
@Roland: найти лучшее нелинейное преобразование в линейной регрессии. Мне нравится создавать категориальную переменную (например, децили) из непрерывных переменных, а затем смотреть на график параметров для каждого дециля и среднего значения в каждом дециле. Это помогает найти нелинейное преобразование переменной. Я хотел бы немного ускорить этот процесс.   -  person Tomek Tarczynski    schedule 05.07.2012


Ответы (3)


Попробуй это. rhs — вектор символов правых сторон, а x и y — данные. Он строит формулу fo для каждого, а затем извлекает параметры и устанавливает для каждого значение 1 в качестве начального значения. Наконец, он запускает nls и возвращает SSE, отсортированные таким образом, что результатом является вектор SSE, названный по правой стороне. Если verbose=TRUE (что по умолчанию), то он также отображает выходные данные каждой подгонки.

sse <- function(rhs, x, y) sort(sapply(rhs, function(rhs, x, y, verbose = TRUE) {
    fo <- as.formula(paste("y", rhs, sep = "~"))
    nms <- setdiff(all.vars(fo), c("x", "y"))
    start <- as.list(setNames(rep(1, length(nms)), nms))
    fm <- nls(fo, data.frame(x, y), start = start)
    if (verbose) { print(fm); cat("---\n") }
    deviance(fm)
}, x = x, y = y))

## test

set.seed(123)
x <- 1:10
y <- rnorm(10, x)

# modify to suit
rhs <- c("a*x+b", "a*x*x+b*x+c")

sse(rhs, x, y)
person G. Grothendieck    schedule 05.07.2012
comment
Забавный результат: полином 2-го порядка всегда соответствует лучшему из первых четырех моделей, использующих несколько наборов данных rnorm. Я предполагаю, что это потому, что rnorm является только псевдослучайным. Но почему это лучше всего аппроксимировать полиномом 2-го порядка каждый раз, странно... - person DWAHL; 06.07.2012
comment
Три параметра подойдут лучше, чем два. Вы можете использовать AIC, если хотите оштрафовать большее количество параметров, или использовать SSE, но просто сравнивать модели с одинаковым количеством параметров. - person G. Grothendieck; 06.07.2012
comment
Спасибо, что показали, как создавать такие функции. Я надеялся, что он уже существует, потому что я могу придумать около 30 простых моделей, что означает, что таких моделей, вероятно, более 100 сотен. @ Г. Гротендик: Как вы думаете, nls даст лучшие результаты, чем nlminb? - person Tomek Tarczynski; 06.07.2012
comment
nls предназначен для нелинейного метода наименьших квадратов, тогда как nlminb предназначен для общей цели, поэтому обычно сначала следует попробовать nls для этого типа задач. - person G. Grothendieck; 06.07.2012

Вы также можете взглянуть на пакеты, предоставляющие функции для вычисления дробных многочленов. На данный момент это mboost (с функцией FP) и mfp (с функцией mfp). Хотя я не пробовал пакеты, теория, лежащая в их основе, соответствует тому, что вам нужно.

Пакет mfp был описан в R-News в 2005 году.

Две ссылки, которые могут представлять интерес:

Ройстон П., Альтман Д. (1994) Регрессия с использованием дробных полиномов непрерывных ковариат. Заявл. стат. 3: 429–467.

Sauerbrei W, Royston P (1999)Построение многомерных прогностических и диагностических моделей: преобразование предикторов с использованием дробных полиномов. Журнал Королевского статистического общества (Серия А) 162: 71–94.

person BenBarnes    schedule 05.07.2012
comment
Спасибо, я никогда не слышал о дробных полиномах. Обязательно прочту об этом! - person Tomek Tarczynski; 06.07.2012

Вы можете подогнать регрессионные сплайны и найти хорошее соответствие, вручную отрегулировав степени свободы несколько раз. Попробуйте следующую функцию:

spline.fit <- function(x, y, df=5) {
  ## INPUT: x, y are two vectors (predictor and response);
  ##        df is the number of spline basis.  Increase "df" to fit more adaptively to the data.
  require(splines) # available as default R Package.
  bx <- bs(x, df)  # B-spline basis matrix as New Predictors (dimension is "length(x)" by "df")
  f <- lm(y ~ bx)  # Linear Regression on Spline Basis (that is, "df" number of new predictors)
  fy <- fitted(f)  # Fitted Response
  plot(x, y); lines(x, fy, col="blue", lwd=2) # Make a plot to show the fit.
  invisible(list(x=bx, y=fy, f=f))    # Return the Basis (new predictors), Fitted Y, Regression
}

if (F) {                                # Unit Test
  spline.fit(1:100, rnorm(100))
  spline.fit(1:100, rnorm(100), df=20)
}
person Feiming Chen    schedule 18.09.2013