R, plyr, 具有复函数

Question

我有以下数据集(CEU):

group  x      y
1     -23     100
1     -0.90   69.62
1     -0.90   72.03
2     -23     100
2      0.69   48.01
2      0.69   45.63

对于 group 的每个值，我想将下面提到的函数应用于 x 和 y 值的每个子集。然后我想组合所有结果并将它们写入表格中以导出。

我不确定如何准确地应用 plyr 函数来做到这一点......如果这确实是正确的行动方案。

x<-c(-23.0000,-0.9031,-0.9031)
y<-c(100,85.72,86.65)

par<-c(16.88,100.28,-.75,4.129)

dcrit<-function(d) { 
    sumsq<-0
    for (i in 1:length(x)){
      sumsq<-sumsq+ (y[i]-(par[1]+(par[2]-par[1])/(1+10^((x[i]-par[3])*d))))^2      
    }
    sumsq
}

S<-function(par) { 
    a<-par[1]
    b<-par[2]
    c<-par[3]
    d<-par[4]
    sumsq<-0
    for (i in 1:length(x)){
      sumsq<-sumsq+ (y[i]-(a+(b-a)/(1+10^((x[i]-c)*d))))^2      
    }
    sumsq
}
optim(par,S)

CEU <- read.csv(file="C:/files/CEU.csv",head=TRUE,sep=",")
CEU

data <- ddply(CEU,.(group),function(xy) 
{
par[1]<-min(y)
par[2]<-100
par[3]<-x[[which.min(abs(y-50))]]
par[4]<-optimize(dcrit,interval=c(-100,100))$minimum

o<-optim(par,S)
par<-o$par

a<-par[1];
b<-par[2];
c<-par[3];
d<-par[4];

k<-(b-a)/(20-a)-1
if (k>0) ec20<-c+1/d*log10(k) else ec20<-NA
ec20

z<-(b-a)/(50-a)-1
 if (z>0) ec50<-c+1/d*log10(z) else ec50<-NA
ec50

j<-(b-a)/(80-a)-1
if (j>0) ec80<-c+1/d*log10(j) else ec80<-NA
ec80

data.frame(ec20, ec50, ec80)

})

data

代码运行没有错误，但仅在由以下设置的原始 x 和 y 值上运行:

 x<-c(-23.0000,-0.9031,-0.9031)
 y<-c(100,85.72,86.65)

ddply 不使用数据集 CEU 中的 x 和 y 值。它们不像处理组值那样以迭代方式替换原始 x 和 y。数据具有适当的组数，并且 ec20/ec50/ec80 值是准确的，但仅适用于原始 x 和 y。

> data
   group       ec20       ec50       ec80
1      1 -0.3652977 -0.6843279 -0.8530892
2      2 -0.3652977 -0.6843279 -0.8530892
3      3 -0.3652977 -0.6843279 -0.8530892
4      4 -0.3652977 -0.6843279 -0.8530892
5      5 -0.3652977 -0.6843279 -0.8530892

Answer 1

看起来你做对了，你只需要产生输出。

我猜这是你的输出在哪里？

k<-(b-a)/(20-a)-1
if (k>0) ec20<-c+1/d*log10(k) else ec20<-NA
ec20

z<-(b-a)/(50-a)-1
 if (z>0) ec50<-c+1/d*log10(z) else ec50<-NA
ec50

j<-(b-a)/(80-a)-1
if (j>0) ec80<-c+1/d*log10(j) else ec80<-NA
ec80

将它们放入 data.frame在函数结束时:

    ...
    data.frame(ec20, ec50, ec80)
}

现在你会得到一个 data.frame与所有这些，三列为 ec20 , ec50和 ec80
对于您的 optim 问题:我认为问题在于

par[3]<-x[which.min(abs(y-50))]

单例 [在 R 中不是常规下标——它得到一个切片——在这种情况下 data.frame列。那条线正在转弯 par从数字向量到 list .添加更多括号:

par[3]<-x[[which.min(abs(y-50))]]