Brownian Motion

If you want to skip all the code and Shiny info you can jump straight to the shiny app: Brownian motion simulator.

In comparative phylogenetics Brownian motion is often either a null model that we compare observations to or is a de facto assumption that our methods depend on. For this reason, I think it can be useful to have nice interactive ways for students to explore this model and develop an intuitive understanding of how it behaves. The R package Shiny offers a great environment to develop something like this. Below you can see just how easy it is to create a slick and responsive web app.

The first step is designing a file called “ui.R” this is the user interface. It consists of two primary code chunks the first is the sidebarLayout this describes the part of the page that the user will interact with. The Shiny package has many implemented widgets built into the package and in this example I use two of these one is sliderInput and the other is the actionButton. SliderInput creates a sliding widget and lets you set the min and max value for the variable as well as an initial starting value. The actionButton widget creates a variable that initially has a value of 1 and this is incremented by 1 each time the button is pressed. I use the actionButton to act as the seed for the stochastic portion of my code. The second part of this simple ui.R file specifies the output to place next to the sidebar. In this example I simply call the object “distPlot” which is created in the server page and I also specify the vertical size to be larger than it would be by default.

 library(shiny)  
 # Define UI for application  
 shinyUI(fluidPage(  
  # Application title  
  titlePanel("Brownian Motion"),  
  # Sidebar for user input  
  sidebarLayout(  
   sidebarPanel(  
    # select how many generations to run  
    sliderInput("gens",  
          "Generations:",  
          min = 10,  
          max = 500,  
          value = 500),  
    # select the number of replicates  
    sliderInput("reps",  
          "Replicates:",  
          min = 10,  
          max = 500,  
          value = 100),  
    # select the part to plot a histogram of  
    sliderInput("mon",  
          "Generation to plot:",  
          min = 10,  
          max = 500,  
          value = 100),  
    # refresh via a seed change  
    actionButton("seed.val", 'Refresh')  
   ),  
   # Show a plot of the simulation result and make it a bit larger with height argument  
   mainPanel(  
    plotOutput("distPlot", height="600px")  
   )  
  )  
 ))

The next component of the Shiny app is the “server.R” file. This is the code that will be run to generate the outpot portion of the webapp. This consists of two components in the example. The first is the reactive component that is the result of the Brownian motion simulation. By using replicate we will create a matrix where the rows represent generations and the columns represent iterations. By placing it inside of a reactive command we insure that it will only be rerun when one of the underlying input variables changes. If you look at the code you can see that this includes: input$seed.val, input$reps, input$gens. In contrast if another variable such as input$mon changes “x” will not be recalculated.

 library(shiny)  
 # I like this color pallete much more than the base options  
 library(viridis)  
 shinyServer(function(input, output) {  
  # this is the actual brownian motion simulation  
  x <- reactive({  
   # insure reproducability  
   set.seed <- input$seed.val  
   replicate(input$reps, cumsum(rnorm(n=input$gens, sd=input$rate)))  
  })  
  # Now we will build the output plot  
  output$distPlot <- renderPlot({   
   # we will have two plots  
   par(mfcol=c(2,1))  
   # need a bit more room   
   par(mar=c(3,2,1,1))  
   # this just keeps our color variation from being hidden when itterations++  
   cols <- sample(viridis(input$reps))  
   # the primary plot  
   matplot(x(),type="l",lty=1, lwd=1.5, col=cols, main="Brownian Motion", ylab="Trait Value", xlab="Generations")  
   # this line will indicate what is being plotted in the histogram below  
   abline(v=input$mon, lwd=3, col="red")  
   # now we plot the distribution switched from a histrogram this looks prettier IMHO  
   plot(density(x()[input$mon,]),   
      xlim=c(min(x()),max(x())),  
      col = 'red',   
      xlab="Trait value",  
      ylab="",  
      lwd=3,  
      main=paste("Trait distribution at generation", input$mon))  
  })  
 })

The second element of the code is the actual plot I want output. Because x is a matrix we can use matplot to plot the results quickly and easily. I choose to use a new color pallete viridis. However, we want to make this interactive by adding a histrogram of any time point on the graph. We do this by using the argument input$mon to select a specific row of x to plot. This determines what data is used for the second graph and adds a vertical line indicating what time point is being plotted.

Here is what the finished product looks like:

Coleopterists Corner

Using Coleoptera to understand evolution; and lots of R

Add a comment