Werden wir Helden für einen Tag

How to rewrite ugly R code #1

Posted on Apr 11, 2011 by Chung-hong Chan

查實用了 R 好耐，已經寫落不少 Code 。最近有點想將這些 Code 整理一下，將常用的東西全部儲在一個叫 csx.R ((別問我 csx 是甚麼，是亂改的。或者可以叫做 CSX: Statistical eXtension 吧。)) 的檔案內。下次要用時，可以用 source() 將所有常用的 function 存入。
上篇文才說過，識 R 只代表識 R 的規則，不代表寫得出好 Code 。我正是只識 R 規則的人。以下是我從 csx 找出來的一段 code 。

# < - crew it.
bmisds <- function (age, gender, bmi){
  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)
    }
  else if (age == 2 && gender == "F")
    {
      bmif <- c(-1.403,15.58,0.07342)
    }
  else if (age == 2 && gender == "M")
    {
      bmif <- c(-2.045,15.861,0.07116)
    }
  else if (age == 3 && gender == "F")
    {
       bmif <- c(-1.612,15.39,0.0741)
    }
  else if (age == 3 && gender == "M")
    {
      bmif <- c(-2.2, 15.59,0.0729)
    }
  else if (age == 4 && gender == "F")
    {
      bmif <- c(-1.735, 15.097, 0.07983)
    }
  else if (age ==4 && gender == "M")
    {
      bmif <- c(-2.287, 15.276,0.07908)
    }
  else if (age == 5 && gender == "F")
    {
      bmif <- c(-1.808, 14.84, 0.08897)
    }
  else if (age == 5 && gender == "M")
    {
      bmif <- c(-2.328,15.072,0.08867)
    }
  else if (age == 6 && gender == "F")
    {
      bmif <- c(-1.851, 14.747,0.09948)
    }
  else if (age == 6 && gender == "M")
    {
      bmif <- c(-2.337,15.026,0.09985)
    }
  else if (age == 7 && gender == "F")
    {
      bmif <- c(-1.874, 14.839, 0.10956)
    }
  else if (age == 7 && gender == "M")
    {
      bmif <- c(-2.323, 15.149, 0.11077)
    }
  else if (age == 8 && gender == "F")
    {
      bmif <- c(-1.882, 15.101, 0.11817)
    }
  else if (age == 8 && gender == "M")
    {
      bmif <- c(-2.296, 15.423, 0.12036)
    }
  else if (age == 9 && gender == "F")
    {
      bmif <- c(-1.88, 15.497, 0.12486)
    }
  else if (age == 9 && gender == "M")
    {
      bmif <- c(-2.26, 15.798, 0.12825)
    }
  else if (age == 10 && gender == "F")
    {
      bmif <- c(-1.872, 15.993, 0.1296)
    }
  else if (age == 10 && gender == "M")
    {
      bmif <- c(-2.221, 16.232, 0.13443)
    }
  else if (age == 11 && gender == "F")
    {
      bmif <- c(-1.86, 16.557, 0.13274)
    }
  else if (age == 11 && gender == "M")
    {
      bmif <- c(-2.18, 16.702,0.13902)
    }
  else if (age == 12 && gender == "F")
    {
      bmif <- c(-1.848,17.155,0.1345)
    }
  else if (age == 12 && gender == "M")
    {
      bmif <- c(-2.139,17.193,0.1423)
    }
  else if (age == 13 && gender == "F")
    {
      bmif <- c(-1.835,17.757,0.13528)
    }
  else if (age == 13 && gender == "M")
    {
      bmif <- c(-2.101,17.697,0.14458)
    }
  else if (age == 14 && gender == "F")
    {
      bmif <- c(-1.822,18.338,0.13544)
    }
  else if (age == 14 && gender == "M")
    {
      bmif <- c(-2.064,18.205,0.14617)
    }
  else if (age == 15 && gender == "F")
    {
      bmif <- c(-1.809,18.884,0.13527)
    }
  else if (age == 15 && gender == "M")
    {
      bmif <- c(-2.03,18.704,0.14733)
    }
  else if (age >= 16 && gender == "F")
    {
      bmif < - c(-1.1797,19.393,0.13497)
    }
  else if (age >= 16 && gender == "M")
    {
      bmif < - c(-1.997,19.186,0.14826)
    }
      (((bmi/bmif[2])**bmif[1] -1)/(bmif[1]*bmif[3]))
}

這個 function 要三個 parameters ，分別是 age, gender 及 BMI ，回報的數值是 BMI Standard Deviation Score (BMI SDS)。這個 Function 幾乎是我每個有份做的研究都會用到的。
為甚麼要用 BMI SDS ？因為兒童的 BMI 標準不像成人那樣，只有單一個標準。兒童在不同成長過程，都有不同的成長標準。故此， BMI 的正常值是隨著年齡及性別轉變。故此比較 16 歲童 BMI = 16 與 8 歲童的 BMI = 16 是沒有意義的。故此，要將 BMI 轉成 SDS ，即是你的 BMI 離開年齡及性別中位數幾多個 SD ，才是有意義的比較。
但是，兒童的 BMI 分佈並不呈常態分佈，故此並不可以用傳統的 z score 計算法，即是

z = (x - mean) / sd

而要用 LMS (Lamda, Mu, Sigma) 方法。 N 年前前中大兒科梁淑芳教授曾經發表香港兒童 BMI 的正常值，就是用到 LMS 方法。每個性別年齡組別都有其 Lamda (Skewness), Mu (Median) 及 Sigma (Coefficient of variation) 值。如圖：

由 LMS 轉成 BMI SDS 的算式是：

由此可見，只是六年級的數學。難就難在要根據年齡及性別從圖表中找回 LMS 三值，再加入原來的 BMI (X) ，就可算出。
從上面的 code 可見，寫得非常爛。真想將它大改造。請留意下期。