Сплайны внутри нелинейных наименьших квадратов в R - PullRequest
5 голосов
/ 04 февраля 2012

Рассмотрим нелинейную модель наименьших квадратов в R, например, следующего вида:

 y ~ theta / ( 1 + exp( -( alpha + beta * x) ) )

(моя настоящая проблема состоит из нескольких переменных, и внешняя функция не логистическая, а немного более сложная; эта задача проще, но я думаю, что если я смогу это сделать, мой случай должен последовать почти сразу)

Я бы хотел заменить термин «альфа + бета * х» на (скажем) естественный кубический сплайн.

вот код для создания примеров данных с нелинейной функцией внутри логистики:

set.seed(438572L)
x <- seq(1,10,by=.25)
y <- 8.6/(1+exp( -(-3+x/4.4+sqrt(x*1.1)*(1.-sin(1.+x/2.9))) )) + rnorm(x, s=0.2 )

Без необходимости в логистике вокруг нее, если бы я был в lm, я мог бы легко заменить линейный член на сплайн-член; поэтому линейная модель примерно такая:

 lm( y ~ x ) 

затем становится

 library("splines")
 lm( y ~ ns( x, df = 5 ) )

Генерация подобранных значений проста и получение прогнозируемых значений с помощью (для пример) пакет rms кажется достаточно простым.

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

Проблема с nls в том, что мне нужно предоставить имена для всех параметров (я очень рад, что назвал их (b1, ..., b5) для одной посадки сплайна (и скажем, c1, ..., c6) для другой переменной - мне нужно будет сделать несколько из них).

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

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

( редактировать для уточнения ) Для этой небольшой проблемы я, конечно, могу сделать это вручную - выписать выражение для внутреннего произведения каждой переменной в матрице, сгенерированной ns, раз вектор параметров. Но затем я должен выписать все это поэлементно снова для каждого сплайна в каждой другой переменной, и снова каждый раз, когда я изменяю df в любом из сплайнов, и снова, если я хочу использовать cs вместо ns. И затем, когда я хочу попытаться сделать некоторый прогноз (/ интерполяцию), мы получаем целый ряд новых проблем, которые необходимо решить. Мне нужно продолжать делать это снова и снова, и, возможно, для существенно большего числа узлов и нескольких переменных, для анализа после анализа - и я подумал, есть ли более аккуратный, простой способ, чем выписывать каждый отдельный термин, без необходимости писать много кода. Я могу видеть довольно простой способ сделать это, который потребует изрядного количества кода для правильной реализации, но, будучи R, я подозреваю, что есть намного более аккуратный способ (или, скорее, 3 или 4 более аккуратных способа), это просто ускользает от меня. Отсюда вопрос.

Я думал, что видел, как кто-то делал что-то подобное в прошлом довольно хорошим способом, но для жизни я не могу найти это сейчас; Я пытался найти его несколько раз.

[В частности, я обычно хотел бы иметь возможность попробовать подгонку любого из нескольких разных сплайнов в каждой переменной - попробовать пару возможностей - чтобы посмотреть, смогу ли я найти простую модель, но все же одну где подгонка подходит для этой цели (шум действительно довольно низкий; некоторый уклон в подгонке вполне приемлем для достижения хорошего плавного результата, но только до определенной точки). Это скорее «найти хорошую, интерпретируемую, но адекватную функцию подбора», чем что-либо, приближающееся к выводу, и анализ данных на самом деле не является проблемой для этой проблемы.]

В качестве альтернативы, если бы это было намного проще, скажем, в gnm или ASSIST или в одном из других пакетов, это было бы полезным знанием, но тогда некоторые указатели о том, как действовать с игрушкой, описанной выше, могли бы помочь.

Ответы [ 2 ]

9 голосов
/ 04 февраля 2012

ns фактически генерирует матрицу предикторов. Что вы можете сделать, так это разбить эту матрицу на отдельные переменные и передать их в nls.

m <- ns(x, df=5)
df <- data.frame(y, m)  # X-variables will be named X1, ... X5
# starting values should be set as appropriate for your data
nls(y ~ theta * plogis(alpha + b1*X1 + b2*X2 + b3*X3 + b4*X4 + b5*X5), data=df,
        start=list(theta=1, alpha=0, b1=1, b2=1, b3=1, b4=1, b5=1))

ETA: вот процесс автоматизации этого для различных значений df. Это создает формулу с использованием анализа текста, а затем использует do.call для вызова nls. Предостережение: не проверено.

my.nls <- function(x, y, df)
{
    m <- ns(x, df=df)
    xn <- colnames(m)
    b <- paste("b", seq_along(xn), sep="")
    fm <- formula(paste("y ~ theta * plogis(1 + alpha + ", paste(b, xn, sep="*",
          collapse=" + "), ")", sep=""))
    start <- c(1, 1, rep(1, length=length(b)))
    names(start) <- c("theta", "alpha", b)
    do.call(nls, list(fm, data=data.frame(y, m), start=start))
}
2 голосов
/ 05 февраля 2012

Осознание, к которому я пришел, когда разъяснял свой собственный вопрос, заставило меня понять, что есть менее неуклюжий путь, чем я видел раньше.

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

Уловка Хонг Ооя с использованием data.frame в матрице, сгенерированной ns для автоматического именования столбцов, довольно милои я использовал это ниже.Скорее всего, я буду использовать paste для их построения, потому что у меня есть несколько переменных, с которыми можно поиграть.

Предполагая, что набор данных, приведенный в вопросе -

lin.expr <- function(p,xn) {
  pn<-paste(p, 1:length(xn), sep = "")
  paste(paste(pn,xn,sep=" * "),collapse=" + ")
  }


m <- ns(x, df=3)
mydf <- data.frame(y, m)  # X-variables will be named X1, X2, ... 
xn <- names(mydf)[2:dim(mydf)[2]]

nspb <- lin.expr("b",xn)

c.form <- paste("y ~ theta * plogis( a + ",nspb,")",sep="")
stl <- list(theta=2, a=-5,b1=10, b2=10, b3=10)
nls( c.form, data=mydf, start= stl)

Моя фактическая формула будет иметь несколько терминов, таких как nspb.Значительные улучшения приветствуются;Я бы предпочел не выбирать свой собственный ответ, но я думаю, что выберу его, если в течение дня или двух ничего не изменится.

edit: добавление Hong Ooi (которое было опубликовано, когда я набирал свой ииспользует аналогичные идеи, но добавьте пару приятных дополнений) в значительной степени это делает;это приемлемый ответ, поэтому я проверил его.

...