使用purrr package學functional programming的觀念(ㄧ):簡介

purrr package 是Hadley Wickham和Lionel Henry所開發的,補強R base中原本的fucntional programming tool,帶入其他語言的特性,整體的風格是使用underscore.js, lodashlazy.js的語法。

R base語法中本身就有很多函數帶有部分Functional Programming意味的風格,如any,which,do.call,apply, sapply, vapply,lapply等等,只是沒有好好的統一,Hadley Wickham恰好把這部分補足了, 且裡頭的函數都是所謂的type-stable pure function導向,正如其package的名稱purrrr。

functional programming 其實是奠基在lambda calulus的觀念上,專注在四個基石上:資料(numbers, strings)、變數(函數的參數)、函數、函數應用(如何用函數修飾函數),概念上有點像是把Object-oriented programming中,function變成一個object,且專注於function功能的實現,用function去處理另一個function,去減少loop傳統寫法的使用,和重複性argument使用的麻煩。

基本上,purrr package中的函數便是依照許多其他functional programming導向的語言如scala、erlang、lisp等等特性來做的,簡單有這幾類的函數:

Category    Related Function    
Map      map(), map_chr(), map_int(), map_dbl(), map_df(), walk()
Reduce      reduce(), reduce_right()
Search      contains(), detect(), detect_index()
Filter   keep(), discrad(), every(), some()

而觀念上的一些使用如recursion寫法、expression的應用等,都可以簡化和提高代碼的效果。

這整個項目的特性:
1. 將函數輸入輸出的特性統一且分開,相對於原生的函數,常常因為輸入的資料格式不同,也會輸出不同的結果,這某些方面的“聰明”,在purrr中都簡潔化,但相對的就會產生很多相同功能,但對應不同資料格式的函數如map_chr(), map_dbl()等。
2. 將pipe %>% 的適用性導入purrr中
3.  改善function error handling的機制(使用safely())

這邊收集相關的發佈資料,能更理解這purrr項目相關進展的脈絡!(把purrr package稱為項目,是覺得這package的開發很有野心,想要把R語言在往可以適應於分散式運算系統,且整個代碼也是在一整個tidyverse下的邏輯開發的,很有一貫性)
purrr 0.1.0 2015/9/29

purrr 0.2.0 2016/1/06

畫3D scatter plot的工具:plotly, scatterplot3d

3D scatterplot的使用可以用來顯示如Principle Component Analysis的結果,目前在R中可以畫3D圖形的package其實不多,scatterplot3dplotly是唯二可以做這件事情的,scatterplot3d可以做出傳統使用R lattice所展現出來的繪圖功能,基本功能完善,但無法即時動態的去旋轉圖,這其實是畫3D Scatterplot最想做的事情,而這部分plotly提供非常好的R接口可以完成,plotly畫出來的圖形是可以即時去調整他的角度和一些顯示的細節,在探索一些三維下有無cluster現象時,真的很方便!

這邊附上官網的範例,非常簡單上手

library(plotly)

mtcars$am[which(mtcars$am == 0)] <- 'Automatic'
mtcars$am[which(mtcars$am == 1)] <- 'Manual'
mtcars$am <- as.factor(mtcars$am)

p <- plot_ly(mtcars, x = ~wt, y = ~hp, z = ~qsec, color = ~am, colors = c('#BF382A', '#0C4B8E')) %>%
  add_markers() %>%
  layout(scene = list(xaxis = list(title = 'Weight'),
                     yaxis = list(title = 'Gross horsepower'),
                     zaxis = list(title = '1/4 mile time')))

幾乎語法上跟ggplot2類似,只是在做輸出時候的layout會嚴謹一點
screenshot.png

使用Rstudio IDE將代碼分段:Code folding and Section

有時候當R代碼越來越長的時候,上下往返,變得很“髒亂”,當然除了本身代碼的簡潔度可以提昇外,可以藉由RStudio IDE的功能,來把一些代碼折疊起來,整個版面會變得清爽許多。

RStudio在處理代碼折疊顯示上,分為 1).自動判斷的折疊, 2).手動分段 這兩種方法,我覺得搭配RStudio IDE code outline的功能,可以讓工作起來更“reproducible”和“tidy”。

在“自動折疊”這塊,RStudio會將下列幾種顯示成可折疊的樣子:
1. 函數的區域內的代碼(function braced region)
2. 在R markdown中代碼區塊,且可以自己標注
3. 在R markdown中非代碼區塊,markdown的標題

另一方面,在可以自己標註代碼區塊的部分,則可以用下面的註釋方式,讓RStudio IDE來辨識,這下面的三種寫法,RStudio都可以自動辨識

# Section One ---------------------------------
# Section Two =================================
###Section Three #############################

 

最後搭配RStudio IDE的outline功能,可以把每次很複雜的分析流程,用outline架起來,顯得比較整齊,和高效!

screenshot.png

 圖左邊,下面即是把代碼折疊起來,整體清爽乾淨,右邊則是panel上的outline功能,所以可以清晰的看清楚整個代碼的架構!

Medscape One-On-One 醫療領域大咖的對話集

Medscape One-on-One是由Eric J. Topol 這位美國知名心臟內科醫師科學家所主持的對話集,Eric Topol是美國Scripps Translational Science Institute的負責人,最近剛獲得NIH一筆百萬美元的研究經費,在美國轉意醫療中佔有非常重要的角色,他的twitter帳號目前也是生命醫學領域影響力前一百名的,總是會轉貼最新的重要醫學發表。

目前Medscape One-on-One每集大概12-18分鐘左右,都會邀請各領域重量級跟醫療相關的人來闡述他們的過去經驗和看法,Eric Topol所邀來的有教授、醫師、創業家等等,每個人都對於改變醫療領域風貌有所貢獻,也都各各個很不一樣。可以直接在網站上看視頻,或是用ipad訂閱podcast,都是很棒的資訊來源。
screenshot.png

目前看到他訪問的來賓名單,每個人的經歷都很不一樣(他挑選的眼光似乎不是所謂golden 履歷,而是特殊貢獻)
Svante Paabo 分子考古學大師,追尋尼安德塔人這本書就是在講他的故事
Nazneen Rahman 英國 the Cancer Genetics Clinical Unit at The Royal Marsden NHS Foundation Trus 的負責人,專精於癌症研究

Malcolm Gladwell 知名暢銷書作者,<異數>是他其中一本膾炙人口的作品,跟他談論關於醫療費用問題,很精彩的對話

screenshot.png

Eric Schadt 計算生物學家

……. 23andme創辦人、2020 cancer moonshot 等等

將R代碼提速:使用foreach package做併行運算(parellel computing)

當使用R做計算,常常以五萬次為單位的時候,就會開始有點不耐煩,在思考如何提高自己代碼的速度,通常這種計算都會“慢”在某個loop中,通常就是做個t.test或是wilcoxon.test之類的。此時就可以考慮把代碼稍微修改,加入並行(parellel)的概念於其中,這時候可以參考使用foreach這個apckage,(在計算量小的時候,就不用牛刀小試,反而會變慢,因為一開始程序會先分推等等並行前的處置,且最後你還會需要寫代碼把資料合起來成自己想要的格式,通常都是把list變成一個data.frame)。
當然傳統基本的apply, lapply, sapply, eapply, mapply, rapply等都可以用來將這類的任務解決,但本質上這些函數依舊是單核單線程在運算的,而這邊介紹的foreach package則可以支援多核以及在叢集電腦上的多node計算。
在foreach裡面的函數,多為binary operator的形式,最簡單的使用方式,是直接將其替代原本的loop寫法:
#最基本的用法就如同for loop
x <- foreach ( i = 1:3) %do% sqrt(i)</div>
#其支援多個variable同時來迭代
x <- foreach(a = 1:3, b=rep(10, 3)) %do% ( a + b )
#可以使用curly brace, 裡頭可以包更複雜的運算
x <- foreach(a= 1:3, b= rep(10 ,3)) %do% {</div>
a + b
}

#使用.combine 選項,可以將運算的複雜度提高
x <- foreach(i = 1:3, .combine=‘c’) %do% exp(i)
x <- foreach( i = 1:3, .combine=‘cbind’) %do% rnorm(4)</div>
#cfun為使用自訂函數
x <- foreach( i = 1:3, .combine=‘cfun’, .multicombine=TRUE) %do% rnorm(4)

基本上使用%do%可以在loop上更多元的操作,而foreach最強大的其實是他並行的運算指令,其實就把%do%改成%dopar%就可以了,不過其實foreach本身是有一些缺欠的,比如他內部跟操作系統中線程的處理,是使用fork而非exec,所以在並行中在指令中間暫時性的variable時,就會出現一些問題,這部分會有另一個package來幫其補齊不足(doParellel package)

PS:當combine pipe使用時,有可能會造成ram需求大升,造成問題!

經典閱讀:Computer Age Statistical Inference: Algorithm, Evidence, And Data Science

screenshot.png

越來越發現複雜的學科,其實很難單純有選修一兩門課來建立清晰明白的思路,反而是必須依靠自己的廣泛閱讀,而像統計學這門“很接地氣”的學科,常常會因所解決領域的問題有些許不同,更是容易在自我學習中迷失方向,所以後來找一兩本非常經典的書來閱讀,才會是系統性建立知識的好習慣,這本書Computer Age Statistical Inference: Algorithms, Evidence, And Data Science熱騰疼剛完成,由劍橋出版社出版,將會是我這半年希望可以努力啃完的好書,從第一面開始就能感受到那種字裡行間的智慧,很令人享受,尤其是讀這種天才教授們所精雕的作品,細膩的科學真的不是文鄒鄒的,讀起來不好理解其實是作者的問題,希望可以有這種功力將複雜的事物細細撥開的能力,並且講述給別人理解!

這是一本由史丹佛統計學家Bradley EfronTrevor Hastie所寫的,兩位統計學家都是如神級般的存在,很多目前常用的統計方法是由他們所提出的,看這兩位大師的履歷,佩服他們努力地把統計方法一步步應用到各式生活裡的資料,並且不斷改善。在這本書中,他們回朔一百年來統計上各種觀念的改變,然後點出電腦計算將如何的改變統計方法,且開啟新的紀元(機器學習、人工智能等等)
“…the role of electronic computation is central to our story. This doesn’t mean that every advance was computer-related. A land bridge had opened up to a new continent but not all were eager to cross.”
這段話很棒,資訊科技的不斷進步,某種新大陸似乎可以藉由搭者這座“計算機”之橋通往,但不是所有人都願意往前走!
書中的內容很系統性的(不是傳統教科書那種系統性,是很“簡潔”且有歷史意義的方式,並且想要傳達出各種方法是屬於何種脈絡下的產物)整理前半世紀所發展的方法:Empirical Bayes, James-Stein estimation. bootsrap, poportional hazrd model, large-scale hypothesis testing, machine learning algorithms (the result of “cross bridge”),另外書中也展開一場很細膩的探索,在frequentist , Bayesian, fisherian三個不同系統中的inference概念,比較其間各自的優點和缺點。
本書主要分成三個大方向:
  1. 傳統經典的統計檢定
  2. 早期電腦時代所衍生的方法
  3. 新時代大資料量下的檢定方式
本書大部分內容皆涵蓋非常高知識密度的東西,蠻期待閱讀的!
話說感覺史丹佛大學統計系非常可怕,擠滿者大師,最近R世界裡面開疆闢土出很多經典package的Hadley Wickham 也要去史丹佛大學教書了,很令人興奮的一個地方。
screenshot.png

Rstudio IDE 1.0可以提升效率的功能

這個月Rstudio剛發布其1.0版本,除了在一些大家正視的熱點如專案管理、封包搭建和notebook互動式coding等等外,其實Rstudio IDE還有做了許多貼心可以改善寫代碼或是分析時候的功能

Tearable Panes

可以拖移的代碼控制板,之前不能的時候,就無法發揮雙螢幕的功效,但現在可以多個代碼穿就顯得非常無痛和舒服了!要是有可以上下或是左右split的功能,就太棒了!tip_tearable.gif

Command History

代碼歷史記錄回朔功能,這部分之前必須用鍵盤上鍵來實現,而Cmd + ⇑ 就可以一次閱覽多次的歷史輸入

tip_console.gif

 

History Panel

代碼歷史控制板,提供對於歷史代碼輸入的搜尋功能tip_history.gif

 

Rename in Scope

提供一次更改整個代碼中某一變數的能力(Rename in Scope),這功能非常強大,也是….一直很痛的問題,還在想說有什麼方法可以改善
tip_rename_in_scope.gif

Jump to Function Definition

F2快捷鍵直接跳進function中的代碼,這功能很棒,很多時候就是…不停的在修各種統計函數的細節
tip_goto_func.gif

Code Snippets

代碼的snippets設定,這功能在大部分其他語言IDE也有,算是可以大大提升速度的功能tip_code_snippet.gif

這邊有 Rstudio release note!其實他們做的gif真的蠻清楚的,一眼就知道可以幹麻!

多重比較的檢定處理 Correction for Multiple Comparisons

在處理基因表現差異的時候,通常都是要比較兩個生物樣品間,幾萬個基因有哪些是有表現差異的,此時就會同時做多次檢定,此時就會遇到一個很有趣的問題!

因為當我們在針對單一個基因在實驗設計下,是否有表現差異時,通常使用\alpha來判斷,也就是控制所謂的Type I Error(不應該reject NULL hypothesis但卻reject了),而我們常常取所謂的p value < 0.05 % 來代表統計上有顯著意義,其實是說因為“搞錯的機率小於5%,所以應該沒問題的概念,簡化為猜一百次,可能有五次猜錯”,而當我們在同一個樣本中,同時比了上萬次,那就代表我們必定會挑出很多p value有顯著,但其實是“搞錯了”的基因,此類問題其實就是所謂的Multiple Comparison 問題

那我們如何評估經過多重比較後,我們必定會搞錯一次的機率呢?

此時就可以使用Family-wise error rate,概念上其實很簡單,把每個檢體中基因間顯不顯著,當作是獨立事件,而其p值,就代表者“搞錯的機率”,假設我們有三個基因A,B,C,同時做檢定,分別得到0.03, 0.04, 0.03,以A基因為例(假設null hypothesis為兩個基因沒有表現差異),當我們要推翻這個假設時,可能有0.03的機率可能無法推翻,換句話說,我們有1-0.03的機率篤定他是沒問題,可以大方推翻檢定的,所以當我們把每個基因間當作獨立事件來看,我們一定沒錯的機率大概是:
(1-0.03)(1-0.04)(1-0.03)

反之,我們一定至少有一個基因,其實沒有顯著差異的機率(在這也就是Family-wise error rate):

$latex FWER = 1 –   (1-0.03)(1-0.04)(1-0.03) $

所以當我們可以評估多重比較後,出錯的機率後,下一步當來就是來調整原本單個基因是否表現差異,只用p = 0.05就顯著的策略,改成隨者我們同時比較的基因數量越多,我們要越嚴格一些,根據 FWER的概念,所衍生出來的調整策略,主要有下面三種:

1.Sid´ak procedure
此方法是由比較直覺的方式來處理多重檢定的問題


\alpha_{e} = 1 - (1 - \alpha_{c} )^{R}
$latex\alpha_{c} = 1 – \sqrt[R]{1-\alpha_{e}} $
$latex\alpha_{c},the probability of a Type I error at the gene level$
$latex\alpha_{e},the probability of a Type I error at the experiment error$

2.Bonferroni
Bonferron的方式主要是建立在Sidak correction上,只是其簡化原本sidak correction中的(1-p)^{R} 的展開:
\alpha_{e}= 1 - (1- \alpha_{c})^{R} =1 - (1-R*\alpha_{c} + ... )  \approx R*\alpha_{c}
因此可以簡化為下面的算式
\alpha_{c} = 1 - \frac{\alpha_{e}}{R}
3.Holm’s method(Stepwise Algorithm)
\widetilde{p}_{(1)} \leq \widetilde{p}_{(2)} \leq ... \leq \widetilde{p}_{(N)}
\widetilde{p}_{(i)} = \max \begin{Bmatrix}1 - ( 1 - p_{(j)})^{N - j + 1}\end{Bmatrix}
但可以發現,這種調整方式,有點暴力,大量基因多重檢定的時候,要顯著的p值會順間變得很嚴格。
相對於FWER,因為其是假設所有基因間是獨立關係而來的數值,但實際上生物體內基因間,多少都會有相關連,所以目前另一種用來作為多重比較的指標是False Discovery Rate(FDR)

跟其衍生而來的方法則有:
1.Benjamini & Hochberg
2.Benjamini & Yekutieli

在R裡面可以簡單使用

p.adjust(p, method = p.adjust.methods, n = length(p))

p.adjust.methods
# c("holm", "hochberg", "hommel", "bonferroni", "BH", "BY",
# "fdr", "none")

在處理像是pathway analysis或是microarray analysis時,使用Bonferroni, Sidak and Holm’s corrections顯得太過保守且沒考慮到生物上的意涵。目前認為使用FDR(False Discovery Rate)的方式相對恰當,以此觀念出發其實也發展出很多FDR的變化,每個都有些微不同的假設,包括Benjamini-Hochberg correction、Benjamini-Yekutieli correction, positive FDR、adaptive Benjamini-Hochberg correction、significant analysis of microarrays, Storey’s q-value, resampling based approach.這之中以Benjamini-Yekutieli方法最為保守。

關於這方面的知識,推薦這本書Bradley Efron, Large-Scale Inference: Empirical Bayes Methods for Estimation, Testing and Predition, Stanford University

Rstudio 開發的shiny基本架構

首先, shiny的基本架構就是由:server和ui兩部分的代碼組成。
ui部分的代碼,shiny會把R的代碼渲染成html,而server端的代碼則會處理運算和reactivity的處理。

在一個最簡單app的template,如下:

library(shiny)

#UI部分,會生成html代碼
ui <- fluidPage()

#server部分,後端進行運算和reactivity處理
server <- function(input, output){

}
shinyApp(ui = ui, server = server)

在ui部分的代碼最主要的概念是分為input和output兩類,這兩類在shiny中分別都有相對應的函數可以使用

Input相關的函數,分別對應其所輸入的資料類型和方式

screenshot.png

Output相關的函數,分別對應其所處理的資料種類screenshot.png

在server端,則有三個原則要遵守:

– 每個object使用output$..的方式來指定
– 每個object的創建使用reader*的函數
– 在object創建過程中,要使用從ui輸入的input時,可以用input$nameID來使用

在撰寫一個shiny app 便可以依據上面的架構,依序完成screenshot.png

從別人的錯誤中學習:PubPeer/Retraction Watch

從失敗或是錯誤中學習通常是比從優秀的例子中學習來得實際,畢竟在當下商業或是科學新聞媒體的風潮下(社會新聞除外,腥羶色是這類新聞的篩選標準),報導“優秀事蹟”才是王道,那些夕陽西下,或是失敗的例子往往不會花太多篇幅去追蹤。

在科學界裡面,論文的發表通常代表者勝利,往往被追捧和談論,而論文發表最後被“退回”,這類事件卻總是在檯面下發生,但這類失敗的故事最具有學習意義,畢竟成功無法複製,但是失敗可以盡量避免。

Retraction WatchPubPeer這兩個網站的出現成為一個很好的“爭議或是退回”論文案件的資料庫。

screenshot.png

Retraction Watch在2010開始由 Ivan OranskyAdam Marcus 建立,其認為科學界在“退稿(retraction)”的過程太過的無聲無息,應該使其透明,因為這些資訊都是很有價值的。

 

 

 

 

screenshot.png

PubPeer則是由Brandon Stell所發起的,目的也是為了建立一個可以作為post-publication peer review的plateform,在裡面可以看到許多令人驚心動魄的“評論”過程,是一個可以學習的地方!