0

我正在使用 dplyr 和 broom 组合,并尝试根据数据组内部的条件拟合回归模型。最后我想提取每组的回归系数。

到目前为止,我对所有组都得到了相同的拟合结果(每个组用字母分隔a:f)。这是主要问题。

library(dplyr)
library(minpack.lm)
library(broom)

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

df <-   data.frame(direc,V,DQ0,DQ1,DQ2,DQ3,No)

拟合条件可以描述如下; direc=North和 if V<J1do 与方程拟合exp((-t_pw)/f0*exp(-del1*(1-V/J1)^2))else ifdirec=SouthV>J2 do 与相同的方程拟合。在这两种情况下,如果V<J1&V>J2不满足返回1每种情况。

更新 我发现通过此链接中的建议,条件nls可以是可能的conditional-formula-for-nls 。

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=1,J1=15,del2=1,J2=-15),trace=T) 

nls_fit

Nonlinear regression model
  model: 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
   del1      J1    del2      J2 
  1.133  23.541   1.079 -20.528 
 residual sum-of-squares: 16.93

Number of iterations to convergence: 4 
Achieved convergence tolerance: 1.49e-08

另一方面,当我尝试拟合其他列时,例如 DQ1、DQ2 和 DQ3;

我试过 nls_fit=nlsLM(df[,3:6]~ifelse(.....

nls.lm(par = start, fn = FCT, jac = jac, control = control, lower = lower, : fn 函数的评估返回不合理的值!

现在问题归结为多柱拟合。我怎样才能适应多列DQ0:DQ3?我检查了如何从数据框中简洁地编写包含许多变量的公式?但找不到在我的数据框中使用的解决方案。


此外,当我DQ0对其组内的列进行拟合时,您可以从输出中看到,为所有组生成相同的 Del 和 J 参数a:f

df_new<- df%>%
  group_by(No)%>%
  do(data.frame(model=tidy()))>%
  ungroup

df_new

Source: local data frame [24 x 6]

   No model.term model.estimate model.std.error model.statistic model.p.value
1   a       del1       1.132546        9024.255    1.255002e-04     0.9999000
2   a         J1      23.540764      984311.373    2.391597e-05     0.9999809
3   a       del2       1.079182       27177.895    3.970809e-05     0.9999684
4   a         J2     -20.527520     2362268.839   -8.689748e-06     0.9999931
5   b       del1       1.132546        9024.255    1.255002e-04     0.9999000
6   b         J1      23.540764      984311.373    2.391597e-05     0.9999809
7   b       del2       1.079182       27177.895    3.970809e-05     0.9999684
8   b         J2     -20.527520     2362268.839   -8.689748e-06     0.9999931
9   c       del1       1.132546        9024.255    1.255002e-04     0.9999000
10  c         J1      23.540764      984311.373    2.391597e-05     0.9999809
.. ..        ...            ...             ...             ...           ...
4

0 回答 0