|
| 1 | +## ----setup, include = FALSE--------------------------------------------------- |
| 2 | +knitr::opts_chunk$set( |
| 3 | + message = TRUE, |
| 4 | + warning = TRUE, |
| 5 | + collapse = TRUE, |
| 6 | + comment = "#>" |
| 7 | +) |
| 8 | + |
| 9 | +## ----------------------------------------------------------------------------- |
| 10 | +sim <- function() { |
| 11 | + fake <- data.frame((x <- rnorm(100, 100, 20)), (y <- rnorm(100, 80, 20))) |
| 12 | + intervals <- t.test(x = x, y = y, data = fake, conf.level = .95)$conf.int[] |
| 13 | +} |
| 14 | + |
| 15 | +set.seed(1031) |
| 16 | + |
| 17 | +z <- replicate(100, sim(), simplify = FALSE) |
| 18 | + |
| 19 | +df <- data.frame(do.call(rbind, z)) |
| 20 | +df$studynumber <- (1:length(z)) |
| 21 | +intrvl.limit <- c("lower.limit", "upper.limit", "studynumber") |
| 22 | +colnames(df) <- intrvl.limit |
| 23 | +df$point <- ((df$lower.limit + df$upper.limit) / 2) |
| 24 | +df$covered <- (df$lower.limit <= 20 & 20 <= df$upper.limit) |
| 25 | +df$coverageprob <- ((as.numeric(table(df$covered)[2]) / nrow(df) * 100)) |
| 26 | + |
| 27 | +library(ggplot2) |
| 28 | + |
| 29 | + |
| 30 | +ggplot(data = df, aes(x = studynumber, y = point, ymin = lower.limit, ymax = upper.limit)) + |
| 31 | + geom_pointrange(mapping = aes(color = covered), size = .40) + |
| 32 | + geom_hline(yintercept = 20, lty = 1, color = "red", alpha = 0.5) + |
| 33 | + coord_flip() + |
| 34 | + labs( |
| 35 | + title = "Simulated 95% Intervals", |
| 36 | + x = "Study Number", |
| 37 | + y = "Estimate", |
| 38 | + subtitle = "Population Parameter is 20" |
| 39 | + ) + |
| 40 | + theme_bw() + # use a white background |
| 41 | + theme(legend.position = "none") + |
| 42 | + annotate( |
| 43 | + geom = "text", x = 102, y = 30, |
| 44 | + label = "Coverage (%) =", size = 2.5, color = "black" |
| 45 | + ) + |
| 46 | + annotate( |
| 47 | + geom = "text", x = 102, y = 35, |
| 48 | + label = df$coverageprob, size = 2.5, color = "black" |
| 49 | + ) |
| 50 | + |
| 51 | +## ----echo=TRUE---------------------------------------------------------------- |
| 52 | +library(concurve) |
| 53 | +library(rstan) |
| 54 | +library(rstanarm) |
| 55 | +library(ggplot2) |
| 56 | +library(cowplot) |
| 57 | +library(bayesplot) |
| 58 | +library(scales) |
| 59 | + |
| 60 | +## ----echo=TRUE---------------------------------------------------------------- |
| 61 | +GroupA <- rnorm(50) |
| 62 | +GroupB <- rnorm(50) |
| 63 | +RandomData <- data.frame(GroupA, GroupB) |
| 64 | +model_freq <- lm(GroupA ~ GroupB, data = RandomData) |
| 65 | + |
| 66 | +## ----results = 'hide', message = FALSE---------------------------------------- |
| 67 | +rstan_options(auto_write = TRUE) |
| 68 | + |
| 69 | +# Using flat prior |
| 70 | +model_bayes <- stan_lm(GroupA ~ GroupB, |
| 71 | + data = RandomData, prior = NULL, |
| 72 | + iter = 5000, warmup = 1000, chains = 4 |
| 73 | +) |
| 74 | + |
| 75 | +## ----echo=TRUE---------------------------------------------------------------- |
| 76 | +randomframe <- curve_gen(model_freq, "GroupB", steps = 10000) |
| 77 | + |
| 78 | +(function1 <- ggcurve(type = "c", randomframe[[1]], nullvalue = TRUE)) |
| 79 | + |
| 80 | +color_scheme_set("teal") |
| 81 | + |
| 82 | +function2 <- mcmc_dens(model_bayes, pars = "GroupB") + |
| 83 | + ggtitle("Posterior Distribution") + |
| 84 | + labs(subtitle = "Function Displays the Full Posterior Distribution", x = "Range of Values", y = "Posterior Probability") + |
| 85 | + scale_y_continuous(breaks = c(0, 0.30, 0.60, 0.90, 1.20, 1.50, 1.80, 2.10, 2.40, 2.70, 3.0)) |
| 86 | + |
| 87 | + |
| 88 | +(breaks1 <- c(0, 0.30, 0.60, 0.90, 1.20, 1.50, 1.80, 2.10, 2.40, 2.70, 3.0)) |
| 89 | + |
| 90 | +(adjustment <- function(x) { |
| 91 | + x / 3 |
| 92 | +}) |
| 93 | + |
| 94 | +(labels <- adjustment(breaks1)) |
| 95 | + |
| 96 | +breaks <- labels |
| 97 | +labels1 <- labels |
| 98 | + |
| 99 | +(function3 <- mcmc_dens(model_bayes, pars = "GroupB") + |
| 100 | + ggtitle("Posterior Distribution") + |
| 101 | + labs(subtitle = "Function Displays the Full Posterior Distribution", x = "Range of Values", y = "Posterior Probability") + |
| 102 | + scale_x_continuous(expand = c(0, 0), breaks = scales::pretty_breaks(n = 10)) + |
| 103 | + scale_y_continuous(expand = c(0, 0), breaks = waiver(), labels = waiver(), n.breaks = 10, limits = c(0, 3.25)) + |
| 104 | + yaxis_text(on = TRUE) + |
| 105 | + yaxis_ticks(on = TRUE) + |
| 106 | + annotate("segment", |
| 107 | + x = 0, xend = 0, y = 0, yend = 3, |
| 108 | + color = "#990000", alpha = 0.4, size = .75, linetype = 3 |
| 109 | + )) |
| 110 | + |
| 111 | +## ----echo=TRUE, fig.height=9, fig.width=7------------------------------------- |
| 112 | +plot_grid(function1, function3, ncol = 1, align = "v") |
| 113 | + |
| 114 | +## ----echo=TRUE---------------------------------------------------------------- |
| 115 | + |
| 116 | +GroupA <- rnorm(500, mean = 2) |
| 117 | +GroupB <- rnorm(500, mean = 1) |
| 118 | +RandomData <- data.frame(GroupA, GroupB) |
| 119 | +model_freq <- lm(GroupA ~ GroupB, data = RandomData) |
| 120 | + |
| 121 | +## ----results = 'hide', message = FALSE---------------------------------------- |
| 122 | + |
| 123 | +# Using flat prior |
| 124 | +model_bayes <- stan_lm(GroupA ~ GroupB, |
| 125 | + data = RandomData, prior = NULL, |
| 126 | + iter = 5000, warmup = 1000, chains = 4 |
| 127 | +) |
| 128 | + |
| 129 | +## ----echo=TRUE---------------------------------------------------------------- |
| 130 | + |
| 131 | +randomframe <- curve_gen(model_freq, "GroupB", steps = 10000) |
| 132 | + |
| 133 | +(function1 <- ggcurve(type = "c", randomframe[[1]], nullvalue = TRUE)) |
| 134 | + |
| 135 | +color_scheme_set("teal") |
| 136 | + |
| 137 | +function2 <- mcmc_dens(model_bayes, pars = "GroupB") + |
| 138 | + ggtitle("Posterior Distribution") + |
| 139 | + labs(subtitle = "Function Displays the Full Posterior Distribution", x = "Range of Values", y = "Posterior Probability") + |
| 140 | + scale_y_continuous(breaks = c(0, 0.30, 0.60, 0.90, 1.20, 1.50, 1.80, 2.10, 2.40, 2.70, 3.0)) |
| 141 | + |
| 142 | + |
| 143 | +(breaks1 <- c(0, 0.30, 0.60, 0.90, 1.20, 1.50, 1.80, 2.10, 2.40, 2.70, 3.0)) |
| 144 | + |
| 145 | +(adjustment <- function(x) { |
| 146 | + x / 3 |
| 147 | +}) |
| 148 | + |
| 149 | +(labels <- adjustment(breaks1)) |
| 150 | + |
| 151 | +breaks <- labels |
| 152 | +labels1 <- labels |
| 153 | + |
| 154 | +(function3 <- mcmc_dens(model_bayes, pars = "GroupB") + |
| 155 | + ggtitle("Posterior Distribution") + |
| 156 | + labs(subtitle = "Function Displays the Full Posterior Distribution", x = "Range of Values", y = "Posterior Probability") + |
| 157 | + scale_x_continuous(expand = c(0, 0), breaks = scales::pretty_breaks(n = 10)) + |
| 158 | + scale_y_continuous(expand = c(0, 0), breaks = waiver(), labels = waiver(), n.breaks = 10, limits = c(0, 9)) + |
| 159 | + yaxis_text(on = TRUE) + |
| 160 | + yaxis_ticks(on = TRUE) + |
| 161 | + annotate("segment", |
| 162 | + x = 0, xend = 0, y = 0, yend = 9, |
| 163 | + color = "#990000", alpha = 0.4, size = .75, linetype = 3 |
| 164 | + )) |
| 165 | + |
| 166 | +## ----echo=TRUE, fig.height=9, fig.width=7------------------------------------- |
| 167 | +plot_grid(function1, function3, ncol = 1, align = "v") |
| 168 | + |
| 169 | +## ----results = 'hide', message = FALSE---------------------------------------- |
| 170 | + |
| 171 | +data(kidiq) |
| 172 | + |
| 173 | +# flat prior |
| 174 | + |
| 175 | +post1 <- stan_lm(kid_score ~ mom_hs, |
| 176 | + data = kidiq, prior = NULL, |
| 177 | + seed = 12345 |
| 178 | +) |
| 179 | + |
| 180 | +## ----------------------------------------------------------------------------- |
| 181 | +post2 <- lm(kid_score ~ mom_hs, data = kidiq) |
| 182 | + |
| 183 | +df3 <- curve_gen(post2, "mom_hs") |
| 184 | + |
| 185 | +(function99 <- ggcurve(df3[[1]])) |
| 186 | + |
| 187 | +summary(post1) |
| 188 | + |
| 189 | +color_scheme_set("teal") |
| 190 | + |
| 191 | +(function101 <- mcmc_areas(post1, pars = "mom_hs", point_est = "none", prob = 1, prob_outer = 1, area_method = "equal height") + |
| 192 | + ggtitle("Posterior Distribution") + |
| 193 | + labs(subtitle = "Function Displays the Full Posterior Distribution", x = "Range of Values", y = "Posterior Probability") + |
| 194 | + yaxis_text(on = TRUE) + |
| 195 | + yaxis_ticks(on = TRUE)) |
| 196 | + |
| 197 | +## ----echo=TRUE, fig.height=9, fig.width=7------------------------------------- |
| 198 | +cowplot::plot_grid(function99, function101, ncol = 1, align = "v") |
| 199 | + |
| 200 | +## ----------------------------------------------------------------------------- |
| 201 | +citation("concurve") |
| 202 | +citation("ggplot2") |
| 203 | +citation("rstan") |
| 204 | +citation("rstanarm") |
| 205 | +citation("cowplot") |
| 206 | +citation("scales") |
| 207 | +citation("bayesplot") |
| 208 | + |
0 commit comments