Skip to content

Instantly share code, notes, and snippets.

@rubenarslan
rubenarslan / thin_brmsfit.R
Last active July 27, 2025 20:12
thin brmsfit (brms fit objects)
#' Thin a brmsfit object post-hoc
#'
#' This function thins a brmsfit object after fitting. It checks the backend
#' used to fit the model ('rstan' or 'cmdstanr') and applies the correct
#' thinning method, as they store samples in different internal structures.
#'
#' @param fit A brmsfit object.
#' @param thin_by An integer factor by which to thin the post-warmup chains.
#' @return A new, thinned brmsfit object.
thin_brmsfit <- function(fit, thin_by) {
@rubenarslan
rubenarslan / partner-preference-sim.R
Created October 15, 2024 20:41
simulation code with memory leak
# Load necessary libraries for data manipulation, visualization, and parallel processing
library(tidyverse) # Provides a collection of R packages for data manipulation
library(furrr) # Parallel processing using purrr
library(dtplyr)
# Set up parallel processing plan with 6 workers using multisession
plan(multisession(workers = 20))
# Check current memory usage
pryr::mem_used()
## dumbed-down simulation of selection bias with noise (aka selecting on fielding of scale before fielding)
# set seed for reproducibility
set.seed(123)
# parameters
n_studies <- 100000
noise_dist <- seq(0.01, 0.07, 0.01)
## draw distribution of true alpha values
@rubenarslan
rubenarslan / reliability_multilevel.R
Created October 1, 2023 16:07
estimating reliability from varying intercepts
library(tidyverse)
# simulate unobserved trait and items with error
N <- 100
dat <- tibble(
id = 1:N,
trait = rnorm(N),
item1 = trait + rnorm(N),
item2 = trait + rnorm(N),
item3 = trait + rnorm(N),
item4 = trait + rnorm(N),
@rubenarslan
rubenarslan / estimand_sex_diff.R
Created August 16, 2023 09:33
misspecified estimand explains result
library(dplyr)
N = 10000
people <- tibble(
male = sample(0:1, N, TRUE),
trait = rnorm(N) + male)
true_means <- people %>%
group_by(male) %>%
summarise(sex_mean = mean(trait))
people <- people %>% left_join(true_means) %>%
left_join(true_means %>% rename(opposite_sex_mean = sex_mean) %>% mutate(male = 1-male)) %>%
@rubenarslan
rubenarslan / evb_strength.R
Created June 19, 2023 17:26
EVB strength DGM
N <- 5000
sim <- tibble(
group = rep(0:1, each = N/2),
evb_strength = 0.1 + group * 0.9,
# four orthogonal traits
trait1 = rnorm(N),
trait2 = rnorm(N),
trait3 = rnorm(N),
trait4 = rnorm(N),
# evaluative bias
@rubenarslan
rubenarslan / gpt_cov.md
Last active February 10, 2023 17:02
gpt covariance matrix
options(width=200)
library(tidyverse)
gpt <- rio::import("~/Downloads/tmp3.csv")

items <- tibble(var=c(str_c("ipipc", 1:10),  str_c("grit", 1:10)), 
                label = names(gpt)[-1])
knitr::kable(items)
n <- 1000
people <- tibble(
children = rpois(n, 1.4),
# for childless people, sat with childcare is not applicable
sat_with_child_care = if_else(children > 0, rnorm(n), NA_real_),
sat_with_housing = rnorm(n),
happiness = rnorm(n) + sat_with_housing +
if_else(children > 0, sat_with_child_care, 0),
)
mean_se_cluster <- function (x, mult = 1, cluster = NULL)
{
x_na <- is.na(x)
x <- x[!x_na]
cluster <- cluster[!x_na]
stopifnot(!is.null(cluster))
mod <- lme4::lmer(x ~ 1 + (1 | cluster))
intercept <- broom.mixed::tidy(mod, effects = "fixed")
se <- mult * intercept$std.error
mean <- intercept$estimate
@rubenarslan
rubenarslan / paper_rejection_cycle.R
Created June 16, 2020 19:44
paper rejection cycle
library(tidyverse)
n_papers <- 1000
papers <- tibble(
paper = 1:n_papers,
quality = rnorm(n_papers),
reviews = 0,
published = FALSE,
journal = NA_real_,
most_recent_assessment = NA_real_