rmcelreath

# prior - likelihood conflict

library(rethinking)

yobs <- 0

mtt <- ulam(
    alist(
        y ~ dstudent(2,mu,1),
        mu ~ dstudent(2,10,1)

rmcelreath

# distribution of null p-values in binomial glm

library(rethinking)

# t test
# pval has uniform distribution under null

f <- function(N=10) {
    y1 <- rnorm(N)
    y2 <- rnorm(N)

rmcelreath

# simulation of simple Kopernican model showing epicycle

library(rethinking)
library(animation)

my_circle <- function(x=0,y=0,r=1,angle=0,...) {
    a <- seq(angle, angle + 2 * pi, length = 360)
    lines( r*cos(a)+x , r*sin(a)+y , ... )
}
# plot(NULL,xlim=c(-1,1),ylim=c(-1,1)); my_circle(0,0,0.5,lty=2)

rmcelreath

library(rethinking)
library(animation)

ani.saveqz <- function (list , dpi=150) {   
    if (missing(list)) 
        list = animation:::.ani.env$.images
    lapply(1:length(list), function(x) {
        dev.hold()
        replayPlot( list[[x]] )
        ani.pause()

rmcelreath

# version that marginalizes out the missing data
flbi_plus <- ulam(
    alist(
        c(M,D) ~ multi_normal( c(mu,nu) , Rho , Sigma ),
        mu <- a1 + b*B1,
        nu <- a2 + b*B2 + m*M,
        c(a1,a2,b,m) ~ normal( 0 , 0.5 ),
        Rho ~ lkj_corr( 2 ),
        Sigma ~ exponential( 1 )
    ), data=dat , chains=4 , cores=4 , cmdstan=TRUE )

rmcelreath

# more exotic example - no instrument (B1 -> D), but have two measures of U
set.seed(1908)
N <- 200 # number of pairs
U <- rnorm(N,0,1) # simulate confound
V <- rnorm(N,U,1)
W <- rnorm(N,U,1)
# birth order and family sizes
B1 <- rbinom(N,size=1,prob=0.5) # 50% first borns
M <- rnorm( N , 2*B1 + U )
B2 <- rbinom(N,size=1,prob=0.5)

rmcelreath

# simulate interventions

# B1 -> D (should be b*m)
post <- extract.samples(flbi)
quantile( with(post,b*m) )

# now do the same with simulation

# B1 = 0
# B1 -> M

rmcelreath

# full-luxury bayesian inference
library(rethinking)
library(cmdstanr)
dat <- list(N=N,M=M,D=D,B1=B1,B2=B2)
set.seed(1908)
flbi <- ulam(
    alist(
        # mom model
            M ~ normal( mu , sigma ),
            mu <- a1 + b*B1 + k*U[i],

rmcelreath

set.seed(1908)
N <- 200 # number of pairs
U <- rnorm(N,0,1) # simulate confounds
# birth order and family sizes
B1 <- rbinom(N,size=1,prob=0.5) # 50% first borns
M <- rnorm( N , 2*B1 + U )
B2 <- rbinom(N,size=1,prob=0.5)
D <- rnorm( N , 2*B2 + U + 0*M )

f <- function(data,indices) with(data[indices,], cov(B1,D) / cov(B1,M) )

rmcelreath

# fit the two regression models
summary( lm( D ~ M ) )
summary( lm( D ~ M + B1 + B2 ) )

# compare the models with AIC
AIC( lm( D ~ M ) )
AIC( lm( D ~ M + B1 + B2 ) )