From 3f2d53f6127c7c9589c61d0df9e0ab47842cfd1a Mon Sep 17 00:00:00 2001
From: flozanoisla <flozanoisla@gmail.com>
Date: Mon, 8 Jul 2024 19:54:12 -0500
Subject: [PATCH] updt

---
 NAMESPACE              |   1 +
 R/plot_diagnostic.R    | 114 +++++++++++++++++++++++++++++++++++++++++
 man/plot_diagnostic.Rd |  33 ++++++++++++
 3 files changed, 148 insertions(+)
 create mode 100644 R/plot_diagnostic.R
 create mode 100644 man/plot_diagnostic.Rd

diff --git a/NAMESPACE b/NAMESPACE
index 7038029..c9b58e8 100644
--- a/NAMESPACE
+++ b/NAMESPACE
@@ -16,6 +16,7 @@ export(mean_comparison)
 export(metamorphosis)
 export(outliers_remove)
 export(plot_diag)
+export(plot_diagnostic)
 export(plot_raw)
 export(plot_smr)
 export(remove_outliers)
diff --git a/R/plot_diagnostic.R b/R/plot_diagnostic.R
new file mode 100644
index 0000000..4cf4114
--- /dev/null
+++ b/R/plot_diagnostic.R
@@ -0,0 +1,114 @@
+#' Diagnostic plots
+#'
+#' Function to plot the diagnostic of models
+#'
+#' @param data Experimental design data frame with the factors and traits.
+#' @param formula Mixed model formula
+#' @param title Plot title
+#'
+#' @return plots
+#' 
+#' @importFrom stats density
+#' 
+#' @export
+#' 
+#' @examples
+#' 
+#' \dontrun{
+#' 
+#' library(inti)
+#' 
+#' plot_diagnostic(data = potato
+#'                 , formula = stemdw ~ (1|bloque) + geno*treat)
+#' 
+#' }
+#' 
+
+plot_diagnostic <- function(data, formula, title = NA) {
+  
+  # data <- inti::potato; formula = stemdw ~ (1|bloque) + geno*treat
+  # formula <- stemdw ~ bloque + geno*treat; title = NA
+  
+  formula <- as.formula(formula)
+  
+  model  <-  if( length(lme4::findbars(formula))>0 ) {
+    
+    data %>% lme4::lmer(formula, data = .)
+    
+  } else {  data %>% aov(formula, data = .) }
+  
+  
+  stopifnot(class(model) %in% c("lm", "aov", "lmerMod"))
+  
+  title <- if(is.null(title) || is.na(title) ) {""} else {title}
+  
+  dt <- if ( inherits(model, "lm") || inherits(model, "aov") ) {
+    
+    ggplot2::fortify(model, data)
+    
+  } else if ( inherits(model, "lmerMod") ) {
+    
+    lme4::fortify.merMod(model, data)
+    
+  }
+  
+  # Histogram
+  p1 <- ggplot(dt , aes(.data$.resid)) +
+    geom_histogram(aes(y = ggplot2::after_stat(density)), 
+                   color = "black", fill = "grey82", bins = 30) + 
+    stat_function(fun = stats::dnorm, color = "blue", args = list(mean = 0, sd = sd(dt$.resid))) +
+    theme_bw() +
+    labs(x = "Residuals"
+         , y = "Frequency"
+         , title = title
+         , subtitle = "Histogram"
+         )
+  
+  # Q-Q plot
+  p2 <- ggplot(dt, aes(sample = .data$.resid)) +
+    stat_qq() +
+    stat_qq_line() +
+    theme_bw() +
+    labs(x = "Theorical Quantiles"
+         , y = "Sample Quantiles"
+         , title = title
+         , subtitle = "Q-Q plot"
+         )
+  
+  # Residual plot
+  p3 <- ggplot(dt, aes(.data$.fitted, .data$.resid)) +
+    geom_point() +
+    stat_smooth(method="loess", formula = 'y ~ x') +
+    geom_hline(yintercept=0, col="red", linetype="dashed") +
+    theme_bw() +
+    labs(x = "Fitted values"
+         , y = "Residuals"
+         , title = title
+         , subtitle = "Residual plot"
+         )
+  
+  # Scale-location | Homoscedasticity
+  p4 <- ggplot(dt, aes(.data$.fitted, sqrt(abs(.data$.resid)))) +
+    geom_point(na.rm=TRUE) +
+    stat_smooth(method="loess", na.rm = TRUE, formula = 'y ~ x') +
+    theme_bw() +
+    labs(x = "Fitted values"
+         , y = expression(sqrt("|Standardized residuals|"))
+         , title = title
+         , subtitle = "Homoscedasticity" 
+         )
+  
+  p5 <- ggplot(dt, aes(stats::hatvalues(model), .data$.resid)) + 
+    stat_smooth(method="loess", na.rm=TRUE) +
+    xlab("Leverage")+ylab("Standardized Residuals") +
+    ggtitle("Residual vs Leverage Plot") + 
+    scale_size_continuous("Cook's Distance", range=c(1,5)) + 
+    theme_bw() + theme(legend.position="bottom")
+  
+  return(list(histogram = p1
+              , qqplot = p2
+              , residual = p3
+              , homoscedasticity = p4
+              ))
+
+}
diff --git a/man/plot_diagnostic.Rd b/man/plot_diagnostic.Rd
new file mode 100644
index 0000000..1d5e577
--- /dev/null
+++ b/man/plot_diagnostic.Rd
@@ -0,0 +1,33 @@
+% Generated by roxygen2: do not edit by hand
+% Please edit documentation in R/plot_diagnostic.R
+\name{plot_diagnostic}
+\alias{plot_diagnostic}
+\title{Diagnostic plots}
+\usage{
+plot_diagnostic(data, formula, title = NA)
+}
+\arguments{
+\item{data}{Experimental design data frame with the factors and traits.}
+
+\item{formula}{Mixed model formula}
+
+\item{title}{Plot title}
+}
+\value{
+plots
+}
+\description{
+Function to plot the diagnostic of models
+}
+\examples{
+
+\dontrun{
+
+library(inti)
+
+plot_diagnostic(data = potato
+                , formula = stemdw ~ (1|bloque) + geno*treat)
+
+}
+
+}