bbolker

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)

bbolker

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

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

bbolker

```{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

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

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

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

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

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)