Ben Bolker bbolker

bbolker / softmax_regression.R

Created November 7, 2024 17:01

fitting a regression with coefficients constrained to sum to 1

	## https://fediscience.org/@[email protected]/113439584737810431

	library(MASS) # for eqscplot
	library(RTMB)

	set.seed(101)
	n <- 1000
	X <- cbind(
	x1 = rnorm(n , 10, 2),
	x2 = rnorm(n, 7, 1),

bbolker / constrained_cor.R

Created August 24, 2024 22:59

attempt to impose arbitrary equality constraints on correlations in a glmmTMB fit

	## simulate example (NOT satisfying constraints, but large and close enough that we
	## probably won't hit singular fits etc. etc.)

	library(glmmTMB)
	set.seed(1001)
	nx <- 3
	N <- 1000
	ng <- 20
	dd <- replicate(nx, rnorm(N), simplify = FALSE) \|> as.data.frame() \|> setNames(paste0("x", 1:nx))
	dd$g <- factor(sample(1:ng, replace = TRUE, size = N))

bbolker / brightstar_scrape.R

Created July 23, 2024 22:34

code to scrape the list of bright stars from Wikipedia and count by hemisphere

	## https://ivelasq.rbind.io/blog/politely-scraping/index.html
	## To clean data
	library(tidyverse)
	library(lubridate)
	library(janitor)
	# To scrape data
	library(rvest)
	library(httr)
	library(polite)

bbolker / lme4_groupRE.R

Created June 23, 2024 23:32

reproduce example for CV question

	## https://stats.stackexchange.com/questions/649052/selecting-the-correct-r-lme4-syntax-for-nested-models

	library(lme4)
	df <- data.frame(Subject = rep(factor(1:120), each = 200),
	Group = rep(LETTERS[1:3], each = 40*200))
	with(df, table(Group))
	with(df, table(table(Subject)))

	df$Accuracy <- simulate(~ Group + (1\|Subject),
	newdata = df,

bbolker / xmeta.R

Last active December 21, 2023 20:01

attempt at matching metafor results with glmmTMB

	library(brms)
	library(metafor)
	library(broom)
	library(broom.mixed)

	set.seed(101)
	y <- rnorm(10)
	sei <- rgamma(10, shape = 1)

	## metafor: fixed-effect meta-analysis

bbolker / find_data_arrays.R

Last active November 27, 2023 16:38

search installed packages for data of specified types

	i1 <- installed.packages()
	null <- matrix(nrow=0, ncol = 2,
	dimnames=list(NULL, c("Item", "type")))
	get_data <- function(pkg) {
	cat(pkg,"\n")
	dd <- data(package=pkg)$results
	## library(pkg, character.only = TRUE)
	## on.exit(try(detach(paste0("package:", pkg))))
	if (nrow(dd) == 0) return(NULL)
	FUN <- function(x) {

bbolker / bayes_tails.R

Last active September 12, 2023 22:20

bayes tails examples

	## original by Richard McElreath at https://twitter.com/rlmcelreath/status/1701165075493470644,
	## https://gist.github.com/rmcelreath/39dd410fc6bb758e54d79249b11eeb2f
	## originally based on https://doi.org/10.1006/jmva.1999.1820
	## with improvements from https://gist.github.com/murphyk/94205bcf335837108ff0e6d51331785a

	post_prior <- function(
	m_pr = 10, m_lik = 0, ## mean values for prior and likelihood
	sd_pr = 1, sd_lik = 1, ## standard deviations
	df_pr = 1000, df_lik = 1000, ## dfs (df = 1000 is effectively Gaussian)
	xlim = c(-5, 15), n = 1001, ## range and delta for x-vector

bbolker / binom_p.R

Created July 8, 2023 17:53

distribution of p-values under the null for a binomial test via GLM

	## adapted from Richard McElreath's code at
	## https://gist.github.com/rmcelreath/4f7010e8d5688c69bbeb7008f0aabe65

	binom_p <- function(N=10, M=5, p=0.25, b=0,
	test = c("Wald", "LRT")) {
	test <- match.arg(test)
	x <- rep(0:1, each = N/2)
	p <- plogis(qlogis(p)+b*x)
	y <- rbinom(N, size=M, prob=p)
	z <- glm( cbind(y, M-y) ~ x , family=binomial )

bbolker / analemma.R

Created June 22, 2023 22:59

compute sunrise/sunset differences

	start_date <- "2023-06-21"
	end_date <- "2023-07-01"
	refdate <- "2023-06-21"
	season <- "summer"

	library(ggplot2); theme_set(theme_bw())
	library(colorspace)
	library(purrr)
	library(dplyr)
	## need archived package for sunrise/sunset calcs

bbolker / cbrm.R

Created May 22, 2023 20:36 — forked from derekpowell/cbrm.R

Wrapper for brm() that supports caching of BRMS models