Нелинейная модель с пятью параметрами (w / nls R) - PullRequest
Новая
       13

Нелинейная модель с пятью параметрами (w / nls R)

0 голосов
/ 19 ноября 2018

Это мой первый вопрос, пожалуйста, дайте мне знать, если я делаю что-то не так.У нас есть df с двумя переменными, и мы хотим смоделировать EPR (скорость производства яиц) как функцию температуры.

Соответствующие пакеты согласно странице nls: 

install.packages("tidyverse")
install.packages("nls.multstart")
install.packages("nlstools")
library(tidyverse)
library(nls.multstart) 
library(nlstools)

Соответствующие переменные от большего df: 

temp=c(9.2,9.9,12.7,12.8,14.3,14.5,16.3,16.5,18,18,19.6,19.6,19.9,19.9,22,22.4,23.2,23.4,25.3,25.6,27,27.3,28.5,30.3,20.9)
EPR=c(1.5,0,0,0,1.27,0.56,3.08,0.575,2.7,3.09,2,6.3,2,3.76,3.7,1.65,7.1,18.9,7.07,3.77,13.79,0,0,0.47,0)
df<-data.frame(temp,EPR)

Здесь я пишу формулу с пятью параметрамиДля оценки (k1, a, b, k2, c) temp будет значением x.Пока все хорошо. 

formula<-function(k1,a,b,k2,c,temp) {
 modelEPR<-k1*1/(1+exp(-a*(temp-b)))-k2*exp(c*temp)
 return(modelEPR)
}

Вот где я застрял;Я уже использую довольно узкий start_lower и upper, так как теперь я знаю параметры, используя Excel Solver несколько успешно.Значения, которые я получаю с помощью этого метода, получат мне модель, хотя и довольно неточную.Да, я дал старту нижний и верхний гораздо больший диапазон в начале, но это не дало лучших результатов. 

fit <- nls_multstart(EPR ~ formula(k1,a,b,k2,c,temp),
                 data = df,
                 iter = 100,
                 start_lower = c(k1 = 14, a = 0.3, b = 20, k2 = 0.02, c = 0.15),
                 start_upper = c(k1 = 15, a = 0.5, b = 21, k2 = 0.08, c = 0.24),
                 supp_errors = 'Y',
                 na.action = na.omit)

fit

Как уже упоминалось, я использовал решатель Excel для успешного создания модели иЯ получил оценки параметров, а затем попытался просто вручную вставить их здесь в R, что делает модель намного лучше.

model<-df %>%
 mutate(pred=(14.69/(1+exp(-0.41*(temp-20.52)))-0.05*exp(0.19 *temp))) %>% 
 ggplot()+
 xlab("Temperature (°C)")+
 ylab("EPR (Eggs per female per day")+
 geom_point(aes(temp,EPR))+
 geom_line(aes(temp,pred),col="red")

model

В конечном итоге у меня есть два вопроса;а) Что я делаю не так?Или это просто странные данные?Кажется, лучше работать с Excel ?!б) Как мне кодировать мост между посадкой и модельюfit даст 5 параметров, но как мне вставить их непосредственно в функцию модели?Могу ли я использовать мутацию как-то здесь?

Буду признателен за любую помощь!

1 Ответ

0 голосов
/ 20 ноября 2018

A. Начальные значения и примерка модели

Чтобы получить начальные значения:

  1. Если k1 = 0, то мы можем изменить формулу следующим образом, а затем использовать результат подбора этой линейной модели в качестве начального значения для c. 

    log(EPR) ~ log(k2) + c * temp

  2. b - это смещение в temp, а a - это масштабирование, поэтому выберите b = mean(temp) и a = 1/sd(temp)

  3. Мы можем использовать algorithm = "plinear", чтобы избежать необходимости указывать начальные значения для линейных параметров, то есть для k1 и k2. При использовании plinear правая часть формулы должна быть такой, чтобы k1 умножить на первый столбец плюс k2 умножить на второй столбец, чтобы получить предсказанное значение EPR.

Это дает следующее. Обратите внимание, что k1 и k2 будут представлены .lin1 и .lin2 в выводе nls. 

fm1 <- lm(log(EPR) ~ temp, df, subset = EPR > 0)
st2 <- list(c = coef(fm1)[[2]], a = 1/sd(df$temp), b = mean(df$temp))
fo2 <- EPR ~ cbind(1/(1+exp(-a*(temp-b))), -exp(c*temp))
fm2 <- nls(fo2, df, start = st2, algorithm = "plinear", 
  control = list(maxiter = 200))
deviance(fm2) # residual sum of squares
## [1] 333.6

Обратите внимание, что это представляет меньшую (лучшую) остаточную сумму квадратов, чем подгонка, показанная в вопросе: 

sum((df$EPR - pred)^2) # residual sum of squares for fit shown in question
## [1] 339.7

Пакеты не использовались.

Мы можем изобразить два совпадения, где подгонка вопроса синего цвета, а подгонка, выполненная здесь, красного цвета. Из графика возникает вопрос, являются ли два больших значения EFR выбросами и должны ли они быть исключены. 

plot(EPR ~ temp, df)
lines(fitted(fm2) ~ temp, df, subset = order(temp), col = "red")
lines(pred ~ temp, df, subset = order(temp), col = "blue")

[продолжение после скриншота]

screenshot

B. Оценка модели по заданным параметрам

Для данной модели, выраженной в формуле, мы можем оценить ее по заданным параметрам, используя пакет nls2. nls2 принимает те же аргументы, что и nls, но если начальное значение представляет собой фрейм данных с одной строкой, а алгоритм - "brute", то он просто возвращает значение правой части, оцененное по начальным значениям. См. ?nls для получения дополнительной информации. 

library(nls2)

fo <- EPR ~ k1*1/(1+exp(-a*(temp-b)))-k2*exp(c*temp)
st <- list(k1 = 14.69, a = 0.41, b = 20.52, k2 = 0.05, c = 0.19)
fm <- nls2(fo, df, start = data.frame(st), algorithm = "brute")

deviance(fm)
## [1] 339.7

fitted(fm) # predictions at parameter values given in st

или с точки зрения функции: 

rhs <- function(a, b, c, k1, k2, temp) k1*1/(1+exp(-a*(temp-b)))-k2*exp(c*temp)
p <- do.call("rhs", c(st, list(temp = df$temp)))
all.equal(p, pred)
## [1] TRUE
...