1

我正在使用一个for循环来创建一系列lattice直方图，比较来自样本数据的长度频率分布，并希望每个直方图显示一个垂直线，abline显示每个分布的平均值并text指示样本数。

我有以下数据：

> head(hist.data,20)
               Scientific_name Count Length..cm. Method
3  Pristipomoides filamentosus     1          60 BotCam
5           Etelis carbunculus     1          43 BotCam
6             Etelis coruscans     1          40 BotCam
12 Pristipomoides filamentosus     1          55 BotCam
16           Aphareus rutilans     1          67 BotCam
17           Aphareus rutilans     1          77 BotCam
20          Etelis carbunculus     1          46 BotCam
21    Pristipomoides sieboldii     1          35 BotCam
23    Pristipomoides sieboldii     1          33 Fishing
25          Etelis carbunculus     1          53 Fishing
26 Pristipomoides filamentosus     1          45 Fishing
27 Pristipomoides filamentosus     1          43 Fishing
28 Pristipomoides filamentosus     1          58 Fishing
29 Pristipomoides filamentosus     1          55 Fishing
30    Pristipomoides sieboldii     1          29 Fishing

我的代码如下：

#create a list of species
sp <- c("Etelis coruscans","Etelis carbunculus","Pristipomoides sieboldii","Pristipomoides filamentosus","Pristipomoides zonatus","Epinephelus quernus","Aphareus rutilans")

#Calculate sample# and mean length by species by method
n <- with(hist.data, tapply(Scientific_name, Method, function(x) count(x)))
mean.length <- aggregate(Length..cm. ~ Scientific_name + Method, data = hist.data, FUN= "mean")

#plot hisotgrams for each spp in 1cm bins
for (i in sp){
    BIN_WIDTH <- 1 #desired bin width
    print(histogram(~ Length..cm. | Method, #create and print histogram
    data = hist.data[hist.data$Scientific_name == i,], 
    nint = (max(hist.data$Length..cm.) - min(hist.data$Length..cm.)+1)/BIN_WIDTH,
    layout = c(1,2),
    type = "density",
    main = substitute(expr = expression(paste("Length-Frequency of ", italic(i), " by Gear")), env = list(i=i)),
    xlab = "Length (cm)",
    panel = function(x, ...){
        #panel.abline(v = 60, col = "red", lty = 2)
        #panel.text(lab = paste("Sample #: ",n$BotCam[1,2]), 90, 100)
        panel.histogram(x,...)
        panel.mathdensity(dmath = dnorm, col = "black",
                      args = list(mean = mean(x), sd= sd(x)), ...)     
    }
    ))


    #save histogram as PDF file
    quartz.save(paste("Length-Frequency of", i, "by method.pdf", sep = " "), type = "pdf")
    dev.off() #close the graphics diver
}

我可以生成以下数组：

n <- with(hist.data, tapply(Scientific_name, Method, function(x) count(x)))
n
$BotCam
                            x freq
1           Aphareus rutilans   16
2          Etelis carbunculus   35
3            Etelis coruscans   20
4 Pristipomoides filamentosus  179
5    Pristipomoides sieboldii  125
6      Pristipomoides zonatus    2

$Fishing
                            x freq
1         Epinephelus quernus    2
2          Etelis carbunculus   68
3            Etelis coruscans   30
4 Pristipomoides filamentosus   24
5    Pristipomoides sieboldii   80
6      Pristipomoides zonatus    5

和

mean.length <- aggregate(Length..cm. ~ Scientific_name + Method, data = hist.data, FUN= "mean")
> mean.length
               Scientific_name  Method Length..cm.
1            Aphareus rutilans  BotCam    58.81250
2           Etelis carbunculus  BotCam    43.65714
3             Etelis coruscans  BotCam    46.55000
4  Pristipomoides filamentosus  BotCam    53.22346
5     Pristipomoides sieboldii  BotCam    35.52000
6       Pristipomoides zonatus  BotCam    35.00000
7          Epinephelus quernus Fishing    74.00000
8           Etelis carbunculus Fishing    42.98529
9             Etelis coruscans Fishing    49.96667
10 Pristipomoides filamentosus Fishing    59.58333
11    Pristipomoides sieboldii Fishing    37.25000
12      Pristipomoides zonatus Fishing    30.80000

我想更换：

#panel.abline(v = 60, col = "red", lty = 2)
#panel.text(lab = paste("Sample #: ",n$BotCam[1,2]), 90, 100)

使用abline基于物种分布的平均值text生成的代码，以及生成ablinefrommean.length和“Sample #:” from的代码n。上图和下图以及loop. 有没有办法做到这一点？我不喜欢lattice，这只是我有一些经验。如果类似的东西ggplot对此更好，请告诉我。

4

2 回答 2

1

在四处寻找之后，我想出了以下似乎可以解决问题的代码。现在我只需要找出panel.text放置输出的最佳方法，因为我的轴会随着循环的每次迭代而变化。

#plot hisotgrams for each spp in 1cm bins
for (i in sp){
    BIN_WIDTH <- 1 #desired bin width
    print(histogram(~ Length..cm. | Method, #create and print histogram
    data = hist.data[hist.data$Scientific_name == i,], 
    nint = (max(hist.data$Length..cm.) - min(hist.data$Length..cm.)+1)/BIN_WIDTH,
    layout = c(1,2),
    type = "density",
    main = substitute(expr = expression(paste("Length-Frequency of ", italic(i), " by Gear")), env = list(i=i)),
    xlab = "Length (cm)",
    panel = function(x, ...){
        mean.values <- mean(x)
        panel.abline(v=mean.values, col.line="red", lty = 2)
        sample.n <- length(x)
        panel.text(lab = paste("Sample size = ", sample.n), mean.values*1.2, .26)
        panel.text(lab = paste("Mean = ", round(mean.values, 1)), mean.values*1.2, .23)
        panel.histogram(x,...)
        panel.mathdensity(dmath = dnorm, col = "black",
                      args = list(mean = mean(x), sd= sd(x)), ...)     
    }
    ))

    #save histogram to PDF file
    quartz.save(paste("Length-Frequency of", i, "by method.png", sep = " "), type = "png")
    dev.off() #close graphics diver
}
于 2013-05-25T00:09:01.610 回答
0

你可以上传你的代码的副本吗？听起来您没有正确定义（x）。

于 2013-07-10T03:58:23.643 回答