使用线性回归计算数组中的趋势

Question

我有一个维度数组c(54,71,360)，其中包含气候数据。前两个维度描述了区域的网格，而第三个维度是时间维度。所以在这种情况下，有 360 个时间步（月）。

这是生成示例数组的代码：

set.seed(5)
my_array <- array(sample(rnorm(100), 600, replace=T), dim= c(54,71,360))

现在我想计算每个网格单元的趋势。趋势等于线性回归方程的斜率。这就是为什么需要计算每个网格单元随时间的线性回归的原因。这正是我正在努力解决的问题。

为了清楚地展示我想要做什么，这里有一个网格单元的示例，它是从数组中获取的，作为长度为 360 的向量：

grid_cell <- my_array[1,1,]

需要计算这个向量随时间的线性回归。为此，我们创建了一个简单的时间向量：

time_vec <- 1:360

由于我只对斜率系数感兴趣，所以可以这样：

trend <- lm(grid_cell ~ time_vec)$coefficients[2]

在这种情况下，这导致1.347029e-05 的值。

我想对数组的每个网格单元执行此操作，以便输出是维度为c(54,71) 的矩阵，这意味着每个网格单元都有一个趋势值。

我尝试了以下方法，但没有成功：

trend_mat <- apply(my_array, 1:2, lm(my_array ~ time_vec)$coefficients[2])

我收到错误消息：Error in model.frame.default: variable lengths differ.

这有点令人惊讶，因为数组的第三维和time_vec 的长度都是 360。

有人知道如何实现这一目标吗？

当然，我也愿意接受其他可能完全不同的解决方案，只要它们导致相同的结果。

Answer 1

问题中代码的问题在于

• apply 的第三个参数应该是一个函数，并且问题的代码提供了一个表达式而不是一个函数。
• 它多次应用lm。我们展示了如何只应用一次lm，而在第二种选择中我们根本不使用lm。如下面的“性能”部分所示，这提供了一个和两个数量级的加速。

如果我们使用更小的数据，如最后的注释所示，则更容易说明。要在您的示例中使用它，只需将 dims 替换为注释中注释掉的行中显示的行。

1) 首先我们将数组reshape成矩阵，执行lm，然后reshape回来。这会调用 lm 一次，而不是调用它 prod(dims[1:2]) 次。

m <- t(matrix(a,,dim(a)[3]))
array(coef(lm(m ~ timevec))[2, ], dim(a)[1:2])
##            [,1]      [,2]      [,3]
## [1,]  0.2636792 0.5682025 -0.255538
## [2,] -0.4453307 0.2338086  0.254682

# check
coef(lm(a[1,1,] ~ timevec))[[2]]
## [1] 0.2636792
coef(lm(a[2,1,] ~ timevec))[[2]]
## [1] -0.4453307
coef(lm(a[1,2,] ~ timevec))[[2]]
## [1] 0.5682025
coef(lm(a[2,2,] ~ timevec))[[2]]
## [1] 0.2338086
coef(lm(a[1,3,] ~ timevec))[[2]]
## [1] -0.255538
coef(lm(a[2,3,] ~ timevec))[[2]]
## [1] 0.254682

2) 或者，我们可以使用斜率系数公式完全删除lm，如下所示：

m <- t(matrix(a,,dim(a)[3]))
array(cov(m, timevec) / var(timevec), dims[1:2])
##            [,1]      [,2]      [,3]
## [1,]  0.2636792 0.5682025 -0.255538
## [2,] -0.4453307 0.2338086  0.254682

性能

我们看到单个 lm 的运行速度比 apply 快约 8 倍，消除 lm 的运行速度比 apply 快约 230 倍。因为apply 在我的笔记本电脑上速度非常慢，所以我只使用了 3 次复制，但如果你有更快的机器或更有耐心，你可以增加它。不过，主要结论不太可能发生太大变化。

library(microbenchmark)

set.seed(5)

dims <- c(54,71,360)
a <- array(rnorm(prod(dims)), dims)
timevec <- seq_len(dim(a)[3])


microbenchmark(times = 3L,
  apply = apply(a, 1:2, function(x) coef(lm(x ~ timevec))[2]),
  lm = {  m <- t(matrix(a,,dim(a)[3]))
          array(coef(lm(m ~ timevec))[2, ], dim(a)[1:2]) 
  },
  cov =  {  m <- t(matrix(a,,dim(a)[3]))
            array(cov(m, timevec) / var(timevec), dims[1:2])
  })

给予：

Unit: milliseconds
  expr        min         lq        mean     median         uq        max neval cld
 apply 13446.7953 13523.6016 13605.25037 13600.4079 13684.4779 13768.5479     3   b
    lm   264.5883   275.7611   476.82077   286.9338   582.9370   878.9402     3  a 
   cov    56.9120    57.8830    58.71573    58.8540    59.6176    60.3812     3  a 

注意

测试数据。

set.seed(5)

# dims <- c(54,71,360)
dims <- 2:4
a <- array(rnorm(prod(dims)), dims)
timevec <- seq_len(dim(a)[3])

Answer 2

问题的回归代码中缺少匿名函数。在这里，我将使用 R 4.1.0 中引入的新 lambda。
我也使用推荐的提取器coef。

set.seed(5)
my_array <- array(sample(rnorm(100), 600, replace=T), dim= c(54,71,360))

time_vec <- 1:360
trend_mat <- apply(my_array, 1:2, \(x) coef(lm(x ~ time_vec))[2])