自动化回归由行

我有一个data.frame自动化回归由行

set.seed(100) 
exp <- data.frame(exp = c(rep(LETTERS[1:2], each = 10)), re = c(rep(seq(1, 10, 1), 2)), age1 = seq(10, 29, 1), age2 = seq(30, 49, 1), 
        h = c(runif(20, 10, 40)), h2 = c(40 + runif(20, 4, 9)))

我想打在数据各行lm设置(h and h2 ~ age1 and age2) 我通过loop

exp$modelh <- 0 

for (i in 1:length(exp$exp)){ 
    age = c(exp$age1[i], exp$age2[i]) 
    h = c(exp$h[i], exp$h2[i]) 
    model = lm(age ~ h) 
    exp$modelh[i] = coef(model)[1] + 100 * coef(model)[2] 

}

和它的作品做但对于非常大的文件需要一些时间。将感谢f.ex更快的解决方案。 dplyr

来源

2015-08-28 Mateusz1981

不，它是行的，即使这些组存在 – Mateusz1981

我很抱歉，你是否试图创建一个只有1个自由度的回归？我可能会建议你重新考虑你的行动计划...... – Jason

@Jason，这只是一个更大的问题的例子 – Mateusz1981

使用dplyr，我们可以试用rowwise()和do。在do的内部，我们将'age1'，'age2'连接（c）以创建'age'，同样，我们可以创建'h'，应用lm，提取coef以创建'modelh'列。

library(dplyr) 
exp %>% 
    rowwise() %>% 
    do({ 
     age <- c(.$age1, .$age2) 
     h <- c(.$h, .$h2) 
     model <- lm(age ~ h) 
     data.frame(., modelh = coef(model)[1] + 100*coef(model)[2]) 
    })

使输出

# exp re age1 age2  h  h2 modelh 
#1 A 1 10 30 19.23298 46.67906 68.85506 
#2 A 2 11 31 17.73018 47.55402 66.17050 
#3 A 3 12 32 26.56967 46.69174 84.98486 
#4 A 4 13 33 11.69149 47.74486 61.98766 
#5 A 5 14 34 24.05648 46.10051 82.90167 
#6 A 6 15 35 24.51312 44.85710 89.21053 
#7 A 7 16 36 34.37208 47.85151 113.37492 
#8 A 8 17 37 21.10962 48.40977 74.79483 
#9 A 9 18 38 26.39676 46.74548 90.34187 
#10 A 10 19 39 15.10786 45.38862 75.07002 
#11 B 1 20 40 28.74989 46.44153 100.54666 
#12 B 2 21 41 36.46497 48.64253 125.34773 
#13 B 3 22 42 18.41062 45.74346 81.70062 
#14 B 4 23 43 21.95464 48.77079 81.20773 
#15 B 5 24 44 32.87653 47.47637 115.95097 
#16 B 6 25 45 30.07065 48.44727 101.10688 
#17 B 7 26 46 16.13836 44.90204 84.31080 
#18 B 8 27 47 20.72575 47.14695 87.00805 
#19 B 9 28 48 20.78425 48.94782 84.25406 
#20 B 10 29 49 30.70872 44.65144 128.39415

我们可以与devel版本的data.table即v1.9.5做到这一点。说明安装devel版本是here。

我们将'data.frame'转换为'data.table'（setDT），使用选项keep.rownames=TRUE创建列'rn'。我们melt通过指定measure中的patterns将数据集从'wide'转换为'long'格式。按照'rn'分组，我们执行lm并获得coef。这可以通过分配（:=）到NULL来分配原始数据集（'exp'）中的新列，同时删除不需要的'rn'列。

library(data.table)#v1.9.5+ 
modelh <- melt(setDT(exp, keep.rownames=TRUE), measure=patterns('^age', '^h'), 
    value.name=c('age', 'h'))[, {model <- lm(age ~h) 
     coef(model)[1] + 100 * coef(model)[2]},rn]$V1 

exp[, modelh:= modelh][, rn := NULL] 
exp 
# exp re age1 age2  h  h2 modelh 
# 1: A 1 10 30 19.23298 46.67906 68.85506 
# 2: A 2 11 31 17.73018 47.55402 66.17050 
# 3: A 3 12 32 26.56967 46.69174 84.98486 
# 4: A 4 13 33 11.69149 47.74486 61.98766 
# 5: A 5 14 34 24.05648 46.10051 82.90167 
# 6: A 6 15 35 24.51312 44.85710 89.21053 
# 7: A 7 16 36 34.37208 47.85151 113.37492 
# 8: A 8 17 37 21.10962 48.40977 74.79483 
# 9: A 9 18 38 26.39676 46.74548 90.34187 
#10: A 10 19 39 15.10786 45.38862 75.07002 
#11: B 1 20 40 28.74989 46.44153 100.54666 
#12: B 2 21 41 36.46497 48.64253 125.34773 
#13: B 3 22 42 18.41062 45.74346 81.70062 
#14: B 4 23 43 21.95464 48.77079 81.20773 
#15: B 5 24 44 32.87653 47.47637 115.95097 
#16: B 6 25 45 30.07065 48.44727 101.10688 
#17: B 7 26 46 16.13836 44.90204 84.31080 
#18: B 8 27 47 20.72575 47.14695 87.00805 
#19: B 9 28 48 20.78425 48.94782 84.25406 
#20: B 10 29 49 30.70872 44.65144 128.39415

来源

2015-08-28 11:21:10 akrun

我不知道结果是否可以通过'apply'获得？ – Mateusz1981

@ Mateusz1981它可以通过'apply'获得，但'for'和'apply'的速度可能没有太大的差别。 – akrun

到目前为止，第一个解决方案是完美的。与2，有安装问题 – Mateusz1981

来自@akrun的好（双）回答。

就像您提到的“这是一个更大的问题的例子”，您的未来分析只是一个建议。很明显，如果你真的对建立模型有兴趣，那么随着年龄和观测值的增加，你会创建越来越多的列。如果你得到N个观测值，那么你只能使用2N个列来表示这2个变量。

我建议使用长数据格式来增加行数而不是列数。

是这样的：如果“大问题”指的是别的东西，这答案是不相关的

exp[1,] # how your first row (model building info) looks like 

# exp re age1 age2  h  h2 
# 1 A 1 10 30 19.23298 46.67906 


reshape(exp[1,],         # how your model building info is transformed 
     varying = list(c("age1","age2"), 
           c("h","h2")), 
     v.names = c("age_value","h_value"), 
     direction = "long") 

#  exp re time age_value h_value id 
# 1.1 A 1 1  10 19.23298 1 
# 1.2 A 1 2  30 46.67906 1

道歉。

来源

2015-08-28 11:54:49 AntoniosK

我知道有些重塑是一种解决方法，但没有找到正确的答案，谢谢 – Mateusz1981

我很高兴它是有用的。现在，建模不是按行，而是按组来划分。你以前的独特行标识符是什么，现在它是你的组标识符。在这个例子中，你的行标识符是变量“exp”和“re”{A，1}的组合，所以在新格式中你的分组仍然是{A，1}，但它现在对应于2行。 – AntoniosK

使用base R，函数sprintf可以帮助我们创建公式。并lapply进行计算。

strings <- sprintf("c(%f,%f) ~ c(%f,%f)", exp$age1, exp$age2, exp$h, exp$h2) 
lst <- lapply(strings, function(x) {model <- lm(as.formula(x));coef(model)[1] + 100 * coef(model)[2]}) 
exp$modelh <- unlist(lst) 
exp 
# exp re age1 age2  h  h2 modelh 
# 1 A 1 10 30 19.23298 46.67906 68.85506 
# 2 A 2 11 31 17.73018 47.55402 66.17050 
# 3 A 3 12 32 26.56967 46.69174 84.98486 
# 4 A 4 13 33 11.69149 47.74486 61.98766 
# 5 A 5 14 34 24.05648 46.10051 82.90167 
# 6 A 6 15 35 24.51312 44.85710 89.21053 
# 7 A 7 16 36 34.37208 47.85151 113.37493 
# 8 A 8 17 37 21.10962 48.40977 74.79483 
# 9 A 9 18 38 26.39676 46.74548 90.34187 
# 10 A 10 19 39 15.10786 45.38862 75.07002 
# 11 B 1 20 40 28.74989 46.44153 100.54666 
# 12 B 2 21 41 36.46497 48.64253 125.34773 
# 13 B 3 22 42 18.41062 45.74346 81.70062 
# 14 B 4 23 43 21.95464 48.77079 81.20773 
# 15 B 5 24 44 32.87653 47.47637 115.95097 
# 16 B 6 25 45 30.07065 48.44727 101.10688 
# 17 B 7 26 46 16.13836 44.90204 84.31080 
# 18 B 8 27 47 20.72575 47.14695 87.00805 
# 19 B 9 28 48 20.78425 48.94782 84.25406 
# 20 B 10 29 49 30.70872 44.65144 128.39416

在lapply函数的表达式as.formula(x)就是在第一行中创建的公式转换成由lm功能可用的格式。

基准

library(dplyr) 
library(microbenchmark) 
set.seed(100) 
big.exp <- data.frame(age1=sample(30, 1e4, T), 
         age2=sample(30:50, 1e4, T), 
         h=runif(1e4, 10, 40), 
         h2= 40 + runif(1e4,4,9)) 

microbenchmark(
    plafort = {strings <- sprintf("c(%f,%f) ~ c(%f,%f)", big.exp$age1, big.exp$age2, big.exp$h, big.exp$h2) 
      lst <- lapply(strings, function(x) {model <- lm(as.formula(x));coef(model)[1] + 100 * coef(model)[2]}) 
      big.exp$modelh <- unlist(lst)}, 

    akdplyr = {big.exp %>% 
    rowwise() %>% 
    do({ 
     age <- c(.$age1, .$age2) 
     h <- c(.$h, .$h2) 
     model <- lm(age ~ h) 
     data.frame(., modelh = coef(model)[1] + 100*coef(model)[2]) 
    })} 

,times=5) 
t: seconds 
    expr  min  lq  mean median  uq  max neval cld 
plafort 13.00605 13.41113 13.92165 13.56927 14.53814 15.08366  5 a 
akdplyr 26.95064 27.64240 29.40892 27.86258 31.02955 33.55940  5 b

（注：我下载今天data.table的最新1.9.5开发人员版本，而是继续尝试着去测试它时收到错误结果也各不相同分数（1.93×10^-8）。舍入可能占的差异。）

all.equal(pl, ak) 
[1] "Attributes: < Component “class”: Lengths (1, 3) differ (string compare on first 1) >" 
[2] "Attributes: < Component “class”: 1 string mismatch >"         
[3] "Component “modelh”: Mean relative difference: 1.933893e-08"

结论

与dplyr相比，lapply方法在速度方面似乎表现良好，但它的5位舍入可能是个问题。改进可能是可能的。转换为矩阵后可能使用apply以提高速度和效率。

来源

2015-08-28 13:41:32

自动化回归由行

回答

相关问题