.Rhistory

# overdominance always results in 1 stable attractor at intermediate allele freq
} else if(w_genotypes["Aa"] >= w_genotypes["AA"]){
output <- list(outcome = "p intermediate")
}
# for h=0.2, it's tricky to draw the arrows bc the stable attractor is near 0.95
return(output)
}
outcome <- get_steady_state(w_gent, delta_p)
outcome
runApp()
runApp()
runApp()
delta_p
w_gent
which(near(delta_p$delta_p, 0, tol = 10^-3))
w_gent <- fitness_genotypes(sel_coef = 0.2, dom_coef = 0.3)
delta_p <- get_p_byp(w_gent)
which(near(delta_p$delta_p, 0, tol = 10^-3))
temp_list <- vector(mode='list', length=8)
h <- seq(from=0.3, to=1.0, by=0.1)
for(i in 1:8){
w_gent <- fitness_genotypes(sel_coef = 0.2, dom_coef = h[i])
temp_list[[i]] <- get_p_byp(w_gent)
}
h[8]
h[1]
for (i in 1:8){
print(which(near(temp_list[[i]]$delta_p, 0, tol = 10^-3)))
}
for (i in 1:8){
print(which(near(temp_list[[i]]$delta_p[5:980], 0, tol = 10^-3)))
}
w_gent <- fitness_genotypes(sel_coef = 0.2, dom_coef = 0.2)
delta_p <- get_p_byp(w_gent)
which(near(delta_p$delta_p, 0, tol = 10^-3))
which(near(delta_p$delta_p[10:1000], 0, tol = 10^-3))
delta_p$p[which(near(delta_p$delta_p[10:1000], 0, tol = 10^-3))]
delta_p$delta_p[which(near(delta_p$delta_p[10:1000], 0, tol = 10^-3))]
w_gent
delta_p$delta_p[990:1000]
runApp()
runApp()
outcome
w_gent
delta_p
min(delta_p$delta_p)
max(delta_p$delta_p)
runApp()
runApp()
runApp()
runApp()
runApp()
delta_p$delta_p
max(delta_p$delta_p)
max(delta_p$delta_p)/2
w_gent
runApp()
runApp()
runApp()
runApp()
runApp()
runApp()
runApp()
runApp()
runApp()
runApp()
runApp()
runApp()
runApp()
runApp()
w_gent
w_gent <- fitness_genotypes(sel_coef = 0.21, dom_coef = 0.8)
w_gent
delta_p <- get_p_byp(w_gent)
outcome <- get_steady_state(w_gent, delta_p)
outcome
# a function indicating the long-term steady state
get_steady_state <- function(w_genotypes, deltap.df){
# some presets for plotting the arrows
min_arrow_length = 0.07
arrow_vertical_whitespace = 0.1
# initialize output variable with NA values:
# depending on outcome, only some may be replaced with non-NA values
output <- list(outcome = NA, stable_pts = NA, unstable_pt = NA, arrow1 = NA, arrow2 = NA)
# underdominance always results in bistability
if(w_genotypes["Aa"] < w_genotypes["aa"]) {
# a hacky way to approximate where the unstable equilibrium point is:
unstable_pt <- mean(deltap.df$p[which(near(deltap.df$delta_p[12:990], 0, tol = 10^-3))])
# draw the arrows symmetrically about x=0 at halfway of the smallest max delta_p values
x_val = min(abs(min(deltap.df$delta_p)), abs(max(deltap.df$delta_p)))/2
# y values for bottom arrow (1): vertically in middle of the 2 equil pts
arrow1_length <- max(unstable_pt-2*arrow_vertical_whitespace, min_arrow_length)
# y values for top arrow (2): vertically in middle of the 2 equil pts
arrow2_length <- max(1-unstable_pt-2*arrow_vertical_whitespace, min_arrow_length)
output <- list(outcome = "bistable",
stable_pts = c(0, 1),
unstable_pt = unstable_pt,
arrow1 = c(x0 = -x_val, y0 = unstable_pt - arrow_vertical_whitespace,
x1 = -x_val, y1 = unstable_pt - arrow_vertical_whitespace - arrow1_length),
arrow2 = c(x0 = x_val, y0 = unstable_pt + arrow_vertical_whitespace,
x1 = x_val, y1 = unstable_pt + arrow_vertical_whitespace + arrow2_length)
)
# additivity / partial dominance always results in A going to fixation (1 stable attractor)
} else if(w_genotypes["Aa"] >= w_genotypes["aa"] &
w_genotypes["Aa"] < w_genotypes["AA"]) {
x_val <- 0.005 # this x-value looks nice
# increase the whitespace of the arrow
arrow_vertical_whitespace <- 3*arrow_vertical_whitespace
arrow1_length <- 1-2*arrow_vertical_whitespace
output <- list(outcome = "p=1",
stable_pts = 1,
arrow1 = c(x0 = x_val, y0 = arrow_vertical_whitespace,
x1 = x_val, y1 = arrow_vertical_whitespace + arrow1_length)
)
# overdominance always results in 1 stable attractor at intermediate allele freq
} else if(w_genotypes["Aa"] > w_genotypes["AA"]){
# a hacky way to approximate where the stable equilibrium point is:
stable_pt <- mean(deltap.df$p[which(near(deltap.df$delta_p[5:980], 0, tol = 10^-3))])
# draw the arrows symmetrically about x=0 at halfway of the smallest max delta_p values
x_val = min(abs(min(deltap.df$delta_p)), abs(max(deltap.df$delta_p)))/2
# y values for bottom arrow (1): vertically in middle of the 2 equil pts
arrow1_length <- max(stable_pt-2*arrow_vertical_whitespace, min_arrow_length)
# y values for top arrow (2): vertically in middle of the 2 equil pts
arrow2_length <- max(1-stable_pt-2*arrow_vertical_whitespace, min_arrow_length)
output <- list(outcome = "p intermediate",
stable_pts = stable_pt,
arrow1 = c(x0 = x_val, y0 = stable_pt - arrow_vertical_whitespace - arrow1_length,
x1 = x_val, y1 = stable_pt - arrow_vertical_whitespace)
)
}
return(output)
}
outcome <- get_steady_state(w_gent, delta_p)
outcome
plot_p_byp <- function(deltap.df, steady_state) {
# change the graphing settings
par(cex.lab = 2, # increase size of axis labels
cex.axis = 1.5, # increase size of tick mark labels
mar = c(4.2, 4.4, 0.7, 0.85)) # decrease the borders
# plot the points
plot(x = deltap.df$delta_p, y = deltap.df$p,
pch = 16, col = "black",
xlab = expression("Change in allele frequency (" * Delta * "p)"),
ylab = "Allele frequency", xlim = range(deltap.df$delta_p),
ylim = c(-0.01, 1.01), yaxs = "i")
# plot the arrows
if(steady_state$outcome == "bistable"){
arrows(x0 = steady_state$arrow1["x0"],
y0 = steady_state$arrow1["y0"],
x1 = steady_state$arrow1["x1"],
y1 = steady_state$arrow1["y1"])
arrows(x0 = steady_state$arrow2["x0"],
y0 = steady_state$arrow2["y0"],
x1 = steady_state$arrow2["x1"],
y1 = steady_state$arrow2["y1"])
} else if(steady_state$outcome == "p=1"){
arrows(x0 = steady_state$arrow1["x0"],
y0 = steady_state$arrow1["y0"],
x1 = steady_state$arrow1["x1"],
y1 = steady_state$arrow1["y1"])
} else if(steady_state$outcome == "p intermediate"){
arrows(x0 = steady_state$arrow1["x0"],
y0 = steady_state$arrow1["y0"],
x1 = steady_state$arrow1["x1"],
y1 = steady_state$arrow1["y1"])
arrows(x0 = steady_state$arrow2["x0"],
y0 = steady_state$arrow2["y0"],
x1 = steady_state$arrow2["x1"],
y1 = steady_state$arrow2["y1"])
}
# Add a dotted vertical line at x = 0
abline(v = 0, lty = 3, lwd = 0.6, col = "black")
}
plot_p_byp(delta_p, outcome)
runApp()
runApp()
runApp()
# a function indicating the long-term steady state
get_steady_state <- function(w_genotypes, deltap.df){
# some presets for plotting the arrows
min_arrow_length = 0.07
arrow_vertical_whitespace = 0.1
# initialize output variable with NA values:
# depending on outcome, only some may be replaced with non-NA values
output <- list(outcome = NA, stable_pts = NA, unstable_pt = NA, arrow1 = NA, arrow2 = NA)
# underdominance always results in bistability
if(w_genotypes["Aa"] < w_genotypes["aa"]) {
# a hacky way to approximate where the unstable equilibrium point is:
unstable_pt <- mean(deltap.df$p[which(near(deltap.df$delta_p[12:990], 0, tol = 10^-3))])
# draw the arrows symmetrically about x=0 at halfway of the smallest max delta_p values
x_val = min(abs(min(deltap.df$delta_p)), abs(max(deltap.df$delta_p)))/2
# y values for bottom arrow (1): vertically in middle of the 2 equil pts
arrow1_length <- max(unstable_pt-2*arrow_vertical_whitespace, min_arrow_length)
# y values for top arrow (2): vertically in middle of the 2 equil pts
arrow2_length <- max(1-unstable_pt-2*arrow_vertical_whitespace, min_arrow_length)
output <- list(outcome = "bistable",
stable_pts = c(0, 1),
unstable_pt = unstable_pt,
arrow1 = c(x0 = -x_val, y0 = unstable_pt - arrow_vertical_whitespace,
x1 = -x_val, y1 = unstable_pt - arrow_vertical_whitespace - arrow1_length),
arrow2 = c(x0 = x_val, y0 = unstable_pt + arrow_vertical_whitespace,
x1 = x_val, y1 = unstable_pt + arrow_vertical_whitespace + arrow2_length)
)
# additivity / partial dominance always results in A going to fixation (1 stable attractor)
} else if(w_genotypes["Aa"] >= w_genotypes["aa"] &
w_genotypes["Aa"] < w_genotypes["AA"]) {
x_val <- 0.005 # this x-value looks nice
# increase the whitespace of the arrow
arrow_vertical_whitespace <- 3*arrow_vertical_whitespace
arrow1_length <- 1-2*arrow_vertical_whitespace
output <- list(outcome = "p=1",
stable_pts = 1,
arrow1 = c(x0 = x_val, y0 = arrow_vertical_whitespace,
x1 = x_val, y1 = arrow_vertical_whitespace + arrow1_length)
)
# overdominance always results in 1 stable attractor at intermediate allele freq
} else if(w_genotypes["Aa"] > w_genotypes["AA"]){
# a hacky way to approximate where the stable equilibrium point is:
stable_pt <- mean(deltap.df$p[which(near(deltap.df$delta_p[5:980], 0, tol = 10^-3))])
# draw the arrows symmetrically about x=0 at one third of the smallest max delta_p values
x_val = min(abs(min(deltap.df$delta_p)), abs(max(deltap.df$delta_p)))/3
# y values for bottom arrow (1): vertically in middle of the 2 equil pts
arrow1_length <- max(stable_pt-2*arrow_vertical_whitespace, min_arrow_length)
# y values for top arrow (2): vertically in middle of the 2 equil pts
arrow2_length <- max(1-stable_pt-2*arrow_vertical_whitespace, min_arrow_length)
output <- list(outcome = "p intermediate",
stable_pts = stable_pt,
arrow1 = c(x0 = x_val, y0 = stable_pt - arrow_vertical_whitespace - arrow1_length,
x1 = x_val, y1 = stable_pt - arrow_vertical_whitespace),
arrow2 = c(x0 = -x_val, y0 = stable_pt + arrow_vertical_whitespace + arrow2_length,
x1 = -x_val, y1 = stable_pt + arrow_vertical_whitespace)
)
}
return(output)
}
outcome <- get_steady_state(w_gent, delta_p)
outcome
runApp()
?arrows
runApp()
runApp()
runApp()
runApp()
runApp()
runApp()
runApp()
runApp()
exp(-1)
exp(-10)
exp(-100)
runApp()
runApp()
runApp()
runApp()
?withMathJax
runApp()
runApp()
runApp()
runApp()
runApp()
runApp()
runApp()
runApp()
runApp()
runApp()
runApp()
runApp()
runApp()
1.25/2
(1.25+1)/2
runApp('C:/Users/hermina/Downloads/for_shiny_server/for_shiny_server/underdom')
version(shiny)
packageVersion(shiny)
packageVersion("shiny")
library(shiny)
packageVersion("shiny")
packageVersion("ggplot2")
ip <- as.data.frame(installed.packages()[,c(1,3:4)])
rownames(ip) <- NULL
ip <- ip[is.na(ip$Priority),1:2,drop=FALSE]
print(ip, row.names=FALSE)
# load the environment
load("./app.R")
# Chunk 1: setup
# load the environment
source("./app.R")
# Chunk 2: the_app
# user interface with fluidPage
ui <- fluidPage(
# title
titlePanel("Diploid model of selection with underdominance"),
# top row is split into 3 columns of equal size (with 3 sliders)
fluidRow(
column(4,
sliderInput(inputId = "h",
label = "Underdominance coefficient (h)",
value = 0.2,
min = -1, max = 0.5,
step = 0.1)
),
column(4,
sliderInput(inputId = "p_init",
label = "Initial allele frequency of A (p)",
value = 0.32,
min = 0.02, max = 0.98,
step = 0.02)
),
column(4,
sliderInput(inputId = "gen",
label = "Number of generations",
value = 50,
min = 5, max = 950,
step = 5)
)
),
# middle row is 2 columns of equal size
fluidRow(
# interactively display the model parameters
column(6,
uiOutput("dynamEq"),
plotOutput("fitLand")),
# display the genotype frequencies over time
column(6,
plotOutput("genoTime")
)
),
# bottom row is 2 columns of equal size
fluidRow(
# display delta p by p
column(6,
plotOutput("delta_p")
),
# display the allele frequency over time
column(6,
plotOutput("alleleTime")
)
)
)
# Server logic
server <- function(input, output) {
# fix the selection coefficient
# be careful to choose a value that does *NOT* result in neutral system
s <- 0.305
######################
# expressions
######################
# a reactive expression to get the fitness of the genotypes
w_gent <- reactive(fitness_genotypes(sel_coef = s,
dom_coef = input$h))
# a reactive expression to get delta p as a function of allele frequecy (p)
delta_p <- reactive(get_p_byp(w_gent()))
# a reactive expression to get the long-term model outcome
outcome <- reactive(get_steady_state(w_genotypes = w_gent(), deltap.df = delta_p()))
# a reactive expression to simulate evolution of the system over time
sims <- reactive(sim_forward_time(num_gens = input$gen,
freq_init = input$p_init,
w_gent()))
######################
# outputs
######################
output$dynamEq <- renderUI({
temp_genotypes <- w_gent()
w_AA <- temp_genotypes["AA"]
w_Aa <- temp_genotypes["Aa"]
w_aa <- temp_genotypes["aa"]
withMathJax(sprintf("$${\\tiny p_{t+1} = \\frac{p_t^2 \\cdot (1+s) + p_t(1-p_t) \\cdot (1+s+h)}{\\overline{w}} = \\frac{p_t^2 \\cdot %.01f + p_t(1-p_t) \\cdot %.01f}{p_t^2 \\cdot %.01f + 2p_t(1-p_t) \\cdot %.01f + (1-p_t)^2 \\cdot %.01f}}$$",
w_AA, w_Aa, w_AA, w_Aa, w_aa))
})
output$fitLand <- renderPlot({plot_schematic(w_gent())})
output$genoTime <- renderPlot({plot_t_genotypes(sims())})
output$delta_p <- renderPlot({plot_p_byp(delta_p(), outcome())})
output$alleleTime <- renderPlot({plot_t_allele(sims())})
}
# Complete app with UI and server components
shinyApp(ui, server)
seq(from=-1, to=0.5, by=0.3)
plot_schematic(fitness_genotypes(sel_coef = 0.305, dom_coef = 0.15))
plot_schematic <- function(w_genotypes) {
# change the graphing settings
par(cex.lab = 2, # increase size of axis labels
cex.axis = 1.5, # increase size of tick mark labels
mar = c(4.2, 4.4, 1, 0.1)) # decrease the borders
# make the barplot in graphics (i.e., base R)
barplot(
w_genotypes,
col = colours_genotypes[names(w_genotypes)],
border = "black",
ylim = c(0, 1.6),
ylab = "Fitness",
xlab = "Genotype"
)
# add a box around the whole plot
box()
}
plot_schematic(fitness_genotypes(sel_coef = 0.305, dom_coef = 0.15))
plot_schematic(fitness_genotypes(sel_coef = 0.305, dom_coef = -0.5))
plot_schematic(fitness_genotypes(sel_coef = 0.305, dom_coef = 0.5))
plot_schematic(fitness_genotypes(sel_coef = 0.305, dom_coef = 0.55))
plot_schematic(fitness_genotypes(sel_coef = 0.305, dom_coef = 0.3))
runApp('temp_example_app.R')
runApp('temp_example_app.R')
runApp('temp_example_app.R')
runApp('temp_example_app.R')
runApp('temp_example_app.R')
runApp('stats_students_trial.R')
runApp('stats_students_trial.R')
runApp('stats_students_trial.R')
runApp('stats_students_trial.R')
runApp('stats_students_trial.R')
runApp('stats_students_trial.R')
runApp('stats_students_trial.R')
runApp('stats_students_trial.R')
runApp('stats_students_trial.R')
runApp('stats_students_trial.R')
runApp('stats_students_trial.R')
runApp('stats_students_trial.R')
runApp('stats_students_trial.R')
colours_genotypes
runApp('stats_students_trial.R')
runApp('stats_students_trial.R')
runApp('stats_students_trial.R')
runApp('stats_students_trial.R')
runApp('stats_students_trial.R')
runApp('stats_students_trial.R')
?bs_theme
runApp('stats_students_trial.R')
runApp('stats_students_trial.R')
runApp('stats_students_trial.R')
runApp('stats_students_trial.R')
runApp('stats_students_trial.R')
?layout_columns
runApp('stats_students_trial.R')
runApp('stats_students_trial.R')
runApp('stats_students_trial.R')
runApp('stats_students_trial.R')
runApp('stats_students_trial.R')
runApp('stats_students_trial.R')
runApp('stats_students_trial.R')
runApp('stats_students_trial.R')
?par
runApp('stats_students_trial.R')
runApp('stats_students_trial.R')
runApp('stats_students_trial.R')
runApp('stats_students_trial.R')
runApp('stats_students_trial.R')
runApp('stats_students_trial.R')
runApp('stats_students_trial.R')
runApp('stats_students_trial.R')
runApp('stats_students_trial.R')
runApp('stats_students_trial.R')
runApp('stats_students_trial.R')
runApp('stats_students_trial.R')
runApp('stats_students_trial.R')
runApp('stats_students_trial.R')
runApp('stats_students_trial.R')
runApp('stats_students_trial.R')
runApp('stats_students_trial.R')
runApp('stats_students_trial.R')
runApp('stats_students_trial.R')
runApp('stats_students_trial.R')
runApp('stats_students_trial.R')
runApp('stats_students_trial.R')
runApp('stats_students_trial.R')
runApp('stats_students_trial.R')
runApp('stats_students_trial.R')
runApp('stats_students_trial.R')
runApp('stats_students_trial.R')
runApp('stats_students_trial.R')
runApp('stats_students_trial.R')
runApp('stats_students_trial.R')
runApp('stats_students_trial.R')
runApp('stats_students_trial.R')
runApp('stats_students_trial.R')
runApp('stats_students_trial.R')
runApp('stats_students_trial.R')
runApp('stats_students_trial.R')
runApp('stats_students_trial.R')
runApp('stats_students_trial.R')
runApp('stats_students_trial.R')
runApp('stats_students_trial.R')
runApp('stats_students_trial.R')
runApp('stats_students_trial.R')
runApp('stats_students_trial.R')
runApp('stats_students_trial.R')
runApp('stats_students_trial.R')
runApp('stats_students_trial.R')
runApp('stats_students_trial.R')
runApp('stats_students_trial.R')
1+1.5
2.5/3
runApp('stats_students_trial.R')
1.5/3
0.75*2
runApp('stats_students_trial.R')
runApp('stats_students_trial.R')
runApp('stats_students_trial.R')
runApp('temp_example_app.R')
runApp('temp_example_app.R')
runApp('temp_example_app.R')
runApp('temp_example_app.R')
runApp('temp_example_app.R')
runApp('temp_example_app.R')
runApp('temp_example_app.R')
runApp('temp_example_app.R')
runApp('temp_example_app.R')
runApp('temp_example_app.R')
?page_navbar
runApp('temp_example_app.R')
runApp('temp_example_app.R')
runApp('temp_example_app.R')
update.packages("bslib")
update.packages("shiny")