在 R 中的变量列表上按组运行线性模型

Question

我有一个数据框，我需要为每个组“站点”运行 6 个 2 变量线性模型。然后，我需要将结果转换为数据框。线性模型中的第二个变量发生了变化。我使用lapply() 将那部分记录下来，但我不知道如何按组运行。我在 SO 上找到了答案，可以回答我的部分问题，但我不知道如何将它们放在一起。

这是一些数据：

structure(list(SiteName = structure(c(1L, 1L, 1L, 1L, 1L, 2L, 
2L, 2L, 2L, 2L, 2L), .Label = c("bp10", "bp12"), class = "factor"), 
    DMWT = c(13.9697916666667, 13.9125, 14.2152083333333, 14.7810416666667, 
    15.1541666666667, 15.7535416666667, 17.3254166666667, 18.4872916666667, 
    20.0564583333333, 21.0595833333333, 21.3925), DMAT = c(16.6714631359947, 
    18.474493439025, 20.9517661662977, 23.7017661662978, 25.5957055602372, 
    20.9688840743375, 23.7188840743375, 25.6128234682769, 27.5143386197921, 
    27.6279749834285, 26.1355507410042), ADD = c(0, 0, 0, 1.90367965367967, 
    5.70129870129876, 0, 1.90367965367967, 5.70129870129876, 
    11.400432900433, 17.2132034632037, 21.53354978355), Air200 = c(7.3229782875097, 
    7.40616010569152, 7.50025101478243, 7.63384949963092, 7.78642525720668, 
    7.51736892282216, 7.65096740767065, 7.80354316524641, 7.97854316524641, 
    8.14729316524641, 8.29592952888278), Air100 = c(15.2711601056916, 
    15.362599499631, 15.512902529934, 15.727296469328, 15.9717661662977, 
    15.5300204379738, 15.7444143773677, 15.9888840743374, 16.2306264985798, 
    16.4472174076707, 16.6433537713071), Air75 = c(16.8986348531664, 
    17.0426752572068, 17.1927762673078, 17.3687358632674, 17.5567156612472, 
    17.2098941753475, 17.3858537713071, 17.5738335692869, 17.7820153874687, 
    18.0100961955496, 18.2532275086809), Air50 = c(19.5072207117523, 
    19.6340388935705, 19.7382813178129, 19.8887358632674, 20.1060085905402, 
    19.7553992258526, 19.9058537713072, 20.1231264985799, 20.4400961955496, 
    20.7669143773678, 20.9841871046405), Air10 = c(21.9214631359947, 
    21.5850994996311, 21.2563116208432, 21.1714631359947, 21.4502510147826, 
    21.2734295288829, 21.1885810440344, 21.4673689228223, 21.9696416500951, 
    22.3779749834284, 22.5476719531254)), .Names = c("SiteName", 
"DMWT", "DMAT", "ADD", "Air200", "Air100", "Air75", "Air50", 
"Air10"), row.names = c(547L, 548L, 549L, 550L, 551L, 1593L, 
1594L, 1595L, 1596L, 1597L, 1598L), class = "data.frame")

这是在模型中使用每个变量的代码。我如何使用这些网站？：

siteslist <- unique(d$SiteName) 
varlist <- names(d)[4:9]
models <- lapply(varlist, function(x) {  # apply the modeling function to our list of air variables
  lm(substitute(DMWT ~ DMAT + i, list(i = as.name(x))), data = d)  # linear model with air variable substituted
})

然后获取模型结果并转换为数据框：

library(relaimpo)
sumfun <- function(x) c(coef(x),
                        summary(x)$adj.r.squared,
                        sqrt(mean(resid(x)^2,na.rm=TRUE)),
                        calc.relimp(x,type="betasq")$betasq[1],
                        calc.relimp(x,type="betasq")$betasq[2],
                        calc.relimp(x,type="pratt")$pratt[1],
                        calc.relimp(x,type="pratt")$pratt[2])
mod.df <- as.data.frame(t(sapply(models,sumfun)))

还尝试结合变量和站点来做这样的事情：

siteslist <- unique(d$SiteName)                              
varlist <- names(d)[4:9]
sets <- expand.grid(SiteName = siteslist, Var = varlist)
models <- lapply(1:nrow(sets), function(x) {  # apply the modeling function to our list of air variables
  lm(substitute(DMWT ~ DMAT + i, list(i = as.name(sets$Var[x]))), data = d[d$SiteName ==  sets$SiteName[x],])  # linear model with air variable substituted
})

...但我收到一个错误"Error in eval(expr, envir, enclos) : object '1' not found"

Answer 1

我会这样做。请注意，这是未经测试的，因为我还没有安装 relaimpo。我真的只是在重新打包你的代码。

一般的方法是 1. 开发一个适用于一组的功能 2. 使用split 将您的数据分组 3.使用lapply将功能应用到每个组 4.（如果需要）将结果组合在一起

对单一站点数据子集的第一次测试：

我所做的唯一更改是 (a) 为一个站点提取数据子集并将其命名为 one_site。 (b) 在您的建模代码中使用one_site。 (c) 我更喜欢将公式作为字符串粘贴在一起而不是使用substitute，所以我进行了更改。 (d) 用于可读性的空格和格式（主要使用 RStudio 的“重新格式化代码”）。

## set up
varlist <- names(d)[4:9]
library(relaimpo)
sumfun <- function(x) {
    c(
        coef(x),
        summary(x)$adj.r.squared,
        sqrt(mean(resid(x) ^ 2, na.rm = TRUE)),
        calc.relimp(x, type = "betasq")$betasq[1],
        calc.relimp(x, type = "betasq")$betasq[2],
        calc.relimp(x, type = "pratt")$pratt[1],
        calc.relimp(x, type = "pratt")$pratt[2]
    )
}

## Testing: this works for one_site
one_site <- subset(d, SiteName == "bp10")

models <- lapply(varlist, function(x) {  # apply the modeling function to our list of air variables
    form <- as.formula(sprintf("DMWT ~ DMAT + %s", x))
    lm(form, data = one_site)  # linear model with air variable substituted
})

## desired result
mod.df <- as.data.frame(t(sapply(models, sumfun)))

把它变成一个函数

一旦您拥有适用于单个网站的代码，我们就会将其转化为函数。唯一的输入似乎是一个站点的数据和varlist 中的变量。我们没有在底部分配结果，而是return它：

fit_one_site = function(one_site, varlist) {
    models <- lapply(varlist, function(x) {
            # apply the modeling function to our list of air variables
            form = as.formula(sprintf("DMWT ~ DMAT + %s", x))
            lm(form, data = one_site)  # linear model with air variable substituted
    })
    return(as.data.frame(t(sapply(models, sumfun))))
}

现在我们可以使用split 将您的数据拆分为SiteName，并使用lapply 将fit_one_site 函数应用于每个部分。

results = lapply(split(d, d$SiteName), FUN = fit_one_site, varlist = names(d)[4:9])

结果应该是数据框列表，每个站点一个。如果您想将它们组合成一个数据框，请参阅我的答案at the list of data frames R-FAQ的相关部分。

Answer 2

我不确定这是否正是您想要做的，但 data.table plyr 包允许您运行由多个变量拆分的模型。下面是一个示例，var1 和 var2 仅表示两个变量，您希望对每个值组合分别建模。

#load packages
library(data.table)
library(plyr)

#break up by variables, then fit the model to each piece
models <- dlply(data, c("var1","var2"),
              function(data)
                lm(DV ~ 
                     IV1 + IV2
                   , data = data, weights = weights))
#apply coef to eah model and return a df
models_coef <- ldply(models, coef)
#print summary
l_ply(models_coef, summary, .print = T)