Werden wir Helden für einen Tag
Home | About | ArchiveHow to rewrite ugly R code #2
研究一下那段 code ，是很簡單的結構，但卻寫得超長的。
這個 code 我想有六七年歷史，最初代版本是 SPSS Syntax ，後來根據 logic 在 R 重寫。這 code 不停重覆的結構是
# < -
if (age ==1 && gender =="F")
{
bmif <- c(-1.013,16.133,0.07656)
}
else if (age ==1 && gender == "M")
{
bmif <- c(-1.79,16.421,0.07149)
}
今時不同住日，現在會就會想試下去改，當係教自己怎樣寫 program 也好。當年寫這些 code 應該用了不少次 copy and paste ，現在看應該要 generalize 那個 pattern 。
結構上的另一個問題，是數據和邏輯混在一起，要想個方法將兩者分拆。功能上也有問題，例如年齡只有等於整數先可以獲得 LMS 數據。兒童年齡常常有 1.5 甚至 2.3 歲。性別也不一定只輸入為 M/F ，可以是 0/1 ，甚至 Male/Female 。
有這麼多的問題要處理，不如就一步一步的行， divide and conquer 。
先是做一個實驗。假如我有一個 data frame 是這個樣的：
age.min;age.max;male;l.value;m.value;s.value
1;1.5;TRUE;-1.013;16.133;0.07656
我另設兩個 variable ，一個是 x = 1.2 ， y = TRUE 我可不可以用 x >= age.min & x < age.max & male == y 來作 query ，查出這個 data frame 所藏的 LMS 數值呢？
實驗如下：
# < -
lmsdata <- data.frame(age.min = 1, age.max = 1.5, male = TRUE, l = -1.013, m=16.133, s=0.07656)
x <- 1.2
y <- TRUE
lmsdata[x >= lmsdata$age.min & x < lmsdata$age.max & y ==lmsdata$male,]
lmsdata$l[x >= lmsdata$age.min & x < lmsdata$age.max & y ==lmsdata$male]
似乎是 work 的。那就易辨了。
首先，是要將上次梁淑芳教授 ( 1998 年) 那堆 LMS 數據轉成 data frame 。人太懶，最易的方法，是用 excel 之類的 spreadsheet 入。今次我是用 Gnumeric ，用甚麼軟件其實不重要。
將數據存為 csv 格式，再 read 入 R 。
# < -
lmsdata <- read.csv("/home/chainsaw/lmsdata_leung.csv")
str(lmsdata)
x <- 1.2
y <- TRUE
lmsdata[x >= lmsdata$age.min & x < lmsdata$age.max & y ==lmsdata$male,]
lmsdata$l[x >= lmsdata$age.min & x < lmsdata$age.max & y ==lmsdata$male]
似乎可以和實驗一樣從 data frame 找回 lms 數據。可以試下透過修改 x 及 y 值 query 其他的年齡和性別組別。似乎最後版本的 function 不再需要原版那一層又一層的 if-else if 的 pattern 了。即開始改寫原來的 function 。因為某些原因，我並不想用回原來的名稱 bmisds() ，我想改名為 sds.cal()
#< -
sds.cal <- function (value, age, male, lmsdata) {
output <- c()
for(i in 1:length(value)) {
lms.list <- c(lmsdata$l[age[i] >= lmsdata$age.min & age[i] < lmsdata$age.max & male[i] ==lmsdata$male], lmsdata$m[age[i] >= lmsdata$age.min & age[i] < lmsdata$age.max & male[i] ==lmsdata$male], lmsdata$s[age[i] >= lmsdata$age.min & age[i] < lmsdata$age.max & male[i] ==lmsdata$male])
output[i]<- ((((value[i]/lms.list[2])**lms.list[1] -1)/(lms.list[1]*lms.list[3])))
}
return(output)
}
用剛才的 lmsdata ，這個 function 好像可以生產出一些結果。
# < -
sds.cal(value=13,age=12,male=TRUE,lmsdata=lmsdata)
但是這個結果可信嗎？
要做個 validation test ，先要隨機找幾個 test cases 。最好的方法是用 R 隨機生產出十個出來了。
# < -
tc.bmi <- runif(10, 1.0, 40.0)
tc.age <- runif(10, 0.0, 20.0)
tc.male <- sample(c(TRUE,FALSE),10, replace=TRUE)
我生產出來的數值是：
> tc.bmi
[1] 28.794667 6.351955 38.040112 35.271790 29.722005 19.566648 11.579569
[8] 35.816091 27.741749 39.451059
> tc.age
[1] 14.176405 7.741322 19.344101 7.509921 13.919720 7.284554 15.643284
[8] 11.755625 16.166327 5.418958
> tc.male
[1] TRUE TRUE TRUE TRUE FALSE TRUE FALSE TRUE FALSE TRUE
這十個 test cases 用手計 ((講大話的，我是用 gnumeric 計的。)) 的 BMI SDS 是：
2.03, -24.64, 2.45, 3.20, 2.42, 1.74, -6.06, 2.64, 1.96, 4.33
可以用這組數字作「黃金標準」，作一點點的 validation 。
# < -
test.results <- sds.cal(tc.bmi, tc.age, tc.male, lmsdata=lmsdata)
epsilon <- 0.01
gold.standard <- c(2.03, -24.64, 2.45, 3.20, 2.42, 1.74, -6.06, 2.64, 1.96, 4.33)
diff.results <- test.results - gold.standard
acceptable <- abs(diff.results) > epsilon
if (sum(acceptable) == 0) {cat("all passed!\n")} else {cat("something wrong!")}
epsilon 是 tolerance 的意思，即是容許 sds.cal 算值與「黃金標準」的差異。這是因為兩個浮點數是很難互相對等的。故此 test.results == gold.standard 多數情況都會回報 FALSE 。 acceptable 是 logical vector ，如果全部都是 FALSE ，而你又去 sum 它，應該是 0 的。今次測驗證實全部答案正確。
最後是小小的整理。
我想 sds.cal 的 default 是用梁淑芳教授的數據作 lmsdata ，只要改少少就可以了。
#< -
sds.cal <- function (value, age, male, lmsdata=leung1998) {
output <- c()
for(i in 1:length(value)) {
lms.list <- c(lmsdata$l[age[i] >= lmsdata$age.min & age[i] < lmsdata$age.max & male[i] ==lmsdata$male], lmsdata$m[age[i] >= lmsdata$age.min & age[i] < lmsdata$age.max & male[i] ==lmsdata$male], lmsdata$s[age[i] >= lmsdata$age.min & age[i] < lmsdata$age.max & male[i] ==lmsdata$male])
output[i]<- ((((value[i]/lms.list[2])**lms.list[1] -1)/(lms.list[1]*lms.list[3])))
}
return(output)
}
leung1998 從何而來？當然可以用原來的方法，每次存入 csv ：
#< -
leung1998 <- read.csv("/home/chainsaw/lmsdata_leung.csv")
但我並不想這樣。我想將數據都儲存在 R source 之內。上網問人，知道有個指令叫 dput() ，可以多加利用。
#< -
dput(lmsdata)
就會 dump 出 lmsdata 的組成方法。 sds.cal 的最終樣貌是這樣的：
# < -
# sds.cal by Chainsaw Riot, released under GPL v3.
# http://www.gnu.org/licenses/gpl.html
sds.cal <- function (value, age, male, lmsdata=leung1998) {
output <- c()
for(i in 1:length(value)) {
lms.list <- c(lmsdata$l[age[i] >= lmsdata$age.min & age[i] < lmsdata$age.max & male[i] ==lmsdata$male], lmsdata$m[age[i] >= lmsdata$age.min & age[i] < lmsdata$age.max & male[i] ==lmsdata$male], lmsdata$s[age[i] >= lmsdata$age.min & age[i] < lmsdata$age.max & male[i] ==lmsdata$male])
output[i]<- ((((value[i]/lms.list[2])**lms.list[1] -1)/(lms.list[1]*lms.list[3])))
}
return(output)
}
leung1998 <- structure(list(age.min = c(0, 0.25, 0.5, 1, 1.5, 2, 2.5, 3, 3.5,
4, 4.5, 5, 5.5, 6, 6.5, 7, 7.5, 8, 8.5, 9, 9.5, 10, 10.5, 11,
11.5, 12, 12.5, 13, 13.5, 14, 14.5, 15, 15.5, 16, 16.5, 17, 17.5,
18, 0, 0.25, 0.5, 1, 1.5, 2, 2.5, 3, 3.5, 4, 4.5, 5, 5.5, 6,
6.5, 7, 7.5, 8, 8.5, 9, 9.5, 10, 10.5, 11, 11.5, 12, 12.5, 13,
13.5, 14, 14.5, 15, 15.5, 16, 16.5, 17, 17.5, 18), age.max = c(0.25,
0.5, 1, 1.5, 2, 2.5, 3, 3.5, 4, 4.5, 5, 5.5, 6, 6.5, 7, 7.5,
8, 8.5, 9, 9.5, 10, 10.5, 11, 11.5, 12, 12.5, 13, 13.5, 14, 14.5,
15, 15.5, 16, 16.5, 17, 17.5, 18, 999, 0.25, 0.5, 1, 1.5, 2,
2.5, 3, 3.5, 4, 4.5, 5, 5.5, 6, 6.5, 7, 7.5, 8, 8.5, 9, 9.5,
10, 10.5, 11, 11.5, 12, 12.5, 13, 13.5, 14, 14.5, 15, 15.5, 16,
16.5, 17, 17.5, 18, 999), male = c(TRUE, TRUE, TRUE, TRUE, TRUE,
TRUE, TRUE, TRUE, TRUE, TRUE, TRUE, TRUE, TRUE, TRUE, TRUE, TRUE,
TRUE, TRUE, TRUE, TRUE, TRUE, TRUE, TRUE, TRUE, TRUE, TRUE, TRUE,
TRUE, TRUE, TRUE, TRUE, TRUE, TRUE, TRUE, TRUE, TRUE, TRUE, TRUE,
FALSE, FALSE, FALSE, FALSE, FALSE, FALSE, FALSE, FALSE, FALSE,
FALSE, FALSE, FALSE, FALSE, FALSE, FALSE, FALSE, FALSE, FALSE,
FALSE, FALSE, FALSE, FALSE, FALSE, FALSE, FALSE, FALSE, FALSE,
FALSE, FALSE, FALSE, FALSE, FALSE, FALSE, FALSE, FALSE, FALSE,
FALSE, FALSE), l = c(-1.868, -1.663, -1.66, -1.79, -1.931, -2.045,
-2.133, -2.2, -2.251, -2.287, -2.313, -2.328, -2.336, -2.337,
-2.332, -2.323, -2.311, -2.296, -2.279, -2.26, -2.241, -2.221,
-2.2, -2.18, -2.159, -2.139, -2.12, -2.101, -2.082, -2.064, -2.047,
-2.03, -2.013, -1.997, -1.982, -1.967, -1.952, -1.938, 0.346,
-0.314, -0.641, -1.013, -1.243, -1.403, -1.522, -1.612, -1.681,
-1.735, -1.776, -1.808, -1.833, -1.851, -1.868, -1.874, -1.879,
-1.882, -1.882, -1.88, -1.877, -1.872, -1.866, -1.86, -1.854,
-1.848, -1.841, -1.835, -1.828, -1.822, -1.816, -1.809, -1.803,
-1.797, -1.791, -1.786, -1.781, -1.775), m = c(13.393, 16.538,
17.055, 16.421, 16.037, 15.861, 15.741, 15.59, 15.426, 15.276,
15.156, 15.072, 15.029, 15.026, 15.066, 15.149, 15.271, 15.423,
15.601, 15.798, 16.009, 16.232, 16.463, 16.702, 16.945, 17.193,
17.444, 17.697, 17.952, 18.205, 18.456, 18.704, 18.947, 19.186,
19.418, 19.644, 19.865, 20.079, 13.125, 15.951, 16.562, 16.133,
15.76, 15.58, 15.485, 15.39, 15.258, 15.097, 14.949, 14.84, 14.772,
14.747, 14.77, 14.839, 14.951, 15.101, 15.284, 15.497, 15.735,
15.993, 16.268, 16.557, 16.854, 17.155, 17.458, 17.757, 18.052,
18.338, 18.646, 18.884, 19.143, 19.393, 19.635, 18.869, 20.096,
20.317), s = c(0.10203, 0.07855, 0.0735, 0.07145, 0.07127, 0.07116,
0.07153, 0.0729, 0.07546, 0.07908, 0.08355, 0.08867, 0.09419,
0.09985, 0.10543, 0.11077, 0.11576, 0.12036, 0.12452, 0.12825,
0.13155, 0.13443, 0.13691, 0.13902, 0.14081, 0.1423, 0.14355,
0.14458, 0.14544, 0.14617, 0.14679, 0.14733, 0.14782, 0.14826,
0.14868, 0.14907, 0.14945, 0.14982, 0.08773, 0.08044, 0.07914,
0.07656, 0.07467, 0.07342, 0.07306, 0.0741, 0.07643, 0.07983,
0.08409, 0.08897, 0.09419, 0.09948, 0.10465, 0.10956, 0.11409,
0.11817, 0.12176, 0.12486, 0.12748, 0.12964, 0.13138, 0.13274,
0.13376, 0.1346, 0.13499, 0.13528, 0.13542, 0.13544, 0.13538,
0.13527, 0.13513, 0.13497, 0.13479, 0.13462, 0.13444, 0.13427
)), .Names = c("age.min", "age.max", "male", "l", "m", "s"), class = "data.frame", row.names = c(NA,
-76L))
至於為何要將原來的名稱 bmisds() 轉名為 sds.cal() ，是有原因的。事實上在兒科研究並不是只有 BMI 要計 SDS ，血壓也要用同樣的方法計算 SDS 。另外，這些 standard 也常常轉變的。現在只要再入這些 LMS 標準便可以換算不同的 SDS 了。