如何将 nlsLM 函数与 R 中的应用族函数之一一起使用
How to Use nlsLM function together with one of the apply family function in R
我需要有关如何按列进行多元回归的指南。
我有一个数据框，我想在其中分别获取每一列的拟合系数。到目前为止，我只能获得一列的结果。
到目前为止我尝试了什么
也许将结果分配给一个新变量
(model.out1 <- lm(y1~x1))
(model.out2 <- lm(y2~x2))
可能有用，但我不想写几个拟合方程让我们每次说 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)}
给出一个我们通常知道的错误
Error in nlsModel(formula, mf, start, wts) : singular gradient
matrix at initial parameter estimates
也许将这些列分开并拟合它们中的每一列并组合拟合系数可能会起作用。但我不知道该怎么做。
这是可重现的数据，供您检查有效性
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 的拟合结果 table。也许这样的东西更可取
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 部分时，我得到了这个错误
可能我遗漏了一些部分。你能对此给出一些明确的答案吗_？不，我们将不胜感激。
首先，我对你的做法深表怀疑。您可能知道，非线性回归是一个迭代过程，其成功在很大程度上取决于起始估计值的选择。不仅如此，您还必须考虑局部最小值的可能性，当然您还需要评估拟合优度，例如通过查看参数的 p 值和检验残差的正态性。您的模型功能相当复杂，因此尝试使这样的过程自动化根本不可能产生结果，即使产生结果也不能保证结果有意义。至少您需要绘制所有情况下的数据与模型函数的关系图。像这样生成一个table是自找麻烦。
其次，你的例子有几个问题。您的模型函数取决于 t_pw 和 f0，AFAICT 您没有在任何地方定义，并且 nlsLM(...) 在包 minpack.lm 中，而不是 minpack.pl（我有无法在任何地方找到后者）。
说了那么多，我可以看出您在制定这个问题和基本问题上付出了很多努力：如何 运行 对任意响应列表进行非线性回归，将数据集按组拆分，很有趣。这是使用 mtcars 数据集执行此操作的一种方法。在此示例中，分组变量为 cyl（气缸数），响应变量为 mpg、qsec 和 hp，（非常简单的）模型函数为: y ~ a * wt / (b + wt)，参数为 a 和 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，运行 s 回归和 returns 系数向量。
行 coefs <- ... 完成了大部分工作。内部 lapply(...) 依次将每一列响应传递给 get.coefs(...)，并将结果作为列表传递给 returns。 do.call(cbind,...) 将列表元素组装成一个系数矩阵，系数在行中，响应变量在列中。外部 lapply(...) 按组（在本例中为圆柱体）拆分原始 data.frame，并将每个分组的子集提交到上述过程。所有这一切的结果，coefs 是一个矩阵列表，每组一个。
最后一行：result <- ... 只是将 coefs 列表重新格式化为您想要的 table。
我需要有关如何按列进行多元回归的指南。 我有一个数据框，我想在其中分别获取每一列的拟合系数。到目前为止，我只能获得一列的结果。
到目前为止我尝试了什么
也许将结果分配给一个新变量
(model.out1 <- lm(y1~x1)) (model.out2 <- lm(y2~x2))
可能有用，但我不想写几个拟合方程让我们每次说 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)}
给出一个我们通常知道的错误
Error in nlsModel(formula, mf, start, wts) : singular gradient matrix at initial parameter estimates
也许将这些列分开并拟合它们中的每一列并组合拟合系数可能会起作用。但我不知道该怎么做。
这是可重现的数据，供您检查有效性
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 的拟合结果 table。也许这样的东西更可取
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 部分时，我得到了这个错误
可能我遗漏了一些部分。你能对此给出一些明确的答案吗_？不，我们将不胜感激。
首先，我对你的做法深表怀疑。您可能知道，非线性回归是一个迭代过程，其成功在很大程度上取决于起始估计值的选择。不仅如此，您还必须考虑局部最小值的可能性，当然您还需要评估拟合优度，例如通过查看参数的 p 值和检验残差的正态性。您的模型功能相当复杂，因此尝试使这样的过程自动化根本不可能产生结果，即使产生结果也不能保证结果有意义。至少您需要绘制所有情况下的数据与模型函数的关系图。像这样生成一个table是自找麻烦。
其次，你的例子有几个问题。您的模型函数取决于 t_pw 和 f0，AFAICT 您没有在任何地方定义，并且 nlsLM(...) 在包 minpack.lm 中，而不是 minpack.pl（我有无法在任何地方找到后者）。
说了那么多，我可以看出您在制定这个问题和基本问题上付出了很多努力：如何 运行 对任意响应列表进行非线性回归，将数据集按组拆分，很有趣。这是使用 mtcars 数据集执行此操作的一种方法。在此示例中，分组变量为 cyl（气缸数），响应变量为 mpg、qsec 和 hp，（非常简单的）模型函数为: y ~ a * wt / (b + wt)，参数为 a 和 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，运行 s 回归和 returns 系数向量。
行 coefs <- ... 完成了大部分工作。内部 lapply(...) 依次将每一列响应传递给 get.coefs(...)，并将结果作为列表传递给 returns。 do.call(cbind,...) 将列表元素组装成一个系数矩阵，系数在行中，响应变量在列中。外部 lapply(...) 按组（在本例中为圆柱体）拆分原始 data.frame，并将每个分组的子集提交到上述过程。所有这一切的结果，coefs 是一个矩阵列表，每组一个。
最后一行：result <- ... 只是将 coefs 列表重新格式化为您想要的 table。