如何将 nlsLM 函数与 R 中的应用族函数之一一起使用

Question

我需要有关如何按列进行多重回归的指南。 我有一个数据框，我想在其中分别获取每列拟合系数。到目前为止，我只能为一列获得这些结果。

到目前为止我尝试了什么

1. 也许将结果分配给一个新变量

   (model.out1

可能有效，但我不想写几个拟合方程，每次都说大约 15 和列名。这不是优雅的解决方案。

 2. using `apply` function 

aa <- apply(df[4:8],2,fit_function)


    fit_function <- function(x){nlsLM(x~ifelse(df$direc=="North"&V<J1, exp((-t_pw)/f0*exp(-del1*(1-V/J1)^2)),1)*ifelse(df$direc=="South"&V>J2, exp((-t_pw)/f0*exp(-del2*(1-V/J2)^2)),1)
    ,data=df,start=c(del1=5,J1=15,del2=1,J2=-5),trace=T)}

给出一个我们通常知道的错误

nlsModel(formula, mf, start, wts) 中的错误：奇异梯度 初始参数估计矩阵

也许将这些列分开并拟合它们中的每一列并结合拟合系数可能会起作用。但我不知道该怎么做。

这里是您检查有效性的可重复数据
df

direc <- rep(rep(c("North","South"),each=10),times=6)
  V <- rep(c(seq(2,40,length.out=10),seq(-2,-40,length.out=10)),times=1)
  DQ0 = c(replicate(2, sort(runif(10,0.001,1))))
  DQ1 = c(replicate(2, sort(runif(10,0.001,1))))
  DQ2 = c(replicate(2, sort(runif(10,0.001,1))))
  DQ3 = c(replicate(2, sort(runif(10,0.001,1))))
  DQ4 = c(replicate(2, sort(runif(10,0.001,1))))
  group  =  c(replicate(1,rep(letters[1:6],each=20)))  
df <- data.frame(group,direc,V,DQ0,DQ1,DQ2,DQ3,DQ4)

library(minpack.pl)

因为我想对所有列 DQ0、DQ1、DQ2、DQ3、DQ4 进行拟合，所以我写下了这个函数。

拟合函数

 f0<-1e-9
 t_pw<-3e-8

nls_fit=nlsLM(DQ0~ifelse(df$direc=="North"&V<J1, exp((-t_pw)/f0*exp(-del1*(1-V/J1)^2)),1)*ifelse(df$direc=="South"&V>J2, exp((-t_pw)/f0*exp(-del2*(1-V/J2)^2)),1)
        ,data=df,start=c(del1=5,J1=15,del2=1,J2=-5),trace=T)

并在每个组内获得拟合结果。

df_new<- df%>%
  group_by(group)%>%
  do(data.frame(model=tidy(nls_fit)))%>%
  select_("delta"="model.term","value"= "model.estimate")

如何获得 DQ1、DQ2、DQ3 和 DQ4 的拟合结果作为表格。也许这样的东西更可取

    group delta  value_DQ0  value_DQ1   value_DQ2  value_DQ3 value_DQ4 
1      a  del1   4.962564       *           *          *         * 
2      a    J1  14.666667       *           *          *         *
3      a  del2   3.496986       *           *          *         *
4      a    J2 -14.468551
5      b  del1   4.962564
6      b    J1  14.666667
7      b  del2   3.496986
8      b    J2 -14.468551
9      c  del1   4.962564
10     c    J1  14.666667
..   ...   ...        ...

编辑 我找到了这个Help with lm and multiple linear regression 也许我可以这样做

dat <- data.frame(x=1:10,y=rnorm(10),z=10:1)
lm(x~., data=dat)

但是当我像上面那样用 DQ0 替换 if else 部分时，我得到了这个错误

可能我错过了某些部分。关于this_你能给出一个明确的答案吗？任何帮助将不胜感激。

Answer 1

首先，我对您的方法深表怀疑。您可能知道，非线性回归是一个迭代过程，其成功很大程度上取决于初始估计的选择。不仅如此，您还必须考虑局部最小值的可能性，当然您还需要评估拟合优度，例如通过查看参数的 p 值并测试残差的正态性。您的模型函数相当复杂，因此尝试自动化这样的过程根本不可能产生结果，即使它确实产生了，您也不能保证结果将是有意义的。至少您需要绘制所有情况下的数据与模型函数。只是生成这样的表是自找麻烦。

其次，您的示例有几个问题。您的模型函数取决于 t_pw 和 f0，您没有在任何地方定义 AFAICT，nlsLM(...) 在包 minpack.lm 中，而不是 minpack.pl（我无法在任何地方找到后者）。

说了这么多，我可以看到您已经付出了很多努力来制定这个问题，以及基本问题：如何使用数据集拆分组对任意响应列表运行非线性回归 -明智，很有趣。这是使用mtcars 数据集的一种方法。在本例中，分组变量为cyl（气缸数），响应变量为mpg、qsec和hp，（非常简单的）模型函数为：y ~ a * wt / (b + wt)，带参数@ 987654332@ 和b。因此，对于每个气缸类别（4、6 和 8），我们将 mpg、qsec 和 hp 建模为 wt 的函数，并确定 a 和 b。

df   <- mtcars                # safer to make a copy
resp <- c("mpg","qsec","hp")  # response variable names
library(minpack.lm)           # for nlsLM(...)
get.coefs <- function(y,df) {
  fit <- nlsLM(y~a*wt/(b+wt), data=data.frame(y=y,df), start=c(a=1,b=-1))
  coef(fit)
}
coefs  <- lapply(split(df,df$cyl),function(df) {do.call(cbind,lapply(df[resp],get.coefs,df))})
result <- do.call(rbind,lapply(names(coefs),function(x) {
  data.frame(group=x, var=rownames(coefs[[x]]), coefs[[x]])
}))
result
#    group var        mpg      qsec           hp
# a      4   a 18.2436308 24.517564  98.80184109
# b      4   b -0.6655570  0.615073   0.42670565
# a1     6   a 14.2066098 62.179060  83.26572253
# b1     6   b -0.8599662  7.639224  -0.97768640
# a2     8   a  9.2212533 21.977518 204.59139171
# b2     8   b -1.4931033  1.213256  -0.08582505

在上面的代码中，函数get.coefs(...) 接受一个包含响应变量的向量y 和一个包含数据集的data.frame df，运行回归并返回一个系数向量。

coefs <- ... 行完成了大部分工作。内部lapply(...) 将响应的每一列依次传递给get.coefs(...)，并将结果作为列表返回。 do.call(cbind,...) 将列表元素组合成一个系数矩阵，系数以行为单位，响应变量以列为单位。外部lapply(...) 按组（在本例中为圆柱体）拆分原始data.frame，并将每个分组子集提交给上述过程。所有这一切的结果，coefs 是一个矩阵列表，每个组一个。

最后一行：result <- ... 只是将coefs 列表重新格式化为您想要的表格。