2

我有一个包含多个响应变量和三种处理的数据集。治疗 2 嵌套在治疗 1 中，治疗 3 嵌套在治疗 2 中。为了简单起见，我只显示了三个响应变量。我想运行这个超过 22 个响应变量，其中 3 个显示在演示表中。

我的目标：

  1. 可视化响应变量如何根据治疗组合变化。我创建了一个脚本来对一个响应变量执行此操作。我正在复制粘贴此代码以运行其他列，这对我来说是一种非常粗暴的方法。这导致了我的第二个目标。
  2. 自动化或修改以下脚本，以便它可以自动循环遍历列并生成所需的表格和图形。

演示数据： demo.table

这是我的脚本：

library(doBy)
length2 <- function (x, na.rm=FALSE) {
    if (na.rm) sum(!is.na(x))
    else       length(x)
}
 attach (demo)
cdataNA <- summaryBy(tyr ~ spp + wat + ins, data=demo, FUN=c(length2,mean,sd), na.rm=TRUE)
# Rename column change.length to just N
names(cdataNA)[names(cdataNA)=="tyr.length2"] <- "N"
# Calculate standard error of the mean
cdataNA$tyr.SE <- cdataNA$tyr.sd / sqrt(cdataNA$N)
cdataNA
# Now create a barplot using ggplot2
library(ggplot2)
a <- ggplot(cdataNA, aes(x = wat, y = tyr.mean, fill = ins))
b <- a + geom_bar(stat = "identity", position = "dodge") + facet_grid (~ spp)
# Now put errorbars.
c <- b + geom_errorbar(aes(ymin=tyr.mean-tyr.SE, ymax=tyr.mean+tyr.SE), 
                  width=.2,                    # Width of the error bars
                  position=position_dodge(.9)) + 
xlab ("wat") + 
ylab ("tyr (PA/PA std)")
c

## esc
library(doBy)
length2 <- function (x, na.rm=FALSE) {
    if (na.rm) sum(!is.na(x))
    else       length(x)
}
cdataNA1 <- summaryBy(esc ~ spp + wat + ins, data=demo, FUN=c(length2,mean,sd), na.rm=TRUE)
# Rename column change.length to just N
names(cdataNA1)[names(cdataNA1)=="esc.length2"] <- "N"
# Calculate standard error of the mean
cdataNA1$esc.SE <- cdataNA1$esc.sd / sqrt(cdataNA1$N)
cdataNA1
# Now create a barplot using ggplot2
library(ggplot2)
a1 <- ggplot(cdataNA1, aes(x = wat, y = esc.mean, fill = ins))
b1 <- a1 + geom_bar(stat = "identity", position = "dodge") + facet_grid (~ spp)
# Now put errorbars.
c1 <- b1 + geom_errorbar(aes(ymin=esc.mean-esc.SE, ymax=esc.mean+esc.SE), 
                  width=.2,                    # Width of the error bars
                  position=position_dodge(.9)) + 
xlab ("wat") + 
ylab ("esc (PA/PA std)")
c1

tyr 的结果表：

  spp  wat ins N tyr.mean      tyr.sd      tyr.SE
1  Bl High  No 4 0.305325 0.034102041 0.017051020
2  Bl High Yes 5 0.186140 0.045165894 0.020198802
3  Bl  Low  No 5 0.310540 0.061810096 0.027642315
4  Bl  Low Yes 5 0.202840 0.029034944 0.012984822
5 Man High  No 4 0.122725 0.075867005 0.037933503
6 Man High Yes 5 0.081800 0.013463469 0.006021046
7 Man  Low  No 5 0.079880 0.009569587 0.004279650
8 Man  Low Yes 4 0.083550 0.018431947 0.009215973

esc 的结果图：esc 的 演示图

所以整个事情都有效，但仍然需要大量的体力劳动，这会阻碍工作流程。实现自动化会很棒。

提前致谢。

4

1 回答 1

0

您可以只用两行来组织数据：

melt.dta <- melt(dta, id.vars = c("spp", "wat", "ins"), measure.vars = "tyr")
cast(melt.dta, spp + wat + ins ~ ., 
     function (x) c("N" = sum(!is.na(x)), 
                    "mean" = mean(x, na.rm = TRUE),                   
                    "sd" = sd(x, na.rm = TRUE),
                    "se" = sd(x, na.rm = TRUE)/sqrt(sum(!is.na(x)))))

它返回：

  spp  wat ins N   mean      sd      se
1  Bl High  No 4 0.3053 0.03410 0.01705
2  Bl High Yes 5 0.1861 0.04517 0.02020
3  Bl  Low  No 5 0.3105 0.06181 0.02764
4  Bl  Low Yes 5 0.2028 0.02903 0.01298
5 Man High  No 4 0.1227 0.07587 0.03793
6 Man High Yes 5 0.0818 0.01346 0.00602
7 Man  Low  No 5 0.0799 0.00957 0.00428
8 Man  Low Yes 4 0.0835 0.01843 0.00922
于 2012-10-20T02:20:03.233 回答