研究一下那段 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 了。