3

我正在使用两个 RasterStack 对象，每个对象由代表单个时间步长的十层组成。

# Mock data
pred.rst.stck <- do.call("stack", lapply(seq(10), function(i) {
  pred.rst <- raster(nrows = 15, ncols = 15, xmn= 0, xmx = 10, ymn = 0, ymx = 10)
  pred.rst[] <- rnorm(225, 50, 10)
  return(pred.rst)
})
resp.rst.stck <- do.call("stack", lapply(seq(10), function(i) {
  resp.rst <- raster(nrows = 10, ncols = 10, xmn = 0, xmx = 10, ymn = 0, ymx = 10)
  resp.rst[] <- rnorm(100, 50, 10)
  return(resp.rst)
})

pred.rst.stck用作一组预测变量和resp.rst.stck一组响应变量。对于预测器 RasterStack 的每个单个单元格，我想在响应 RasterStack 的每个单元格上拟合一个线性模型，为每个拟合模型提取相应的 R 平方并将它们相加。长话短说，这是我迄今为止使用 Rparallel包的最快方法：

# Parallelization
library(parallel)
n.cores <- detectCores()
clstr <- makePSOCKcluster(n.cores)

# Extract cell values from RasterStack objects
pred.vals <- getValues(pred)
resp.vals <- getValues(resp)
clusterExport(clstr, c("pred.vals", "resp.vals"))

# Loop through all predictor cells
rsq.sums <- parLapply(clstr, seq(nrow(pred.vals)), function(i) {
  # For each predictor cell, loop through all response cells, 
  # fit linear model, extract and sum up single R-squared 
  do.call("sum", lapply(seq(nrow(resp.vals)), function(j) {
    summary(lm(resp.vals[j, ] ~ pred.vals[i, ]))$r.squared
  }))
})

虽然parLapply与普通相比性能要好得多lapply，但我想知道是否有一种优雅的方式可以加快整个过程。有什么建议么？

干杯，
弗洛里安

4

1 回答 1

5

您可以尝试一些技巧。我不太明白您创建线性模型的方式，但是r.squared您从线性模型计算的结果相当于 Pearson 相关系数的平方（cor在 R 中使用默认参数），它的计算速度比线性模型快得多模型。

使用您的数据比较这两个函数：

# Finding r.squared using a lm
f1 <- function(n){summary(lm(resp.vals[n,] ~pred.vals[n,]))$r.squared}
# Finding r.squared using Pearson's
f2 <- function(n){ cor(resp.vals[n,],pred.vals[n,])^2}

# Both give the same result
f1(1)
[1] 0.0009032986
f2(1)
[1] 0.0009032986

在计时这些函数的单个调用方面：

require(microbenchmark)
microbenchmark( f1(1) , f2(1) )
Unit: microseconds
  expr      min        lq   median        uq      max neval
 f1(1) 2009.328 2037.0680 2071.045 2124.9225 4102.079   100
 f2(2)   73.075   77.7365   84.690   89.7215 5485.368   100

因此，您应该能够通过使用cor而不是lm.

交换原始函数以供使用的快速系统时间比较cov()^2表明情况如下：

#Using cov()
   user  system elapsed 
  0.013   0.002   1.085

#Using lm()
   user  system elapsed 
  0.159   0.028  26.334
于 2013-04-16T18:15:46.400 回答