将函数重复应用于多个主题

Question

我有一个数据框如下，

> mydata
date  station  treatment  subject   par
A       a         0         R1      1.3    
A       a         0         R1      1.4    
A       a         1         R2      1.4   
A       a         1         R2      1.1    
A       b         0         R1      1.5    
A       b         0         R1      1.8     
A       b         1         R2      2.5     
A       b         1         R2      9.5    
B       a         0         R1      0.3    
B       a         0         R1      8.2    
B       a         1         R2      7.3    
B       a         1         R2      0.2    
B       b         0         R1      9.4    
B       b         0         R1      3.2    
B       b         1         R2      3.5    
B       b         1         R2      2.4 
....

地点：

date 是一个具有 2 个级别 A/B 的因子； station 是具有 2 个级别 a/b 的因子； treatment 是具有 2 个级别 0/1 的因子；

subject 是分配给处理的重复 R1 到 R20（10 到 treatment 0 和 10 到处理 1）；

和 par是我的参数，是在每个日期和站点重复测量每个受试者的粒径

我需要做的是： 将标准分成 10 个相等的箱子，并计算每个箱子中的数量。这必须在由日期站和主题组合定义的mydata 的子集中完成。最终结果必须是 daframe myres，如下所示：

> myres
    date  station  treatment  bin.centre  freq
    A       a         0         1.2        4 
    A       a         0         1.3        3    
    A       a         0         1.4        2 
    A       a         0         1.5        1    
    A       a         1         1.2        4    
    A       a         1         1.3        3    
    A       a         1         1.4        2     
    A       a         1         1.5        1    
    B       b         0         2.3        5   
    B       b         0         2.4        4    
    B       b         0         2.5        3    
    B       b         0         2.6        2   
    B       b         1         2.3        5   
    B       b         1         2.4        4   
    B       b         1         2.5        3   
    B       b         1         2.6        2
    ....

这是我到目前为止所做的：

#define the number of bins
num.bins<-10

#define the width of each bins
bin.width<-(max(par)-min(par))/num.bins

#define the lower and upper boundaries of each bins
bins<-seq(from=min(par), to=max(par), by=bin.width)

#define the centre of each bins
bin.centre<-c(seq(min(bins)+bin.width/2,max(bins)-bin.width/2,by=bin.width))

#create a vector to store the frequency in each bins

  freq<-numeric(length(length(bins-1)))

 # this is the loop that counts the frequency of particles between the lower and upper boundaries
 of each bins and store the result in freq

 for(i in 1:10){
    freq[i]<-length(which(par>=bins[i] &
    par<bins[i+1]))
     }

 #create the data frame with the results
 res<-data.frame(bin.centre,res)

我的第一种方法是手动对 mydata 进行子集化，使用 subset()，针对每个主题站和日期的组合，并对每个子集应用上述命令序列，然后使用 @ 组合每个 res 构建最终数据帧987654334@，但这个过程非常复杂，容易传播错误。 我想做的是自动化上述过程，以便计算每个主题的分箱频率分布。我的直觉是，最好的方法是创建一个函数来估计这个粒子分布，然后通过 for 循环将其应用于每个对象。但是，我不确定如何执行此操作。任何建议将不胜感激。

谢谢 马特奥。

Answer 1

您可以使用plyr 包中的功能通过几个步骤完成此操作。这使您可以将数据拆分为所需的块，对每个块应用统计信息，然后组合结果。

首先我设置了一些虚拟数据：

set.seed(1)
n <- 100
dat <- data.frame(
    date=sample(LETTERS[1:2], n, replace=TRUE),
    station=sample(letters[1:2], n, replace=TRUE),
    treatment=sample(0:1, n, replace=TRUE),
    subject=paste("R", sample(1:2, n, replace=TRUE), sep=""),
    par=runif(n, 0, 5)
)
head(dat)

  date station treatment subject       par
1    A       b         0      R2 3.2943880
2    A       a         0      R1 0.9253498
3    B       a         1      R1 4.7718907
4    B       b         0      R1 4.4892425
5    A       b         0      R1 4.7184853
6    B       a         1      R2 3.6184538

现在我使用 base 中名为 cut 的函数将您的面值分成大小相等的箱子：

dat$bin <- cut(dat$par, breaks=10)

现在是有趣的部分。加载包plyr 并使用函数ddply 进行拆分、应用和合并。因为您需要频率计数，我们可以使用函数length 来计算每个重复出现在该 bin 中的次数：

library(plyr)
res <- ddply(dat, .(date, station, treatment, bin), 
  summarise, freq=length(treatment))
head(res)

  date station treatment             bin freq
1    A       a         0 (0.00422,0.501]    1
2    A       a         0   (0.501,0.998]    2
3    A       a         0      (1.5,1.99]    4
4    A       a         0     (1.99,2.49]    2
5    A       a         0     (2.49,2.99]    2
6    A       a         0     (2.99,3.48]    1