2015-08-28 61 views
2

我有一個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

+0

不,它是行的,即使這些組存在 – Mateusz1981

+0

我很抱歉,你是否試圖創建一個只有1個自由度的迴歸?我可能會建議你重新考慮你的行動計劃...... – Jason

+0

@Jason,這只是一個更大的問題的例子 – Mateusz1981

回答

4

使用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.tablev1.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 
+0

我不知道結果是否可以通過'apply'獲得? – Mateusz1981

+0

@ Mateusz1981它可以通過'apply'獲得,但'for'和'apply'的速度可能沒有太大的差別。 – akrun

+0

到目前爲止,第一個解決方案是完美的。與2,有安裝問題 – Mateusz1981

2

來自@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 

道歉。

+0

我知道有些重塑是一種解決方法,但沒有找到正確的答案,謝謝 – Mateusz1981

+0

我很高興它是有用的。現在,建模不是按行,而是按組來劃分。你以前的獨特行標識符是什麼,現在它是你的組標識符。在這個例子中,你的行標識符是變量「exp」和「re」{A,1}的組合,所以在新格式中你的分組仍然是{A,1},但它現在對應於2行。 – AntoniosK

2

使用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以提高速度和效率。