在 R 中计算数据表中每行的几何平均值

Question

我有一个数据表，想要在几列中为每一行计算几何平均值。其中某些数值为零，因此我需要将其排除。

维基百科上的几何平均数定义为：“几何平均数被定义为n个数字的乘积的n次方根”，因此对于2个数字，它们的几何平均数就是它们的乘积的平方根。

在我的情况下，每行的n次方根将根据其中非零值的数量而变化。

在下面的示例中，结果列的前两行计算如下：

1：(a * c)^(1/2)

2：(a * b * c)^(1/3)

因此，我需要让公式查看a:c列，取非零值的乘积，然后取出非零值的数量的n次方根。

library(data.table)
dt <- data.table(a = c(0.5, 0.3,0,0.6), b = c(0,0.4,0.1,0), 
c = c(0.9,0.5,0.1,0), Result = c(0.67, 0.39, 0.1, 0.6))

- MidnightDataGeek

---

“Result”是预期的输出吗？ - akrun

所有的数字都是非负数吗？ - Sandipan Dey

抱歉，结果是期望的输出，所有值都是非负数。 - MidnightDataGeek

Answer 1

我们可以尝试使用 data.table 方法。

dt[, v1 := Reduce(`+`, lapply(.SD, function(x) x!=0)), .SDcols = 1:3]
dt[, result2 := round((Reduce(`*`, lapply(.SD, function(x) 
    replace(x, x==0, 1))))^(1/v1), 2), .SDcols = 1:3][, v1 := NULL][]
#    a   b   c Result result2
#1: 0.5 0.0 0.9   0.67    0.67
#2: 0.3 0.4 0.5   0.39    0.39
#3: 0.0 0.1 0.1   0.10    0.10
#4: 0.6 0.0 0.0   0.60    0.60

另一个不太高效的选项是按行序分组，然后在每一行上执行操作。

dt[, result2 := {
           u1 <- unlist(.SD)
           round(prod(u1[u1!=0])^(1/sum(u1!=0)), 2)} , 1:nrow(dt), .SDcols = 1:3]
dt
#     a   b   c Result result2
#1: 0.5 0.0 0.9   0.67    0.67
#2: 0.3 0.4 0.5   0.39    0.39
#3: 0.0 0.1 0.1   0.10    0.10
#4: 0.6 0.0 0.0   0.60    0.60

注意：这两种方法都是data.table的方法。

或者另一种选项是由@DavidArenburg贡献的。

dt[, Result := round(Reduce(`*`, replace(.SD, .SD == 0, 1))^(1/rowSums(.SD != 0)), 2)]

另一个向量化的选择是转换为矩阵。

library(matrixStats)
m1 <- as.matrix(setDF(dt)[1:3])
round(rowProds(replace(m1, !m1, 1))^(1/rowSums(m1!=0)), 2)
#[1] 0.67 0.39 0.10 0.60

- akrun

---

代码块 dt[, result2....' 返回错误...在 eval(expr, envir, enclos) 中出错：找不到对象 'v1'。 - MidnightDataGeek

我正在使用1.9.6版本。我已经修改了下面的答案，严格来说它还是“data.table”方法吗？对我来说速度很关键，所以我尽量使用DT完成所有操作。 - MidnightDataGeek

@MidnightDataGeek 好的，这就是原因。在1.10.0版本中，使用.SDcols时也可以访问其他列。是的，前两个是data.table方法。但是，如果您正在使用apply，它可能会转换为matrix，并且使用data.table没有任何好处。 - akrun

也许更简单的表达是：dt[, Result := round(Reduce(`*`, replace(.SD, .SD == 0, 1))^(1/rowSums(.SD != 0)), 2)]。 - David Arenburg

@DavidArenburg 谢谢。你的解决方案快速易读。 真是太神奇了，你们能如此快速地想出解决方案！ - MidnightDataGeek

Answer 2

prod(a)^(1/length(a)) 给出向量 a 的几何平均值

- Ajay Ohri

---

如此陈述，被点赞的Geomean公式是错误的！它应该是1/n次幂。不应该有“sum”。我很惊讶之前没有人点赞这个评论！ - RichardBJ

Answer 3

这也适用，假设所有值均为非负数。

dt$Result <- apply(dt, 1, function(x) (prod(x[x!=0]))^(1/sum(x!=0)))
dt
#     a   b   c    Result
#1: 0.5 0.0 0.9 0.6708204
#2: 0.3 0.4 0.5 0.3914868
#3: 0.0 0.1 0.1 0.1000000
#4: 0.6 0.0 0.0 0.6000000

- Sandipan Dey

---

谢谢。我根据我的需求进行了修改... dt[, Result := apply(.SD, 1, function(x) (prod(x[x != 0])) ^ (1/sum(x!=0))), .SDcols = 1:3] - MidnightDataGeek

Answer 4

其他选项：

m1 <- as.matrix(setDF(dt)[1:3])
exp(rowMeans(log(m1)))

- Tom Wenseleers