bbolker

library(tidyverse)
library(png)
## https://stackoverflow.com/questions/9917049/inserting-an-image-to-ggplot2
## 2023 date is when McMaster's university office announced to recipients
##  (before ResearchNet posting, dates may vary by university)
this_year <- 2024
(df <- read.table(text="
2023 March 29
2022 April 11
2021 April 14

bbolker

N <- 100
set.seed(101)
dd <- data.frame(
    y = rpois(N, exp(rnorm(N))),
    f = factor(seq(N))
)
dd <- dd[order(dd$y), ]
library(lme4)
m <- glmer(y ~ 1 + (1|f), family = poisson, data = dd)

bbolker

## https://forward.com/fast-forward/389673/finally-mathematical-proof-that-dreidel-is-a-terrible-game/
## Feinerman, Robert. “An Ancient Unfair Game.” American Mathematical Monthly 83, no. 8 (1976). https://www.jstor.org/stable/2319887.

library(emdbook)
library(ggplot2)
library(rayshader)
library(colorspace)

pvec <- 2:5
nvec <- 2:20

bbolker

library(lme4)
library(sjPlot)
set.seed(101)
n <- 1e3 ## obs per group
ns <- 4  ## groups
dd <- data.frame(s = rep(1:ns, each= n),
                 f = factor(sample(1:20, size = n*ns, replace = TRUE)),
                 x = rnorm(n*ns))
sims <- lapply(1:ns,
               function(i) simulate(~x + (1|f),

bbolker

## hinges as computed by stat_boxplot
my_hinges <- function(y, type = 7, coef = 1.5) {
    qs <- as.numeric(stats::quantile(y, c(0.25, 0.75), type = type))
    iqr <- diff(qs)
    outliers <- y < qs[1]-coef*iqr | y > qs[2] + coef*iqr
    hinges <- range(c(qs, y[!outliers]))
    return(hinges)
}

f_hinges <- function(x) {

bbolker

library(lme4)
library(glmmTMB)
library(glmulti)
set.seed(101)

## updated 3 March 2024; take random effects as formulas so that e.g. `ar1()` terms work

## FIXME: better data so we don't get convergence weirdness etc.?
## A random vector of count data
vy1 <- round(runif(100, min=1,max=20)*round(runif(100,min=1,max=20)))

bbolker

## Comparison of bias-correction methods

## n.b. both "logit-normal" and "logistic-normal" are used
##  for logit(Y) ~ Gaussian(mu, sigma)
library(logitnorm)
library(emmeans)

##
eta <- -2
ranef_sd <- 1

bbolker

library(brms)
library(lme4)
system.time(b1 <- brm(Reaction ~ Days + (Days|Subject), sleepstudy)) ## 39 sec

scode <- make_stancode(Reaction ~ Days + (Days|Subject), data = sleepstudy)
sdata <- make_standata(Reaction ~ Days + (Days|Subject), data = sleepstudy)

library(microbenchmark)
m1 <- microbenchmark(
    lmer = lmer(Reaction ~ Days + (Days|Subject), sleepstudy),

bbolker

library(ape)
cat("(((Strix_aluco:4.2,Asio_otus:4.2):3.1,",
    "Athene_noctua:7.3):6.3,Tyto_alba:13.5);",
    file = "ex.tre", sep = "\n")
tree.owls <- read.tree("ex.tre")
set.seed(101)
dd <- data.frame(y = rnorm(4),
                 obs = factor(1:4),
                 g = factor(1))
library(lme4)

bbolker

## example of updated refit, now incorporated in code
## * NO convergence testing/warnings/etc
## * won't work on weird binomial input (i.e. anything other than 0/1)
## * results seem only to be identical up to floating-point accuracy
## * doesn't yet do other clever things like resetting starting values

remotes::install_github("glmmTMB/glmmTMB/glmmTMB@smart_refit")
library(glmmTMB)
library(rbenchmark)
data("sleepstudy", package = "lme4")