I want to display an animation of a simulation using two plots.
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.
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)