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

【问题标题】：Linear regression on subsets with dependent variable per column using dlply() in R使用 R 中的 dlply() 对每列具有因变量的子集进行线性回归
【发布时间】：2016-10-09 23:32:41
【问题描述】：

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

我的数据框包括一列时间类别，一列 (slope$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 和 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.

希望有人能帮帮我。

提前非常感谢！

【问题讨论】：

您能说得更具体些吗？我在阅读本文时遇到的一个问题是了解您是否想要对每列进行简单线性回归或对整个子集 t0、t1、...进行多元线性回归。
我很抱歉造成混乱。我正在尝试对每一列进行简单的线性回归。
你只想要斜率、截距等。你到底想要什么？运行lm()函数总结了很多信息。
主要我想要每个线性回归的系数用于进一步计算。但是，不仅提取系数，而且无论如何拥有模型都会非常有帮助。在这种情况下，我想绘制回归图以查看值的拟合程度。

标签： r statistics subset plyr linear-regression

【解决方案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'，它相当于您的“时间点”变量。在这种情况下，变量 'mpg' 是 y 变量，相当于数据中的 'Abs'。

根据我上面对斜率和截距的解释，很明显我们需要三个表/数据集：您的 y 相对于您的 x 的相关数据集每个组，每个变量和组的标准差表，以及每个组和每个变量的均值表。

要获取相关数据集，我们要按'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。

【讨论】：

对不起，我长时间的沉默，我已经尝试过你的答案，但不幸的是我遇到了一些错误，所以我决定改用不同的方法。不过还是谢谢你这么努力地帮助我！
实际上，当您最初发布它时，我已经这样做了，现在我又这样做了。但是它说比声望15的人的投票被记录但没有显示......？所以我想这是我的问题，我是这个社区的新手。

【解决方案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() 函数是可行的。

【讨论】：