基本介紹

• 許多時候我們希望可以將我們的工作成果與其他人分享，同時也希望別人可以簡單的享用這些資源，這時候WebApp就是一個很好的媒介。

• WebApp有下列優點：

1. 只需要連上網路即可使用

2. 透過互動清楚的展示你所想表達的事物

3. 對於使用者來說非常簡單上手

• 以SIR model為例，這個網站能讓使用者在簡單決定重要參數後，畫出隨時間不同的人數變化：

R語言

• R語言可能是所有程式語言中最容易上手的，並且它支援非常多統計方法的使用，除此之外，他還是個免費的資源。

• 要開始使用R，你必須先到CRAN下載檔案，安裝後才可以使用

• 但是若是希望更快速的上手，我推薦大家先至Rstudio官網下載Rstudio，它可以使用圖形介面支援R裡面的一些基礎功能

安裝套件

• 套件(Package)是R的核心，而撰寫Webapp是使用一個叫做shiny的套件，他的介紹可以參考Shiny官網，我們可以在Rstudio裡面用下列的方式安裝這個套件

開始創造一個簡單的App - Step 1

• 在最開始的時候，我們要先建立一個新的資料夾(Myapp)，這個資料夾內以後將裝著兩個R script檔案，分別是ui.R(主管使用者介面)以及server.R(主管伺服器端的處理)。

開始創造一個簡單的App - Step 2

• 接著進入Myapp資料夾後，我們要新增一個R script檔案。

開始創造一個簡單的App - Step 3

• 接著在檔案內貼上下列程式碼，然後儲存為ui.R

library(shiny)

# Define UI for application that plots random distributions 
shinyUI(pageWithSidebar(

  # Application title
  headerPanel("Hello Shiny!"),

  # Sidebar with a slider input for number of observations
  sidebarPanel(
    sliderInput("obs", "Number of observations:", min = 0, max = 1000, value = 500)
  ),

  # Show a plot of the generated distribution
  mainPanel(
    plotOutput("distPlot")
  )
))

開始創造一個簡單的App - Step 4

• 再新增一個R script檔案，並且貼上下列程式碼，然後儲存為server.R

library(shiny)

# Define server logic required to generate and plot a random distribution
shinyServer(function(input, output) {
  
  # Expression that generates a plot of the distribution. The expression is
  # wrapped in a call to renderPlot to indicate that:
  # 
  # 1) It is 'reactive' and therefore should be automatically re-executed
  # when inputs change 2) Its output type is a plot
  output$distPlot = renderPlot({
    
    # generate an rnorm distribution and plot it
    dist = rnorm(input$obs)
    hist(dist)
  })
  
})

開始創造一個簡單的App - Step 5

• 都儲存完成後，接著點選Run App，你就可以完成你寫的App了。

• 完成的檔案如下圖所示：

練習-1

• 重複上述步驟，試著建立出你的第一個WebApp。

– 如果你想複製貼上點不同的東西，可以在shiny-examples內找到更多例子，

• 對於程度較佳的同學，可以試著改變App的標題，或是一些簡單的參數。

簡單介紹Shinyapps

• 前面已經說過，使用shiny(R package)所撰寫的App，他的基本構造是一個包含ui.R(主管使用者介面)以及server.R(主管伺服器端的處理)的資料夾。

• Shiny app的基本運作流程為：

1. 使用者自ui.R中的給定一個參數。

2. 這個參數傳到server.R裡面，使用反應函數進行計算。

3. 反應完成後，再回傳至ui.R輸出反應結果。

簡介ui.R內的參數

• 下列是剛剛你所貼上的程式碼(ui.R)

library(shiny)

# Define UI for application that plots random distributions 
shinyUI(pageWithSidebar(

  # Application title
  headerPanel("Hello Shiny!"),

  # Sidebar with a slider input for number of observations
  sidebarPanel(
    sliderInput("obs", "Number of observations:", min = 0, max = 1000, value = 500)
  ),

  # Show a plot of the generated distribution
  mainPanel(
    plotOutput("distPlot")
  )
))

• 其中是以shinyUI()這個函數為開頭，首先要先指定說需要有控制區域，因此再使用pageWithSidebar()函數，而控制區函數內必須包含3個子函數，分別是：

1. headerPanel()用來定義網頁標題

2. sidebarPanel()用來定義的控制選單內含哪些可控參數，本例中只有一個滑動輸入元件sliderInput()，元件為obs

3. mainPanle()則是用來定義輸出區域的輸出結果，本例中只有一個圖片輸出元件plotOutput()，元件為distPlot

簡介server.R內的參數

• 下列是剛剛你所貼上的程式碼(server.R)

library(shiny)

# Define server logic required to generate and plot a random distribution
shinyServer(function(input, output) {

    # Expression that generates a plot of the distribution. The expression is
    # wrapped in a call to renderPlot to indicate that:
    # 
    # 1) It is 'reactive' and therefore should be automatically re-executed
    # when inputs change 2) Its output type is a plot
    output$distPlot = renderPlot({

        # generate an rnorm distribution and plot it
        dist = rnorm(input$obs)
        hist(dist)
    })

})

• 其中是以shinyServer()這個函數為開頭，裡面包含著一個函數function()要求你指定input及output，由於我們只要做一個反應函數，由於我們想要畫一張圖，所以使用renderPlot()函數，在裡面我們指定它畫圖的過程，並且將結果儲存在output裡面的distPlot。

• 這個畫圖的過程很簡單，就是先指定一個數字（由input$obs提供，也就是剛剛在ui.R中使用者所給定的參數)，然後要求R使用rnorm()隨機產生n個平均數為0，標準差為1的數列，而這個數列儲存在dist元件內。接著在以hist(dist)畫出這個數列的直方圖。

• 因此，在renderPlot()函數內，我們根據使用者指定的參數(input$obs)產生了一張直方圖，而這張直方圖將會儲存在output裡面的distPlot。

• 接著這個物件distPlot就會回到ui.R中，而根據我們在ui.R中所下的指令，他將會使用plotOutput()函數使這張圖形呈現在輸出區上。

Note：所有的輸入元件都存取在input這個List內；而所有的輸出元件都存取在output這個List內。

學習R裡面的基本函數

• 至此為止，我們已經學會了webApp的基本寫法。然而這時候若需要進一步的增加它的應用面，我們需要的是學習更多R裡面的函數，讓我們的App可以擁有更強大的功能。

• 請在R內練習下列這段程式碼的使用，並試圖改變裡面的obs,M,S,Coler的數字，了解他的功能

        obs = 500
        M = 170
        S = 10
        Coler = "blue"
        dist = rnorm(obs,mean=M,sd=S)
        hist(dist,col=Coler)

增加可控參數至剛剛的App中

• 在剛剛的練習中我們已經清楚了在R裡面，改變參數可以改變圖形的常態分布，接著我們可以開始增加這些參數至我們剛剛寫的App內。

• 除了已經存在的obs之外，我們將繼續增加M,S,Coler在我們的控制區域之內。

• ui.R

library(shiny)

# Define UI for application that plots random distributions 
shinyUI(pageWithSidebar(

  # Application title
  headerPanel("Hello Shiny!"),

  # Sidebar with inputs for number of observations, mean, SD, and coler.
  sidebarPanel(
    sliderInput("obs", "Number of observations:", min = 0, max = 1000, value = 500),
    numericInput("M", "Mean of this normal distribution:", min = -200, max = 200, value = 100),
    numericInput("S", "SD of this normal distribution:", min = 0, max = 50, value = 10),
    radioButtons("Coler", "Select the color of histogram:", choices = c("Red" = "red", "Blue" = "blue", "Green" = "green"))
  ),

  # Show a plot of the generated distribution
  mainPanel(
    plotOutput("distPlot")
  )
))

• server.R (注意在server.R之中，所有在ui.R內所出現的參數前面都要加上“input$”，表示他其實是在input這個List之內的東西。)

library(shiny)

# Define server logic required to generate and plot a random distribution
shinyServer(function(input, output) {

    # Expression that generates a plot of the distribution. The expression is
    # wrapped in a call to renderPlot to indicate that:
    # 
    # 1) It is 'reactive' and therefore should be automatically re-executed
    # when inputs change 2) Its output type is a plot
    output$distPlot = renderPlot({

        # generate an rnorm distribution and plot it
        dist = rnorm(input$obs,mean=input$M,sd=input$S)
        hist(dist,col=input$Coler)
    })

})

學習將自己的計算工具包裝成一個自訂函數

• 我們先把R做為一個計算機來使用，並且使用函數功能做一些較複雜的運算。

• 下面這個範例，是試著計算出在進行民調抽樣時的樣本數估計，假定我們可以接受的誤差為3%，我們將可以透過下列公式計算出所需樣本數(X=可接受的誤差值；p=預期的支持率；Z=可接受之信賴區間對應之Z值)

Sig = 0.05  #two-sided significance level
p = 0.5     #The expected rate
X = 0.03    #Acceptable error range
Z = qnorm(1-Sig/2)   #The Z value based on the two-sided significance level 
n = p*(1-p)*Z^2/X^2  #Calculate the needed number of sample
n           #Print the outcome

• 接著，我們可以把它包為一個函數(請試試它的功能)

Samplesize = function(Sig,p,X) {
  Z = qnorm(1-Sig/2)
  n = p*(1-p)*Z^2/X^2
  n
}

Samplesize(Sig = 0.05,p = 0.5, X = 0.03) #Test 1
Samplesize(Sig = 0.05,p = 0.2, X = 0.03) #Test 2
Samplesize(Sig = 0.05,p = 0.2, X = 0.05) #Test 3

創造一個新的App使我們能更方便的計算樣本數

• 在R裡面我們學會函數運算之後，接著我們會希望其他人也可以輕易的使用我們剛剛設定好的函數，因此我們會希望將他寫成像這樣子的App，供大家使用。

• ui.R

library(shiny)

# Define UI for application that calculate the needed sample size 
shinyUI(pageWithSidebar(

  # Application title
  headerPanel("Sample size calculator"),

  # Sidebar with numeric inputs for significance level, expected rate, acceptable error range.
  sidebarPanel(
    numericInput("Sig", "Two-sided significance level:", min = 0, max = 1.000, value = 0.05),
    numericInput("p", "The expected rate:", min = 0, max = 1, value = 0.5),
    numericInput("X", "Acceptable error range:", min = 0, max = 1, value = 0.03)
  ),

  # Show the outcome
  mainPanel(
    h3(textOutput("Size"))
  )
))

• server.R

library(shiny)

# Your function
Samplesize = function(Sig,p,X) {
  Z = qnorm(1-Sig/2)
  n = p*(1-p)*Z^2/X^2
  n
}

# Define server logic required to calculate the needed sample size
shinyServer(function(input, output) {

    output$Size = renderText({
      N = Samplesize(Sig = input$Sig,p = input$p, X = input$X)      
      paste("We need",round(N),"samples for this sampling.")
    })

})

練習-2

• 下列公式是Case control study所需要的樣本數估計公式，如果可以的話，請試著建立函數計算Sample size，並且建立在App內

– 原始參數共有5個，分別是:

1. alpha: 雙尾顯著水準（請在運算前把它轉換為適當的Z value）

2. power: 單尾檢力（請在運算前把它轉換為適當的Z value）

3. OR: 臨床上有意義之最小OR

4. r: Case組人數是Control組人數的幾倍

5. p0: 預期Control組之暴露率

• 如果你真的無法自己寫出上述公式，那請你複製下列這串程式碼

Samplesize_OR = function(alpha,power,OR,r,p0,corr=FALSE) {
  Za=qnorm(1-alpha/2)
  Zb=qnorm(power)
  q0=1-p0
  p1=OR*p0/(OR*p0+q0)
  q1=1-p1
  pbar=(p0+r*p1)/(r+1)
  qbar=1-pbar
  n0=(Za*sqrt((r+1)*pbar*qbar)+Zb*sqrt(r*p0*q0+p1*q1))^2/(r*(p0-p1)^2)
  n1=r*n0
  n0corr=n0/4*(1+sqrt(1+(2*(1+r)/(n1*r*sqrt((p0-p1)^2)))))^2
  n1corr=r*n0corr
  if (corr==FALSE) {return(paste("Control:",round(n0)+1,"; Case:",round(n1)+1,"; Total:",round(n0)+round(n1)+2,sep=""))}
  if (corr==TRUE) {return(paste("Control:",round(n0corr)+1,"; Case:",round(n1corr)+1,"; Total:",round(n0corr)+round(n1corr)+2,sep=""))}
}

Samplesize_OR(alpha=0.05,power=0.8,OR=1.5,r=1,p0=0.2,corr=FALSE)
Samplesize_OR(alpha=0.05,power=0.8,OR=1.5,r=1,p0=0.2,corr=TRUE)

練習-2 答案

• ui.R

library(shiny)

shinyUI(pageWithSidebar(
  
  headerPanel("Sample size for Case Control study"),
  
  sidebarPanel(
    numericInput("alpha", "Two-sided significance level:", min = 0.001, max = 0.999, value = 0.05),
    numericInput("power", "Power:", min = 0.001, max = 0.999, value = 0.8),
    numericInput("OR", "Smallest difference of clinical/biological importance:", min = 0.01, max = 100, value = 1.5),
    numericInput("r", "The ratio of Case/Control:", min = 0.01, max = 100, value = 1),
    numericInput("p0", "Proportion of controls with exposures:", min = 0.001, max = 0.999, value = 0.2)    
  ),
  
  mainPanel(
    h3(textOutput("text1")),
    h4(textOutput("text2")),
    h3(textOutput("text3")),
    h4(textOutput("text4"))
  )
))

• server.R

library(shiny)

Samplesize_OR = function(alpha,power,OR,r,p0,corr=FALSE) {
  Za=qnorm(1-alpha/2)
  Zb=qnorm(power)
  q0=1-p0
  p1=OR*p0/(OR*p0+q0)
  q1=1-p1
  pbar=(p0+r*p1)/(r+1)
  qbar=1-pbar
  n0=(Za*sqrt((r+1)*pbar*qbar)+Zb*sqrt(r*p0*q0+p1*q1))^2/(r*(p0-p1)^2)
  n1=r*n0
  n0corr=n0/4*(1+sqrt(1+(2*(1+r)/(n1*r*sqrt((p0-p1)^2)))))^2
  n1corr=r*n0corr
  if (corr==FALSE) {return(paste("Control:",round(n0)+1,"; Case:",round(n1)+1,"; Total:",round(n0)+round(n1)+2,sep=""))}
  if (corr==TRUE) {return(paste("Control:",round(n0corr)+1,"; Case:",round(n1corr)+1,"; Total:",round(n0corr)+round(n1corr)+2,sep=""))}
}

shinyServer(function(input, output) {
  
  output$text1 = renderText({
    "Sample size without continuity correction:"
  })
  
  output$text2 = renderText({
    Samplesize_OR(alpha=input$alpha,power=input$power,OR=input$OR,r=input$r,p0=input$p0)
  })
  
  output$text3 = renderText({
    "Sample size with continuity correction:"
  })
  
  output$text4 = renderText({
    Samplesize_OR(alpha=input$alpha,power=input$power,OR=input$OR,r=input$r,p0=input$p0,corr=TRUE)
  })
  
})

將兩種Sample size的計算工具合成同個App

• 如果可以的話，我們會更希望兩種計算工具可以合在同個App內，那我們在ui.R中會需要一個新的函數，叫做conditionalPanel()

– conditionalPanel()可以讓我們的控制bar僅在特定條件下出現

• 這邊有個範例，設計一個按鈕選擇我要產生常態分布還是t分佈

– 常態分布必須要有兩個參數(平均數&標準差)決定，而t分布則必須要由一個參數決定(自由度)，因此我們就必須設計當選擇常態分不時，控制區要出現的是平均數&標準差；若選擇了t分布，控制區要出現的則是自由度

• ui.R

library(shiny)
shinyUI(pageWithSidebar(
  
  headerPanel("Generating distribution"),
  
  sidebarPanel(
    selectInput("method", "Choose distribution:", choices=c("Normal"="norm", "Student t"="st")),
    helpText("Setting parameter(s) for distribution model"),
    conditionalPanel(condition="input.method=='norm'",
                     numericInput(inputId="mu", label="mean", value=0),
                     numericInput(inputId="sd", label="standard deviation", value=1, min=0)
    ),
    conditionalPanel(condition="input.method=='st'",
                     numericInput(inputId="df", label="Df", value=10, min=1)
    ),
    sliderInput(inputId="obs", 
                label="Number of observations:", 
                min = 1, max = 1000, value = 500)
  ),
  
  mainPanel(
    plotOutput("distPlot")
  )

))

• server.R

library(shiny)

shinyServer(function(input, output) {
    output$distPlot <- renderPlot({
        if (input$method == "norm") 
            {dist <- rnorm(input$obs, mean = input$mu, sd = input$sd)}
        if (input$method == "st") 
            {dist <- rt(input$obs, df = input$df)}
        hist(dist)
    })
})

練習-3

• 請試著使用conditionalPanel()，將剛剛練習中的兩種樣本數估計的方式寫在同個App之內，並讓使用者自行選擇。

練習-3 答案

• ui.R

library(shiny)

shinyUI(pageWithSidebar(
  
  headerPanel("Sample size calculator"),
  
  sidebarPanel(
    selectInput("method", "Choose a method:", choices=c("Rate"="rate", "Case Control Study"="cc")),
    conditionalPanel(condition="input.method=='rate'",
    numericInput("Sig", "Two-sided significance level:", min = 0, max = 1.000, value = 0.05),
    numericInput("p", "The expected rate:", min = 0, max = 1, value = 0.5),
    numericInput("X", "Acceptable error range:", min = 0, max = 1, value = 0.03)
    ),
    conditionalPanel(condition="input.method=='cc'",
    numericInput("alpha", "Two-sided significance level:", min = 0.001, max = 0.999, value = 0.05),
    numericInput("power", "Power:", min = 0.001, max = 0.999, value = 0.8),
    numericInput("OR", "Smallest difference of clinical/biological importance:", min = 0.01, max = 100, value = 1.5),
    numericInput("r", "The ratio of Case/Control:", min = 0.01, max = 100, value = 1),
    numericInput("p0", "Proportion of controls with exposures:", min = 0.001, max = 0.999, value = 0.2) 
    )  
  ),
  
  mainPanel(
    h3(textOutput("text1")),
    h4(textOutput("text2")),
    h3(textOutput("text3")),
    h4(textOutput("text4"))
  )
))

• server.R

library(shiny)

Samplesize = function(Sig,p,X) {
  Z = qnorm(1-Sig/2)
  n = p*(1-p)*Z^2/X^2
  n
}

Samplesize_OR = function(alpha,power,OR,r,p0,corr=FALSE) {
  Za=qnorm(1-alpha/2)
  Zb=qnorm(power)
  q0=1-p0
  p1=OR*p0/(OR*p0+q0)
  q1=1-p1
  pbar=(p0+r*p1)/(r+1)
  qbar=1-pbar
  n0=(Za*sqrt((r+1)*pbar*qbar)+Zb*sqrt(r*p0*q0+p1*q1))^2/(r*(p0-p1)^2)
  n1=r*n0
  n0corr=n0/4*(1+sqrt(1+(2*(1+r)/(n1*r*sqrt((p0-p1)^2)))))^2
  n1corr=r*n0corr
  if (corr==FALSE) {return(paste("Control:",round(n0)+1,"; Case:",round(n1)+1,"; Total:",round(n0)+round(n1)+2,sep=""))}
  if (corr==TRUE) {return(paste("Control:",round(n0corr)+1,"; Case:",round(n1corr)+1,"; Total:",round(n0corr)+round(n1corr)+2,sep=""))}
}

shinyServer(function(input, output) {
  
  output$text1 = renderText({
    if (input$method=="rate") {
      N = Samplesize(Sig = input$Sig,p = input$p, X = input$X)      
      return(paste("We need",round(N),"samples for this sampling."))
     }
    if (input$method=="cc") {"Sample size without continuity correction:"}
  })
  
  output$text2 = renderText({
    if (input$method=="rate") {return(NULL)}
    if (input$method=="cc") {Samplesize_OR(alpha=input$alpha,power=input$power,OR=input$OR,r=input$r,p0=input$p0)}
  })
  
  output$text3 = renderText({
    if (input$method=="rate") {return(NULL)}
    if (input$method=="cc") {"Sample size with continuity correction:"}
  })
  
  output$text4 = renderText({
    if (input$method=="rate") {return(NULL)}
    if (input$method=="cc") {Samplesize_OR(alpha=input$alpha,power=input$power,OR=input$OR,r=input$r,p0=input$p0,corr=TRUE)}
  })
  
})

建立產生SIR model的App

• SIR model是將整個族群分成三種人，其中S為易感受個案，I為受感染個案(具感染力)，R為已康復個案(具免疫力)。在SIR model中，所有在S狀態的人都會慢慢經由I狀態轉變為R狀態。

– 接著所有S狀態的人會根據參數β，慢慢經由S狀態改變為I狀態。因此參數β可以視為一種傳染係數。

– 而所有I狀態的人則會根據參數λ，慢慢經由I狀態改變為R狀態。因此參數γ可以視為一種康復係數。

• SIR model的關係式如下：

• 由於我們數學不好，我們覺得解SIR model的積分方程太困難了。但我們仍然能夠用暴力把函數圖形給畫出來(使用暴力最大的缺點就是程式運行速度較慢，但只要能先寫出來慢一點點又何妨？)。

在R內使用迴圈功能

• 所謂的暴力其時就是迴圈，因此在寫出SIR model的函數以前我們要先學會迴圈。

– 下列是由1加到100的指令

x=0
for (i in 1:100) {
  x=x+i
}
x

– 下列是每隔3秒鐘讀取CDC的寫出一個檔案

for (i in 1:10) {
  data=readLines("http://nidss.cdc.gov.tw/download/Dengue_Daily_last12m.csv",encoding = "UTF-8", n = 50) #只讀取前50列
  write.table(data,paste0("e:/dengue",i,".csv"),col.names=FALSE,row.names=FALSE,quote=FALSE)
  Sys.sleep(3)
}

使用迴圈建立SIR model

• 我們依照下面的思路撰寫程式碼：

1. 設置m為預測時間範圍，可由使用者指定，在這設為50。

2. 接著令n.S為初始易感受個案人口數(S狀態，假設是999人)。

3. 接著令n.I為初始受感染個案人口數(I狀態，假設是1人)。

4. 接著令n.R為初始受感染個案人口數(R狀態，假設是0人)。

5. 令B為參數β(假設是1.5)，以及令G為參數γ(假設是0.3)

6. 先產生一個物件t(時間)，他是一個從0到m的數列。除此之外，產生對應的S、I、R數列，長度與t相同，預設理面全部都為初始值。

7. 計算總人口數N。

8. 建立迴圈計算各時間點S、I、R的個案數

9. 畫圖。

• 下列是建立SIR model的函數：

#Parameters setting
m=50
n.S=999
n.I=1
n.R=0
B=0.7
G=0.3

#Generate vectors
t=0:m                 #t=seq(0,m,by=0.1)
S=rep(n.S,m+1)        #S=rep(n.S,length(t))
I=rep(n.I,m+1)        #I=rep(n.I,length(t))
R=rep(n.R,m+1)        #R=rep(n.R,length(t))
N=n.S+n.I+n.R

#for loop
for (i in 1:(length(t)-1)) {
  #S->I
  SI=B*S[i]*I[i]/N
  #if time space is not equal to 1, then use follow program:
  #SI=B*S[i]*I[i]/N*(t[i+1]-t[i])
  
  #I->R
  IR=G*I[i]
  #if time space is not equal to 1, then use follow program:
  #IR=G*I[i]*(t[i+1]-t[i])
  
  #To caculate the number of S, I, R in each time point
  S[i+1]=S[i]-SI
  I[i+1]=I[i]+SI-IR
  R[i+1]=R[i]+IR
}

#Type 1 plot

plot(t,S,type="l",col="blue",ylab="Population",xlab="Time",xlim=c(0,m),ylim=c(0,N),lwd=2)
lines(t,I,col="red",lwd=2)
lines(t,R,col="darkgreen",lwd=2)
legend("topright",c("S","I","R"),pch=15,col=c("blue","red","darkgreen")
       ,bg="gray90",text.col=c("blue","red","green"),cex=1.5)

#Type 2 plot

plot(-1,-1,xlim=c(0,m),ylim=c(0,N),type="l",col="blue",ylab="Population",xlab="Time")
polygon(c(t,t[length(t):1]),c(rep(0,length(t)),S[length(S):1]),col="blue",border=NA)
polygon(c(t,t[length(t):1]),c(S,S[length(S):1]+I[length(I):1]),col="red",border=NA)
polygon(c(t,t[length(t):1]),c(S+I,rep(N,length(t))),col="darkgreen",border=NA)
legend("topright",c("S","I","R"),pch=15,col=c("blue","red","darkgreen")
       ,bg="gray90",text.col=c("blue","red","green"),cex=1.5)

整合上述過程為一個自訂函數

SIR = function (m, n.S, n.I, n.R, B, G, type = "1") {
    #Generate vectors
    t=0:m                 #t=seq(0,m,by=0.1)
    S=rep(n.S,m+1)        #S=rep(n.S,length(t))
    I=rep(n.I,m+1)        #I=rep(n.I,length(t))
    R=rep(n.R,m+1)        #R=rep(n.R,length(t))
    N=n.S+n.I+n.R
    
    #for loop
    for (i in 1:(length(t)-1)) {
        #S->I
        SI=B*S[i]*I[i]/N
        #if time space is not equal to 1, then use follow program:
        #SI=B*S[i]*I[i]/N*(t[i+1]-t[i])
        
        #I->R
        IR=G*I[i]
        #if time space is not equal to 1, then use follow program:
        #IR=G*I[i]*(t[i+1]-t[i])
        
        #To caculate the number of S, I, R in each time point
        S[i+1]=S[i]-SI
        I[i+1]=I[i]+SI-IR
        R[i+1]=R[i]+IR
    }
    
    if (type=="1") {
        #Type 1 plot
        plot(t,S,type="l",col="blue",ylab="Population",xlab="Time",xlim=c(0,m),ylim=c(0,N),lwd=2)
        lines(t,I,col="red",lwd=2)
        lines(t,R,col="darkgreen",lwd=2)
        legend("topright",c("S","I","R"),pch=15,col=c("blue","red","darkgreen")
               ,bg="gray90",text.col=c("blue","red","green"),cex=1.5)
    } else {
        #Type 2 plot
        plot(-1,-1,xlim=c(0,m),ylim=c(0,N),type="l",col="blue",ylab="Population",xlab="Time")
        polygon(c(t,t[length(t):1]),c(rep(0,length(t)),S[length(S):1]),col="blue",border=NA)
        polygon(c(t,t[length(t):1]),c(S,S[length(S):1]+I[length(I):1]),col="red",border=NA)
        polygon(c(t,t[length(t):1]),c(S+I,rep(N,length(t))),col="darkgreen",border=NA)
        legend("topright",c("S","I","R"),pch=15,col=c("blue","red","darkgreen")
               ,bg="gray90",text.col=c("blue","red","green"),cex=1.5)
    }
}

#Testing function
SIR(m = 100, n.S = 999, n.I = 1, n.R = 0, B = 0.3, G = 0.1, type = 1)

SIR(m = 100, n.S = 999, n.I = 1, n.R = 0, B = 0.3, G = 0.1, type = 2)

練習-4

• 請試著使用函數SIR，結合WebApp讓使用者能輸入參數產生圖形。

練習-4 解答

• ui.R

library(shiny)

shinyUI(pageWithSidebar(
  
  headerPanel("SIR model calculator"),
  
  sidebarPanel(
    sliderInput("m", "Infection time:", min = 10, max = 200, value = 30, step = 1),
    numericInput("n.S", "Number of initial susceptible cases:", min = 100, max = 1000000, value = 999, step = 1),
    numericInput("n.I", "Number of initial infected cases:", min = 1, max = 10000, value = 1, step = 1),
    numericInput("n.R", "Number of initial recovered cases:", min = 0, max = 1000000, value = 0, step = 1),
    numericInput("B", "The infected parameter (Beta):", min = 0.01, max = 10, value = 0.6, step = 0.01),
    numericInput("G", "The recovered parameter(Gamma):", min = 0.01, max = 1, value = 0.3, step = 0.01),
    radioButtons("type", "The type of plot:",  c("Type 1" = "1", "Type 2" = "2"))
  ),
  
  mainPanel(
    plotOutput("PLOT",width = "800px",height = "600px")
  )
))

• server.R

library(shiny)

SIR = function (m, n.S, n.I, n.R, B, G, type = "1") {
    #Generate vectors
    t=0:m                 #t=seq(0,m,by=0.1)
    S=rep(n.S,m+1)        #S=rep(n.S,length(t))
    I=rep(n.I,m+1)        #I=rep(n.I,length(t))
    R=rep(n.R,m+1)        #R=rep(n.R,length(t))
    N=n.S+n.I+n.R
    
    #for loop
    for (i in 1:(length(t)-1)) {
        #S->I
        SI=B*S[i]*I[i]/N
        #if time space is not equal to 1, then use follow program:
        #SI=B*S[i]*I[i]/N*(t[i+1]-t[i])
        
        #I->R
        IR=G*I[i]
        #if time space is not equal to 1, then use follow program:
        #IR=G*I[i]*(t[i+1]-t[i])
        
        #To caculate the number of S, I, R in each time point
        S[i+1]=S[i]-SI
        I[i+1]=I[i]+SI-IR
        R[i+1]=R[i]+IR
    }
    
    if (type=="1") {
        #Type 1 plot
        plot(t,S,type="l",col="blue",ylab="Population",xlab="Time",xlim=c(0,m),ylim=c(0,N),lwd=2)
        lines(t,I,col="red",lwd=2)
        lines(t,R,col="darkgreen",lwd=2)
        legend("topright",c("S","I","R"),pch=15,col=c("blue","red","darkgreen")
               ,bg="gray90",text.col=c("blue","red","green"),cex=1.5)
    } else {
        #Type 2 plot
        plot(-1,-1,xlim=c(0,m),ylim=c(0,N),type="l",col="blue",ylab="Population",xlab="Time")
        polygon(c(t,t[length(t):1]),c(rep(0,length(t)),S[length(S):1]),col="blue",border=NA)
        polygon(c(t,t[length(t):1]),c(S,S[length(S):1]+I[length(I):1]),col="red",border=NA)
        polygon(c(t,t[length(t):1]),c(S+I,rep(N,length(t))),col="darkgreen",border=NA)
        legend("topright",c("S","I","R"),pch=15,col=c("blue","red","darkgreen")
               ,bg="gray90",text.col=c("blue","red","green"),cex=1.5)
    }
}

shinyServer(function(input, output) {
  
  output$PLOT = renderPlot({
    SIR(m = input$m, n.S = input$n.S, n.I = input$n.I, n.R = input$n.R, B = input$B, G = input$G, type = input$type)
  })
  
})

小結

• 今天的課程教了幾項重點，請大家有空多多練習：

1. 安裝套件
2. 撰寫App(ui.R以及server.R)
3. 在R內自創函數
4. 利用conditionalPanel()增加使用者介面的自由度
5. 使用迴圈