将SAS中的PROC NLIN转换为R

Question

最近我被分配了一项工作，需要将SAS代码翻译成R。到目前为止，我已经成功地完成了80%的工作，但卡在了使用PROC NLIN的部分。根据我的阅读，PROC NLIN用于拟合非线性模型，但我不确定该代码是否实际上正在这样做，因此不知道如何在R中处理它。以下是该代码：

proc nlin data=ds1 outest=estout;
 parms ET= 0 f= 10.68;
  E= f- R*(1-ET*M); 
  L   = E*E;
  model.like = sqrt(E*E);
  by Name ; 
run;

样例数据如下 -

Name    M           R
Anna    0.5456231   4.118197
Anna    0.5359164   4.240243
Anna    0.541881    3.943975
Anna    0.5436047   3.822222
Anna    0.5522962   3.58813
Anna    0.5561487   3.513195
Anna    0.5423374   3.666507
Anna    0.525836    3.715371
Anna    0.5209941   3.805572
Anna    0.5304675   3.750689
Anna    0.5232541   3.788292

当我查看 SAS 帮助中的 PROC NLIN 页面时，其中使用了参数“MODEL”来指定方程，但是这里的代码没有模型方程。Model.like 用于指定似然函数（第4316页-https://support.sas.com/documentation/cdl/en/statugnlin/61811/PDF/default/statugnlin.pdf）。那么这段代码是在做什么？我完全困惑了。起初我认为可以使用 R 中的 nls() 来完成这个任务，我尝试了以下代码——

fit = nls(E~ f - R*(1-eta*M),sample, start=list(eta=0,phi=10.86)
      ,trace=T)

但是我很快意识到这是错误的，因为即使进行了5000次迭代，模型仍然没有收敛。这是因为我的数据集中没有'E'列。那么，SAS是如何做到的呢？

非常感谢您的帮助！

- RHelp

Answer 1

首先让我们弄清楚SAS代码在做什么。PROC NLIN可以被欺骗成执行各种最小化问题，但设置有时很反直觉。您需要定义一个因变量（$y$）和一个基于其他变量和一些参数的预测值（$f(x,\beta)$），它将最小化$\sum_i [y_i - f(x_i,\beta)]^2$。

定义$y$和$f$的关键是以下行：

model.like = sqrt(E*E)

这相当于

model like = sqrt(E*E)

因此，这意味着 $\sum [like - \sqrt{E\cdot E}]^2$ 将被最小化。根据您提供的示例，我认为变量 like 已经在之前被定义，并且已经被设置为常数0。这意味着正在最小化 $\sum [0- \sqrt{E\cdot E}]^2 = \sum E^2$。 E 被定义为 f- R*(1-ET*M)，因此实际上正在最小化 $\sum [f- R*(1-ET*M)]^2$，其中 f 和 ET 是未知参数。我不确定它的含义，但这就是正在发生的事情。

将其重写为 R 确实可以使用 nls，我们可以使用相同的技巧：预测零。

sample <- read.table(textConnection(
"Name    M           R
 Anna    0.5456231   4.118197
 Anna    0.5359164   4.240243
 Anna    0.541881    3.943975
 Anna    0.5436047   3.822222
 Anna    0.5522962   3.58813
 Anna    0.5561487   3.513195
 Anna    0.5423374   3.666507
 Anna    0.525836    3.715371
 Anna    0.5209941   3.805572
 Anna    0.5304675   3.750689
 Anna    0.5232541   3.788292"), header=TRUE)

nls(0 ~ f - R*(1-eta*M), data=sample, start=list(eta=0,f=10.86), trace=T)

输出结果

546.5988 :   0.00 10.86
0.06273518 :  1.7259120 0.2731282
Nonlinear regression model
  model: 0 ~ f - R * (1 - eta * M)
   data: sample
   eta      f 
1.7259 0.2731 
 residual sum-of-squares: 0.06274

Number of iterations to convergence: 1 
Achieved convergence tolerance: 4.345e-07

请注意，SAS代码是按名称执行的，因此您需要确保R代码适用于每个名称的不同模型。

- Aniko

---

非常感谢！这正是我想要理解的！ - RHelp

Answer 2

使用tidyverse的方法

library(tidyverse)
library(broom)

sample <- read.table(textConnection(
  "Name    M           R
  Anna    0.5456231   4.118197
  Anna    0.5359164   4.240243
  Anna    0.541881    3.943975
  Anna    0.5436047   3.822222
  Anna    0.5522962   3.58813
  Anna    0.5561487   3.513195
  Anna    0.5423374   3.666507
  Anna    0.525836    3.715371
  Anna    0.5209941   3.805572
  Anna    0.5304675   3.750689
  Anna    0.5232541   3.788292"), header=TRUE)


x <- sample %>%
  group_by(Name) %>%
  nest() %>%
  mutate(
    model = data %>% map(~nls(0 ~ f - R*(1-eta*M), data= . , start=list(eta=0,f=10.86), trace=T)),
    coef = map(model, tidy),
    quali = map(model, glance),
    resid = map(model, augment)
  )

unnest(select(x, coef))
# A tibble: 2 x 6
    Name  term  estimate std.error statistic      p.value
  <fctr> <chr>     <dbl>     <dbl>     <dbl>        <dbl>
1   Anna   eta 1.7259120 0.2260999 7.6334045 3.213398e-05
2   Anna     f 0.2731282 0.4645288 0.5879683 5.710103e-01

unnest(select(x, quali))
# A tibble: 1 x 8
       sigma isConv       finTol   logLik       AIC       BIC   deviance df.residual
       <dbl>  <lgl>        <dbl>    <dbl>     <dbl>     <dbl>      <dbl>       <int>
1 0.08348998   TRUE 4.345363e-07 12.80868 -19.61736 -18.42368 0.06273518           9

- Italo Cegatta

---

谢谢Italo！我能问一下如何为每个组（在这种情况下是每个名称）传递不同的起始参数吗？ - Ale