在多列数据上使用 rollapply 和 lm答案

【问题标题】：Using rollapply and lm over multiple columns of data在多列数据上使用 rollapply 和 lm
【发布时间】：2016-02-22 04:05:06
【问题描述】：

我有一个类似于以下的数据框，共有 500 列：

Probes <- data.frame(Days=seq(0.01, 4.91, 0.01), B1=5:495,B2=-100:390, B3=10:500,B4=-200:290)

我想计算一个滚动窗口线性回归，其中我的窗口大小为 12 个数据点，每个顺序回归由 6 个数据点分隔。对于每个回归，“天数”始终是模型的 x 分量，y 将是其他每一列（B1，然后是 B2、B3 等）。然后，我想将系数保存为具有现有列标题（B1、B2 等）的数据框。

我认为我的代码很接近，但不能正常工作。我使用了动物园图书馆的 rollapply。

slopedata<-rollapply(zoo(Probes), width=12, function(Probes) { 
 coef(lm(formula=y~Probes$Days, data = Probes))[2]
 }, by = 6, by.column=TRUE, align="right")

如果可能，我还希望将“xmins”保存到向量中以添加到数据框中。这意味着每个回归中使用的最小 x 值（基本上它是“天”列中的每 6 个数字。）感谢您的帮助。

【问题讨论】：

如果我正确理解了您的问题，您将需要在两个方向上“循环”。 zoo 中的 rollapply 将为您提供窗口方向（向下移动行）。但是，如果您要进行多个 y 回归，您还需要遍历每一种可能性（每一列都是另一个回归）
好的，谢谢。你是说除了 rollapply 函数之外我还需要循环它吗？我不太确定放在一起会是什么样子。
你也认为 plyr 的 colwise 函数可以解决这个问题吗？
这取决于你到底想要什么。你想要回归 y1 vs x, y2 vs x, y3 vs x,...等。通过他们自己。或者除了滚动窗口 y1 与 x（仅前 12 即 1-12）、y1 与 x（仅 18-30）、...等
我想要第二个....除了滚动窗口之外的列回归。几乎你写的除了滚动窗口每 6 个单元格移动一次，像这样：Y1[1:12]~x、Y1[6:18]~x、Y1[12:24] 等等。所以我会基本上最终得到一个矩阵，其值是我最初拥有的值的 1/6。将为此制作一个示例矩阵，以便清楚。

标签： r zoo lm rollapply

【解决方案1】：

你也可以使用rollRegres包如下

# setup data
Probes <- data.frame(
  # I changed the days to be intergers
  Days=seq(1L, 491L, 1L),
  B1=5:495, B2=-100:390, B3=10:500 , B4=-200:290)

# setup grp argument
grp_arg <- as.integer((Probes$Days - 1L) %/% 6)

# estimate coefs. width argument is realtive in grp units
library(rollRegres)
X <- cbind(1, Probes$Days / 100)
Ys <- as.matrix(Probes[, 2:5])
out <- lapply(1:ncol(Ys), function(i)
  roll_regres.fit(x = X, y = Ys[, i], width = 2L, grp = grp_arg)$coefs)
out <- do.call(cbind, out)

# only keep the complete.cases and the unique values
colnames(out) <- sapply(1:4, function(i) paste0("B", i, 0:1))
out <- out[c(T, grp_arg[-1] != head(grp_arg, -1)), ]
out <- out[complete.cases(out), ]
head(out)
#R      B10 B11  B20 B21 B30 B31  B40 B41
#R [1,]   4 100 -101 100   9 100 -201 100
#R [2,]   4 100 -101 100   9 100 -201 100
#R [3,]   4 100 -101 100   9 100 -201 100
#R [4,]   4 100 -101 100   9 100 -201 100
#R [5,]   4 100 -101 100   9 100 -201 100
#R [6,]   4 100 -101 100   9 100 -201 100

该解决方案比 zoo 解决方案快很多

library(zoo) coefs <- function(z) c(unlist(as.data.frame(coef(lm(z[,-1] ~ z[,1]))))) microbenchmark::microbenchmark(   rollapply = {
    z <- zoo(Probes, Probes[[1]])
    rz <- rollapply(z, 12, by = 6, coefs, by.column = FALSE, align = "left")   },   roll_regres = {
    grp_arg <- as.integer((Probes$Days - 1L) %/% 6)

    X <- cbind(1, Probes$Days / 100)
    Ys <- as.matrix(Probes[, 2:5])
    out <- lapply(1:ncol(Ys), function(i)
      roll_regres.fit(x = X, y = Ys[, i], width = 2L, grp = grp_arg)$coefs)
    out <- do.call(cbind, out)

    colnames(out) <- sapply(1:4, function(i) paste0("B", i, 0:1))
    out <- out[c(T, grp_arg[-1] != head(grp_arg, -1)), ]
    out <- out[complete.cases(out), ]
    head(out)   } )
#R Unit: microseconds
#R        expr       min        lq      mean     median        uq       max neval
#R   rollapply 53392.614 56330.492 59793.106 58363.2825 60902.938 119206.76   100
#R roll_regres   865.186   920.297  1074.161   983.9015  1047.705   5071.41   100

目前，由于版本0.1.0 的验证出错，您需要从 Github 安装包。因此，运行

devtools::install_github("boennecd/rollRegres", upgrade_dependencies = FALSE,
                         build_vignettes = TRUE)

【讨论】：

【解决方案2】：

试试这个：

# here are the xmin values you wanted
xmins <- Probes$Days[seq(1,nrow(Probes),6)]

# here we build a function that will run regressions across the columns
# y1 vs x, y2 vs x, y3 vs x...
# you enter the window and by (12/6) in order to limit the interval being
# regressed. this is later called in do.call
runreg <- function(Probes,m,window=12,by=6){

  # beg,end are used to specify the interval
  beg <- seq(1,nrow(Probes),by)[m]
  end <- beg+window-1

  # this is used to go through all the columns
  N <- ncol(Probes)-1
  tmp <- numeric(N)
  # go through each column and store the coefficients in tmp
  for(i in 1:N){
     y <- Probes[[i+1]][beg:end]
     x <- Probes$Days[beg:end]
     tmp[i] <- coef(lm(y~x))[2][[1]]
  }
  # put all our column regressions into a dataframe
  res <- rbind('coeff'=tmp)
  colnames(res) <- colnames(Probes)[-1]

  return(res)
}

# now that we've built the function to do the column regressions
# we just need to go through all the window-ed regressions (row regressions)
res <- do.call(rbind,lapply(1:length(xmins),function(m) runreg(Probes,m)))

# these rownames are the index of the xmin values
rownames(res) <- seq(1,nrow(Probes),6)
res <- data.frame(res,xmins)

【讨论】：

是的，这很好用。我在 Excel 中仔细检查一些斜率时得到的值略有不同，但我猜这可能是由于两个不同程序拟合线性模型的方式不同，而不是编码问题。感谢您的帮助。
这也可能是索引的问题。也许我是一个单位？但这很容易解决。一般结构应该有帮助！尝试更改end <- beg+window-1，我相信这将使其正确。索引现在将变为：1:(1+12-1) 而不是之前的 1:(1+12)

【解决方案3】：

1) 定义一个动物园对象z，其数据包含Probes，其索引取自Probes的第一列，即Days。注意lm 允许y 是一个矩阵，定义了一个计算回归系数的coefs 函数。最后rollapply 超过z。请注意，返回对象的索引给出了 xmin。

library(zoo)

z <- zoo(Probes, Probes[[1]])

coefs <- function(z) c(unlist(as.data.frame(coef(lm(z[,-1] ~ z[,1])))))
rz <- rollapply(z, 12, by = 6, coefs, by.column = FALSE, align = "left")

给予：

> head(rz)
     B11 B12  B21 B22 B31 B32  B41 B42
0.01   4 100 -101 100   9 100 -201 100
0.07   4 100 -101 100   9 100 -201 100
0.13   4 100 -101 100   9 100 -201 100
0.19   4 100 -101 100   9 100 -201 100
0.25   4 100 -101 100   9 100 -201 100
0.31   4 100 -101 100   9 100 -201 100

请注意，如果您需要 rz 的数据框表示，则可以使用 DF <- fortify.zoo(rz)。

2)另一种有点类似的方法是在行号上使用rollaplly：

library(zoo)
y <- as.matrix(Probes[-1])
Days <- Probes$Days
n <- nrow(Probes)
coefs <- function(ix) c(unlist(as.data.frame(coef(lm(y ~ Days, subset = ix)))), 
      xmins = Days[ix][1])
r <- rollapply(1:n, 12, by = 6, coefs)

【讨论】：

请注意，这已被修改为更短。
谢谢，这段代码也可以工作，并且值与我在 Excel 中测试一些后得到的值相匹配。我不需要数据框中的截距，所以我猜如果我想消除这些列，我只需在 coef 函数的某处添加一个 [-1] ？
在这种情况下，coefs 函数可以简化为：coefs <- function(z) coef(lm(z[,-1] ~ z[,1]))[2,] in (1)。
@user507 我想知道你是否可以同时运行这两个代码，看看哪一个更高效
在上述两个相似的代码中，第二个（滚动应用行号）对于我的目的来说是两个中更好的，也是我最终使用的那个。下面的代码有效，但我无法解决值略微偏离的扭结。正如您所建议的，我最初也认为这可能是一个索引问题，但尝试多个不同的索引更改并没有解决问题。