在 R 中使用 dlply() 对每列具有因变量的子集进行线性回归

Question

我想分别为每个类别的数据框自动生成线性回归。

我的数据框包括一列时间类别，一列（斜率$Abs）作为因变量，几列应该用作自变量。

head(slope)
   timepoint   Abs      In1      In2      In3     Out1     Out2     Out3 ...
1:        t0 275.0 2.169214 2.169214 2.169214 2.069684 2.069684 2.069684
2:        t0 275.5 2.163937 2.163937 2.163937 2.063853 2.063853 2.063853
3:        t0 276.0 2.153298 2.158632 2.153298 2.052088 2.052088 2.057988
4: ...

总而言之，对于每个时间点，我有 40 个变量，我想对每个组合进行线性回归。如 In1~Abs[t0], In1~Abs[t1] 等每一列。当然我可以手动完成这个，但我想一定有更优雅的方式来完成这项工作。

我做了研究，发现 dlply() 可能是我正在寻找的功能。但是，我的尝试导致错误。

所以我以某种方式尝试结合以前发现的问题的答案： On individual variables per column and on subsets per category

我想出了一个这样的函数：

lm.fun <- function(x) {summary(lm(x ~ slope$Abs, data=slope))}
lm.list <- dlply(.data=slope, .variables=slope$timepoint, .fun=lm.fun )

但是我得到以下错误：

Error in eval.quoted(.variables, data) : 
   envir must be either NULL, a list, or an environment.

希望有人能帮帮我。

非常感谢！

Answer 1

根据我的研究，R 中的 dplyr 包不能很好地将 y~x 形式的公式接受到其函数中。所以另一种选择是人工计算。现在让我首先通知您 slope = cor(x,y)*sd(y)/sd(x)（在此处找到的参考资料：http://faculty.cas.usf.edu/mbrannick/regression/regbas.html）和 intercept = mean(y) - slope*mean(x)。简单线性回归要求我们在找到截距时使用质心作为参考点，因为它是一个无偏估计量。使用单个点只会让您截取该单个点，而不是整体截距。

现在对于这个解释，我将使用 mtcars 数据集。我只想要数据的一个子集，所以我使用变量 c('mpg', 'cyl', 'disp', 'hp', 'drat', 'wt', 'qsec') 来基本上模仿您的数据集。在我的示例中，我的分组变量是 'cyl'，它等同于您的 'timepoint' 变量。在这种情况下，变量 'mpg' 是 y 变量，相当于数据中的 'Abs'。

根据我上面对斜率和截距的解释，很明显我们需要三个 tables/datasets：y 相对于 [=每个组 65=]x，每个变量和组的标准差 table，每个组和每个变量的均值 table。

要获取相关数据集，我们要按 'cyl' 分组并计算的相关系数，您应该使用：

df <- mtcars[c('mpg', 'cyl', 'disp', 'hp', 'drat', 'wt', 'qsec')]
corrs <- data.frame(k1 %>% group_by(cyl) %>% do(head(data.frame(cor(.[,c(1,3:7)])), n = 1)))

由于我的数据集的结构方式，第二个变量 (df[ ,2]) 是 'cyl'。对于你，你应该使用

do(head(data.frame(cor(.[,c(2:40)])), n = 1)))

因为你的第一列是分组变量，它不是数字。本质上，您想要遍历所有数字变量。不使用 head 会产生相关矩阵，但由于您有兴趣找到彼此独立的斜率 x-变量，因此您只需要具有相关系数的行您的 y-变量等于 1 (r_yy = 1)。

要获得每个组、每个变量的标准差和均值，请使用

sds     <- data.frame(k1 %>% group_by(cyl) %>% summarise_each(funs(sd)))
means   <- data.frame(k1 %>% group_by(cyl) %>% summarise_each(funs(mean)))

您的组名将是第一列，因此请确保为每个数据集 corrs、sds 和 means 重命名您的行并删除第 1 列。

rownames(corrs) <- rownames(means) <- rownames(sds) <- corrs[ ,1]
corrs <- corrs[ ,-1]; sds <- sds[ ,-1]; means <- means[ ,-1]

现在我们需要计算sd(y)/sd(x)。我完成并看到它完成的最好方法是使用 apply 附属函数。

sdst <- data.frame(t(apply(sds, 1, function(X) X[1]/X)))

我使用 X[1] 因为 sds 中的第一个变量是我的 y-变量。删除 timepoint 后的第一个变量是 Abs，这是您的 y-变量。所以用那个。

现在剩下的就很简单了。由于所有内容都保存为数据框，要找到坡度，您需要做的就是

slopes    <- sdst*corrs
inter     <- slopes*means
intercept <- data.frame(t(apply(inter, 1, function(x) x[1]-x)))

同样在这里，因为我们的 y-变量在第一列，所以我们使用 x[1]。要检查是否一切正常，y 变量的斜率应为 1，截距应为 0。

Answer 2

我已经用更简单的方法解决了这个问题，所以我想更新答案。

为了让生活更轻松，我转换了数据帧结构，以便使用 reshape 包的 melt() 函数将所有列转换为行。

melt(slope, id = c("Abs", "timepoint"), variable_name = "Sites")

输出的列名默认为 "value"。

然后创建一列，添加两个预测变量 paste()。

slope$FullTreat <- paste(slope$Sites,slope$timepoint, sep="_")

运行通过数据集为每个治疗组合创建单独模型的函数。

models <- dlply(slope, ~ FullTreat, function(df) { 
          lm(value ~ Abs, data = df)
          })

简单地提取系数运行

coefs <- ldply(models, coef)

然后使用 colsplit() 也从 reshape 再次将 FullTreat 列拆分为单独的列。另外，将截距和斜率添加到新数据框：

coefs <- cbind(colsplit(coefs$FullTreat, split="_",
         c("Sites","Timepoint")), coefs[,2:3])

我还没有研究过绘制模型所有回归的函数，但我想这对于 ldply() 函数是可行的。

在 R 中使用 dlply() 对每列具有因变量的子集进行线性回归

Linear regression on subsets with dependent variable per column using dlply() in R

statistics

r

subset

linear-regression

plyr