bbolker

## 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

## 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

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

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

## 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

## 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

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 <- function(formula,
                 data,
                 family = gaussian(),
                 prior = NULL,
                 autocor = NULL,
                 cov_ranef = NULL,
                 sample_prior = c("no", "yes", "only"),
                 sparse = FALSE,
                 knots = NULL,
                 stan_funs = NULL,

bbolker

library(coxme)
library(Matrix)

#' @param fit a fitted \code{coxme} model
#' @param data the data frame for which to evaluate conditional SDs
#' @param grpvar (character) grouping variable (should be able to
#' extract this from the model, but it's a nuisance)
#' @param rform random-effects term (ditto)
#' @param sd.only return only sqrt(diag(V))? (otherwise return full covariance matrix)
coxme.condsd <- function(fit,

bbolker

library(lme4)
library(nlme)
library(glmmTMB)
library(Matrix)

## what do various R packages do when confronted with an
## undefined REML problem (fixed effects *and* random effects
## for every subject?

## Twitter thread: https://twitter.com/ten_photos/status/1657399290166222850