尝试让脚本为每 24 行计算一个值(使用函数)

Trying to make a script calculate a value (using a function) for every 24 rows

我在 Whosebug 上找不到类似问题的解决方案。希望有人能帮忙!

我用的是R环境

我有龟巢的数据。每个巢中有两种类型的小时数据。第一个是每小时温度,它有一个相关的每小时发育("anatomical" 胚胎发育量”)。

我正在计算加权 中位数 。在这种情况下,中位数是温度,它是按发展加权的。

我这里有一个脚本,我用它来计算加权中位数:

weighted.median <- function(x, w, probs=0.5, na.rm=TRUE) {
  x <- as.numeric(as.vector(x))
  w <- as.numeric(as.vector(w))
  if(anyNA(x) || anyNA(w)) {
    ok <- !(is.na(x) | is.na(w))
        x <- x[ok]
    w <- w[ok]
  }
  stopifnot(all(w >= 0))
  if(all(w == 0)) stop("All weights are zero", call.=FALSE)
  #'
  oo <- order(x)
  x <- x[oo]
  w <- w[oo]
  Fx <- cumsum(w)/sum(w)
  #'
  result <- numeric(length(probs))
  for(i in seq_along(result)) {
    p <- probs[i]
    lefties <- which(Fx <= p)
    if(length(lefties) == 0) {
      result[i] <- x[1]
    } else {
      left <- max(lefties)
      result[i] <- x[left]
      if(Fx[left] < p && left < length(x)) {
        right <- left+1
        y <- x[left] + (x[right]-x[left]) * (p-Fx[left])/(Fx[right]-        Fx[left])
        if(is.finite(y)) result[i] <- y
      }
    }
  }
  names(result) <- paste0(format(100 * probs, trim = TRUE), "%")
  return(result)
}

所以从函数中你可以看到我需要两个输入向量,x 和 w(分别是温度和发展)。

我遇到的问题是我的每小时温度曲线持续 5 天到 53 天(即 120 小时到 1272 小时)。

我想计算一个巢内所有日期的每日加权中位数(即,取 x 和 w 的 24 行,计算加权中位数,然后移动到第 25-48 行,依此类推。 ) 因此,输出向量将是长度为 n/24 的每日加权中位数列表(其中 n 是 x 中的总行数)。

换句话说,我想自动分析我的数据,其方式等同于手动执行此操作(nest1 是 Nest 1 的数据表,其中包含两个向量,temp 和 devo(devo 是权重)) :

`weighted.median(nest1$temp[c(1,1:24)],nest1$devo[c(1,1:24)],na.rm=TRUE)`

其次是

weighted.median(nest1$temp[c(1,25:48)],nest1$devo[c(1,25:48)],na.rm=TRUE)

其次是

weighted.median(nest1$temp[c(1,49:72)],nest1$devo[c(1,49:72)],na.rm=TRUE)

一直到

`weighted.median(nest1$temp[c(1,n-23:n)],nest1$devo[c(1,n-23:n)],na.rm=TRUE)`

恐怕我什至不知道从哪里开始。非常感谢任何帮助或线索。

主要思想是为第 1 天、第 2 天、...、第 n/24 天创建一个新列,按天将数据框拆分为子集,并将您的函数应用于每个子集。

首先我创建了一些示例数据:

set.seed(123)
n <- 120 # number of rows
nest1 <- data.frame(temp = rnorm(n), devo = rpois(n, 5))

创建拆分变量:

nest1$day <- rep(1:(nrow(nest1)/24), each = 24)

然后,使用by()函数将nest1除以nest1$day并将该函数应用于每个子集:

out <- by(nest1, nest1$day, function(d) {
  weighted.median(d$temp, d$devo, na.rm = TRUE)
})
data.frame(day = dimnames(out)[[1]], x = as.vector(out))
#   day           x
# 1   1 -0.45244433
# 2   2  0.15337312
# 3   3  0.07071673
# 4   4  0.23873174
# 5   5 -0.27694709

除了使用 by,您还可以使用 dplyr 包中的 group_by + summarise 函数:

library(dplyr)
nest1 %>%
  group_by(day) %>%
  summarise(x = weighted.median(temp, devo, na.rm = TRUE))
# # A tibble: 5 x 2
#     day       x
#   <int>   <dbl>
# 1     1 -0.452 
# 2     2  0.153 
# 3     3  0.0707
# 4     4  0.239 
# 5     5 -0.277