Первые шаги с нелинейной регрессией в R

Проведение прямой линии через облако точек (т. е. выполнение линейной регрессии) - это самый простой анализ данных, который можно сделать. Иногда он хорошо согласуется с данными, но в некоторых (многих) ситуациях отношения между переменными не являются линейными. В этом случае можно использовать три различных способа: (1) попытаться линеаризовать взаимосвязь путем преобразования данных, (2) подогнать полиномиальные или сложные модели сплайнов к данным или (3) подогнать нелинейные функции к данным. Как вы уже догадались из названия, этот пост будет посвящен третьему варианту.

Что такое нелинейная регрессия?

В нелинейной регрессии аналитик определяет функцию с набором параметров, соответствующих данным. Самый простой способ оценить такие параметры - использовать метод нелинейных наименьших квадратов (функция nls в R), который в основном аппроксимирует нелинейную функцию с помощью линейной и итеративно пытается найти лучшие значения параметров (wiki). Приятной особенностью нелинейной регрессии в прикладном контексте является то, что оцениваемые параметры имеют четкую интерпретацию (Vmax в модели Михаэлиса-Ментен - это максимальная скорость), которую, например, будет сложнее получить, используя линейные модели на преобразованных данных.

Подгонка нелинейной функции методом наименьших квадратов

Первый пример с использованием уравнения Михаэлиса-Ментен:

#simulate some data
set.seed(20160227)
x<-seq(0,50,1)
y<-((runif(1,10,20)*x)/(runif(1,0,10)+x))+rnorm(51,0,1)
#for simple models nls find good starting values for the parameters even if it throw a warning
m<-nls(y~a*x/(b+x))
#get some estimation of goodness of fit
cor(y,predict(m))
[1] 0.9496598

Код ниже строит график:

#plot
plot(x,y)
lines(x,predict(m),lty=2,col="red",lwd=3)

Проблема поиска правильных начальных значений

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

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

#simulate some data, this without a priori knowledge of the parameter value
y<-runif(1,5,15)*exp(-runif(1,0.01,0.05)*x)+rnorm(51,0,0.5)
#visually estimate some starting parameter values
plot(x,y)
#from this graph set approximate starting values
a_start<-8 #param a is the y value when x=0
b_start<-2*log(2)/a_start #b is the decay rate
#model
m<-nls(y~a*exp(-b*x),start=list(a=a_start,b=b_start))
#get some estimation of goodness of fit
cor(y,predict(m))
#plot the fit
lines(x,predict(m),col="red",lty=2,lwd=3)
[1] 0.9811831

Использование функции SelfStarting

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

где Nt - количество особей в момент времени t, r - скорость роста популяции, а K - емкость. Мы можем переписать это как дифференциальное уравнение:

library(deSolve)
#simulating some population growth from the logistic equation and estimating the parameters using nls
log_growth <- function(Time, State, Pars) {
  with(as.list(c(State, Pars)), {
    dN <- R*N*(1-N/K)
    return(list(c(dN)))
  })
}
#the parameters for the logisitc growth
pars  <- c(R=0.2,K=1000)
#the initial numbers
N_ini  <- c(N=1)
#the time step to evaluate the ODE
times <- seq(0, 50, by = 1)
#the ODE
out   <- ode(N_ini, times, log_growth, pars)
#add some random variation to it
N_obs<-out[,2]+rnorm(51,0,50)
#numbers cannot go lower than 1
N_obs<-ifelse(N_obs<1,1,N_obs)
#plot
plot(times,N_obs)

Эта часть была просто имитацией некоторых данных со случайной ошибкой, теперь переходим к сложной части для оценки начальных значений.

Теперь в R есть встроенная функция для оценки начальных значений параметра логистического уравнения (SSlogis), но она использует следующее уравнение:

#find the parameters for the equation
SS<-getInitial(N_obs~SSlogis(times,alpha,xmid,scale),data=data.frame(N_obs=N_obs,times=times))

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

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

#we used a different parametrization
K_start<-SS["alpha"]
R_start<-1/SS["scale"]
N0_start<-SS["alpha"]/(exp(SS["xmid"]/SS["scale"])+1)
#the formula for the model
log_formula<-formula(N_obs~K*N0*exp(R*times)/(K+N0*(exp(R*times)-1)))
#fit the model
m<-nls(log_formula,start=list(K=K_start,R=R_start,N0=N0_start))
#estimated parameters
summary(m)
#get some estimation of goodness of fit
cor(N_obs,predict(m))
Formula: N_obs ~ K * N0 * exp(R * times)/(K + N0 * (exp(R * times) - 1))

Parameters:
    Estimate Std. Error t value Pr(>|t|)    
K  1.012e+03  3.446e+01  29.366   <2e-16 ***
R  2.010e-01  1.504e-02  13.360   <2e-16 ***
N0 9.600e-01  4.582e-01   2.095   0.0415 *  
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 49.01 on 48 degrees of freedom

Number of iterations to convergence: 1 
Achieved convergence tolerance: 1.537e-06

[1] 0.9910316

Построим график:

lines(times,predict(m),col="red",lty=2,lwd=3)

Было немного сложно перейти от параметризации SSlogis к нашей собственной, но оно того стоило!

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

Оригинал: First steps with Non-Linear Regression in R