2

我有以下数据集（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

对于组的每个值，我想将下面提到的函数应用于 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
4

1 回答 1

3

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

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

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函数末尾的 a 中：

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

现在你会得到一个data.frame包含所有这些的，三列代表ec20,ec50和ec80

对于 optim 的问题：我认为问题在于

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

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

par[3]<-x[[which.min(abs(y-50))]]
于 2011-08-29T23:02:42.013 回答