minor updates

Author: florianhartig

.gitignore

@@ -17,3 +17,5 @@ DHARMa_0.0.0.9000.tgz
 DHARMa/Rplots.pdf
 
 DHARMa/vignettes/Demo.html
+
+*.gz

DHARMa/DESCRIPTION

@@ -1,14 +1,17 @@
 Package: DHARMa
 Title: Residual Diagnostics for HierArchical (Multi-level / Mixed) Regession Models
-Version: 0.0.0.9000
+Version: 0.0.1.0
+Date: 16.6.2016
 Authors@R: c(person("Florian", "Hartig", email = "florian.hartig@biom.uni-freiburg.de", role = c("aut", "cre")))
 Author: Florian Hartig [aut, cre],
 Maintainer: Florian Hartig <florian.hartig@biom.uni-freiburg.d>
 Description: The package creates easily interpretable, residuals for generalized linear mixed models that are standardized to values between 0 and 1. This is achieved by a simulation-based approach, similar to the Bayesian p-value or the parametric bootstrap: 1) simulate new data from the fitted model 2) calculate the cummulative empirical density function 3) residual is the value of the empirical density function at the value of the observed data.  
 Depends:
     R (>= 3.1.2)
 Imports:
-    gap
+    gap,
+    quantreg,
+    plyr
 Suggests:
     lme4
 License: What license is it under?

DHARMa/R/createData.R

@@ -1,7 +1,7 @@
 #' Function to simulate data
 #' @export
 #' @param overdispersion if this is a numeric value, it will be used as the sd of a random normal variate that is added to the linear predictor. Alternatively, a random function can be provided that takes as input the linear predictor. 
-#' @example /inst/examples/createPrior.R
+#' @example /inst/examples/createDataHelp.R
 createData <- function(replicates=1, sampleSize = 10, intercept = 0, fixedEffects = 1, numGroups = 10, randomEffectVariance = 1, overdispersion = 0.5, family = poisson(), scale = 1, cor = 0, roundPoissonVariance = NULL, quadraticFixedEffects = NULL){
 
   nPredictors = length(fixedEffects)
@@ -39,7 +39,7 @@ createData <- function(replicates=1, sampleSize = 10, intercept = 0, fixedEffect
     else if (family$family == "binomial") observedResponse = rbinom(n = sampleSize, 1, prob = linkResponse)
     else if (family$family == "poisson") {
       if(is.null(roundPoissonVariance)) observedResponse = rpois(n = sampleSize, lambda = linkResponse)
-      else observedResponse = round(rnorm(n = length(linkResponse), mean = linkResponse, sd = sqrt(linkResponse) * roundPoissonVariance))
+      else observedResponse = round(rnorm(n = length(linkResponse), mean = linkResponse, sd = roundPoissonVariance))
     }
 
     else stop("wrong link argument supplied")

DHARMa/R/plotSimulatedResiduals.R

@@ -1,16 +1,68 @@
 #' creates Poisson data overdispersion and random intercept
-plotSimulatedResiduals <- function(simulationOutput){
+plotSimulatedResiduals <- function(simulationOutput, quantreg = F){
+  
+  
+  pred = simulationOutput$fittedPredictedResponse
+  
+  resid = simulationOutput$scaledResiduals
+  
+  dat = data.frame(pred, resid)
+  
+  
   par(mfrow = c(1,2), oma = c(0,1,2,1))
 
-  gap::qqunif(simulationOutput$scaledResiduals,pch=21,bg="blue",bty="n", logscale = F)
+  gap::qqunif(simulationOutput$scaledResiduals,pch=2,bty="n", logscale = F, col = "black", cex = 0.6, main = "QQ plot residuals", cex.main = 1)
   #hist(simulationOutput$scaledResiduals, main = "Distribution of scaled residuals", breaks = 50, freq = F)
   #lines(density(simulationOutput$scaledResiduals, na.rm = T, from = 0, to = 1, kernel = "rectangular", bw = 0.01, cut = 0.01), col = "red")
 
-  plot(simulationOutput$fittedPredictedResponse, simulationOutput$scaledResiduals, xlab = "predicted", ylab = "Residual", main = "Residual vs. predicted")
+  plot(simulationOutput$fittedPredictedResponse, simulationOutput$scaledResiduals, xlab = "Predicted", ylab = "Residual", main = "Residual vs. predicted\n lines should match", cex.main = 1)
+  
+
+  
+
+  if(quantreg == F){
+    
+    lines(smooth.spline(simulationOutput$fittedPredictedResponse, simulationOutput$scaledResiduals, df = 10), lty = 2, lwd = 2, col = "red")
+    
+    abline(h = 0.5, col = "red", lwd = 2)
+
+  }else{
+    
+    require(quantreg)
+    dat <- plyr::arrange(dat,pred)
+    fit<-quantreg::rqss(resid~qss(pred,constraint="N"),tau=0.5,data = dat)
+    lines(unique(dat$pred)[-1],fit$coef[1] + fit$coef[-1], col = "red", lwd = 2)
+    abline(h = 0.5, col = "red", lwd = 2)
+    
+    fit<-quantreg::rqss(resid~qss(pred,constraint="N"),tau=0.25,data = dat)
+    lines(unique(dat$pred)[-1],fit$coef[1] + fit$coef[-1], col = "green", lwd = 2, lty =2)
+    abline(h = 0.25, col = "green", lwd = 2, lty =2)
+    
+    fit<-quantreg::rqss(resid~qss(pred,constraint="N"),tau=0.75,data = dat)
+    lines(unique(dat$pred)[-1],fit$coef[1] + fit$coef[-1], col = "blue", lwd = 2, lty = 2)
+    abline(h = 0.75, col = "blue", lwd = 2, lty =2)   
+  }
+  
+  
+
 
   mtext("DHARMa scaled residual plots", outer = T)
 }
 
+#' Plots a generic residual plot with spline 
+plotResiduals <- function(pred, res){
+  
+  plot(pred, res)
+  
+  lines(smooth.spline(simulationOutput$fittedPredictedResponse, simulationOutput$scaledResiduals, df = 10), lty = 2, lwd = 2, col = "red")
+  
+  abline(h = 0.5, col = "red", lwd = 2)
+  
+}
+
+
+
+#plotSimulatedResiduals(simulationOutput)
 
 #plot(simulationOutput$observedResponse, simulationOutput$scaledResiduals, xlab = "predicted", ylab = "Residual", main = "Residual vs. predicted")
 

DHARMa/man/createData.Rd

@@ -9,7 +9,7 @@ vignette: >
 ---
 
 ```{r global_options, include=FALSE}
-knitr::opts_chunk$set(fig.width=7, fig.height=4, warning=FALSE, message=FALSE, cache = T)
+knitr::opts_chunk$set(fig.width=8, fig.height=5, warning=FALSE, message=FALSE, cache = F)
 ```
 
 ***Summary**: The DHARMa package creates easily interpretable, residuals for generalized linear mixed models that are standardized to values between 0 and 1. This is achieved by a simulation-based approach, similar to the Bayesian p-value or the parametric bootstrap: 1) simulate new data from the fitted model 2) calculate the cummulative empirical density function 3) residual is the value of the empirical density function at the value of the observed data.*  
@@ -61,7 +61,7 @@ The above exampel shows a misspecified model. The reason is that we created data
 For comparison, also "standard residual plots" (deviance, pearso and raw) are shown
 
 ```{r}
-testData = createData(sampleSize = 500, intercept = 0, overdispersion = 0, family = poisson(), randomEffectVariance = 0)
+testData = createData(sampleSize = 250, intercept = 0, overdispersion = 0, family = poisson(), randomEffectVariance = 0)
 fittedModel <- glmer(observedResponse ~ Environment1 + (1|group) , family = "poisson", data = testData)
 
 plotConcentionalResiduals(fittedModel)
@@ -105,9 +105,11 @@ plotSimulatedResiduals(simulationOutput = simulationOutput)
 
 
 ```{r}
-testData = createData(sampleSize = 500, intercept=1, fixedEffects = 2, overdispersion = 0, family = poisson(), roundPoissonVariance = 0.6)
+testData = createData(sampleSize = 500, intercept=0, fixedEffects = 2, overdispersion = 0, family = poisson(), roundPoissonVariance = 0.001, randomEffectVariance = 0)
 fittedModel <- glmer(observedResponse ~ Environment1 + (1|group) , family = "poisson", data = testData)
 
+summary(fittedModel)
+
 plotConcentionalResiduals(fittedModel)
 
 simulationOutput <- simulateResiduals(fittedModel = fittedModel)
@@ -173,8 +175,8 @@ Difficult to see with the overall pattern, but it becomes clear if we plot again
 
 ```{r}
 par(mfrow = c(1,2))
-plot(testData$Environment1,  simulationOutput$scaledResiduals)
-plot(testData$Environment2,  simulationOutput$scaledResiduals)
+plotResiduals(testData$Environment1,  simulationOutput$scaledResiduals)
+plotResiduals(testData$Environment2,  simulationOutput$scaledResiduals)
 ```
 
 

DHARMa/man/plotResiduals.Rd

@@ -1,2 +1,2 @@
-rsync -avz --delete /Users/Florian/Home/Projekte/Papers_in_Progress/15-RegressionDiagnosticsGLMM/DHARMa/DHARMa_0.0.0.9000.tar.gz /Users/Florian/Dropbox/Deploy/DHARMa/DHARMa_0.0.0.9000.tar.gz
+rsync -avz --delete /Users/Florian/Home/Projekte/Papers_in_Progress/15-RegressionDiagnosticsGLMM/DHARMa/DHARMa_0.0.* /Users/Florian/Dropbox/Deploy/DHARMa/DHARMa_0.0.0.9000.tar.gz