Skip to content

Instantly share code, notes, and snippets.

View bbolker's full-sized avatar

Ben Bolker bbolker

View GitHub Profile
library(lme4)
load("areaRep26_glmer.Rdata")
n <- nrow(areaRep26_glmer)
nsp <- length(unique(areaRep26_glmer$especie))
## experimental design: each species appears in exactly 1 'sistema'
## each species appears once in each 'tempo'
with(areaRep26_glmer,table(especie,sistema))
with(areaRep26_glmer,table(especie,tempo))
summary(areaRep26_glmer)
@bbolker
bbolker / simulate_formula.R
Created September 20, 2020 20:27
methods and tests for simulate-formula methods
simulate.formula <- function(object, nsim=1, seed=NULL, ...) {
## utility fun for generating new class
cfun <- function(cc) {
c(paste0("formula_lhs_", cc), "formula_lhs", class(object))
}
if(length(object)==3) { ## two-sided formula
lhs <- object[[2L]]
.Basis <- try(eval(lhs, envir=environment(object),
@bbolker
bbolker / glmmTMB_GH639.R
Created December 5, 2020 00:42
workspace/minimal example
EVCounts <- structure(list(name = c("ANUBIS", "ANUBIS", "ANUBIS", "ANUBIS",
"ANUBIS", "ANUBIS", "ANUBIS", "ANUBIS", "ANUBIS", "ANUBIS", "ANUBIS",
"ANUBIS", "ASHERA", "ASHERA", "ASHERA", "ASHERA", "ATHENA", "ATHENA",
"ATHENA", "ATHENA", "ATHENA", "ATHENA", "ATHENA", "ATUTAHI",
"ATUTAHI", "ATUTAHI", "BASTET", "BASTET", "BASTET", "BASTET",
"BASTET", "BASTET", "BAUBO", "BAUBO", "BAUBO", "BAUBO", "BRAN",
"BRAN", "BRAN", "BRAN", "BRAN", "CHAC", "CHAC", "CHAC", "CHAC",
"CHAC", "CHAC", "CHAC", "COEUS", "COEUS", "COEUS", "COEUS", "COEUS",
"DIONYSUS", "DIONYSUS", "DIONYSUS", "DIONYSUS", "DON", "DON",
"DON", "DON", "DON", "ERIS", "ERIS", "ERIS", "ERIS", "ERIS",
@bbolker
bbolker / wid.csv
Created April 4, 2021 17:24
data from https://wid.world/ (US, pre-tax income)
percentile year inc1 inc2
p20p80 1980 31718.0283 29193.2179
p20p80 1981 31735.9894 29279.6271
p20p80 1982 30590.3481 28139.5886
p20p80 1983 30775.1342 28381.5553
p20p80 1984 32216.2634 29989.8252
p20p80 1985 32851.6509 30679.849
p20p80 1986 33408.1257 31126.8753
p20p80 1987 34021.8476 31658.5276
p20p80 1988 34504.1039 32200.9287
@bbolker
bbolker / gepd_prior.R
Created May 13, 2021 20:48
generalized exponential power priors
## https://cran.r-project.org/web/packages/gnorm/vignettes/gnormUse.html
library(gnorm)
#' @param lwr lower (loose) bound
#' @param upr upper (loose) bound
#' @param tail_prob prob outside {lwr, upr}
#' @param ctr_prob probability of being in the middle 50\% of {lwr, upr}
get_gnorm <- function(lwr=-1, upr=1, tail_prob=2*pnorm(lwr),
ctr_prob=abs(diff(pnorm(c(-1,1)*lwr/2)))) {
## default tail_prob/ctr_prob assume lwr/upr symmetric around 0 ...
## start from Gaussian
@bbolker
bbolker / trunc_phenom.R
Created July 26, 2021 00:00
exploring the properties of truncated distributions with 'simple' data
library(readxl)
my_sum <- function(f) {
C <- contr.sum(length(levels(f)))
colnames(C) <- paste0(".sum_",levels(f)[-length(levels(f))])
return(C)
}
dd <- read_excel("Truncated.compois.problem.xlsx")
dd$pop <- as.factor(dd$pop)
@bbolker
bbolker / trunc_compois_phenom.Rmd
Created July 27, 2021 19:47
exploring inference for 'separated' trunc distribution
```{r}
dd <- read_excel("Truncated.compois.problem.xlsx")
dd <- transform(dd, pop = factor(pop), morph = factor(morph))
m0 <- glmmTMB(mates ~ pop*morph,
family='truncated_compois', data=dd)
m1 <- update(m0, contrasts = list(pop = my_sum(dd$pop), morph = my_sum(dd$morph)))
tt <- purrr::map_dfr(list(treatment=m0,sum=m1),
tidy, effects = "fixed", .id = "contrasts", conf.int = TRUE) %>%
filter(term != "(Intercept)") %>%
@bbolker
bbolker / youngadult_ON_vaxproj.R
Created August 1, 2021 23:58
projection of young adult full vax rates in ON
##https://www.dataquest.io/blog/r-api-tutorial/
## https://tonyelhabr.rbind.io/post/nested-json-to-tidy-data-frame-r/
## https://api.covid19tracker.ca/docs/1.0/vaccine-age-groups
library(httr)
library(jsonlite)
library(tibble)
library(purrr)
library(ggplot2)
library(dplyr)
library(nlme)
library(tidyverse)
library(dotwhisker)
library(broom.mixed)
data <- read_csv("https://raw.githubusercontent.com/mgree013/Rethinking_Biodiversity_Streams/main/data.csv",
col_types = cols())
## add centred/scaled vars
## might not help but can't hurt
data_s <- mutate(data,
@bbolker
bbolker / crossed_gamlss.R
Created September 29, 2021 00:42
unbalanced crossed random effects in lme4, lme4, gamlss
library(lme4)
library(nlme)
set.seed(101)
## want to sample a *mostly* balanced design?
np <- nrow(Penicillin)
pu <- Penicillin[sample(np, size=round(np/2)),]
pu$dummy <- 1
## WILL NOT WORK without a "grouped data" structure
## set up a dummy grouping variable for the top level (necessary)