rmcelreath

# sim confounding of age --> happiness by controlling for marital status, when in truth age + happiness --> marriage
# assume:
# (1) no relationship between happiness and age
# (2) 5 happiest people each year get married (after age 20)

N_years <- 200
N_born_per_year <- 40
N_marry_per_year <- 5
mortality <- 0.05
max_age <- 100

rmcelreath

## code to reproduce Figure 2.5 on page 30 of Statistical Rethinking

library(rethinking)

col1 <- "black"
col2 <- col.alpha( "black" , 0.7 )

# show posterior as data comes in
# 1 indicates 'water'; 0 indicates 'land'
d <- c(1,0,1,1,1,0,1,0,1) # length 9

rmcelreath

###################
# example workflow
# 1 build model
# 2 simulate data
# 3 analyze data

################
# 1 build model

# example context: tradeoff between egg number and egg survival

rmcelreath

// contest model
data{
    int N_choices;
    int N_items;
    int N_raters;
    int y[N_choices];
    int item1[N_choices];
    int item2[N_choices];
    int rater[N_choices];
    real x[N_items];      // item feature

rmcelreath

// simple reinforcement model
data{
    int N;          // total number of observations
    int N_trials;   // num trials per individual
    int N_individuals;
    int trial[N];   // trial number in individual
    int y[N];       // choices
    real P[N,2];       // payoffs to each option
    int id[N];      // individual ID
}

rmcelreath

# script to illustrate collider bias described in Breen 2018 doi:10.1093/esr/jcy037

library(rethinking)

N <- 200  # number of grandparent-parent-child triads
b_gp <- 1 # direct effect of G on P
b_gc <- 0 # direct effect of G on C
b_pc <- 1 # direct effect of P on C
b_U <- 2  # direct effect of U on P and C

rmcelreath

library(rethinking)

HMC2 <- function (U, grad_U, epsilon, L, current_q , ... ) {
  q = current_q
  p = rnorm(length(q),0,1)  # independent standard normal variates
  current_p = p
  # Make a half step for momentum at the beginning
  p = p - epsilon * grad_U( q , ... ) / 2
  # Alternate full steps for position and momentum
  qtraj <- matrix(NA,nrow=L+1,ncol=length(q))

rmcelreath

HMC2 <- function (U, grad_U, epsilon, L, current_q , ... ) {
  q = current_q
  p = rnorm(length(q),0,1)  # independent standard normal variates
  current_p = p
  # Make a half step for momentum at the beginning
  p = p - epsilon * grad_U( q , ... ) / 2
  # Alternate full steps for position and momentum
  qtraj <- matrix(NA,nrow=L+1,ncol=length(q))
  ptraj <- qtraj
  qtraj[1,] <- current_q

rmcelreath

# R script to accompany EHBEA 2019 presentation "Causal Salad in Human Evolution & Ecology"

# to install latest rethinking package, see:
# https://github.com/rmcelreath/statrethinking_winter2019
library(rethinking)

# Example 1: Grandparent -> child
# what is influence of grandparent presence on grandchild welfare?
# confounded by parent condition (wealth), 
#  because grandparents live with wealthier children

rmcelreath

# install rethinking library with instructions here:
# https://github.com/rmcelreath/statrethinking_winter2019

library(rethinking)
library(dagitty)
library(shape)

# the structural model
# W: wages
# Q: quarter of year of birth