bbolker

---
title: "runr test"
date: "`r format(Sys.time(), '%H:%M %d %B %Y ')`"
author: Ben Bolker
---

```{r opts,message=FALSE,echo=FALSE,warning=FALSE}
library("knitr")
library("runr")
py <- proc_python()

bbolker

#' apply a function to a factorial combination of elements of lists
#' returns (if \code{FLATTEN=TRUE}) a flat list (with length equal to the product of the
#' lengths of the input lists) of results, along with a \code{grid} attribute containing
#' a data frame giving the values used for each element in the list
#' @param FUN function to apply
#' @param ... list of vectors to apply to
#' @param FLATTEN
#' @param MoreArgs additional arguments to pass to \code{FUN}
#' @param GlobalVars list of names of variables in global environment needed for parallel runs
#' @param .progress progress bar?

bbolker

library(lme4)
set.seed(101)
dd <- data.frame(x=rnorm(200),f=factor(rep(1:10,each=20)))
dd$y <- simulate(~x+(1|f),
                 newdata=dd,
                 newparams=list(beta=c(10,1),
                                theta=1,
                                sigma=5),
                 family=Gamma(link="identity"))[[1]]
## NB not sure what sigma really is here??? must be shape parameter

bbolker

library(lme4)
set.seed(101)
dd <- data.frame(x=rnorm(200),f=factor(rep(1:10,each=20)))
dd$y <- simulate(~x+(1|f),
                 newdata=dd,
                 newparams=list(beta=c(10,1),
                                theta=1,
                                sigma=5),
                 family=Gamma(link="identity"))[[1]]
## NB not sure what sigma really is here??? must be shape parameter

bbolker

bigcorPar <- function(x, nblocks = 10, verbose = TRUE, ncore="all", ...){
  library(ff, quietly = TRUE)
  require(doMC)
	if(ncore=="all"){
		ncore = multicore:::detectCores()
		registerDoMC(cores = ncore)
	} else{
		registerDoMC(cores = ncore)
	}

bbolker

## adapted from https://rviews.rstudio.com/2020/03/05/covid-19-epidemiology-with-r/

## CC-BY-SA 3.0 https://creativecommons.org/licenses/by-sa/3.0/

library(tidyverse)
library(lubridate)
library(directlabels)
library(colorspace)
library(emmeans)
library(broom)

bbolker

f1 <- expression(eps/(2-x/eps))
f2 <- expression(eps*log(exp(x/eps)+1))
flist <- list(f1,f2)


dfun <- function(x,e,n=0,eps=0.001) {
    alt_less <- switch(n+1,
                  eval(e),
                  eval(D(e,"x")),
                  eval(D(D(e,"x"),"x")))

bbolker

library(nlme)
packageVersion('nlme')

## the values are fixed for all models
rho <- 0.3
fixform <- distance ~ age + factor(Sex)
## remove the grouping and other classes to keep it simple
## full balanced data set
Obal<-as.data.frame(Orthodont)
## create unbalanced Orthodont data set with no singletons

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)

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