0

I want to display an animation of a simulation using two plots.

  1. One is a "population" view of the simulation where points represent individuals. At every step in the simulation, I want to draw random points within a circle keeping those sampled in the next generation. The simulation stops when the frequency of any of the two types of individuals reaches 1 or 0.

  2. The other one is a frequency chart with frequency on the y-axis and generations on the x-axis. Ideally, I want this graph to expand at every generation and stop when the frequency reaches 1/0.

My first problem is that I can't get reactiveTimer() to work as want it to. It does not self-update or if it does it goes back to the starting point without "remembering" previous states.

My second problem is that if I use an if statement for the condition to keep the simulation going it only iterates a single generation after pressing run. Alternatively, if I use a while loop it will just go to directly to the last generation, skipping all the middle parts of the simulation.

My third problem is that I cannot grow a data.frame within a reactive environment so that I can plot the frequencies after each generation.

Code:

library(shiny)
library(ggplot2)

# function to make a circle data.frame
# https://stackoverflow.com/questions/6862742/draw-a-circle-with-ggplot2
circleFun <- function(center=c(0,0), diameter=10, npoints=100){
    r = diameter / 2
    tt = seq(0,2*pi,length.out = npoints)
    xx = center[1] + r * cos(tt)
    yy = center[2] + r * sin(tt)
    return(data.frame(x = xx, y = yy))
}

# ui
ui <- fluidPage(
    titlePanel("Genetic Drift Simulator"),
    sidebarLayout(
        sidebarPanel(
                # Input: select from menu
                numericInput(inputId = "population_size",
                          label = "Population size:",
                          value = 10,
                          step = 10),
                sliderInput(inputId = "initial_frequency",
                          label = "Initial frequency of allele 1:",
                          min = 0,
                          max = 1,
                          step = 0.1,
                          value = 0.5),
                actionButton(inputId = "run",
                            label = "Run simulation"),
                actionButton(inputId = "reset",
                        label = "Reset values")
                ),
        mainPanel(
            fluidRow(
                column(4,
                    verbatimTextOutput("text")
                    )
                ),
            fluidRow(
                column(8,
                    plotOutput("pop_plot")
                    )
                ),
            fluidRow(
                column(8,
                    plotOutput("freq_plot")
                    )
                )
            )
        )
    )

server <- function(input, output, session) {

    waits <- reactiveValues(timer = reactiveTimer(Inf))

    returns <- reactiveValues(
        z=NULL,
        x=NULL,
        y=NULL,
        freq=NULL,
        circle=NULL,
        i=NULL
        )

    frequencies <- reactiveValues(df=NULL)

    observe({
        returns$z=rbinom(input$population_size, 1, input$initial_frequency)
        returns$x=rnorm(input$population_size)
        returns$y=rnorm(input$population_size)
        returns$freq=input$initial_frequency
        returns$circle=circleFun()
        returns$i=0

        frequencies$df = data.frame(x=returns$i, y=returns$i)
    })

    population <- reactive({
        data.frame(x=returns$x, y=returns$y, z=returns$z)
    })

    grow_freq <- function(df, x, y){
        rbind(df, c(x,y))
    }

    grow <- reactive({
        frequencies$df = grow_freq(frequencies$df, returns$i, returns$freq)
    })

    drift <- reactive({
        returns$z = sample(returns$z, replace=T)
        # random locations
        returns$x = rnorm(input$population_size)
        returns$y = rnorm(input$population_size)
        # calculate frequency
        returns$freq = sum(returns$z == 1)/input$population_size
        # increase to next generation
        returns$i = returns$i+1
    })

    observeEvent(input$run, {
        if (returns$freq < 1 & returns$freq > 0){
            # observeEvent(reactiveTimer(200), {
                drift()
                grow()
            # })
        }
        # else {
        #     waits$timer <- reactiveTimer(Inf)
        # }
    })

    observeEvent(input$reset, {
        timer = reactiveTimer(Inf)

        returns$z = rbinom(input$population_size, 1, input$initial_frequency)
        returns$x = rnorm(input$population_size)
        returns$y = rnorm(input$population_size)
        returns$freq = input$initial_frequency
        returns$circle=circleFun(center=c(0,0), diameter=10, npoints=100)
        returns$i = 0

        frequencies$df = data.frame(x=returns$i, y=returns$i)
    })

    output$text <- renderText({
        text = paste("Population size: ",input$population_size,"\n",
                    "Frequency allele 1: ",returns$freq,"\n",
                    "Generation: ",returns$i, sep="")
        print(text)
        })

    output$pop_plot <- renderPlot({
        ggplot(data=population(), aes(x, y)) +
            geom_point(aes(color=factor(z)), size=5, alpha=0.7) +
            geom_path(data=returns$circle, color="black", size=2) +
            scale_color_brewer(type="qual", palette=1, name="allele") +
            theme(axis.title=element_blank(), axis.text=element_blank()) +
            theme(legend.title=element_text(size=16), legend.text=element_text(size=14))
        },
        height = 400, width = 450)

    output$freq_plot <- renderPlot({
        ggplot(data=frequencies$df, aes(x, y)) +
            geom_point() +
            geom_line() +
            ylim(0,1)
        },
        height = 300, width = 500)

}

shinyApp(ui = ui, server = server)

ssanch
  • 389
  • 2
  • 6

2 Answers2

2

The three problems you noted are resolved in the code below. I changed your first observe to observeEvent. Please note that you need to click on reset when the simulation ends.

# ui
ui <- fluidPage(
  titlePanel("Genetic Drift Simulator"),
  sidebarLayout(
    sidebarPanel(
      # Input: select from menu
      numericInput(inputId = "population_size",
                   label = "Population size:",
                   value = 10,
                   step = 10),
      sliderInput(inputId = "initial_frequency",
                  label = "Initial frequency of allele 1:",
                  min = 0,
                  max = 1,
                  step = 0.1,
                  value = 0.5),
      actionButton(inputId = "run",
                   label = "Run simulation"),
      actionButton(inputId = "reset",
                   label = "Reset values")
    ),
    mainPanel(
      fluidRow(
        column(4,
               verbatimTextOutput("text")
        )
      ),
      fluidRow(
        column(8,
               plotOutput("pop_plot")
        )
      ),
      fluidRow(
        column(8,
               plotOutput("freq_plot")
        )
      )
    )
  )
)

server <- function(input, output, session) {

  # Anything that calls autoInvalidate will automatically invalidate every 2 seconds.
  autoInvalidate <- reactiveTimer(2000)
  #waits <- reactiveValues(timer = reactiveTimer(Inf))

  returns <- reactiveValues(
    z=NULL,
    x=NULL,
    y=NULL,
    freq=NULL,
    circle=NULL,
    i=NULL
  )

  frequencies <- reactiveValues(df=NULL)
  #observe({
  observeEvent(list(input$population_size,input$initial_frequency), {
    returns$z=rbinom(input$population_size, 1, input$initial_frequency)
    returns$x=rnorm(input$population_size)
    returns$y=rnorm(input$population_size)
    returns$freq=input$initial_frequency
    returns$circle=circleFun()
    returns$i=0

    frequencies$df = data.frame(x=returns$i, y=returns$i)
  })

  population <- reactive({
    data.frame(x=returns$x, y=returns$y, z=returns$z)
  })

  grow_freq <- function(df, x, y){
    rbind(df, c(x,y))
  }

  grow <- reactive({
    frequencies$df <- grow_freq(frequencies$df, returns$i, returns$freq)
  })

  drift <- reactive({
    returns$z = sample(returns$z, replace=T)
    # random locations
    returns$x = rnorm(input$population_size)
    returns$y = rnorm(input$population_size)
    # calculate frequency
    returns$freq = sum(returns$z == 1)/input$population_size
    # increase to next generation
    returns$i = returns$i+1
  })

  observeEvent(list(input$run,autoInvalidate()), {

    if (returns$freq < 1 & returns$freq > 0){
      # observeEvent(reactiveTimer(200), {
      drift()
      grow()
      # })
    }
    # else {
    #   waits$timer <- reactiveTimer(200)
    # }
  })

  observeEvent(input$reset, {
    #timer = reactiveTimer(Inf)

    returns$z = rbinom(input$population_size, 1, input$initial_frequency)
    returns$x = rnorm(input$population_size)
    returns$y = rnorm(input$population_size)
    returns$freq = input$initial_frequency
    returns$circle=circleFun(center=c(0,0), diameter=10, npoints=100)
    returns$i = 0

    frequencies$df = data.frame(x=returns$i, y=returns$i)
  })

  output$text <- renderText({
    text = paste("Population size: ",input$population_size,"\n",
                 "Frequency allele 1: ",returns$freq,"\n",
                 "Generation: ",returns$i, sep="")
    print(text)
  })

  output$pop_plot <- renderPlot({
    autoInvalidate()
    ggplot(data=population(), aes(x, y)) +
      geom_point(aes(color=factor(z)), size=5, alpha=0.7) +
      geom_path(data=returns$circle, color="black", size=2) +
      scale_color_brewer(type="qual", palette=1, name="allele") +
      theme(axis.title=element_blank(), axis.text=element_blank()) +
      theme(legend.title=element_text(size=16), legend.text=element_text(size=14))
  },
  height = 400, width = 450)

  output$freq_plot <- renderPlot({
    autoInvalidate()
    ggplot(data=frequencies$df, aes(x, y)) +
      geom_point() +
      geom_line() +
      ylim(0,1)
  },
  height = 300, width = 500)

}

shinyApp(ui = ui, server = server)
YBS
  • 19,324
  • 2
  • 9
  • 27
  • @ssanch, please check the updated code. Do let me know if the `reactiveTimer` response is what you expected. – YBS Oct 08 '20 at 13:20
  • That's great, thanks! I'm glad to see I wasn't that far off. Is there a way to activate the timer only after clicking the "simulate" button? I tried a few things but nothing worked. Also, the frequency chart should start with y=0.5 and it starts at y=0. Not sure why that happens. – ssanch Oct 09 '20 at 04:09
  • I tweaked the code and it now works as expected, with the exception of the "reset" bottom. Not sure of the reason. – ssanch Oct 09 '20 at 05:02
  • Final tweak, it now works perfectly. Thanks for the help @YBS ! – ssanch Oct 09 '20 at 05:25
0

My edits to the accepted answers were rejected, so I'm posting the revised code here as an answer. The reason was that I was not addressing the OP's questions. Which is strange because the OP is me.

In short:

  1. I removed reactiveTimer functions from plots.
  2. The counter now starts after pressing the run button and not automatically.
  3. Fixed a bug with initial frequency on 2nd plot which was set to 0 and not 0.5
  4. Removed else statement on the simulation block to allow the reset button to work.

An example of the ShinyApp can be found in here.

# ui
ui <- fluidPage(
    titlePanel("Genetic Drift Simulator"),
    sidebarLayout(
        sidebarPanel(
            # Input: select from menu
            numericInput(inputId = "population_size",
            label = "Population size:",
            value = 10,
            step = 10),
            sliderInput(inputId = "initial_frequency",
            label = "Initial frequency of allele 1:",
            min = 0,
            max = 1,
            step = 0.1,
            value = 0.5),
            actionButton(inputId = "run",
            label = "Run simulation"),
            actionButton(inputId = "reset",
            label = "Reset values")
        ),
        mainPanel(
            fluidRow(
                column(4,
                    verbatimTextOutput("text")
                )
            ),
            fluidRow(
                column(8,
                    plotOutput("pop_plot")
                )
            ),
            fluidRow(
                column(8,
                    plotOutput("freq_plot")
                )
            )
        )
    )
)

server <- function(input, output, session) {

    # Anything that calls autoInvalidate will automatically invalidate.
    autoInvalidate <- reactiveValues(timer=NULL)

    returns <- reactiveValues(
        z=NULL,
        x=NULL,
        y=NULL,
        freq=NULL,
        circle=NULL,
        i=NULL
    )

    frequencies <- reactiveValues(df=NULL)

    observeEvent(list(input$population_size,input$initial_frequency), {
        returns$z=rbinom(input$population_size, 1, input$initial_frequency)
        returns$x=rnorm(input$population_size)
        returns$y=rnorm(input$population_size)
        returns$freq=input$initial_frequency
        returns$circle=circleFun()
        returns$i=0

        frequencies$df = data.frame(x=returns$i, y=returns$freq)

        autoInvalidate$timer = reactiveTimer(Inf)

    })

    population <- reactive({
        data.frame(x=returns$x, y=returns$y, z=returns$z)
    })

    grow_freq <- function(df, x, y){
        rbind(df, c(x,y))
    }

    grow <- reactive({
        frequencies$df <- grow_freq(frequencies$df, returns$i, returns$freq)
    })

    drift <- reactive({
        returns$z = sample(returns$z, replace=T)
        # random locations
        returns$x = rnorm(input$population_size)
        returns$y = rnorm(input$population_size)
        # calculate frequency
        returns$freq = sum(returns$z == 1)/input$population_size
        # increase to next generation
        returns$i = returns$i+1
    })

    observeEvent(input$run, {
        autoInvalidate$timer = reactiveTimer(1000) # changed to 1 second
        drift()
        grow()
    })

    observeEvent(autoInvalidate$timer(), {
        if (returns$freq < 1 & returns$freq > 0 & returns$i != 0){
            autoInvalidate$timer()
            drift()
            grow()
        }
        # else if (returns$freq == 0 | returns$freq == 1) {
        #     autoInvalidate$timer = reactiveTimer(Inf)
        # }
    })

    observeEvent(input$reset, {
        returns$z = rbinom(input$population_size, 1, input$initial_frequency)
        returns$x = rnorm(input$population_size)
        returns$y = rnorm(input$population_size)
        returns$freq = input$initial_frequency
        returns$circle=circleFun(center=c(0,0), diameter=10, npoints=100)
        returns$i = 0

        frequencies$df = data.frame(x=returns$i, y=returns$freq)

        autoInvalidate$timer = reactiveTimer(Inf)
    })

    output$text <- renderText({
        #autoInvalidate$timer()
        text = paste("Population size: ",input$population_size,"\n",
        "Frequency allele 1: ",returns$freq,"\n",
        "Generation: ",returns$i, sep="")
        print(text)
    })

    output$pop_plot <- renderPlot({
        #autoInvalidate$timer()
        ggplot(data=population(), aes(x, y)) +
        geom_point(aes(color=factor(z)), size=5, alpha=0.7) +
        geom_path(data=returns$circle, color="black", size=2) +
        scale_color_brewer(type="qual", palette=1, name="allele") +
        theme(axis.title=element_blank(), axis.text=element_blank()) +
        theme(legend.title=element_text(size=16), legend.text=element_text(size=14))
    },
    height = 400, width = 450)

    output$freq_plot <- renderPlot({
        #autoInvalidate$timer()
        if (dim(frequencies$df)[1] == 1){
            ggplot(data=frequencies$df, aes(x, y)) +
            geom_hline(yintercept=0.5) +
            geom_point() +
            labs(x="generation",y="frequency") +
            ylim(0,1)
        } else {
            ggplot(data=frequencies$df, aes(x, y)) +
            geom_hline(yintercept=0.5) +
            geom_point() +
            geom_line() +
            labs(x="generation",y="frequency") +
            ylim(0,1)
        }
    },
    height = 300, width = 500)

}

shinyApp(ui = ui, server = server)
ssanch
  • 389
  • 2
  • 6