在 R 中为研究中的几个人运行脚本

Question

我是动物行为专业的硕士生，R 语言非常初学者。我正在为我的论文做统计数据，并决定在我求助于 PSPP 之前尝试一下 R，在我上过关于语言的小型入门课程之后，使用 R Studio。

我写了一个非常基本的脚本

1. 从 csv 导入来自 51 个人的数据。文件，并导入一个 csv。用于存储结果的文件。
2. 使用来自一个人的数据创建一个子集，然后根据数据做出假设，运行二项式检验并确定假设是否被拒绝。假设、二项式检验的 p 值和每个个体对假设的拒绝/接受都存储在结果文件中。

当我为每个人运行该脚本时，它是完全正常的。对于代码质量不佳，我提前道歉，但如果您想看一下，这是一个人的脚本：

## Import data and results file
odor_behavior <- read.csv("odor_behavior.csv")
percalf_binomial_odor_results <- read.csv("percalf_binomial_odor_results.csv")

# Create binomial dataset per calf (example here: 958)
pref958 <- subset(odor_behavior, Calf_ID == 958, c("Preference"), drop=FALSE)

# Count the number of testing sessions for this calf
n958 <- nrow(pref958)

# Calculate how many times a suckling behavior appears (Control=1 or Modified=3).
mat958 <- as.matrix(pref958)
fmat958 <- factor(mat958,levels=c(1,3),ordered=TRUE)
fpref958 <- table(fmat958)
data958 <- as.data.frame(fpref958)

# Hypothesize a preference according to the observations and save it. 
# Also, save variables for the binomial test.
## Find the row in which to store the results
row958<- which(percalf_binomial_odor_results$Calf_ID == 958)
## Decide the hypothesis. 
  if (data958[1,2] > data958[2,2]) {
    x958 <- data958[1,2]
    percalf_binomial_odor_results[row958,c("H_preference")] <- "Control"
  } else if (data958[1,2] < data958[2,2]) {
    x958 <- data958[2,2]
    percalf_binomial_odor_results[row958,c("H_preference")] <- "Odor"
  } else if (data958[1,2] == data958[2,2]) {
    x958 <- data958[1,2]
    percalf_binomial_odor_results[row958,c("H_preference")] <- "Either"
  } else {print("Check your data, Maria!")
  }


# Run binomial test
binom.test (x=x958, n=n958, p=0.5, conf.level = 0.95)

# Save the p.value and decide whether the hypothesized preference is real or not
pbin958 <- binom.test (x=x958, n=n958, p=0.5, conf.level = 0.95)$p.value
percalf_binomial_odor_results[row958,c("Bin_odor_pvalue")] <- pbin958
  if (pbin958 <0.05) {
    percalf_binomial_odor_results[row958,c("Real_preference")] <- "Yes"
  } else {
    percalf_binomial_odor_results[row958,c("Real_preference")] <- "No"
  }

如前所述，当我通过复制粘贴代码并手动更改“Calf_ID”五十次来为每个人重复脚本时，一切运行顺利。

问题：脚本是否可以为每个人单独运行？我尝试通过将所有 51 个“Calf_ID”放在一个矩阵中，并将脚本放在一个函数中，该函数使用矩阵中的元素来遍历每个人。

Lapply 返回正确的结果，但仅限于假设是否被接受/拒绝（即脚本的最后一部分），这让我认为该脚本确实适用于所有 51 个人。但是，在此过程中，没有结果存储在相应的文件中（很可能是因为我首先做错了什么）。当我执行手动复制粘贴/编辑过程时，所有结果（假设、p 值和每个人的假设结果）都正确存储在结果文件中。

我可以通过复制粘贴的方法得到我的结果；这只是效率问题，因为现在我有 51 个人，但是当他们有 1000 个人时，一个人会做什么？很高兴知道如何让它发挥作用！

编辑
根据 cmets 的要求，以下是结果（在运行个人 958 的脚本之后）到

dput(head(odor_behavior))：
结构（列表（Calf_ID = c(958L, 958L, 958L, 958L, 958L, 958L), Preference = c(1L, 1L, 3L, 1L, 1L, 3L)), .Names = c("Calf_ID", "Preference "), row.names = c(NA, 6L), class= "data.frame")

和

dput(head(percalf_binomial_odor_results))：
结构（列表（Calf_ID = c（958L，7015L，7017L，959L，7018L，7019L）， H_preference = c("控制", NA, NA, NA, NA, NA), Bin_odor_pvalue = c(0.453125, NA, NA, NA, NA, NA), Real_preference = c("否", NA, NA, NA, NA, NA)), .Names = c("Calf_ID", "H_preference", "Bin_odor_pvalue", "Real_preference"), row.names = c(NA, 6L), class= "data.frame")

更新
这是同一主题的另一个示例，基于为初始脚本提供的答案@alexis_laz。在这种情况下，我正在导入一个包含 3 列（Calf_ID、Control_time、Odor_time）的数据集，其中包括给定小牛 ID 的多个控制和气味时间测量值（如初始示例中的首选项）。该脚本尝试对每个 Calf_ID 的 Control_time - Odor_time 对执行 Wilcoxon 符号秩检验。以下是脚本的外观：

    ### Import data and load library with test
odor_times <- read.csv("odor_times.csv")
library(coin)

### Create a function that examines suckling preference
prefanalysis <- function(Control_time, Odor_time)
{
  sessions <- length(Control_time)

  #Run a Wilcoxon signed rank test. Set y ~ x as control and odor, 
  mypvalue <- pvalue(wilcoxsign_test(Control_time ~ Odor_time))
  myzvalue <- statistic(wilcoxsign_test(Control_time ~ Odor_time))

  #Decide if there is a true difference in suckling times according to p value, and save it
  true_diff <- if(mypvalue < 0.05) "Yes" else "No"

  # Create a data frame with hypothesis, p value and real preference
  data.frame(sessions, myzvalue, mypvalue, true_diff, stringsAsFactors = FALSE)
}

### Run the function for all calves, with preference depending on calf ID. Save results.
results <- as.data.frame(do.call(data.frame, aggregate(cbind(Control_time, Odor_time) ~ Calf_ID, odor_times, prefanalysis)))

在这种情况下，我收到以下错误：

     Error in eval(expr, envir, enclos) : 
  argument "Odor_time" is missing, with no default

当使用 Debug 重新运行时，我可以看到脚本开始正确读取第一个个体的 Control_time（即具有特定 Calf_ID 的行）。然而，对于 Odor_time，它似乎无法做到这一点。我无法（终生）看到我做错了什么。

Answer 1

您似乎将计算应用于按每个“Calf_ID”分组的“偏好”。我们可以利用 R 的 group-by 函数，而不是重复地为每个 id 设置子集 odor_behavior$Preference，首先，构建一个函数，在给定适当输入的情况下执行所有工作。

我试图在一个函数中简化你的脚本——希望我没有错过很难容纳的细节：

ff = function(pref)
{
    fx = factor(pref, levels = c(1L, 3L), labels = c("Control", "Odor"))
    tab = table(fx)

    p = binom.test(max(tab), length(pref), p = 0.5, conf.level = 0.95)$p.value

    real_pref = if(p < 0.05) "yes" else "no"
    H_pref = if(tab["Control"] == tab["Odor"]) "Either" else names(which.max(tab))

    data.frame(p = p, real_pref = real_pref, H_pref = H_pref, stringsAsFactors = FALSE)
}

然后，将其应用于所有“Calf_ID”个人：

do.call(data.frame, aggregate(pref ~ ID, dat, ff))
#   ID     pref.p pref.real_pref pref.H_pref
#1   1 0.07295139             no     Control
#2   2  0.1689778             no     Control
#3   3  0.8919232             no        Odor
#4   4  0.6029232             no        Odor
#5   5          1             no     Control
#6   6  0.7659918             no     Control
#7   7          1             no        Odor
#8   8   0.889884             no     Control
#9   9          1             no      Either
#10 10  0.5758493             no     Control

dat 在哪里：

set.seed(1821)
dat = data.frame(ID = sample(10, 500, TRUE), pref = sample(c(1L, 3L), 500, TRUE))