Экспоненциальная кривая с nls с использованием групп data.table - PullRequest
Новая
       67

Экспоненциальная кривая с nls с использованием групп data.table

0 голосов
/ 26 сентября 2018

Я хотел бы подогнать экспоненциальные кривые к группам 1 и 2 в таблице данных, показанной ниже, и получить новый столбец, содержащий остаточную стандартную ошибку, соответствующую каждой группе.Экспоненциальная кривая должна следовать за y=a*exp(b*x)+c

## Example data table
DT <- data.table(
x = c(1,2,3,4,5,6,7,8,1,2,3,4,5,6,7,8),
y = c(15.4,16,16.4,17.7,20,23,27,35,25.4,26,26.4,27.7,30,33,37,45),
groups = c(1,1,1,1,1,1,1,1,2,2,2,2,2,2,2,2)

Однако я знаю только, как подгонять кривые nls и получать остаточную стандартную ошибку отдельных групп, используя код, приведенный ниже, который оценивает хорошие начальные параметры a , b и c : 

subsetDT <- DT[group == 1]
c.0 <- min(subsetDT[,y]) * 0.5
model.0 <- lm(log(y- c.0) ~ x, data=subsetDT)
start <- list(a=exp(coef(model.0)[1]), b=coef(model.0)[2], c=c.0)
model <- nls(y ~ a * exp(b * x) + c,
         data = subsetDT, start = start, 
         control = nls.control(maxiter=500))
sigma <- summary(model)$sigma

Я не хочу подгруппировать DT группой в цикле для вычисления sigma идругая информация о модели.

Я знаю, что если бы я использовал lm, я мог бы сделать следующее, чтобы получить новые столбцы, содержащие информацию о модели: 

DT[, `:=` (r.squared=summary(lm(log(y)~x))$r.squared,
           int=coef(lm(log(y)~x))[1],
           coeff=coef(lm(log(y)~x))[2]
          ), by=c("groups")]

Как я могуиспользуйте :=, чтобы соответствовать экспоненциальной кривой и включить мои nls параметры a , b и c ?

1 Ответ

0 голосов
/ 26 сентября 2018

Если вы хотите добавить сигма, a, b, c в качестве новых столбцов в исходном наборе данных, вы можете сделать следующее: 

DT[, c("sigma", "a", "b", "c") := {
        c.0 <- min(y) * 0.5
        model.0 <- lm(log(y - c.0) ~ x, data=.SD)
        start <- list(a=exp(coef(model.0)[1]), b=coef(model.0)[2], c=c.0)
        model <- nls(y ~ a * exp(b * x) + c,
            data=.SD, 
            start=start, 
            control=nls.control(maxiter=500))
        c(.(sigma=summary(model)$sigma), as.list(coef(model)))
    },
    by=.(groups)]

output: 

    x    y groups     sigma         a         b        c
 1: 1 15.4      1 0.2986243 0.5265405 0.4565363 14.56728
 2: 2 16.0      1 0.2986243 0.5265405 0.4565363 14.56728
 3: 3 16.4      1 0.2986243 0.5265405 0.4565363 14.56728
 4: 4 17.7      1 0.2986243 0.5265405 0.4565363 14.56728
 5: 5 20.0      1 0.2986243 0.5265405 0.4565363 14.56728
 6: 6 23.0      1 0.2986243 0.5265405 0.4565363 14.56728
 7: 7 27.0      1 0.2986243 0.5265405 0.4565363 14.56728
 8: 8 35.0      1 0.2986243 0.5265405 0.4565363 14.56728
 9: 1 25.4      2 0.2986243 0.5265404 0.4565363 24.56728
10: 2 26.0      2 0.2986243 0.5265404 0.4565363 24.56728
11: 3 26.4      2 0.2986243 0.5265404 0.4565363 24.56728
12: 4 27.7      2 0.2986243 0.5265404 0.4565363 24.56728
13: 5 30.0      2 0.2986243 0.5265404 0.4565363 24.56728
14: 6 33.0      2 0.2986243 0.5265404 0.4565363 24.56728
15: 7 37.0      2 0.2986243 0.5265404 0.4565363 24.56728
16: 8 45.0      2 0.2986243 0.5265404 0.4565363 24.56728
...