Skip to content

Instantly share code, notes, and snippets.

Avatar
🦔

Richard McElreath rmcelreath

🦔
View GitHub Profile
@rmcelreath
rmcelreath / MCpostpred.r
Created Oct 16, 2021
Monte Carlo posterior predictive simulation
View MCpostpred.r
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
rmcelreath / RFDT7.R
Created Jun 29, 2021
Code example for RFDT
View RFDT7.R
# 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
rmcelreath / RFDT6.R
Created Jun 29, 2021
Code eamples for RFDT
View RFDT6.R
# 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
rmcelreath / RFDT5.R
Created Jun 29, 2021
Code examples for RFDT
View RFDT5.R
# 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
rmcelreath / RFDT4.R
Last active Jun 29, 2021
Code examples for RFDT
View RFDT4.R
# 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
rmcelreath / RFDT3.R
Created Jun 18, 2021
Code examples for RFDT
View RFDT3.R
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
rmcelreath / RFDT2.R
Created Jun 15, 2021
Code examples for RFDT
View RFDT2.R
# 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 ) )
@rmcelreath
rmcelreath / RFDT1.R
Last active Jun 18, 2021
Code examples for RFDT
View RFDT1.R
set.seed(1908)
N <- 200 # number of pairs
U <- rnorm(N) # 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 ) # change the 0 to turn on causal influence of mom
@rmcelreath
rmcelreath / data_caring.md
Last active Dec 31, 2020
Data sharing, caring and access conventions
View data_caring.md

Data Sharing, Caring, and Access

Most major scientific journals now require data archiving after publication. Many require it for review. So anthropology must figure out how to responsibly provide access to data. Luckily, other fields have already invested a lot in figured this out.

The important thing to realize is that providing access does not mean complete loss of control, either for the communities that provide the data or for researchers. For example, access can be restricted to researchers affiliated with research institutions, and limits can be placed on rights of reuse. But some level of access is necessary, if for no other reason than to allow other scholars the ability to verify our analyses. At the same time, researchers who want access to data must respect that data are political and continued collaboration with communities that provide data requires privacy and limits on reuse.

A standard data sharing memorandum grants all collaborators non-exclusive use of the data for scientific analysis

@rmcelreath
rmcelreath / nested_effects_ulam.R
Last active Apr 21, 2020
Nested varying effects in ulam example
View nested_effects_ulam.R
# nested varying effects in ulam
# simulate an example
# people in nations, allowing varying among people to vary by nation
N_nations <- 10
N_id_per_nation <- sample( 5:100 , size=N_nations , replace=TRUE )
N_id <- sum(N_id_per_nation)
# variation among individuals in each nation