Блог

Главная - Вопросы по программированию - Как я могу сравнить параметры nls по группам в R

Как я могу сравнить параметры nls по группам в R

#r

Вопрос:

У меня есть такой набор данных

ме образец obs
1.5625 s 0.312
1.5625 s 0.302
3.125 s 0.335
3.125 s 0.333
6.25 s 0.423
6.25 s 0.391
12.5 s 0.562
12.5 s 0.56
25 s 0.84
25 s 0.843
50 s 1.202
50 s 1.185
100 s 1.408
100 s 1.338
200 s 1.42
200 s 1.37
1.5625 t 0.317
1.5625 t 0.313
3.125 t 0.345
3.125 t 0.343
6.25 t 0.413
6.25 t 0.404
12.5 t 0.577
12.5 t 0.557
25 t 0.863
25 t 0.862
50 t 1.22
50 t 1.197
100 t 1.395
100 t 1.364
200 t 1.425
200 t 1.415

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

 proc nlin data＝assay;
　model obs＝D＋(A－D)/(1＋(iu/((cs∗(sample＝“S”)
＋Ct∗(sample＝“T”))))∗∗(B));
　parms　D=1　B=1　Cs=1　Ct=1 A=1;
run;

Поэтому я пишу что-то вроде этого, а затем получаю

 nlm_1 <- nls(obs ~ (a - d) / (1   (iu / c[sample]) ^ b)   d, data = csf_1, start = list(a = 0.3, b = 1.8, c = c(25, 25), d = 1.4))
Error in numericDeriv(form[[3L]], names(ind), env) :
Missing value or an infinity produced when evaluating the model

Но без [sample] этого модель может быть рассчитана

 nlm_1 <- nls(obs ~ (a - d) / (1   (iu / c) ^ b)   d, data = csf_1, start = list(a = 0.3, b = 1.8, c = c(25), d = 1.4))
summary(nlm_1)
Formula: obs ~ (a - d)/(1   (iu/c)^b)   d

Parameters:
  Estimate Std. Error t value Pr(>|t|)
a  0.31590    0.00824   38.34   <2e-16 ***
b  1.83368    0.06962   26.34   <2e-16 ***
c 25.58422    0.55494   46.10   <2e-16 ***
d  1.44777    0.01171  123.63   <2e-16 ***
---
Signif. codes:
0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 0.02049 on 28 degrees of freedom

Number of iterations to convergence: 4
Achieved convergence tolerance: 6.721e-06

Я не понимаю, может ли кто-нибудь сказать мне, что не так с моим кодом и как я могу достичь своей цели с помощью R? Спасибо!

Комментарии:

1. Вы пытаетесь использовать sample категориальный индекс в качестве индекса, основанного на коде

Ответ №1:

Спасибо @akrun. После того, как я перешел csf_1$sample в фактор, я, наконец, получил то, что хотел.

 csf_1[, 2] <- as.factor(c(rep("s", 16), rep("t", 16)))
nlm_1 <- nls(obs ~ (a - d) / (1   (iu / c[sample]) ^ b)   d, data = csf_1, start = list(a = 0.3, b = 1.8, c = c(25, 25), d = 1.4))
summary(nlm_1)
Formula: obs ~ (a - d)/(1   (iu/c[sample])^b)   d

Parameters:
    Estimate Std. Error t value Pr(>|t|)
a   0.315874   0.008102   38.99   <2e-16 ***
b   1.833303   0.068432   26.79   <2e-16 ***
c1 26.075317   0.656779   39.70   <2e-16 ***
c2 25.114050   0.632787   39.69   <2e-16 ***
d   1.447901   0.011518  125.71   <2e-16 ***
---
Signif. codes:
0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 0.02015 on 27 degrees of freedom

Number of iterations to convergence: 4
Achieved convergence tolerance: 6.225e-06
Метки: