Второй шаг с нелинейной регрессией в R: добавление предикторов

Оригинал: Second step with non-linear regression: adding predictors

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

Пример, который я буду использовать в этом посте - это логистическая функция роста, она часто используется в экологии для моделирования роста населения. Например, предположим, вы подсчитываете количество бактериальных клеток в чашке Петри, вначале количество клеток будет расти экспоненциально, но через некоторое время из-за ограничений в ресурсах (будь то пространство или еда) популяция бактерий достигнет равновесия. В результате будет получена классическая S-образная нелинейная логистическая функция роста. Логистическая функция роста имеет три параметра: темп роста, называемый «r», размер популяции в состоянии равновесия, называемый «K», и размер популяции в начале, называемый «n0».

Первый пример: категориальный предиктор

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

#load libraries
library(nlme)
#first try effect of treatment on logistic growth
Ks <- c(100,200,150) 
n0 <- c(5,5,6)
r <- c(0.15,0.2,0.15)
time <- 1:50
#this function returns population dynamics following
#a logistic curves
logF <- function(time,K,n0,r){
  d <- K * n0 * exp(r*time) / (K + n0 * (exp(r*time) - 1))
  return(d)
}
#simulate some data
dat <- data.frame(Treatment=character(),Time=numeric(),
                  Abundance=numeric())
for(i in 1:3){
  Ab <- logF(time = time,K=Ks[i],n0=n0[i],r=r[i])
  tmp <- data.frame(Treatment=paste0("T",i),Time=time,
                    Abundance=Ab+rnorm(time,0,5))
  #note that random deviates were added to the simulated 
  #population density values
  dat <-rbind(dat,tmp)
}
#the formula for the models
lF<-formula(Abundance~K*n0*exp(r*Time)/(K+n0*(exp(r*Time)-1)) | Treatment)

#fit the model
(m <- nlsList(lF,data=dat,start=list(K=150,n0=10,r=0.5)))
Call:
  Model: Abundance ~ K * n0 * exp(r * Time)/(K + n0 * (exp(r * Time) - 1)) | Treatment 
   Data: dat 

Coefficients:
          K       n0         r
T1 105.2532 4.996502 0.1470740
T2 199.5149 4.841359 0.2006150
T3 150.2882 5.065196 0.1595293

Degrees of freedom: 150 total; 141 residual
Residual standard error: 5.085229

Код моделировал значения совокупности с использованием трех наборов параметров (r, K и n0). Затем мы задали формулу нелинейной регрессии, используя вертикальную черту «|» символ, чтобы явно просить подобрать различные параметры для каждого типа пищи. Выходные данные модели дают нам оценочные параметры для каждого типа пищи. Теперь мы можем сопоставить прогнозы модели с реальными данными:

#derive the predicted lines
dat$Pred <-predict(m)
#plot
ggplot(dat,aes(x=Time,y=Abundance,color=Treatment))+
       geom_point()+geom_line(aes(y=Pred))

Второй пример: непрерывный предиктор

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

Для этого я воспользуюсь другим пакетом: bbmle и его функцией mle2:

#load libraries
library(bbmle)
library(reshape2)
#parameter for simulation
K <- 200
n0 <- 5
#the gradient in temperature
Temp <- ceiling(seq(0,20,length=10))
 

#simulate some data
mm <- sapply(Temp,function(x){ 
             rpois(50,logF(time = 1:50,K = K,n0=n0,r = 0.05 + 0.01 * x))})
#sample data from a poisson distribution with lambda parameter equal
#to the expected value, note that we do not provide one value for the
#parameter r but a vector of values representing the effect of temperature
#on r
#some reshaping
datT <- melt(mm)
names(datT) <- c("Time","Temp","Abund")
datT$Temp <- Temp[datT$Temp]

#fit the model
(mll <- mle2(Abund~dpois(K*n0*exp(r*Time)/(K+n0*(exp(r*Time)-1))),
             data=datT,parameters = list(r~Temp),
             start=list(K=100,n0=10,r=0.05)))
Call:
mle2(minuslogl = Abund ~ dpois(K * n0 * exp(r * Time)/(K + n0 * 
    (exp(r * Time) - 1))), start = list(K = 100, n0 = 10, r = 0.05), 
    data = datT, parameters = list(r ~ Temp))

Coefficients:
            K            n0 r.(Intercept)        r.Temp 
 200.13176811    4.60546966    0.05195730    0.01016994 

Log-likelihood: -1744.01 

Первый аргумент для mle2 - это либо функция, возвращающая отрицательное логарифмическое правдоподобие, либо формула вида «ответ ~ некоторое распределение (ожидаемое значение)». В нашем случае мы подбираем распределение Пуассона с ожидаемыми значениями, полученными из логистического уравнения и его трех параметров. Зависимость между скоростью роста и температурой задается с помощью аргумента «parameter». Выходные данные базовой модели дают нам все оценочные параметры, представляющие интерес здесь: «r. (Intercept)», который представляет собой скорость роста, когда температура равна 0, и «r.Temp», которая представляет собой увеличение скорости роста для каждого повышения температуры на 1 градус.

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

datT$pred <- predict(mll)
ggplot(datT,aes(x=Time,y=Abund,color=factor(Temp)))+
       geom_point()+geom_line(aes(y=pred))

Третий пример: больше непрерывных предикторов

Крутая вещь с mle2 заключается в том, что вы можете подбирать любые модели, которые вы можете себе представить, если вы можете записывать функции логарифмического правдоподобия. Мы увидим это с расширением предыдущей модели. В биологии принято считать, что виды должны иметь некоторую оптимальную температуру, поэтому мы можем ожидать колоколообразной связи между равновесием популяции и температурой. Итак, мы добавим к предыдущей модели квадратичную зависимость между равновесием населения («K») и температурой.

#simulate some data
mm <- sapply(Temp,function(x){
 rpois(50,logF(time = 1:50,K = 100+20*x-x**2,n0=n0,r = 0.05 + 0.01 * x))})
#note that this time both K and r are given as vectors to represent their
#relation with temperature

#some reshaping
datT <- melt(mm)
names(datT) <- c("Time","Temp","Abund")
datT$Temp <- Temp[datT$Temp]

#the negative log-likelihood function
LL <- function(k0,k1,k2,n0,r0,r1){
  K <- k0 + k1*Temp + k2*Temp**2
  r <- r0 + r1*Temp
  lbd <- K*n0*exp(r*Time)/(K+n0*(exp(r*Time)-1))
  -sum(dpois(Abund,lbd,log=TRUE))
}

(mll2 <- mle2(LL,data=datT,start=list(k0=100,k1=1,k2=-0.1,n0=10,r0=0.1,r1=0.1)))
Call:
mle2(minuslogl = LL, start = list(k0 = 100, k1 = 1, k2 = -0.1, 
    n0 = 10, r0 = 0.1, r1 = 0.1), data = datT)

Coefficients:
          k0           k1           k2           n0           r0           r1 
100.47610915  20.09677780  -1.00836974   4.84940342   0.04960558   0.01018838 

Log-likelihood: -1700.91 

На этот раз я определил функцию отрицательного логарифмического правдоподобия, чтобы она соответствовала модели. Эта функция принимает в качестве аргумента параметр модели. В этом случае у нас есть 3 параметра для описания квадратичной связи между K и температурой, 1 параметр для постоянного n0 и 2 параметра для описания линейной связи между r и температурой. Важно понимать, что делает последняя строка. В основном, вызов dpois возвращает вероятности получения наблюдаемых значений численности («Abund») на основе ожидаемого среднего значения («lbd»). Функция LL вернет одно отрицательное значение логарифмического правдоподобия (сумму), а задача mle2 - минимизировать это значение. Для этого mle2 будет пробовать много разных комбинаций параметров каждый раз, вызывая функцию LL, и когда она достигнет минимума, она вернет использованный набор параметров. Это оценка максимального правдоподобия (MLE). Как и раньше, выходные данные модели дают нам оценочные параметры, и мы можем использовать их для построения подобранных регрессий:

cc <- coef(mll2)
#sorry for the awful coding, to get the model predicted values we need
#to code the relation by hand which is rather clumsy ...
datT$pred <- melt(sapply(Temp,function(x) {logF(time=1:50,K=cc[1]+cc[2]*x+cc[3]*x**2,n0=cc[4],r=cc[5]+cc[6]*x)}))[,3]
ggplot(datT,aes(x=Time,y=Abund,color=factor(Temp)))+
       geom_point()+geom_line(aes(y=pred))

Заключение

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