2016-08-15 91 views
1

我使用模型來查看作物區域的溫室氣體排放量。爲了測試數據的標準偏差,我試圖通過多次迭代對它進行蒙特卡羅風格分析。R遍歷數據表列中的模型

model parameters 
a <- 0.1474 # Alpha 
b <- 0.0005232 # Beta 
g <- -0.00001518 # Gamma 
d <- 0.000003662 # Delta 
rain <- crm$rain # rainfall value for that location from the col 'rain' 

的數據是在data.table爲以下但在N列從N1-N100運行:

 rn rain Wheat  N1  N2  N3  N4  N5  N6 
# 1: 10007 1049.61 0.1718 0.6363109 0.939479 0.9242736 0.9018818 0.6556216 0.1150655 
# 2: 10018 1114.31 0.1629 0.6363109 0.939479 0.9242736 0.9018818 0.6556216 0.1150655 
# 3: 10023 1361.61 0.1082 0.6363109 0.939479 0.9242736 0.9018818 0.6556216 0.1150655 
# 4: 10024 1407.20 0.0494 0.6363109 0.939479 0.9242736 0.9018818 0.6556216 0.1150655 
# 5: 10025 1499.56 0.0200 0.6363109 0.939479 0.9242736 0.9018818 0.6556216 0.1150655 
# 6: 10026 1654.13 0.0040 0.6363109 0.939479 0.9242736 0.9018818 0.6556216 0.1150655 

所以我的問題是如何申請我下面的模型,每個N列,並添加結果到數據表的結尾?該模型的N值固定,但我正在努力如何從每列中獲取值到模型中。

logN2O <- function (x) {a+(b*rain)+(g*N)+(d*rain*N)} 

非常感謝提前。

編輯

爲了澄清我想先運行與N1的價值模型,並作出新的山坳與結果在最後。然後對N2的值進行相同的操作,直到N列的末尾。

+2

並不很清楚你是什麼之後,但也許你想'lapply(.SD,函數(X){任何需要做的事情,以每列}' – MichaelChirico

+0

@ MichaelChirico我將編輯問題得到更清楚的說明,本質上我想針對每行的N(1:6)值對該行的降雨值運行模型,並將模型輸出放在數據表的末尾 –

回答

1

我認爲這應該工作:

n <- 1:6 
cols <- paste0("N",n,"_res") 
dt[,(cols) := lapply(.SD ,function (x) {a + (b*dt[,rain]) + (g*x) + (d*dt[,rain]*x)}), .SDcols = paste0("N",n)] 

基本上你只指定「N」是你要通過環(在這種情況下,N1 - N6),然後將其存儲附帶「_res」的結果 - 例如「N1_res」。

數據:

dt <- structure(list(rn = c(10007L, 10018L, 10023L, 10024L, 10025L, 
10026L), rain = c(1049.61, 1114.31, 1361.61, 1407.2, 1499.56, 
1654.13), Wheat = c(0.1718, 0.1629, 0.1082, 0.0494, 0.02, 0.004 
), N1 = c(0.6363109, 0.6363109, 0.6363109, 0.6363109, 0.6363109, 
0.6363109), N2 = c(0.939479, 0.939479, 0.939479, 0.939479, 0.939479, 
0.939479), N3 = c(0.9242736, 0.9242736, 0.9242736, 0.9242736, 
0.9242736, 0.9242736), N4 = c(0.9018818, 0.9018818, 0.9018818, 
0.9018818, 0.9018818, 0.9018818), N5 = c(0.6556216, 0.6556216, 
0.6556216, 0.6556216, 0.6556216, 0.6556216), N6 = c(0.1150655, 
0.1150655, 0.1150655, 0.1150655, 0.1150655, 0.1150655)), .Names = c("rn", 
"rain", "Wheat", "N1", "N2", "N3", "N4", "N5", "N6"), class = c("data.table", 
"data.frame"), row.names = c(NA, -6L)) 
+0

完美的工作,非常感謝你的幫助!我的谷歌搜索目前正在使我失望--SD做的是什麼? –

+1

.SD代表** D ** ata的** S ** ubset,實際上是它自己的單獨'data.table'。它通常用於每個'by'組,但在這種情況下,我們只用它來訪問我們感興趣的列(在'.SDcols'中指定的那些列)。 –