Merge pull request #61 from florianhartig/parallel

Remodeled simulations / parallel options

Code/Experimental/arm.R

@@ -0,0 +1,121 @@
+library(lme4)
+library(arm)
+
+data(VerbAgg, package = 'lme4')
+
+verb_mod <- glmer(r2 ~ (Anger + Gender + btype + situ)^2 + (1|id) + (1|item), family = binomial, data = VerbAgg, control = glmerControl(optimizer = "bobyqa", optCtrl = list(maxfun = 100000)))
+summary(verb_mod)
+
+par(mfcol=c(1, 2))
+binnedplot(predict(verb_mod, type="response", re.form=NULL), resid(verb_mod, type="response"), nclass=40, main='With random effects')
+binnedplot(predict(verb_mod, type="response", re.form=NA), resid(verb_mod, type="response"), nclass=40, main='Without random effects')
+
+plot(predict(verb_mod, type="response", re.form=NULL), predict(verb_mod, type="response", re.form=NA))
+
+
+library(DHARMa)
+
+dat = createData(sampleSize = 500, family = binomial(), randomEffectVariance = 1, fixedEffects = c(1,1), quadraticFixedEffects = c(5,5))
+fit = glmer(observedResponse ~ Environment1 + (1|group), family = binomial, data = dat)
+par(mfcol=c(1, 2))
+binnedplot(predict(fit, type="response", re.form=NULL), resid(fit, type="response"), nclass=40, main='With random effects')
+binnedplot(predict(fit, type="response", re.form=NA), resid(fit, type="response"), nclass=40, main='Without random effects')
+
+
+res <- simulateResiduals(fit)
+plot(res)
+plotResiduals(dat$Environment1, res$scaledResiduals)
+plotResiduals(dat$Environment2, res$scaledResiduals)
+
+
+res <- simulateResiduals(verb_mod)
+plot(res)
+
+par(mfrow = c(2,2))
+plotResiduals(VerbAgg$Anger, res$scaledResiduals, main = "Anger")
+plotResiduals(VerbAgg$Gender, res$scaledResiduals, main = "Gender")
+plotResiduals(VerbAgg$btype, res$scaledResiduals, main = "btype")
+plotResiduals(VerbAgg$situ, res$scaledResiduals, main = "situ")
+
+ranEst = ranef(verb_mod)
+plotResiduals(predict(verb_mod, type="response", re.form=NULL), res$scaledResiduals)
+plotResiduals(ranEst$id[,1][VerbAgg$id], res$scaledResiduals)
+
+
+
+# Create new data based on the fitted model and refit
+VerbAgg$newResponse =  simulate(verb_mod)$sim_1
+verb_mod <- glmer(newResponse ~ (Anger + Gender + btype + situ)^2 + (1|id) + (1|item), family = binomial, data = VerbAgg)
+
+# Check residuals
+res <- simulateResiduals(verb_mod)
+ranEst = ranef(verb_mod)
+plotResiduals(predict(verb_mod, type="response", re.form=NULL), res$scaledResiduals)
+plotResiduals(ranEst$id[,1][VerbAgg$id], res$scaledResiduals)
+
+
+# ==========================================
+
+
+#for sim function
+library(arm)
+
+## Annette Dobson (1990) "An Introduction to Generalized Linear Models".
+## Page 9: Plant Weight Data.
+ctl <- c(4.17,5.58,5.18,6.11,4.50,4.61,5.17,4.53,5.33,5.14)
+trt <- c(4.81,4.17,4.41,3.59,5.87,3.83,6.03,4.89,4.32,4.69)
+group <- gl(2, 10, 20, labels = c("Ctl","Trt"))
+weight <- c(ctl, trt)
+lm.D9 <- lm(weight ~ group)
+lm.null <- lm(weight ~ 1)
+
+#from rethinking library for numerically stable log sums
+log_sum_exp <- function (x) {
+  xmax <- max(x)
+  xsum <- sum(exp(x - xmax))
+  xmax + log(xsum)
+}
+
+
+#function for WAIC from an LM
+waic.lm <- function(mod, n.sims=1e3){
+  mod_sims <- sim(mod, n.sims=10)
+  mod_X <- model.matrix(mod)
+  mod_Y <- mod$fitted.values+mod$residuals
+
+  #generate distribution of observations
+  pred_sims <- apply(mod_sims@coef, 1, function(b) mod_X %*% b )
+  # pred_sims_err <- sapply(1:nrow(pred_sims), 
+  #                     function(i) rnorm(ncol(pred_sims), pred_sims[i,], mod_sims@sigma[i]))
+
+  #I hate nested loops, but it's expedient
+  ll <- sapply(1:ncol(pred_sims), 
+               function(j){
+                 sapply(1:nrow(pred_sims),
+                        function(i) dnorm(pred_sims[i,j], mod_Y[i], mod_sims@sigma[i], log=TRUE))
+               })
+
+  #get the things we'll need...
+  lppd <- apply(ll, 1, function(arow) log_sum_exp(arow) - log(length(arow)))
+  pWAIC <- apply(ll, 1, var)
+
+  -2*sum(lppd) + 2*sum(pWAIC)             
+}
+
+waic.lm(lm.D9)
+waic.lm(lm.null)
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+

Code/Experimental/glmmtmb-correlationStructures.R

@@ -0,0 +1,23 @@
+
+
+# mildly long-tailed
+plot()
+
+# An ARIMA simulation
+x = runif(100)
+time = 1:100
+id = 1:100
+
+ts.sim <- arima.sim(n = 100, list(ar = c(0.8897, -0.4858), ma = c(-0.2279, 0.2488)),
+                    rand.gen = function(n, ...) sqrt(0.1796) * rt(n, df = 5)) + 5*x 
+ts.plot(ts.sim)
+
+m1 <- glmmTMB(ts.sim ~ x + ar1(time -1 |id) )
+summary(m1)
+
+library(nlme)
+library(glmmTMB)
+data(sleepstudy,package="lme4")
+sleepstudy$row <- factor(1:180)
+gt_min <- glmmTMB(Reaction ~ (1|Subject) + ar1(row + 0 | Subject), sleepstudy)
+summary(gt_min)

DHARMa/DESCRIPTION

@@ -21,17 +21,18 @@ Imports:
     graphics,
     utils,
     grDevices,
+    parallel,
+    doParallel,
+    foreach,
     gap,
     qrnn,
     lmtest,
     ape,
     sfsmisc,
     MASS,
-    doParallel,
-    foreach,
     lme4,
     mgcv,
-    glmmTMB
+    glmmTMB (> 0.2)
 Suggests:
     knitr,
     testthat

DHARMa/NAMESPACE

@@ -5,15 +5,14 @@ S3method(plot,DHARMa)
 S3method(print,DHARMa)
 S3method(refit,glmmTMB)
 S3method(refit,lm)
-export(benchmarkOverdispersion)
-export(benchmarkP)
-export(benchmarkUniformity)
 export(createDHARMa)
 export(createData)
+export(getRandomState)
 export(plotConventionalResiduals)
 export(plotQQunif)
 export(plotResiduals)
 export(plotSimulatedResiduals)
+export(runBenchmarks)
 export(simulateResiduals)
 export(testOverdispersion)
 export(testOverdispersionParametric)