Skip to content

Instantly share code, notes, and snippets.

@mrecos
Created August 23, 2016 01:33
  • Star 0 You must be signed in to star a gist
  • Fork 0 You must be signed in to fork a gist
Star You must be signed in to star a gist
Save mrecos/1d577157c132398763a05dcbfdd5b1b9 to your computer and use it in GitHub Desktop.
Code for Pokemon post-evolution CP model blog post
########## FUNCTIONS ###################
sign_flip <- function(x){
flip <- rbinom(1,1,0.5)
x <- ifelse(flip == 0, x, -x)
return(x)
}
sim_CP_new <- function(dat, spec, shape1, shape2, reps = 100){
spec_dat <- filter(dat1, species == eval(spec))
b_spec <- rbeta(10000, shape1, shape2)
b_spec <- sapply(b_spec, FUN = sign_flip)
b_spec <- b_spec + spec_dat$spec_mean_pdcp[1]
sim_spec <- replicate(reps,spec_dat$cp * (1 + sample(b_spec, nrow(spec_dat))))
return(sim_spec)
}
prep_sim_data_for_plot <- function(sim_spec, spec_dat){
sim_spec <- data.frame(sim_spec,
cp_new = spec_dat$cp_new)
spec_melt <- reshape2::melt(sim_spec)
spec_melt$model <- rep(c("sim","data"),
times = c(((ncol(sim_spec)-1)*nrow(sim_spec))
,nrow(sim_spec)))
spec_melt$species <- as.character(unique(spec_dat$species))
return(spec_melt)
}
sim_beta_RMSE <- function(dat, melt, draws_from_beta){
melt <- filter(melt, model == "sim")
melt$cp_new <- rep(dat$cp_new, times = draws_from_beta)
melt$id <- rep(1:nrow(dat), times = draws_from_beta)
beta_rmse <- group_by(melt, id) %>%
summarise(means = mean(value),
sd = sd(value),
rmse = sqrt(mean((cp_new - value)^2)),
mae = mean(abs(cp_new - value)),
species = eval(unique(as.character(dat$species))),
cp_new = mean(cp_new)) %>%
data.frame()
return(beta_rmse)
}
sim_CP_new_by_param <- function(dat, spec, param_dist, reps = 100){
shape_sample <- param_dist[sample(nrow(param_dist),1),]
shape1 <- shape_sample[1,1]
shape2 <- shape_sample[1,2]
spec_dat <- filter(dat1, species == eval(spec))
f_spec <- fitdist(abs(spec_dat$cent_spec_mean_pdcp),"beta")
b_spec <- rbeta(10000, shape1, shape2)
b_spec <- sapply(b_spec, FUN = sign_flip)
b_spec <- b_spec + spec_dat$spec_mean_pdcp[1]
# simulate cp_new of known pokes under [reps] different beta draws
sim_spec <- replicate(reps,spec_dat$cp * (1 + sample(b_spec, nrow(spec_dat))),
simplify = TRUE)
return(sim_spec)
}
sim_param_RMSE <- function(dat, melt, param_draws){
melt$cp_new <- rep(dat$cp_new, param_draws)
sim_rmse <- group_by(melt, Var1) %>%
summarise(means = mean(value),
sd = sd(value),
rmse = sqrt(mean((cp_new - value)^2)),
mae = mean(abs(cp_new - value)),
species = eval(unique(as.character(dat$species))),
cp_new = mean(cp_new)) %>%
data.frame()
}
prep_psim_for_plot <- function(psim_spec_dat, spec_dat){
psim_spec <- matrix(psim_spec_dat, nrow = nrow(spec_dat),
ncol = param_draws * beta_draws,
byrow = FALSE)
psim_spec_melt <- melt(psim_spec)
psim_spec_means <- group_by(psim_spec_melt, Var1) %>%
summarise(means = mean(value)) %>%
arrange(means)
psim_spec_melt$ord <- factor(psim_spec_melt$Var1, levels = psim_spec_means$Var1)
spec_dat$id <- seq(1:nrow(spec_dat))
spec_dat$ord <- factor(spec_dat$id, levels = psim_spec_means$Var1)
return(list(psim_spec_melt, spec_dat))
}
pedict_simulated_cp_new <- function(sim_dat, shape1, shape2, reps = 100){
# replicate a bunch of beta values from parametrized distribution
b_spec <- rbeta(10000, shape1, shape2)
# use sign_flip to randomly flip sign (p = 0.5)
b_spec <- sapply(b_spec, FUN = sign_flip)
# add the pos/neg beta (variation) to the mean pdcp for that species
b_spec <- b_spec + sim_dat$spec_mean_pdcp[1]
# replicate evolutions scenerios sampled from b_spec on simulated pokes
sim_spec <- replicate(reps,sim_dat$cp * (1 + sample(b_spec, nrow(sim_dat))))
return(sim_spec)
}
mround <- function(x,base){
base*round(x/base)
}
sim_new_species <- function(spec_ranges, spec_dat){
spec <- eval(unique(spec_dat$specie))
spec_cp_range <- data.frame(species = eval(unique(spec_dat$species)),
cp = mround(seq(spec_ranges[which(spec_ranges$species == spec),"min_cp"],
spec_ranges[which(spec_ranges$species == spec),"max_cp"],
length.out = length_out), base),
spec_mean_pdcp = spec_dat$spec_mean_pdcp[1])
return(spec_cp_range)
}
make_pred_cp_plot_dat <- function(predicted_spec, spec_cp_range){
spec_cp_new_ranges <- group_by(melt(predicted_spec), Var1) %>%
summarise(
cp_new_min = floor(min(value)),
cp_new_5 = floor(quantile(value, 0.05)),
cp_new_mean = floor(mean(value)),
cp_new_95 = ceiling(quantile(value, 0.95)),
cp_new_max = ceiling(max(value))) %>%
mutate(cp = spec_cp_range$cp) %>%
data.frame()
return(spec_cp_new_ranges)
}
########## END FUNCTIONS ######################
########## Require Libraries ##################
library("dplyr") # for data wrangling
library("corrplot") # for correlation plot
library("fitdistrplus") # for distribution GoF
library("ggplot2") # for plotting
library("ggalt") # for added plotting functions
library("broom") # for tidy
library("tidyr") # for gather()
library("tibble") # for display of data
library("xtable") # for output of tables for blog
library("ggthemes") # for color-blind safe plot colors
library("reshape2") # for melt function
########## END Require Libraries ##############
setwd("~/Dropbox/R/pokemon/")
pokedat <- as_tibble(read.csv("pokemon.csv"))
dat1 <- mutate(pokedat, delta_cp = cp_new - cp) %>% # net change in cp
# the % change in cp from old to new
mutate(pcnt_d_cp = delta_cp / cp) %>%
# group by species to calculate additonal varaibles
group_by(species) %>%
# species grouped mean percent change
mutate(spec_mean_pdcp = mean(pcnt_d_cp)) %>%
# species grouped std dev of % changes from old to new cp
mutate(spec_mean_pdcp_sd = sd(pcnt_d_cp)) %>%
# difference in % delta cp (pdcp) change from species group mean
mutate(cent_spec_mean_pdcp = pcnt_d_cp - spec_mean_pdcp) %>%
# z score for pdcp change within species group
mutate(z_spec_mean_pdcp = cent_spec_mean_pdcp / spec_mean_pdcp_sd) %>%
data.frame()
# Z score to answer the more appropriate question:
# what leads to an above average gain in CP per evolution within a given species?
print(xtable(table(dat1$species)), type = "html")
# plot the full 2k by 2k base pairs() plot for fun
# corr plot
c2 <- cor(dat1[,c(3:8,11,14:20,23,26,28:33)])
corrplot(c2, method = "square", order = "hclust",
tl.col = "black", tl.cex = 1.5)
# mean CP & HP Table
dat1 %>%
group_by(species) %>%
summarize(Mean_CP = round(mean(cp),2),
Mean_HP = round(mean(hp),2)) %>%
xtable(.) %>%
print(., type = "html", include.rownames=FALSE)
# R and R^2 table
dat1 %>%
group_by(species) %>%
summarize(r = round(cor(cp, hp),3)) %>%
mutate(r2 = round(r^2,3)) %>%
xtable(.) %>%
print(., type = "html", include.rownames=FALSE)
# plot linear correlation of CP & HP by species
ggplot(data = dat1, aes(x = cp, y = hp, color = species)) +
geom_point() +
geom_smooth(method = "lm", se = FALSE) +
# scale_color_pokemon("infernape") +
scale_color_colorblind() +
theme_bw() +
labs(title="Correlation of Pokémon CP and HP",
caption="Data: https://www.openintro.org/stat/data/?data=pokemon",
x = "CP",
y = "HP") +
theme(
axis.text.x = element_text(size = 10, family = "Trebuchet MS"),
axis.text.y = element_text(size = 10, family = "Trebuchet MS"),
axis.title = element_text(size = 12, family = "Trebuchet MS"),
plot.caption = element_text(size = 7, hjust=0, margin=margin(t=5),
family = "Trebuchet MS"),
plot.title=element_text(family="TrebuchetMS-Bold"),
panel.grid.major = element_blank(),
panel.grid.minor = element_blank()
)
ggsave("CP_HP_corr.png", width = 6, height = 5)
# plot distriubtion of centered PDCP
ggplot(dat1, aes(x = cent_spec_mean_pdcp, fill = species)) +
geom_density() +
geom_vline(xintercept = 0, color = "gray60", linetype = 2) +
# scale_fill_pokemon("infernape") +
scale_fill_colorblind() +
facet_wrap(~species, nrow = 1) +
theme_bw() +
labs(title="Distribution or Percent Change in CP from Species Mean",
caption="Data: https://www.openintro.org/stat/data/?data=pokemon",
x = "Centered PDCP",
y = "Density") +
theme(
strip.background = element_rect(colour = "white", fill = "white"),
strip.text.x = element_text(colour = "black", size = 10, face = "bold",
family = "Trebuchet MS"),
panel.margin = unit(0, "lines"),
panel.border = element_rect(colour = "gray90"),
axis.text.x = element_text(size = 8, family = "Trebuchet MS"),
axis.text.y = element_text(size = 8, family = "Trebuchet MS"),
axis.title = element_text(size = 10, family = "Trebuchet MS"),
plot.caption = element_text(size = 7, hjust=0, margin=margin(t=5),
family = "Trebuchet MS"),
plot.title=element_text(family="TrebuchetMS-Bold"),
legend.position = "none",
panel.grid.major = element_blank(),
panel.grid.minor = element_blank()
)
ggsave("centered_PDCP.png", width = 6, height = 5)
# resdiuals of centered PDCP
ggplot(data = dat1, aes(x = cp, y = cent_spec_mean_pdcp, color = species)) +
geom_hline(yintercept = 0, color = "gray60", linetype = 2) +
geom_point() +
scale_color_colorblind() +
facet_wrap(~species, nrow = 2) +
theme_bw() +
labs(title="Centered Percent Change in CP from Species Mean",
caption="Data: https://www.openintro.org/stat/data/?data=pokemon",
x = "CP",
y = "Density") +
theme(
strip.background = element_rect(colour = "white", fill = "white"),
strip.text.x = element_text(colour = "black", size = 10, face = "bold",
family = "Trebuchet MS"),
panel.margin = unit(0, "lines"),
panel.border = element_rect(colour = "gray90"),
axis.text.x = element_text(size = 8, family = "Trebuchet MS"),
axis.text.y = element_text(size = 8, family = "Trebuchet MS"),
axis.title = element_text(size = 10, family = "Trebuchet MS"),
plot.caption = element_text(size = 7, hjust=0, margin=margin(t=5),
family = "Trebuchet MS"),
plot.title=element_text(family="TrebuchetMS-Bold"),
legend.position = "none",
panel.grid.major = element_blank(),
panel.grid.minor = element_blank()
)
ggsave("centered_PDCP_over_CP.png", width = 6, height = 5)
# trend of mean species PDCP and variation in PDCP
ggplot(data = dat1, aes(x = spec_mean_pdcp, y = spec_mean_pdcp_sd, color = species)) +
geom_line(color = "gray65") +
geom_point(size = 4) +
scale_color_colorblind() +
theme_bw() +
labs(title="Relationship if Mean and Variation in PDCP by Species",
caption="Data: https://www.openintro.org/stat/data/?data=pokemon",
x = "Mean PDCP",
y = "Std Dev of PDCP") +
theme(
axis.text.x = element_text(size = 8, family = "Trebuchet MS"),
axis.text.y = element_text(size = 8, family = "Trebuchet MS"),
axis.title = element_text(size = 10, family = "Trebuchet MS"),
plot.caption = element_text(size = 7, hjust=0, margin=margin(t=5),
family = "Trebuchet MS"),
plot.title=element_text(family="TrebuchetMS-Bold")
)
ggsave("mean_vs_stdDev_PDCP.png", width = 6, height = 5)
# look for universal deviation from mean
quantile(dat1$cent_spec_mean_pdcp) # difference from mean
# distrobution of Centered PDCP for all species
ggplot(dat1, aes(x = cent_spec_mean_pdcp)) +
geom_vline(xintercept = 0, color = "gray60", linetype = 2) +
geom_density(fill = "skyblue") +
scale_x_continuous(breaks = c(seq(-0.45,0.15,0.05))) +
geom_vline(xintercept = 0, color = "gray60", linetype = 2) +
theme_bw() +
labs(title="Distribution of Deviations from Centered PDCP",
caption="Data: https://www.openintro.org/stat/data/?data=pokemon",
x = "Deviation from Centered PDCP",
y = "Density") +
theme(
axis.text.x = element_text(size = 8, family = "Trebuchet MS"),
axis.text.y = element_text(size = 8, family = "Trebuchet MS"),
axis.title = element_text(size = 10, family = "Trebuchet MS"),
plot.caption = element_text(size = 7, hjust=0, margin=margin(t=5),
family = "Trebuchet MS"),
plot.title=element_text(family="TrebuchetMS-Bold")
)
ggsave("dist_dev_from_centered_PDCP.png", width = 6, height = 5)
# distrobution of absolute centered PDCP for all species
ggplot(dat1, aes(x = abs(cent_spec_mean_pdcp))) +
geom_vline(xintercept = 0, color = "gray60", linetype = 2) +
geom_density(fill = "skyblue") +
scale_x_continuous(breaks = c(seq(0,0.5,0.05))) +
geom_vline(xintercept = 0, color = "gray60", linetype = 2) +
theme_bw() +
labs(title="Distribution of Absolute Deviations from Centered PDCP",
caption="Data: https://www.openintro.org/stat/data/?data=pokemon",
x = "Deviation from Centered PDCP",
y = "Density") +
theme(
axis.text.x = element_text(size = 8, family = "Trebuchet MS"),
axis.text.y = element_text(size = 8, family = "Trebuchet MS"),
axis.title = element_text(size = 10, family = "Trebuchet MS"),
plot.caption = element_text(size = 7, hjust=0, margin=margin(t=5),
family = "Trebuchet MS"),
plot.title=element_text(family="TrebuchetMS-Bold")
)
ggsave("ABS_dist_dev_from_centered_PDCP.png", width = 6, height = 5)
### Fit beta parameters to data via MLE
# look at dist. being parametric about this
f3 <- fitdist(abs(dat1$cent_spec_mean_pdcp),"beta")
# plot boostrapped common distributions reatlive to data for comparison
descdist(abs(dat1$cent_spec_mean_pdcp), boot = 1001)
# save via quartz, may not work on Windows
quartz.save(file = "Cullen_Frey_graph.png", height = 5, width = 6,
type = "png", device = dev.cur(), dpi = 300)
# looking at MLE fit of beta to data
summary(f3)
gofstat(f3) # Goodness-of-Fit statistics
plot(f3) # CDF, QQ< PP, Density plot
# sample form beta distribution based on the MLE estimate of [f3]
b <- rbeta(10000,f3$estimate[1], f3$estimate[2])
# plot empirical data versus samples from parameterized beta distribution
ggplot() +
geom_density(data = data.frame(beta_samp = b),
aes(x = beta_samp), fill = "orange") +
geom_density(data = dat1, aes(x = abs(cent_spec_mean_pdcp)), fill = "skyblue3", alpha = 0.55) +
theme_bw() +
scale_x_continuous(breaks = c(seq(0,0.5,0.05))) +
geom_vline(xintercept = 0, color = "gray60", linetype = 2) +
theme_bw() +
labs(title="Comparirson of Empirical and Simulated Beta Distribution",
subtitle="Orange = Beta(0.57, 15.68), Blue = empirical",
caption="Data: https://www.openintro.org/stat/data/?data=pokemon",
x = "Absolute Deviation from Centered PDCP",
y = "Density") +
theme(
axis.text.x = element_text(size = 8, family = "Trebuchet MS"),
axis.text.y = element_text(size = 8, family = "Trebuchet MS"),
axis.title = element_text(size = 10, family = "Trebuchet MS"),
plot.caption = element_text(size = 7, hjust=0, margin=margin(t=5),
family = "Trebuchet MS"),
plot.title=element_text(family="TrebuchetMS-Bold"),
plot.subtitle=element_text(family="TrebuchetMS-Bold"),
panel.grid.major = element_blank(),
panel.grid.minor = element_blank()
)
ggsave("Beta_vs_empirical_dist_centered_PDCP.png", width = 6, height = 5)
###### Simulate plausible evolution paths for each observed Pokémon in sample
### Sequence to sample shape parameters from parametrized beta
### Apply these beta draws to each pokémon species through
### data generating model
draws_from_beta = 100
shape1 = f3$estimate[1]
shape2 = f3$estimate[2]
# Simulate plausible evolution scenarios based on our model
# the melt data for plotting purposes
# sim_pidgey
sim_pidgey <- sim_CP_new(dat1, "Pidgey", shape1, shape2, draws_from_beta)
pidgey_dat <- filter(dat1, species == "Pidgey")
pidgey_melt <- prep_sim_data_for_plot(sim_pidgey, pidgey_dat)
# sim_eevee
sim_eevee <- sim_CP_new(dat1, "Eevee", shape1, shape2, draws_from_beta)
eevee_dat <- filter(dat1, species == "Eevee")
eevee_melt <- prep_sim_data_for_plot(sim_eevee, eevee_dat)
# sim_caterpie
sim_caterpie <- sim_CP_new(dat1, "Caterpie",shape1, shape2, draws_from_beta)
caterpie_dat <- filter(dat1, species == "Caterpie")
caterpie_melt <- prep_sim_data_for_plot(sim_caterpie, caterpie_dat)
# sim_weedle
sim_weedle <- sim_CP_new(dat1, "Weedle",shape1, shape2, draws_from_beta)
weedle_dat <- filter(dat1, species == "Weedle")
weedle_melt <- prep_sim_data_for_plot(sim_weedle, weedle_dat)
# combine all
sim_pokemon <- rbind(pidgey_melt, eevee_melt, caterpie_melt, weedle_melt)
## plot density of empirical post-evolution CP (red line) over
## density of models simulated evolution scenarios (blue lines)
## scenarios based on model and universal beta shape parameter estimates from [f3]
ggplot(data = sim_pokemon, aes(x = value,group = variable, color = model)) +
geom_line(stat="density", size=0.75, alpha=0.25) +
geom_line(data = filter(sim_pokemon, model == "data"),
aes(x = value), stat="density", color = "red", size = 0.25) +
guides(color = guide_legend(override.aes = list(alpha=1))) +
scale_color_manual(values = c("red", "skyblue")) +
facet_grid(species~., scales = "free") +
theme_bw() +
labs(title="Density of Simulated and Empirical Post-Evolution CP per Sepcies",
caption="Data: https://www.openintro.org/stat/data/?data=pokemon",
x = "Post-Evolution CP",
y = "Density") +
theme(
strip.background = element_rect(colour = "white", fill = "white"),
strip.text.x = element_text(colour = "black", size = 10, face = "bold",
family = "Trebuchet MS"),
strip.text.y = element_text(colour = "black", size = 10, face = "bold",
family = "Trebuchet MS"),
panel.margin = unit(0, "lines"),
panel.border = element_rect(colour = "gray90"),
axis.text.x = element_text(size = 8, family = "Trebuchet MS"),
axis.text.y = element_text(size = 8, family = "Trebuchet MS"),
axis.title = element_text(size = 10, family = "Trebuchet MS"),
plot.caption = element_text(size = 7, hjust=0, margin=margin(t=5),
family = "Trebuchet MS"),
plot.title=element_text(family="TrebuchetMS-Bold"),
legend.position = "bottom",
panel.grid.major = element_blank(),
panel.grid.minor = element_blank()
)
ggsave("Species_simulated_CP_New.png", width = 8, height = 8)
# Check errors via RMSE of simulated evolution scenerios vs data
pidgey_sim_beta_rmse <- sim_beta_RMSE(pidgey_dat, pidgey_melt)
eevee_sim_beta_rmse <- sim_beta_RMSE(eevee_dat, eevee_melt)
caterpie_sim_beta_rmse <- sim_beta_RMSE(caterpie_dat, caterpie_melt)
weedle_sim_beta_rmse <- sim_beta_RMSE(weedle_dat, weedle_melt)
# combine together
species_sim_beta_metrics <- rbind(pidgey_sim_beta_rmse, eevee_sim_beta_rmse,
caterpie_sim_beta_rmse, weedle_sim_beta_rmse)
# plot jittered distribution of RMSE for the mean post-evolution CP
# for each observed pokémon by species
ggplot(species_sim_beta_metrics, aes(x = species, y = rmse,
group = species, color = species)) +
geom_jitter(width = 0.3) +
theme_bw() +
scale_color_colorblind() +
labs(title="RMSE Error for Mean Beta Sample Post-evolution CP ",
subtitle = "Mean of n = 100 beta samples for each species observation",
caption="Data: https://www.openintro.org/stat/data/?data=pokemon",
x = "Species",
y = "RMSE") +
theme(
axis.text.x = element_text(size = 8, family = "Trebuchet MS"),
axis.text.y = element_text(size = 8, family = "Trebuchet MS"),
axis.title = element_text(size = 10, family = "Trebuchet MS"),
plot.caption = element_text(size = 7, hjust=0, margin=margin(t=5),
family = "Trebuchet MS"),
plot.title=element_text(family="TrebuchetMS-Bold"),
legend.position = "right"
)
ggsave("Species_beta_CP_New_mean_RMSE.png", width = 6, height = 3)
# table of RMSE results
sim_beta_rmse_table <- group_by(species_sim_beta_metrics, species) %>%
summarise(cp_new_mean = round(mean(cp_new),2),
sim_cp_new_mean = round(mean(means),2),
sd_mean = round(mean(sd),2),
rmse_mean = round(mean(rmse),2),
mae_mean = round(mean(mae),2)) %>%
mutate(pcnt_mae = round(mae_mean / cp_new_mean,3) * 100) %>%
data.frame()
# xtable for blog
print(xtable(sim_beta_rmse_table), type = "html", include.rownames=FALSE)
##### Look at pluasibe variation in beta shape parameters
### instead of one universal MLE estimate, this samples many plausible
### shape paramers via bootstrap, then estiamtes evolution scenerios for
### a range of models that are coherent with the given data
# Boostrap 1001 different pluasible parameters pairs for beta distribution
beta_boot <- bootdist(f3, niter = 1001)
# extract estimates
boot_est <- beta_boot$estim
# get qunatiles of boostrapped parameter estimates and put in table for blog
shp1_boot_est <- round(quantile(boot_est[,1], c(0.025, 0.25, 0.5, 0.75, 0.95)),2)
shp2_boot_est <- round(quantile(boot_est[,2], c(0.025, 0.25, 0.5, 0.75, 0.95)),2)
boot_est_m <- rbind(shp1_boot_est, shp2_boot_est)
print(xtable(boot_est_m), type = "html")
# plot distribution and density of plausible shape parameter pairs
ggplot(beta_boot$estim, aes(x = shape1, y = shape2)) +
geom_point(alpha = 0.95, color = "skyblue3") +
geom_density_2d(color = "gray50", alpha = 0.85) +
theme_bw() +
labs(title="Distribution of Shape Parameters for Beta Distribution",
subtitle="Boostraps, n = 1001",
caption="Data: https://www.openintro.org/stat/data/?data=pokemon",
x = "Shape 1",
y = "Shape 2") +
theme(
axis.text.x = element_text(size = 8, family = "Trebuchet MS"),
axis.text.y = element_text(size = 8, family = "Trebuchet MS"),
axis.title = element_text(size = 10, family = "Trebuchet MS"),
plot.caption = element_text(size = 7, hjust=0, margin=margin(t=5),
family = "Trebuchet MS"),
plot.title=element_text(family="TrebuchetMS-Bold"),
plot.subtitle=element_text(family="TrebuchetMS-Bold")
)
ggsave("Boostrapped_shape_parameters.png", width = 6, height = 5)
# samples of 100 different pidgeys (e.g. PDCPs), under different parameters consistent with the data
# sim_pidgey with params
param_draws = 100
beta_draws = 1
param_dist <- beta_boot$estim
# Simulate beta_draws post-evolution scenerios from param_draws different
# beta distirbution parameterizations (e.g. other plausbile models that incorprate
# uncertainty in limited number ofobservations)
psim_pidgey <- replicate(param_draws, sim_CP_new_by_param(dat1, "Pidgey", param_dist, beta_draws))
psim_eevee <- replicate(param_draws, sim_CP_new_by_param(dat1, "Eevee", param_dist, beta_draws))
psim_weedle <- replicate(param_draws, sim_CP_new_by_param(dat1, "Weedle", param_dist, beta_draws))
psim_caterpie <- replicate(param_draws, sim_CP_new_by_param(dat1, "Caterpie", param_dist, beta_draws))
# prepare psim data for plotting
psim_pidgey_list <- prep_psim_for_plot(psim_pidgey, pidgey_dat)
psim_eevee_list <- prep_psim_for_plot(psim_eevee, eevee_dat)
psim_weedle_list <- prep_psim_for_plot(psim_weedle, weedle_dat)
psim_caterpie_list <- prep_psim_for_plot(psim_caterpie, caterpie_dat)
# Pidgey parameter Sim plot
ggplot() +
geom_point(data = psim_pidgey_list[[1]], aes(x = ord, y = value),
alpha = 0.50, color = "skyblue3") +
geom_point(data = psim_pidgey_list[[2]], aes(x = ord, y = cp_new), color = "red") +
theme_bw() +
labs(title="Range of Simulated Post Evolution CP - Pidgey",
subtitle="Simulated Beta Parameters, n = 100",
caption="Data: https://www.openintro.org/stat/data/?data=pokemon",
x = "Pokémon Ordered by Mean Simulated cp_new",
y = "Post Evolution CP") +
theme(
axis.text.x = element_blank(),
axis.text.y = element_text(size = 8, family = "Trebuchet MS"),
axis.title = element_text(size = 10, family = "Trebuchet MS"),
plot.caption = element_text(size = 7, hjust=0, margin=margin(t=5),
family = "Trebuchet MS"),
plot.title=element_text(family="TrebuchetMS-Bold"),
plot.subtitle=element_text(family="TrebuchetMS-Bold")
)
ggsave("Simulated_shape_parameters_pidgey.png", width = 6, height = 5)
# Eevee parameter Sim plot
ggplot() +
geom_point(data = psim_eevee_list[[1]], aes(x = ord, y = value),
alpha = 0.50, color = "skyblue3") +
geom_point(data = psim_eevee_list[[2]], aes(x = ord, y = cp_new), color = "red") +
theme_bw() +
labs(title="Range of Simulated Post Evolution CP = Eevee",
subtitle="Simulated Beta Parameters, n = 100",
caption="Data: https://www.openintro.org/stat/data/?data=pokemon",
x = "Pokémon Ordered by Mean Simulated cp_new",
y = "Post Evolution CP") +
theme(
axis.text.x = element_blank(),
axis.text.y = element_text(size = 8, family = "Trebuchet MS"),
axis.title = element_text(size = 10, family = "Trebuchet MS"),
plot.caption = element_text(size = 7, hjust=0, margin=margin(t=5),
family = "Trebuchet MS"),
plot.title=element_text(family="TrebuchetMS-Bold"),
plot.subtitle=element_text(family="TrebuchetMS-Bold")
)
ggsave("Simulated_shape_parameters_eevee.png", width = 6, height = 5)
# Weedle parameter Sim plot
ggplot() +
geom_point(data = psim_weedle_list[[1]], aes(x = ord, y = value),
alpha = 0.50, color = "skyblue3") +
geom_point(data = psim_weedle_list[[2]], aes(x = ord, y = cp_new), color = "red") +
theme_bw() +
labs(title="Range of Simulated Post Evolution CP = Weedle",
subtitle="Simulated Beta Parameters, n = 100",
caption="Data: https://www.openintro.org/stat/data/?data=pokemon",
x = "Pokémon Ordered by Mean Simulated cp_new",
y = "Post Evolution CP") +
theme(
axis.text.x = element_blank(),
axis.text.y = element_text(size = 8, family = "Trebuchet MS"),
axis.title = element_text(size = 10, family = "Trebuchet MS"),
plot.caption = element_text(size = 7, hjust=0, margin=margin(t=5),
family = "Trebuchet MS"),
plot.title=element_text(family="TrebuchetMS-Bold"),
plot.subtitle=element_text(family="TrebuchetMS-Bold")
)
ggsave("Simulated_shape_parameters_weedle.png", width = 6, height = 5)
# Caterpie parameter Sim plot
ggplot() +
geom_point(data = psim_caterpie_list[[1]], aes(x = ord, y = value),
alpha = 0.50, color = "skyblue3") +
geom_point(data = psim_caterpie_list[[2]], aes(x = ord, y = cp_new), color = "red") +
theme_bw() +
labs(title="Range of Simulated Post Evolution CP = Caterpie",
subtitle="Simulated Beta Parameters, n = 100",
caption="Data: https://www.openintro.org/stat/data/?data=pokemon",
x = "Pokémon Ordered by Mean Simulated cp_new",
y = "Post Evolution CP") +
theme(
axis.text.x = element_blank(),
axis.text.y = element_text(size = 8, family = "Trebuchet MS"),
axis.title = element_text(size = 10, family = "Trebuchet MS"),
plot.caption = element_text(size = 7, hjust=0, margin=margin(t=5),
family = "Trebuchet MS"),
plot.title=element_text(family="TrebuchetMS-Bold"),
plot.subtitle=element_text(family="TrebuchetMS-Bold")
)
ggsave("Simulated_shape_parameters_caterpie.png", width = 6, height = 5)
# RMSE results from simualted beta parameters
pidgey_sim_rmse <- sim_param_RMSE(pidgey_dat, psim_pidgey_list[[1]])
eevee_sim_rmse <- sim_param_RMSE(eevee_dat, psim_eevee_list[[1]])
caterpie_sim_rmse <- sim_param_RMSE(caterpie_dat, psim_caterpie_list[[1]])
weedle_sim_rmse <- sim_param_RMSE(weedle_dat, psim_weedle_list[[1]])
# combine RMSE
species_sim_metrics <- rbind(pidgey_sim_rmse, eevee_sim_rmse,
caterpie_sim_rmse, weedle_sim_rmse)
# plot RMSE of simulated models parameters
ggplot(species_sim_metrics, aes(x = species, y = rmse,
group = species, color = species)) +
geom_jitter(width = 0.3) +
theme_bw() +
scale_color_colorblind() +
labs(title="RMSE Error for Mean Simulated Post-evolution CP ",
subtitle = "Mean of n = 100 simulations for each species observation",
caption="Data: https://www.openintro.org/stat/data/?data=pokemon",
x = "Species",
y = "RMSE") +
theme(
axis.text.x = element_text(size = 8, family = "Trebuchet MS"),
axis.text.y = element_text(size = 8, family = "Trebuchet MS"),
axis.title = element_text(size = 10, family = "Trebuchet MS"),
plot.caption = element_text(size = 7, hjust=0, margin=margin(t=5),
family = "Trebuchet MS"),
plot.title=element_text(family="TrebuchetMS-Bold"),
legend.position = "right"
)
ggsave("Species_simulated_CP_New_mean_RMSE.png", width = 6, height = 3)
# put metris into a table
sim_rmse_table <- group_by(species_sim_metrics, species) %>%
summarise(cp_new_mean = round(mean(cp_new),2),
sim_cp_new_mean = round(mean(means),2),
sd_mean = round(mean(sd),2),
rmse_mean = round(mean(rmse),2),
mae_mean = round(mean(mae),2)) %>%
mutate(pcnt_mae = round(mae_mean / cp_new_mean,3) * 100) %>%
data.frame()
# print table HTML for blog
print(xtable(sim_rmse_table), type = "html", include.rownames=FALSE)
#### Model parameters in Full Bayes??!!! Are you crazy?!
## estimation using JAGS
library("rjags") # if you dare!!!
# extract abs_cent_spec_mean_pdcp for jags estimation
y <- abs(dat1$cent_spec_mean_pdcp)
# build model text for JAGE
model_string <- "model{
for(i in 1:length(y)) {
y[i] ~ dbeta(alpha, beta) # y distributed as a beta
}
alpha ~ dt(0,1,1)T(0,) # half cauchy prior
beta ~ dt(0,1,1)T(0,) # half cauchy prior
}"
# initialize model with data
model <- jags.model(textConnection(model_string), data = list(y = y), n.chains = 3, n.adapt= 10000)
update(model, 10000); # Burnin for 10000 samples
# sample alpha (shape1) and beta (shape2) from posterior. Can take a while
mcmc_samples <- coda.samples(model, variable.names=c("alpha", "beta"), n.iter=20000)
# plot(mcmc_samples) # plot chain convergence and parameter density
# model parameter distribution stats
summary(mcmc_samples)
# table of stats and print HTML for blog
jags_sum <- summary(mcmc_samples)$quantiles
print(xtable(jags_sum), type = "html")
# plot shape1 and shape 2 parameters
alpha <- mcmc_samples[[1]][,1]
beta <- mcmc_samples[[1]][,2]
# plot correlated parameters with MLE from earlier param estimation and MAP from posterior
ggplot(data = data.frame(alpha = as.numeric(alpha), beta = as.numeric(beta)),
aes(x = alpha, y = beta)) +
geom_point(alpha = 0.15, color = "skyblue3") +
geom_density_2d(color = "gray50", alpha = 0.85) +
geom_point(aes(x = 0.57, y = 15.68), color = "yellow", size = 2) + # MLE
geom_point(aes(x = mean(alpha), y = mean(beta)), color = "red", size = 2) + #MAP
theme_bw() +
labs(title="Posterior of Beta Distribution Shape Parameters",
subtitle="n = 20K samples, red = MAP, yellow = MLE",
caption="Data: https://www.openintro.org/stat/data/?data=pokemon",
x = "Shape 1",
y = "Shape 2") +
theme(
axis.text.x = element_text(size = 8, family = "Trebuchet MS"),
axis.text.y = element_text(size = 8, family = "Trebuchet MS"),
axis.title = element_text(size = 10, family = "Trebuchet MS"),
plot.caption = element_text(size = 7, hjust=0, margin=margin(t=5),
family = "Trebuchet MS"),
plot.title=element_text(family="TrebuchetMS-Bold"),
plot.subtitle=element_text(family="TrebuchetMS-Bold")
)
ggsave("Bayes_simulated_shape_parameters.png", width = 6, height = 5)
# prior, post, plot
post_plot <- data.frame(value = c(alpha, beta),
parameter = rep(c("alpha", "beta"), each = length(alpha)),
model = "posterior")
alpha_half_cauchy <- rt(nrow(post_plot),1,1)
alpha_half_cauchy <- alpha_half_cauchy[alpha_half_cauchy > 0]
beta_half_cauchy <- rt(nrow(post_plot),1,1)
beta_half_cauchy <- beta_half_cauchy[beta_half_cauchy > 0]
prior_plot <- data.frame(value = c(alpha_half_cauchy,beta_half_cauchy),
parameter = rep(c("alpha", "beta"),
times = c(length(alpha_half_cauchy), length(beta_half_cauchy))),
model = "prior"
)
jags_plot_dat <- rbind(post_plot, prior_plot)
ggplot(data = filter(jags_plot_dat, parameter == "alpha"),
aes(x = value, group = model, color = model)) +
geom_density() +
scale_x_continuous(limits = c(0,10), breaks = seq(0,10,0.5)) +
scale_color_manual(values = c("orange", "purple")) +
theme_bw() +
labs(title="Prior and Posterior of Shape1 Parameter of Beta Distribution",
subtitle="Half Cauchy as T distribution (1,1) > 0",
caption="20K Samples from JAGS model",
x = "Shape1 Parameter",
y = "Density") +
theme(
axis.text.x = element_text(size = 8, family = "Trebuchet MS"),
axis.text.y = element_text(size = 8, family = "Trebuchet MS"),
axis.title = element_text(size = 10, family = "Trebuchet MS"),
plot.caption = element_text(size = 7, hjust=0, margin=margin(t=5),
family = "Trebuchet MS"),
plot.title=element_text(family="TrebuchetMS-Bold"),
plot.subtitle=element_text(family="TrebuchetMS-Bold")
)
ggsave("alpha_plot.png", width = 6, height = 5)
ggplot(data = filter(jags_plot_dat, parameter == "beta"),
aes(x = value, group = model, color = model)) +
geom_density() +
scale_x_continuous(limits = c(0,100), breaks = seq(0,100,5)) +
scale_color_manual(values = c("orange", "purple")) +
theme_bw() +
labs(title="Prior and Posterior of Shape2 Parameter of Beta Distribution",
subtitle="Half Cauchy as T distribution (1,1) > 0",
caption="20K Samples from JAGS model",
x = "Shape2 Parameter",
y = "Density") +
theme(
axis.text.x = element_text(size = 8, family = "Trebuchet MS"),
axis.text.y = element_text(size = 8, family = "Trebuchet MS"),
axis.title = element_text(size = 10, family = "Trebuchet MS"),
plot.caption = element_text(size = 7, hjust=0, margin=margin(t=5),
family = "Trebuchet MS"),
plot.title=element_text(family="TrebuchetMS-Bold"),
plot.subtitle=element_text(family="TrebuchetMS-Bold")
)
ggsave("beta_plot.png", width = 6, height = 5)
###### Generalized predictions
base = 10
sim_reps = 10000
length_out = 10
# don't need all of this for here, but good summary for all species
spec_ranges <- group_by(dat1, species) %>%
summarise(min_cp = min(cp),
max_cp = max(cp),
mean_cp = mean(cp),
mean_cp_new = mean(cp_new),
mean_pdcp = mean(pcnt_d_cp),
mean_abs_cent_spec_mean_pdcp = mean(abs(cent_spec_mean_pdcp))) %>%
data.frame()
print(xtable(spec_ranges), type = "html", include.rownames=FALSE)
## Predict siualted Pokémon species based on empirical range
# Pidgey
pidgey_cp_range <- sim_new_species(spec_ranges, pidgey_dat)
predicted_pidgey <- pedict_simulated_cp_new(pidgey_cp_range, f3$estimate[1], f3$estimate[2], reps = sim_reps)
pidgey_cp_new_ranges <- make_pred_cp_plot_dat(predicted_pidgey, pidgey_cp_range)
# Eevee
eevee_cp_range <- sim_new_species(spec_ranges, eevee_dat)
predicted_eevee <- pedict_simulated_cp_new(eevee_cp_range, f3$estimate[1], f3$estimate[2], reps = sim_reps)
eevee_cp_new_ranges <- make_pred_cp_plot_dat(predicted_eevee, eevee_cp_range)
# Weedle
weedle_cp_range <- sim_new_species(spec_ranges, weedle_dat)
predicted_weedle <- pedict_simulated_cp_new(weedle_cp_range, f3$estimate[1], f3$estimate[2], reps = sim_reps)
weedle_cp_new_ranges <- make_pred_cp_plot_dat(predicted_weedle, weedle_cp_range)
# Caterpie
caterpie_cp_range <- sim_new_species(spec_ranges, caterpie_dat)
predicted_caterpie <- pedict_simulated_cp_new(caterpie_cp_range, f3$estimate[1], f3$estimate[2], reps = sim_reps)
caterpie_cp_new_ranges <- make_pred_cp_plot_dat(predicted_caterpie, caterpie_cp_range)
# line plot with error ribbons, ok, but not too flashy
ggplot(data = pidgey_cp_new_ranges, aes(x = cp, y = cp_new_mean)) +
geom_ribbon(aes(ymin = cp_new_min, ymax = cp_new_max),alpha=0.4, fill = "skyblue") +
geom_ribbon(aes(ymin = cp_new_5, ymax = cp_new_95),alpha=0.75, fill = "skyblue3") +
geom_point(data = pidgey_dat, aes(x = cp, y = cp_new), color = "gray30", alpha = 0.60) +
geom_line(color = "firebrick", size = 0.75) +
scale_x_continuous(breaks = seq(min(pidgey_cp_new_ranges$cp),
max(pidgey_cp_new_ranges$cp),by = 20)) +
scale_y_continuous(breaks = seq(min(pidgey_cp_new_ranges$cp_new_min),
max(pidgey_cp_new_ranges$cp_new_max),by = 50)) +
theme_bw() +
labs(title="Range of Model Predictions for CP_new - Pidgey",
subtitle="Mean, 95% prediction interval & min/max prediction",
caption="Data: https://www.openintro.org/stat/data/?data=pokemon",
x = "Pre-Evolution CP",
y = "Modeled Post-Evolution CP") +
theme(
axis.text.x = element_text(size = 8, family = "Trebuchet MS"),
axis.text.y = element_text(size = 8, family = "Trebuchet MS"),
axis.title = element_text(size = 10, family = "Trebuchet MS"),
plot.caption = element_text(size = 7, hjust=0, margin=margin(t=5),
family = "Trebuchet MS"),
plot.title=element_text(family="TrebuchetMS-Bold"),
plot.subtitle=element_text(family="TrebuchetMS-Bold")
)
ggsave("Predicted_cp_new_pidgey.png", width = 6, height = 5)
# eevee plot
ggplot(data = eevee_cp_new_ranges, aes(x = cp, y = cp_new_mean)) +
geom_ribbon(aes(ymin = cp_new_min, ymax = cp_new_max),alpha=0.4, fill = "skyblue") +
geom_ribbon(aes(ymin = cp_new_5, ymax = cp_new_95),alpha=0.75, fill = "skyblue3") +
geom_point(data = eevee_dat, aes(x = cp, y = cp_new), color = "gray30", alpha = 0.50) +
geom_line(color = "firebrick", size = 0.75) +
scale_x_continuous(breaks = seq(min(eevee_cp_new_ranges$cp),
max(eevee_cp_new_ranges$cp),by = 20)) +
scale_y_continuous(breaks = seq(min(eevee_cp_new_ranges$cp_new_min),
max(eevee_cp_new_ranges$cp_new_max),by = 50)) +
theme_bw() +
labs(title="Range of Model Predictions for CP_new - Eevee",
subtitle="Mean, 95% prediction interval & min/max prediction",
caption="Data: https://www.openintro.org/stat/data/?data=pokemon",
x = "Pre-Evolution CP",
y = "Modeled Post-Evolution CP") +
theme(
axis.text.x = element_text(size = 8, family = "Trebuchet MS"),
axis.text.y = element_text(size = 8, family = "Trebuchet MS"),
axis.title = element_text(size = 10, family = "Trebuchet MS"),
plot.caption = element_text(size = 7, hjust=0, margin=margin(t=5),
family = "Trebuchet MS"),
plot.title=element_text(family="TrebuchetMS-Bold"),
plot.subtitle=element_text(family="TrebuchetMS-Bold")
)
ggsave("Predicted_cp_new_eevee.png", width = 6, height = 5)
# weedle plot
ggplot(data = weedle_cp_new_ranges, aes(x = cp, y = cp_new_mean)) +
geom_ribbon(aes(ymin = cp_new_min, ymax = cp_new_max),alpha=0.4, fill = "skyblue") +
geom_ribbon(aes(ymin = cp_new_5, ymax = cp_new_95),alpha=0.75, fill = "skyblue3") +
geom_point(data = weedle_dat, aes(x = cp, y = cp_new), color = "gray30", alpha = 0.50) +
geom_line(color = "firebrick", size = 0.75) +
scale_x_continuous(breaks = seq(min(weedle_cp_new_ranges$cp),
max(weedle_cp_new_ranges$cp),by = 20)) +
scale_y_continuous(breaks = seq(min(weedle_cp_new_ranges$cp_new_min),
max(weedle_cp_new_ranges$cp_new_max),by = 50)) +
theme_bw() +
labs(title="Range of Model Predictions for CP_new - weedle",
subtitle="Mean, 95% prediction interval & min/max prediction",
caption="Data: https://www.openintro.org/stat/data/?data=pokemon",
x = "Pre-Evolution CP",
y = "Modeled Post-Evolution CP") +
theme(
axis.text.x = element_text(size = 8, family = "Trebuchet MS"),
axis.text.y = element_text(size = 8, family = "Trebuchet MS"),
axis.title = element_text(size = 10, family = "Trebuchet MS"),
plot.caption = element_text(size = 7, hjust=0, margin=margin(t=5),
family = "Trebuchet MS"),
plot.title=element_text(family="TrebuchetMS-Bold"),
plot.subtitle=element_text(family="TrebuchetMS-Bold")
)
ggsave("Predicted_cp_new_weedle.png", width = 6, height = 5)
# caterpie plot
ggplot(data = caterpie_cp_new_ranges, aes(x = cp, y = cp_new_mean)) +
geom_ribbon(aes(ymin = cp_new_min, ymax = cp_new_max),alpha=0.4, fill = "skyblue") +
geom_ribbon(aes(ymin = cp_new_5, ymax = cp_new_95),alpha=0.75, fill = "skyblue3") +
geom_point(data = caterpie_dat, aes(x = cp, y = cp_new), color = "gray30", alpha = 0.50) +
geom_line(color = "firebrick", size = 0.75) +
scale_x_continuous(breaks = seq(min(caterpie_cp_new_ranges$cp),
max(caterpie_cp_new_ranges$cp),by = 20)) +
scale_y_continuous(breaks = seq(min(caterpie_cp_new_ranges$cp_new_min),
max(caterpie_cp_new_ranges$cp_new_max),by = 50)) +
theme_bw() +
labs(title="Range of Model Predictions for CP_new - Caterpie",
subtitle="Mean, 95% prediction interval & min/max prediction",
caption="Data: https://www.openintro.org/stat/data/?data=pokemon",
x = "Pre-Evolution CP",
y = "Modeled Post-Evolution CP") +
theme(
axis.text.x = element_text(size = 8, family = "Trebuchet MS"),
axis.text.y = element_text(size = 8, family = "Trebuchet MS"),
axis.title = element_text(size = 10, family = "Trebuchet MS"),
plot.caption = element_text(size = 7, hjust=0, margin=margin(t=5),
family = "Trebuchet MS"),
plot.title=element_text(family="TrebuchetMS-Bold"),
plot.subtitle=element_text(family="TrebuchetMS-Bold")
)
ggsave("Predicted_cp_new_caterpie.png", width = 6, height = 5)
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment