第一節：進階驗證篩選(1)

上週，我們發現這樣還是沒有辦法「精準的」找到錯誤的筆數，也許整份資料其實只要能夠刪掉幾筆就能成立。

– 我們再使用CKD門診衛教計劃的範例資料，請在這裡下載範例資料。

dat = read.csv("validated_example.csv", header = TRUE, fileEncoding = 'CP950') 
head(dat)

##   Patient       Date MDRD.GFR Stage   WBC  RBC  HB  Hct  MCV Urea.Nitrogen
## 1    1472  2011/5/31     15.9     4 10.10 2.76 8.7 27.0 98.0            75
## 2    1472   2011/8/8     24.0     4  5.79 1.80 5.7 17.9 99.4            41
## 3    1472 2011/10/17     22.7     4  5.20 2.44 7.8 23.0 94.3            49
## 4    1472  2012/1/10     18.7     4 12.99 2.95 9.4 28.2 95.6            46
## 5    1472  2013/5/14     30.4     3  6.78 3.02 9.7 30.0 99.3            48
## 6    1472  2013/8/22     26.8     4  6.88 3.73 8.6 26.0 95.2            82
##   Creatinine Uric.Acid  Na   K Albumin
## 1        3.0       3.9 139 3.3     4.2
## 2        2.1       4.8 143 4.6     4.2
## 3        2.2       4.7 146 4.8     4.5
## 4        2.6       4.0 141 4.2     4.6
## 5        1.7       3.9 141 4.4     4.1
## 6        1.9       3.0 138 4.1     4.3

– 我們先進行一個前處理，刪掉日期不符的變項

dat$Date = as.Date(dat[,"Date"])
dat = dat[!dat$Date %in% NA,]
levels.Patient = levels(as.factor(dat$Patient))

第一節：進階驗證篩選(2)

– 首先，我們還是先創造一個新變項

dat$Wrong.Date_interval = NA

接著，我們先處理一個人。先叫出資料來，但其實我們並不一定要先把這個人的資料存出來，可以透過索引函數直接作業，這樣也可以直接把Wrong.Date_interval填完

– 我們選擇第5個人

i = 5
dat[dat$Patient==levels.Patient[i],c("Patient", "Date", "Wrong.Date_interval")]

##    Patient       Date Wrong.Date_interval
## 41     566 2008-03-21                  NA
## 42     566 2008-05-17                  NA
## 43     566 2008-06-10                  NA
## 44     566 2008-06-16                  NA
## 46     566 2008-07-12                  NA

在現在的規則中，第一筆絕對是正確的，故直接先將第一筆填入FALSE，剩下的再想辦法

i = 5
dat[dat$Patient==levels.Patient[i],"Wrong.Date_interval"][1] = FALSE
dat[dat$Patient==levels.Patient[i],c("Patient", "Date", "Wrong.Date_interval")]

##    Patient       Date Wrong.Date_interval
## 41     566 2008-03-21               FALSE
## 42     566 2008-05-17                  NA
## 43     566 2008-06-10                  NA
## 44     566 2008-06-16                  NA
## 46     566 2008-07-12                  NA

第一節：進階驗證篩選(3)

透過迴圈，我們能夠從第二筆開始檢查他是否間格不夠。

– 值得注意的是，這位病患在『2008-03-21』、『2008-05-17』、『2008-06-10』這三天分別被申報，其間格分別是27與65，雖然『2008-05-17』必須被核刪，但考慮到刪除這天後，『2008-03-21』與『2008-06-10』就相距81天，這是一個可以接受的日期，所以我們必須想一下，該怎麼解決這件事。

我們可以在迴圈進行時，找尋『Wrong.Date_interval』為FALSE的最後一筆出來，之後並以他為標記做相減

– 記得避免遺漏值！

i = 5
dat[dat$Patient==levels.Patient[i],"Wrong.Date_interval"][1] = FALSE
n.date = length(dat[dat$Patient==levels.Patient[i],"Date"])

k = 2
false.dates = dat[dat$Patient==levels.Patient[i] & dat$Wrong.Date_interval == FALSE & !is.na(dat$Wrong.Date_interval),"Date"]
last.date = false.dates[length(false.dates)]
dif = dat[dat$Patient==levels.Patient[i],"Date"][k] - last.date
dat[dat$Patient==levels.Patient[i],"Wrong.Date_interval"][k] = dif < 60

dat[dat$Patient==levels.Patient[i],c("Patient", "Date", "Wrong.Date_interval")]

##    Patient       Date Wrong.Date_interval
## 41     566 2008-03-21               FALSE
## 42     566 2008-05-17                TRUE
## 43     566 2008-06-10                  NA
## 44     566 2008-06-16                  NA
## 46     566 2008-07-12                  NA

第一節：進階驗證篩選(4)

接著我們只要將i以及k打包成迴圈就完成了

– 需要注意的是，如果n.date>1才需要做後續檢查

levels.Patient = levels(as.factor(dat$Patient))
n.Patient = length(levels.Patient)
dat$Wrong.Date_interval = NA

pb = txtProgressBar(max = n.Patient, style=3)

for (i in 1:n.Patient) {
  dat[dat$Patient==levels.Patient[i],"Wrong.Date_interval"][1] = FALSE
  n.date = length(dat[dat$Patient==levels.Patient[i],"Date"])
  if (n.date>1) {
    for (k in 2:n.date) {
      false.dates = dat[dat$Patient==levels.Patient[i] & dat$Wrong.Date_interval == FALSE & !is.na(dat$Wrong.Date_interval),"Date"]
      last.date = false.dates[length(false.dates)]
      dif = dat[dat$Patient==levels.Patient[i],"Date"][k] - last.date
      dat[dat$Patient==levels.Patient[i],"Wrong.Date_interval"][k] = dif < 60
    }
  }
  setTxtProgressBar(pb, i)
}

close(pb)

練習1：把剛剛學到的推廣到GFR的變化量

GFR在短時間內有劇烈變化很有可能是因為住院的關係，因此該筆紀錄是可以被刪除的。

– 因此請透過類似的方法找出GFR的異常值，如果出現短時間內變化過大的值請把它找出來（單月斜率超過2以上），之後在檢查時請忽略它！

一定要使用進度條函數才能掌握你的進度喔！

練習1答案

這題沒有想像中簡單，需要注意錯誤訊息：

#Read data
dat = read.csv("validated_example.csv", header = TRUE, fileEncoding = 'CP950') 

#Rule 1: check date-format
dat$Date = as.Date(dat[,"Date"])
dat$Wrong.Date = is.na(dat$Date)
dat = dat[dat$Wrong.Date == FALSE,]

#Rule 2: check eGFR change
levels.Patient = levels(as.factor(dat$Patient))
n.Patient = length(levels.Patient)
dat$Wrong.eGFR_change = NA

pb = txtProgressBar(max = n.Patient, style=3)

for (i in 1:n.Patient) {
  dat[dat$Patient==levels.Patient[i], "Wrong.eGFR_change"][1] = FALSE
  n.date = length(dat[dat$Patient==levels.Patient[i],"Date"])
  if (n.date>1) {
    for (k in 2:n.date) {
      false.dates = dat[dat$Patient==levels.Patient[i] & dat$Wrong.eGFR_change == FALSE & !is.na(dat$Wrong.eGFR_change),"Date"]
      last.date = false.dates[length(false.dates)]
      diff_date = dat[dat$Patient==levels.Patient[i],"Date"][k] - last.date
      diff_date = as.numeric(diff_date, units = 'days') # 注意要轉成數字，否則無法相除

      false.eGFRs = dat[dat$Patient==levels.Patient[i] & dat$Wrong.eGFR_change == FALSE & !is.na(dat$Wrong.eGFR_change),"MDRD.GFR"]
      last.eGFR = false.eGFRs[length(false.eGFRs)]
      diff_eGFR = dat[dat$Patient==levels.Patient[i],"MDRD.GFR"][k] - last.eGFR
      
      slope = diff_eGFR / diff_date * 30
      
      dat[dat$Patient==levels.Patient[i],"Wrong.eGFR_change"][k] = abs(slope) > 2
    }
  }
  setTxtProgressBar(pb, i)
}

close(pb)

第二節：快速讀寫檔案(1)

在這節課之前，我們使用的函數種類相當有限，但我們也在之前的範例中發現了不同程式碼執行同樣的事情在速度上有顯著的不同。

– 不同的主因其實牽涉到程式碼執行過程，但目前我們的知識很有限，做出來就不容易了更別說還要求速度快，因此比較可行的方式是去「抄」程式碼！

在讀取檔案的部分，我們這裡先介紹一個特別的套件：data.table，它裡面有眾多協助我們加速處理大量資料的函數

第二節：快速讀寫檔案(2)

– 至於安裝套件的方法，假設你是用Rstudio，可以看到右下角有個Packages分頁，點選後你能看到Install按鍵，透過這種方式就能安裝套件了

F01

假設你是使用一般的R程式，可以在Console地方打上這個指令：

install.packages("data.table")

安裝完成後，未來我們要使用這個套件只要打上這串指令：

library(data.table)

第二節：快速讀寫檔案(3)

讓我們來體驗一下data.table套件內一個神奇的函數：「fread」，它可以以數倍的速度讀取一個csv檔案：

– 至於大Data在哪，我們可以使用第六節課用到的大檔案，請在這裡下載。

這是使用「read.csv」讀取的速度：

t0 = Sys.time()
dat1 = read.csv('laboratory_2.csv', header = TRUE, fileEncoding = 'CP950')
Sys.time() - t0

## Time difference of 1.707349 secs

這是使用「fread」讀取的速度：

t0 = Sys.time()
dat2 = fread('laboratory_2.csv', header = TRUE)
Sys.time() - t0

## Time difference of 0.07944632 secs

第二節：快速讀寫檔案(4)

比較討厭的地方在使用「fread」讀取進來的物件格式是「data.table」，這是一個新的物件格式。

class(dat1)

## [1] "data.frame"

class(dat2)

## [1] "data.table" "data.frame"

解決這個問題的方式是在使用「fread」的時候多給一個參數：

t0 = Sys.time()
dat2 = fread('laboratory_2.csv', header = TRUE, data.table = FALSE)
Sys.time() - t0

## Time difference of 0.1181285 secs

class(dat2)

## [1] "data.frame"

第二節：快速讀寫檔案(5)

但現在兩者還不完全一樣，問題出在哪呢?

all.equal(dat1, dat2)

## [1] "Component \"COLLECTIONDATE\": 'current' is not a factor"
## [2] "Component \"TESTNAME\": 'current' is not a factor"      
## [3] "Component \"UNITS\": 'current' is not a factor"

原來是使用「read.csv」有另一個參數「stringsAsFactors」，預設是TRUE，但使用「fread」的時候預設卻是FALSE，所以我們再改變一下：

dat1 = read.csv('laboratory_2.csv', header = TRUE, stringsAsFactors = FALSE, fileEncoding = 'CP950')

除此之外，還有encoding的問題，但這個我們必須透過手動修正「laboratory_2.csv」編碼解決。

all.equal(dat1, dat2)

## [1] "Component \"COLLECTIONDATE\": 181483 string mismatches"

第二節：快速讀寫檔案(6)

除了讀取檔案之外，寫出檔案也可以使用「fwrite」函數進行加速：

– 這是使用「write.csv」寫出的速度：

t0 = Sys.time()
write.csv(dat1, 'laboratory_test.csv', row.names = FALSE, quote = TRUE)
Sys.time() - t0

## Time difference of 0.9308543 secs

– 這是使用「fwrite」寫出的速度：

t0 = Sys.time()
fwrite(dat2, 'laboratory_test.csv', row.names = FALSE, quote = TRUE)
Sys.time() - t0

## Time difference of 0.0409317 secs

第二節：快速讀寫檔案(7)

另外還有一種儲存、載入物件的方式，使用R內建的「save」與「load」函數甚至可以完全保留物件原有的所有屬性，並且能夠把任何物件存出，也具有較小的儲存空間：

– 這是儲存

t0 = Sys.time()
save(dat1, file = 'laboratory_test.RData')
Sys.time() - t0

## Time difference of 0.6009369 secs

– 這是載入

t0 = Sys.time()
load('laboratory_test.RData')
Sys.time() - t0

## Time difference of 0.204529 secs

練習2：讀取圖片及顯示圖片

學會運用套件後，我們了解到R可以透過套件的擴充它的性能，在此之前我們完全沒辦法想像要怎樣把圖片讀進R裡面，現在請你隨便找一張圖像檔案，並試著從網路上找找看要怎樣讀檔以及顯示圖片！

– 透過Google搜尋「R display image」後，你將可以找到這個頁面，其中第一個連結進去後你會發現這裡已經有人發問和回答了：

F02

請試著利用它所提供的方式讀取及顯示圖片吧！

練習2答案

你應該會找到，關鍵的套件是「jpeg」，而裡面的函數「readJPEG」可以把圖讀進來：

library("jpeg")
img <- readJPEG(system.file("img", "Rlogo.jpg", package = "jpeg"))
plot(0:1, 0:1, type = "n", ann = FALSE, axes = FALSE)
rasterImage(img, -0.04, -0.04, 1.04, 1.04)

第三節：增加程式可讀性(1)

在R裡面我們經常會使用到組合函數，像是假使我們想要知道生化值的數量共有幾項，我們可以透過下面這種方法得知：

length(levels(factor(dat1$TESTNAME)))

## [1] 25

而這樣程式寫起來會很複雜，比較好的方式應該是這樣：

factorized_TESTNAME = factor(dat1$TESTNAME)
lvl_TESTNAME = levels(factorized_TESTNAME)
length(lvl_TESTNAME)

## [1] 25

但這樣在R裡面會額外儲存很多垃圾物件，並且會拖慢運算速度，因此如何在程式可讀性與執行速度上做平衡呢?

第三節：增加程式可讀性(2)

這裡我們介紹另一個套件「magrittr」，它擁有一個特殊的運算符號「%>%」，功能是把從左到右依序執行任務

– 後面函數的「.」代表上一步的結果

library(magrittr)

n.TESTNAME = dat1$TESTNAME %>% factor %>% levels %>% length
n.TESTNAME

## [1] 25

n.TESTNAME = dat1$TESTNAME %>% factor() %>% levels() %>% length()
n.TESTNAME

## [1] 25

n.TESTNAME = dat1$TESTNAME %>% factor(.) %>% levels(.) %>% length(.)
n.TESTNAME

## [1] 25

為什麼需要「.」這種方式呢，這是因為假使你的函數需要指令參數，那可以透過這種方式來寫出：

f = function(x, a, b) {a*x^2 + b}
1:5 %>% f(., 2, 5)

## [1]  7 13 23 37 55

1:5 %>% f(2, ., 5)

## [1]  9 13 17 21 25

1:5 %>% f(2, 5, .)

## [1] 21 22 23 24 25

第三節：增加程式可讀性(3)

我們再多介紹幾個特殊的運算符號(這兩個比較不常用)：

– 「%<>%」：不要顯示結果，而是改變物件內容

a = 1
a %<>% add(1)
a

## [1] 2

– 「%$%」：指定物件內的索引格式

n.TESTNAME = dat1 %$% TESTNAME %>% factor %>% levels %>% length
n.TESTNAME

## [1] 25

第三節：增加程式可讀性(4)

學習特殊運算符號的目標除了是增加自己程式的可讀性之外，更重要的是會增加及執行速度！
讓我們看看結合了眾多改變後，再回頭看看第六課的練習一這個任務要花多久：

t0 = Sys.time()

dat1$COLLECTIONDATE = dat1[,3] %>% as.Date

levels.TESTNAME = dat1[,4] %>% factor %>% levels
n.TESTNAME = levels.TESTNAME %>% length
levels.PATNUMBER = dat1[,1] %>% factor %>% levels
n.PATNUMBER = levels.PATNUMBER %>% length

dat_list = list()

for (i in 1:n.PATNUMBER) {
  
  subdat = dat1[dat1[,1]==levels.PATNUMBER[i],]
  levels.COLLECTIONDATE = subdat[,3] %>% factor %>% levels
  n.COLLECTIONDATE = levels.COLLECTIONDATE %>% length

  submatrix = matrix(NA, nrow = n.COLLECTIONDATE, ncol = n.TESTNAME + 2)
  colnames(submatrix) = c("PATNUMBER", "COLLECTIONDATE", levels.TESTNAME)

  submatrix[,1] = levels.PATNUMBER[i]
  submatrix[,2] = levels.COLLECTIONDATE

  for (j in 1:n.COLLECTIONDATE) {
    subsubdat = subdat[subdat[,3]==levels.COLLECTIONDATE[j],]
    for (k in 1:nrow(subsubdat)) {
      NAME = subsubdat[k,4]
      position = which(NAME == levels.TESTNAME) + 2
      submatrix[j, position] = subsubdat[k,5]
    }
  }
  
  dat_list[[i]] = submatrix
}

final.data = do.call("rbind", dat_list)

Sys.time() - t0

## Time difference of 1.657264 mins

head(final.data)

##      PATNUMBER COLLECTIONDATE Albumin Albumin body fluid AST BUN BUN Fluid
## [1,] "26"      "2011-05-12"   NA      NA                 NA  NA  NA       
## [2,] "26"      "2011-05-30"   NA      NA                 NA  NA  NA       
## [3,] "26"      "2011-05-31"   NA      NA                 NA  NA  NA       
## [4,] "26"      "2011-06-01"   NA      NA                 NA  NA  NA       
## [5,] "26"      "2011-06-02"   NA      NA                 NA  NA  NA       
## [6,] "26"      "2011-06-06"   NA      NA                 NA  NA  NA       
##      Cholesterol Fluid Creatinine Creatinine Fluid GLU(AC) HDL-Cholesterol
## [1,] NA                "1.8"      NA               NA      NA             
## [2,] NA                "3"        NA               NA      NA             
## [3,] NA                "2.9"      NA               NA      NA             
## [4,] NA                "2.9"      NA               NA      NA             
## [5,] NA                "2.4"      NA               NA      NA             
## [6,] NA                "1.9"      NA               NA      NA             
##      IP    K  LDL-Cholesterol Na    NA Fluid Total Calcium Total Cholesterol
## [1,] NA    NA NA              "140" NA       NA            NA               
## [2,] NA    NA NA              "139" NA       NA            NA               
## [3,] "4.6" NA NA              "145" NA       "7.8"         "134"            
## [4,] NA    NA NA              "144" NA       "6.1"         NA               
## [5,] NA    NA NA              "138" NA       "7.1"         NA               
## [6,] NA    NA NA              "134" NA       "8"           NA               
##      Triglyceride Triglycerol Fluid Uric Acid urine Calcium urine Phosphorus
## [1,] NA           NA                NA        NA            NA              
## [2,] NA           NA                NA        NA            NA              
## [3,] "131"        NA                NA        "5.6"         "41.4"          
## [4,] NA           NA                NA        NA            NA              
## [5,] NA           NA                NA        NA            NA              
## [6,] NA           NA                NA        NA            NA              
##      urine Potassium urine Sodium urine Uric Acid
## [1,] NA              NA           NA             
## [2,] NA              NA           NA             
## [3,] "29"            "61"         NA             
## [4,] NA              NA           NA             
## [5,] NA              NA           NA             
## [6,] NA              NA           NA

居然1分多鐘就做完了，配上讀取/寫出檔案的時間加起來也不超過1分半，看來大檔案的處理也並不怎麼可怕！

final.data = as.data.frame(final.data, stringsAsFactors = FALSE)

fwrite(final.data, 'final_data.csv', row.names = FALSE, quote = TRUE)

R語言程式設計導論