week 2: linear model and causal inference

geocentric models

Annoucements and such

  • We will start using brms today! Install this package now, if you haven’t already.
install.packages(c("brms","tidybayes"))

Workspace setup

library(here)
library(tidyverse)
library(cowplot)
library(brms)
library(tidybayes)
library(patchwork)

Workflow

  1. State a clear question.

  2. Sketch your causal assumptions.

  3. Use the sketch to define a generative model.

  4. Use the generative model to build an estimator.

  5. Profit.

Model “recipes”

  1. Recognize a set of variables to work with. (Data and parameters.)
  2. Define each variable either in terms of the other variables OR in terms of a probability distribution.
  3. The combination of variables and their probability distributions defines a joint generative model that can be used to simulate hypothetical observations and analyze real ones.

Here’s an example:

\[\begin{align*} y_i &\sim \text{Normal}(\mu_i,\sigma) \\ \mu_i &= \beta x_i \\ \beta &\sim \text{Normal}(0,10) \\ \sigma &\sim \text{Exponential}(1) \\ x_i &\sim \text{Normal}(0,1) \\ \end{align*}\]

Model for globe-tossing

Here’s the model for last week’s globe-tossing experiment:

\[\begin{align*} W &\sim \text{Binomial}(N,p) \\ p &\sim \text{Uniform}(0,1) \\ \end{align*}\]

  • \(W\) is the observed count of water.
  • \(N\) is the total number of tosses.
  • \(p\) is the proportion of water on the globe.

The whole model can be read as:

The count \(W\) is distributed binomially with sample size \(N\) and probability \(p\). The prior for \(p\) is assumed to be uniform between 0 and 1.

Model for globe-tossing

Here’s the model for last week’s globe-tossing experiment:

\[\begin{align*} W &\sim \text{Binomial}(N,p) \\ p &\sim \text{Uniform}(0,1) \\ \end{align*}\]

Estimating the posterior using brms

Last week, we used grid approximation to estimate the posterior distribution. In the video you watched for today, McElreath moves on to something called QUADRATIC APPROXIMATION. It’s good to understand what that’s doing, but you and I are moving right along to MCMC. We won’t go into the details of what’s happening for a few weeks, but let’s start with the code to estimate the model.

Let’s say we tossed the globe 9 times and observed 6 waters:

m1 <-
  brm(data = list(w = 6),
      family = binomial(link = "identity"),
      w | trials(9) ~ 0 + Intercept,
      prior(uniform(0, 1), class = b),
      iter = 5000, warmup = 1000, seed = 3, chains=1,
      file = here("files/models/m21.1"))
1
Data can be a data frame or a list. Just make sure the variable names match what’s in your formula.
2
How you assume your outcome variable is distributed.
3
The formula for your outcome.
4
Priors for every parameter in your model. Here, we only have one parameter, so we only need 1 prior. This happens to be a flat prior between 0 and 1.
5
Some choices about how we want our model to run. We’ll go more into this later.
6
These models can take a long time to run. You can automatically save the output to a file; when you do this, the next time you run this code, it won’t actually estimate the model, but will instead pull the output from your stated file. Be WARNED: if you change the data or the model code, it will NOT restimate your model until you delete the file.

Sampling from the posterior

Grid approximation gave us the calculated probability of each possible value of our parameter, \(p\). But our method of conducting bayes will no longer give us such a neat solution. Here’s how you get the posterior distribution for \(p\):

samples_from_post = as_draws_df(m1)
samples_from_post
# A draws_df: 4000 iterations, 1 chains, and 3 variables
   b_Intercept lprior lp__
1         0.46      0 -2.1
2         0.49      0 -1.9
3         0.56      0 -1.5
4         0.66      0 -1.3
5         0.71      0 -1.3
6         0.62      0 -1.3
7         0.68      0 -1.3
8         0.55      0 -1.5
9         0.63      0 -1.3
10        0.67      0 -1.3
# ... with 3990 more draws
# ... hidden reserved variables {'.chain', '.iteration', '.draw'}
samples_from_post %>%  
  ggplot(aes(x=b_Intercept)) +
  geom_density(fill = "grey", color = "white") +
  labs(x="Proportion water")

Posterior predictive distribution

ppd = posterior_predict(m1)
dim(ppd)
[1] 4000    1
ppd
        [,1]
   [1,]    4
   [2,]    4
   [3,]    8
   [4,]    4
   [5,]    8
   [6,]    4
   [7,]    6
   [8,]    3
   [9,]    5
  [10,]    7
  [11,]    4
  [12,]    6
  [13,]    6
  [14,]    9
  [15,]    6
  [16,]    8
  [17,]    6
  [18,]    7
  [19,]    6
  [20,]    9
  [21,]    9
  [22,]    7
  [23,]    2
  [24,]    6
  [25,]    4
  [26,]    5
  [27,]    6
  [28,]    5
  [29,]    3
  [30,]    5
  [31,]    7
  [32,]    4
  [33,]    4
  [34,]    6
  [35,]    3
  [36,]    5
  [37,]    8
  [38,]    2
  [39,]    4
  [40,]    6
  [41,]    7
  [42,]    7
  [43,]    7
  [44,]    4
  [45,]    4
  [46,]    7
  [47,]    7
  [48,]    3
  [49,]    7
  [50,]    8
  [51,]    6
  [52,]    5
  [53,]    6
  [54,]    5
  [55,]    4
  [56,]    5
  [57,]    5
  [58,]    2
  [59,]    5
  [60,]    7
  [61,]    4
  [62,]    3
  [63,]    6
  [64,]    6
  [65,]    7
  [66,]    6
  [67,]    6
  [68,]    5
  [69,]    6
  [70,]    8
  [71,]    6
  [72,]    5
  [73,]    7
  [74,]    4
  [75,]    7
  [76,]    6
  [77,]    6
  [78,]    7
  [79,]    9
  [80,]    3
  [81,]    4
  [82,]    6
  [83,]    4
  [84,]    6
  [85,]    7
  [86,]    4
  [87,]    5
  [88,]    7
  [89,]    4
  [90,]    7
  [91,]    7
  [92,]    4
  [93,]    1
  [94,]    4
  [95,]    7
  [96,]    7
  [97,]    7
  [98,]    8
  [99,]    6
 [100,]    9
 [101,]    5
 [102,]    6
 [103,]    7
 [104,]    8
 [105,]    3
 [106,]    6
 [107,]    6
 [108,]    3
 [109,]    4
 [110,]    5
 [111,]    7
 [112,]    5
 [113,]    5
 [114,]    6
 [115,]    5
 [116,]    4
 [117,]    6
 [118,]    5
 [119,]    7
 [120,]    4
 [121,]    7
 [122,]    4
 [123,]    5
 [124,]    7
 [125,]    6
 [126,]    8
 [127,]    6
 [128,]    6
 [129,]    5
 [130,]    9
 [131,]    7
 [132,]    8
 [133,]    6
 [134,]    6
 [135,]    5
 [136,]    8
 [137,]    6
 [138,]    4
 [139,]    8
 [140,]    6
 [141,]    8
 [142,]    7
 [143,]    3
 [144,]    6
 [145,]    4
 [146,]    9
 [147,]    7
 [148,]    4
 [149,]    6
 [150,]    5
 [151,]    4
 [152,]    7
 [153,]    7
 [154,]    7
 [155,]    6
 [156,]    3
 [157,]    5
 [158,]    8
 [159,]    5
 [160,]    3
 [161,]    7
 [162,]    3
 [163,]    7
 [164,]    6
 [165,]    3
 [166,]    5
 [167,]    7
 [168,]    5
 [169,]    6
 [170,]    7
 [171,]    8
 [172,]    4
 [173,]    4
 [174,]    5
 [175,]    5
 [176,]    5
 [177,]    7
 [178,]    9
 [179,]    5
 [180,]    8
 [181,]    6
 [182,]    8
 [183,]    9
 [184,]    6
 [185,]    6
 [186,]    7
 [187,]    7
 [188,]    7
 [189,]    1
 [190,]    3
 [191,]    7
 [192,]    2
 [193,]    5
 [194,]    8
 [195,]    5
 [196,]    6
 [197,]    7
 [198,]    5
 [199,]    8
 [200,]    8
 [201,]    8
 [202,]    5
 [203,]    7
 [204,]    9
 [205,]    5
 [206,]    6
 [207,]    5
 [208,]    4
 [209,]    5
 [210,]    6
 [211,]    6
 [212,]    6
 [213,]    9
 [214,]    7
 [215,]    6
 [216,]    8
 [217,]    8
 [218,]    7
 [219,]    5
 [220,]    2
 [221,]    7
 [222,]    6
 [223,]    3
 [224,]    9
 [225,]    6
 [226,]    4
 [227,]    7
 [228,]    7
 [229,]    6
 [230,]    7
 [231,]    7
 [232,]    4
 [233,]    8
 [234,]    8
 [235,]    7
 [236,]    9
 [237,]    4
 [238,]    6
 [239,]    3
 [240,]    3
 [241,]    8
 [242,]    5
 [243,]    7
 [244,]    7
 [245,]    7
 [246,]    8
 [247,]    4
 [248,]    6
 [249,]    2
 [250,]    6
 [251,]    5
 [252,]    9
 [253,]    8
 [254,]    8
 [255,]    8
 [256,]    6
 [257,]    7
 [258,]    7
 [259,]    6
 [260,]    6
 [261,]    7
 [262,]    7
 [263,]    7
 [264,]    7
 [265,]    8
 [266,]    5
 [267,]    6
 [268,]    6
 [269,]    6
 [270,]    6
 [271,]    6
 [272,]    7
 [273,]    5
 [274,]    8
 [275,]    6
 [276,]    8
 [277,]    7
 [278,]    9
 [279,]    9
 [280,]    6
 [281,]    5
 [282,]    5
 [283,]    1
 [284,]    1
 [285,]    5
 [286,]    2
 [287,]    1
 [288,]    9
 [289,]    7
 [290,]    8
 [291,]    6
 [292,]    8
 [293,]    8
 [294,]    3
 [295,]    2
 [296,]    4
 [297,]    8
 [298,]    6
 [299,]    3
 [300,]    4
 [301,]    8
 [302,]    5
 [303,]    7
 [304,]    4
 [305,]    8
 [306,]    7
 [307,]    4
 [308,]    7
 [309,]    6
 [310,]    8
 [311,]    2
 [312,]    5
 [313,]    6
 [314,]    5
 [315,]    3
 [316,]    7
 [317,]    7
 [318,]    7
 [319,]    7
 [320,]    6
 [321,]    4
 [322,]    4
 [323,]    5
 [324,]    3
 [325,]    7
 [326,]    8
 [327,]    7
 [328,]    6
 [329,]    6
 [330,]    6
 [331,]    5
 [332,]    8
 [333,]    7
 [334,]    8
 [335,]    8
 [336,]    8
 [337,]    9
 [338,]    4
 [339,]    5
 [340,]    3
 [341,]    3
 [342,]    5
 [343,]    6
 [344,]    8
 [345,]    4
 [346,]    5
 [347,]    4
 [348,]    8
 [349,]    7
 [350,]    7
 [351,]    9
 [352,]    5
 [353,]    4
 [354,]    3
 [355,]    4
 [356,]    7
 [357,]    6
 [358,]    5
 [359,]    4
 [360,]    4
 [361,]    3
 [362,]    5
 [363,]    7
 [364,]    6
 [365,]    5
 [366,]    7
 [367,]    4
 [368,]    3
 [369,]    9
 [370,]    9
 [371,]    4
 [372,]    0
 [373,]    2
 [374,]    3
 [375,]    1
 [376,]    4
 [377,]    3
 [378,]    4
 [379,]    3
 [380,]    6
 [381,]    4
 [382,]    7
 [383,]    6
 [384,]    9
 [385,]    7
 [386,]    6
 [387,]    9
 [388,]    7
 [389,]    6
 [390,]    6
 [391,]    9
 [392,]    7
 [393,]    6
 [394,]    8
 [395,]    5
 [396,]    7
 [397,]    6
 [398,]    7
 [399,]    7
 [400,]    7
 [401,]    5
 [402,]    5
 [403,]    7
 [404,]    5
 [405,]    6
 [406,]    8
 [407,]    6
 [408,]    7
 [409,]    7
 [410,]    9
 [411,]    6
 [412,]    8
 [413,]    5
 [414,]    4
 [415,]    7
 [416,]    4
 [417,]    6
 [418,]    7
 [419,]    7
 [420,]    8
 [421,]    8
 [422,]    6
 [423,]    4
 [424,]    6
 [425,]    7
 [426,]    3
 [427,]    2
 [428,]    9
 [429,]    5
 [430,]    6
 [431,]    5
 [432,]    6
 [433,]    5
 [434,]    6
 [435,]    9
 [436,]    8
 [437,]    6
 [438,]    7
 [439,]    5
 [440,]    4
 [441,]    6
 [442,]    5
 [443,]    8
 [444,]    4
 [445,]    6
 [446,]    1
 [447,]    4
 [448,]    5
 [449,]    7
 [450,]    3
 [451,]    5
 [452,]    6
 [453,]    7
 [454,]    6
 [455,]    6
 [456,]    3
 [457,]    8
 [458,]    8
 [459,]    6
 [460,]    6
 [461,]    7
 [462,]    9
 [463,]    8
 [464,]    5
 [465,]    7
 [466,]    2
 [467,]    5
 [468,]    6
 [469,]    8
 [470,]    4
 [471,]    5
 [472,]    7
 [473,]    7
 [474,]    7
 [475,]    6
 [476,]    7
 [477,]    5
 [478,]    6
 [479,]    7
 [480,]    8
 [481,]    6
 [482,]    5
 [483,]    7
 [484,]    5
 [485,]    4
 [486,]    5
 [487,]    7
 [488,]    5
 [489,]    6
 [490,]    5
 [491,]    7
 [492,]    7
 [493,]    1
 [494,]    9
 [495,]    4
 [496,]    3
 [497,]    7
 [498,]    8
 [499,]    7
 [500,]    6
 [501,]    8
 [502,]    6
 [503,]    7
 [504,]    7
 [505,]    6
 [506,]    5
 [507,]    5
 [508,]    7
 [509,]    6
 [510,]    5
 [511,]    4
 [512,]    4
 [513,]    3
 [514,]    6
 [515,]    7
 [516,]    5
 [517,]    6
 [518,]    8
 [519,]    6
 [520,]    2
 [521,]    1
 [522,]    3
 [523,]    3
 [524,]    2
 [525,]    4
 [526,]    9
 [527,]    6
 [528,]    3
 [529,]    4
 [530,]    3
 [531,]    8
 [532,]    7
 [533,]    7
 [534,]    6
 [535,]    5
 [536,]    2
 [537,]    8
 [538,]    9
 [539,]    8
 [540,]    5
 [541,]    5
 [542,]    5
 [543,]    7
 [544,]    8
 [545,]    6
 [546,]    8
 [547,]    4
 [548,]    6
 [549,]    7
 [550,]    2
 [551,]    7
 [552,]    9
 [553,]    4
 [554,]    5
 [555,]    4
 [556,]    5
 [557,]    2
 [558,]    6
 [559,]    3
 [560,]    6
 [561,]    5
 [562,]    4
 [563,]    2
 [564,]    8
 [565,]    5
 [566,]    5
 [567,]    4
 [568,]    6
 [569,]    8
 [570,]    7
 [571,]    4
 [572,]    5
 [573,]    5
 [574,]    2
 [575,]    7
 [576,]    4
 [577,]    7
 [578,]    6
 [579,]    5
 [580,]    8
 [581,]    6
 [582,]    7
 [583,]    5
 [584,]    7
 [585,]    2
 [586,]    1
 [587,]    2
 [588,]    3
 [589,]    0
 [590,]    3
 [591,]    5
 [592,]    3
 [593,]    4
 [594,]    4
 [595,]    5
 [596,]    7
 [597,]    3
 [598,]    5
 [599,]    5
 [600,]    4
 [601,]    6
 [602,]    7
 [603,]    1
 [604,]    8
 [605,]    8
 [606,]    8
 [607,]    6
 [608,]    3
 [609,]    4
 [610,]    6
 [611,]    4
 [612,]    6
 [613,]    4
 [614,]    6
 [615,]    6
 [616,]    2
 [617,]    6
 [618,]    7
 [619,]    9
 [620,]    7
 [621,]    4
 [622,]    4
 [623,]    7
 [624,]    6
 [625,]    3
 [626,]    6
 [627,]    6
 [628,]    5
 [629,]    5
 [630,]    1
 [631,]    7
 [632,]    6
 [633,]    5
 [634,]    5
 [635,]    4
 [636,]    6
 [637,]    8
 [638,]    3
 [639,]    2
 [640,]    2
 [641,]    6
 [642,]    3
 [643,]    8
 [644,]    2
 [645,]    5
 [646,]    6
 [647,]    6
 [648,]    3
 [649,]    3
 [650,]    6
 [651,]    8
 [652,]    8
 [653,]    6
 [654,]    6
 [655,]    4
 [656,]    8
 [657,]    8
 [658,]    6
 [659,]    8
 [660,]    6
 [661,]    6
 [662,]    3
 [663,]    3
 [664,]    3
 [665,]    8
 [666,]    4
 [667,]    8
 [668,]    6
 [669,]    6
 [670,]    7
 [671,]    9
 [672,]    5
 [673,]    6
 [674,]    3
 [675,]    2
 [676,]    6
 [677,]    6
 [678,]    5
 [679,]    8
 [680,]    6
 [681,]    8
 [682,]    7
 [683,]    3
 [684,]    5
 [685,]    2
 [686,]    6
 [687,]    6
 [688,]    5
 [689,]    6
 [690,]    8
 [691,]    6
 [692,]    2
 [693,]    3
 [694,]    5
 [695,]    6
 [696,]    6
 [697,]    5
 [698,]    2
 [699,]    3
 [700,]    4
 [701,]    4
 [702,]    6
 [703,]    1
 [704,]    9
 [705,]    5
 [706,]    7
 [707,]    5
 [708,]    7
 [709,]    8
 [710,]    5
 [711,]    3
 [712,]    8
 [713,]    6
 [714,]    5
 [715,]    5
 [716,]    5
 [717,]    2
 [718,]    4
 [719,]    6
 [720,]    7
 [721,]    8
 [722,]    5
 [723,]    4
 [724,]    6
 [725,]    6
 [726,]    7
 [727,]    4
 [728,]    6
 [729,]    6
 [730,]    5
 [731,]    5
 [732,]    6
 [733,]    5
 [734,]    4
 [735,]    4
 [736,]    7
 [737,]    8
 [738,]    7
 [739,]    6
 [740,]    6
 [741,]    8
 [742,]    7
 [743,]    5
 [744,]    3
 [745,]    7
 [746,]    8
 [747,]    5
 [748,]    5
 [749,]    6
 [750,]    5
 [751,]    5
 [752,]    8
 [753,]    5
 [754,]    6
 [755,]    4
 [756,]    6
 [757,]    7
 [758,]    7
 [759,]    3
 [760,]    7
 [761,]    4
 [762,]    5
 [763,]    7
 [764,]    8
 [765,]    8
 [766,]    9
 [767,]    7
 [768,]    6
 [769,]    7
 [770,]    7
 [771,]    8
 [772,]    4
 [773,]    6
 [774,]    6
 [775,]    9
 [776,]    6
 [777,]    5
 [778,]    9
 [779,]    4
 [780,]    7
 [781,]    6
 [782,]    6
 [783,]    4
 [784,]    5
 [785,]    2
 [786,]    4
 [787,]    4
 [788,]    4
 [789,]    5
 [790,]    8
 [791,]    8
 [792,]    8
 [793,]    8
 [794,]    9
 [795,]    8
 [796,]    7
 [797,]    7
 [798,]    7
 [799,]    9
 [800,]    9
 [801,]    9
 [802,]    8
 [803,]    6
 [804,]    4
 [805,]    6
 [806,]    4
 [807,]    5
 [808,]    5
 [809,]    7
 [810,]    5
 [811,]    6
 [812,]    7
 [813,]    5
 [814,]    7
 [815,]    8
 [816,]    4
 [817,]    7
 [818,]    6
 [819,]    9
 [820,]    7
 [821,]    5
 [822,]    6
 [823,]    5
 [824,]    5
 [825,]    5
 [826,]    7
 [827,]    6
 [828,]    5
 [829,]    6
 [830,]    6
 [831,]    2
 [832,]    7
 [833,]    5
 [834,]    6
 [835,]    4
 [836,]    6
 [837,]    5
 [838,]    4
 [839,]    8
 [840,]    6
 [841,]    7
 [842,]    8
 [843,]    8
 [844,]    6
 [845,]    6
 [846,]    3
 [847,]    3
 [848,]    3
 [849,]    3
 [850,]    6
 [851,]    8
 [852,]    9
 [853,]    5
 [854,]    2
 [855,]    8
 [856,]    6
 [857,]    5
 [858,]    7
 [859,]    7
 [860,]    6
 [861,]    7
 [862,]    3
 [863,]    5
 [864,]    8
 [865,]    8
 [866,]    9
 [867,]    8
 [868,]    6
 [869,]    5
 [870,]    6
 [871,]    4
 [872,]    5
 [873,]    7
 [874,]    4
 [875,]    1
 [876,]    1
 [877,]    4
 [878,]    3
 [879,]    3
 [880,]    3
 [881,]    3
 [882,]    7
 [883,]    7
 [884,]    2
 [885,]    2
 [886,]    5
 [887,]    8
 [888,]    7
 [889,]    5
 [890,]    4
 [891,]    4
 [892,]    5
 [893,]    5
 [894,]    7
 [895,]    7
 [896,]    8
 [897,]    5
 [898,]    7
 [899,]    5
 [900,]    5
 [901,]    7
 [902,]    6
 [903,]    7
 [904,]    5
 [905,]    6
 [906,]    6
 [907,]    7
 [908,]    9
 [909,]    8
 [910,]    8
 [911,]    9
 [912,]    8
 [913,]    7
 [914,]    7
 [915,]    6
 [916,]    6
 [917,]    3
 [918,]    8
 [919,]    5
 [920,]    6
 [921,]    6
 [922,]    4
 [923,]    3
 [924,]    7
 [925,]    6
 [926,]    8
 [927,]    2
 [928,]    1
 [929,]    4
 [930,]    5
 [931,]    6
 [932,]    8
 [933,]    7
 [934,]    9
 [935,]    4
 [936,]    2
 [937,]    4
 [938,]    9
 [939,]    6
 [940,]    8
 [941,]    5
 [942,]    6
 [943,]    6
 [944,]    8
 [945,]    6
 [946,]    7
 [947,]    4
 [948,]    4
 [949,]    5
 [950,]    4
 [951,]    3
 [952,]    4
 [953,]    6
 [954,]    5
 [955,]    3
 [956,]    5
 [957,]    5
 [958,]    4
 [959,]    9
 [960,]    9
 [961,]    7
 [962,]    7
 [963,]    3
 [964,]    2
 [965,]    2
 [966,]    7
 [967,]    4
 [968,]    3
 [969,]    6
 [970,]    3
 [971,]    7
 [972,]    8
 [973,]    5
 [974,]    4
 [975,]    4
 [976,]    4
 [977,]    6
 [978,]    8
 [979,]    9
 [980,]    7
 [981,]    6
 [982,]    9
 [983,]    5
 [984,]    4
 [985,]    3
 [986,]    5
 [987,]    4
 [988,]    6
 [989,]    4
 [990,]    2
 [991,]    2
 [992,]    9
 [993,]    7
 [994,]    6
 [995,]    6
 [996,]    7
 [997,]    3
 [998,]    7
 [999,]    4
[1000,]    3
 [ reached 'max' / getOption("max.print") -- omitted 3000 rows ]
data.frame(obs = ppd) %>% 
  ggplot(aes(x=obs)) +
  geom_histogram() +
  labs(x="Observed water (out of 9)")

Sampling parameters from the prior

Oops, we’ve jumped ahead of ourselves! Best practice is to simulate values from your prior first and check to see if those priors are reasonable.

m1p <-
  brm(data = list(w = 6),                            
      family = binomial(link = "identity"),          
      w | trials(9) ~ 0 + Intercept,                 
      prior(beta(1, 1), class = b, lb = 0, ub = 1),  
      iter = 5000, warmup = 1000, seed = 3, chains=1,          
      sample_prior = "only")
samples_from_prior = as_draws_df(m1p)
samples_from_prior %>% 
  ggplot(aes(x=b_Intercept)) +
  geom_density(fill = "grey", color = "white") +
  labs(x="Proportion water", title="Prior")

Simulating observations from the prior

We may also want the PRIOR PREDICTIVE DISTRIBUTION which is the expected observiations given our prior.

prior_pd = posterior_predict(m1p)
data.frame(obs = prior_pd) %>% 
  ggplot(aes(x=obs)) +
  geom_histogram() +
  labs(x="Observed water (out of 9)")

Simulating from your priors – prior predictive simulation – is an essential part of modeling. This allows you to see what your choices imply about the data. You’ll be able to diagnose bad choices.

an aside about learning in R

At this point in the course, I’m going to start throwing a lot of code at you. Do I expect you to memorize this code? Of course not.

Do you need to understand every single thing that’s happening in the code? Nope.

But, you’ll learn a lot by taking the time to figure out what’s happening in a code chunk. Class time will frequently include exercises where I ask you to adapt code I’ve shared in the slides to a new dataset or to answer a new problem. When doing so, go back through the old code and figure out what’s going on. Run the code one line at a time. Always observe the output and take some time to look at the object that was created or modified. Here are some functions that will be extremely useful:

str() # what kind of object is this? what is its structure?
dim() # what are the dimensions (rows/columns) of this object
head() # give me the first bit of this object
str(prior_pd)
 int [1:4000, 1] 5 8 7 5 6 8 7 2 3 2 ...
 - attr(*, "dimnames")=List of 2
  ..$ : NULL
  ..$ : NULL
dim(prior_pd)
[1] 4000    1
head(prior_pd)
     [,1]
[1,]    5
[2,]    8
[3,]    7
[4,]    5
[5,]    6
[6,]    8

Continous outcomes

The globe tossing example is cute and easy to work with, but let’s move towards the kinds of variables we more frequently work with. Let’s create a model for some outcome, \(y\) that is continuous.

\[\begin{align*} y_i &\sim \text{Normal}(\mu, \sigma) \\ \mu &\sim \text{Normal}(0, 10) \\ \sigma &\sim \text{Uniform}(0, 5) \end{align*}\]

set.seed(9)
y = rnorm(n = 31, mean = 4, sd = .5)
m2 = brm(
  data = list(y=y),
  family = gaussian,
  y ~ 1,
  prior = c(prior( normal(0,10), class=Intercept),
            prior( uniform(0,5), class=sigma)),  
      iter = 5000, warmup = 1000, seed = 3, chains=1,
  file = here("files/models/m21.2")
)

An example: weight and height

Using the Howell data (don’t load the rethinking package because it can interfere with brms).

data("Howell1", package = "rethinking")
d <- Howell1
str(d)
'data.frame':   544 obs. of  4 variables:
 $ height: num  152 140 137 157 145 ...
 $ weight: num  47.8 36.5 31.9 53 41.3 ...
 $ age   : num  63 63 65 41 51 35 32 27 19 54 ...
 $ male  : int  1 0 0 1 0 1 0 1 0 1 ...
library(measurements)
d$height <- conv_unit(d$height, from = "cm", to = "feet")
d$weight <- conv_unit(d$weight, from = "kg", to = "lbs")
rethinking::precis(d)
             mean         sd      5.5%    94.5%      histogram
height  4.5362072  0.9055921  2.661042   5.4375      ▁▁▁▁▂▂▇▇▁
weight 78.5079631 32.4502290 20.636856 120.1583 ▁▂▃▂▂▁▁▃▅▇▇▃▂▁
age    29.3443934 20.7468882  1.000000  66.1350      ▇▅▅▃▅▂▂▁▁
male    0.4724265  0.4996986  0.000000   1.0000     ▇▁▁▁▁▁▁▁▁▇
d2 <- d[ d$age >= 18, ]

exercise

Write a mathematical model for the weights in this data set. (Don’t worry about predicting from other variables yet.)

solution

\[\begin{align*} w &\sim \text{Normal}(\mu, \sigma) \\ \mu &\sim \text{Normal}(130, 20) \\ \sigma &\sim \text{Uniform}(0, 25) \\ \end{align*}\]

\[\begin{align*} w &\sim \text{Normal}(\mu, \sigma) \\ \mu &\sim \text{Normal}(130, 20) \\ \sigma &\sim \text{Uniform}(0, 25) \\ \end{align*}\]

exercise

Simulate from your priors (parameters values and prior predictive values).

solution

Sample from your priors:

m3p = brm(
  data = d2,
  family = gaussian,
  weight ~ 1,
  prior = c(prior( normal(130,20), class=Intercept),
            prior( uniform(0,25), class=sigma, lb=0, ub=25)),  
      iter = 5000, warmup = 1000, seed = 3, chains=1,
  sample_prior = "only")

Sampling parameter estimates for the Intercept.

pairs(m3p)

This is a different (shorter) way to plot your posterior. Good things: it automatically includes the scatterplot so you can see the implications of how these parameters correlate. Bad things: not customizable and not useable when you have a lot of parameters.

Simulate values of weight.

prior_pd = posterior_predict(m3p)
dim(prior_pd)
[1] 4000  352
as.data.frame(prior_pd) %>% 
  pivot_longer(everything()) %>% 
  ggplot(aes(x=value)) +
  geom_histogram() +
  labs(x="Expected observed weights (based on prior)")

Another shorter way:

pp_check(m3p)

Fit the model

m3 = brm(
  data = d2,
  family = gaussian,
  weight ~ 1,
  prior = c(prior( normal(130,20), class=Intercept),
            prior( uniform(0,25), class=sigma, lb=0, ub=25)),  
      iter = 5000, warmup = 1000, seed = 3, chains=1,
  file = here("files/models/m21.3"))
posterior_summary(m3)
                Estimate  Est.Error         Q2.5        Q97.5
b_Intercept    99.221909 0.75720356    97.751074   100.737488
sigma          14.291987 0.55198322    13.254363    15.428183
Intercept      99.221909 0.75720356    97.751074   100.737488
lprior         -8.318377 0.05827127    -8.433538    -8.203915
lp__        -1441.294276 1.02695551 -1444.235605 -1440.292326
pairs(m3)
Code 
posterior_predict(m3) %>% 
  as.data.frame() %>% 
  pivot_longer(everything()) %>% 
  ggplot(aes(x=value)) +
  geom_density(fill = "grey", color = "white") +
  geom_density( aes(x = weight), data=d2, inherit.aes = F) 
pp_check(m3)

Adding in a linear component

We might assume that height and weight are associated with each other. Indeed, within our sample:

plot(d2$weight ~ d2$height)

exercise

Update your mathematical model to incorporate height.

\[\begin{align*} w_i &\sim \text{Normal}(\mu_i, \sigma) \\ \mu_i &= \alpha + \beta h_i \\ \alpha &\sim \text{Normal}(??, ??) \\ \beta &\sim \text{Normal}(0, 25) \\ \sigma &\sim \text{Uniform}(0, 25) \\ \end{align*}\]

exercise

Update your mathematical model to incorporate height.

\[\begin{align*} w_i &\sim \text{Normal}(\mu_i, \sigma) \\ \mu_i &= \alpha + \beta (h_i - \bar{h}) \\ \alpha &\sim \text{Normal}(130, 20) \\ \beta &\sim \text{Normal}(0, 25) \\ \sigma &\sim \text{Uniform}(0, 25) \\ \end{align*}\]

\[\begin{align*} w_i &\sim \text{Normal}(\mu_i, \sigma) \\ \mu_i &= \alpha + \beta (h_i - \bar{h}) \\ \alpha &\sim \text{Normal}(130, 20) \\ \beta &\sim \text{Normal}(0, 25) \\ \sigma &\sim \text{Uniform}(0, 25) \\ \end{align*}\]

To update our brms code:

d2$height_c = d2$height - mean(d2$height)
  
m4p = brm(
  data = d2,
  family = gaussian,
  weight ~ 1 + height_c,
  prior = c(prior( normal(130,20), class=Intercept),
            prior( normal(0,25),   class=b),
            prior( uniform(0,25),  class=sigma, lb=0, ub=25)),  
      iter = 5000, warmup = 1000, seed = 3, chains=1,
  sample_prior = "only")
set.seed(9)
samples_from_prior = as_draws_df(m4p)
str(samples_from_prior)
draws_df [4,000 × 9] (S3: draws_df/draws/tbl_df/tbl/data.frame)
 $ b_Intercept: num [1:4000] 126 142 148 118 115 ...
 $ b_height_c : num [1:4000] -16.652 -9.498 -12.363 0.248 -10.576 ...
 $ sigma      : num [1:4000] 2.63 19.37 16.55 4.3 12.92 ...
 $ Intercept  : num [1:4000] 126 142 148 118 115 ...
 $ lprior     : num [1:4000] -11.5 -11.5 -11.8 -11.5 -11.7 ...
 $ lp__       : num [1:4000] -10.66 -10.04 -10.08 -10.18 -9.83 ...
 $ .chain     : int [1:4000] 1 1 1 1 1 1 1 1 1 1 ...
 $ .iteration : int [1:4000] 1 2 3 4 5 6 7 8 9 10 ...
 $ .draw      : int [1:4000] 1 2 3 4 5 6 7 8 9 10 ...
Code 
d2 %>% 
  ggplot(aes(x=height_c, y=weight)) +
  geom_blank() +
  geom_abline( aes(intercept=b_Intercept, slope=b_height_c), 
               data=samples_from_prior[1:50, ],
               alpha=.3) +
  scale_x_continuous(name = "height(feet)", 
                     breaks=seq(4,6,by=.5)-mean(d2$height),
                     labels=seq(4,6,by=.5))

Describe in words what’s wrong with our priors.

Slope should not be negative. How can we fix this?

Could use a uniform distribution bounded by 0.

m4p = brm(
  data = d2,
  family = gaussian,
  weight ~ 1 + height_c,
  prior = c(prior( normal(130,20), class=Intercept),
            prior( uniform(0,25),   class=b),
            prior( uniform(0,25),  class=sigma, lb=0, ub=25)),  
      iter = 5000, warmup = 1000, seed = 3, chains=1,
  sample_prior = "only")

exercise

Fit the new weight model to the data.

solution

m4 = brm(
  data = d2,
  family = gaussian,
  weight ~ 1 + height_c,
  prior = c(prior( normal(130,20), class=Intercept),
            prior( uniform(0,25),   class=b),
            prior( uniform(0,25),  class=sigma, lb=0, ub=25)),  
      iter = 5000, warmup = 1000, seed = 3, chains=1,
  file=here("files/models/m21.4"))
posterior_summary(m4)
               Estimate Est.Error        Q2.5       Q97.5
b_Intercept    99.21508 0.5154139    98.17317   100.20955
b_height_c     24.74661 0.2609169    24.05150    24.99355
sigma          10.36090 0.3898531     9.65234    11.12228
Intercept      99.21508 0.5154139    98.17317   100.20955
lprior        -11.53739 0.0396881   -11.61861   -11.46176
lp__        -1332.15882 1.3936194 -1335.90639 -1330.52002

exercise

Draw lines from the posterior distribution and plot with the data.

solution

Code 
set.seed(9)
samples_from_posterior = as_draws_df(m4)
d2%>% 
  ggplot(aes(x=height_c, y=weight)) +
  geom_point(size=.5) +
  geom_abline( aes(intercept=b_Intercept, slope=b_height_c), 
               data=samples_from_posterior[1:50, ],
               alpha=.3,
               color="#1c5253") +
  scale_x_continuous(name = "height(feet)", 
                     breaks=seq(4,6,by=.5)-mean(d2$height),
                     labels=seq(4,6,by=.5))

A side note: a major concern or critique of Bayesian analysis is that the subjectivity of the priors allow for nefarious behavior. “Putting our thumbs on the scale,” so to speak. But priors are quickly overwhelmed by data. Case in point:

m4e = brm(
  data = d2,
  family = gaussian,
  weight ~ 1 + height_c,
  prior = c(prior( normal(130,20), class=Intercept),
            prior( normal(-5,5),   class=b),
            prior( uniform(0,25),  class=sigma, lb=0, ub=25)),  
      iter = 5000, warmup = 1000, seed = 3, chains=1,
  file=here("files/models/m21.4e"))
posterior_summary(m4e)
                Estimate Est.Error         Q2.5       Q97.5
b_Intercept    99.197733 0.5035653    98.214092   100.15932
b_height_c     35.775818 1.8832462    32.177572    39.50070
sigma           9.529429 0.3687178     8.839379    10.27233
Intercept      99.197733 0.5035653    98.214092   100.15932
lprior        -44.172476 3.0738968   -50.456409   -38.49994
lp__        -1334.629214 1.1849503 -1337.645046 -1333.23509

You’ll only really get into trouble with uniform priors that have a boundary, if true population parameter is outside your boundary. A good rule of thumb is to avoid the uniform distribution. We’ll cover other options for priors for \(\sigma\) in future lectures, but as a preview, the exponential distribution works very well for this!