GAMLSS Tutorial for Psychological Data (HQFigures).Rmd

---
title: "A GAMLSS Tutorial for Assessing Potential Heteroscedasticity in Psychological Data"
author: "Raydonal Ospina & Juan C. Correa"
date: "`r Sys.Date()`"
output:
  prettydoc::html_pretty:
    theme: cayman  # Keep the 'architect' theme
    highlight: github
    css: styles.css   # Add a separate CSS file for customization
header-includes:    
    geometry: margin = 1.0in
---

```{r, echo = FALSE, warning = FALSE}
# Packages
suppressPackageStartupMessages(library(tidyverse)) # Data wrangling
suppressPackageStartupMessages(library(scales)) # Data wrangling
suppressPackageStartupMessages(library(dplyr)) # Data wrangling
suppressPackageStartupMessages(library(gridExtra)) # arrange plots of ggplot 
suppressPackageStartupMessages(library(ggpubr)) # arrange plots of ggplot 
suppressPackageStartupMessages(library(MASS)) # for shapiro test
suppressPackageStartupMessages(library(car)) # for ncvTest    
suppressPackageStartupMessages(library(lmtest)) # dor Breusch-Pagan Test
suppressPackageStartupMessages(library(kableExtra)) # tables
suppressPackageStartupMessages(library(ggmosaic)) # tables
suppressPackageStartupMessages(library(knitr))
suppressPackageStartupMessages(library(GGally)) # gamlss distributions
suppressPackageStartupMessages(library(gamlss)) # gamlss fits
suppressPackageStartupMessages(library(gamlss.dist)) # gamlss distributions
suppressPackageStartupMessages(library(gamlss.util)) # gamlss fits
suppressPackageStartupMessages(library(broom)) # Convert Statistical Objects into Tidy Tibbles
suppressPackageStartupMessages(library(ggthemes)) # ggplot themes
suppressPackageStartupMessages(library(cowplot)) # ggplot arrange
suppressPackageStartupMessages(library(gtable)) # ggplot table
suppressPackageStartupMessages(library(RColorBrewer)) # cool colors
suppressPackageStartupMessages(library(colorspace)) # cool colors
suppressPackageStartupMessages(library(gvlma)) # Top-level function for Global Validation of Linear Models Assumptions.


# Auxiliar - name of response in linear regression model
responseName <- function (model, ...) deparse(attr(terms(model), "variables")[[2]])

# knitr::opts_chunk$set(fig.path="figures", comment=NA,  fig.width=13, echo = FALSE)

#Get data file in place
# setwd("~/MEGA/Pesquisas/Marmolejo-Ramos/Dunnin-Kruger-effect")

# Load data - data in same RMD directory
data<-read.csv("Dunning-Kruger_INTELL_UPLOAD.csv")
```

```{r setup, include=FALSE}
knitr::opts_chunk$set(warning = FALSE, message = FALSE) 
```

## Background of the data

This tutorial relies on the data set originally described and analysed by Gignac, G. E. & Zajenkowski, M. (2020). The Dunning-Kruger Effect is (mostly) a statistical artefact: Valid approaches to testing the hypothesis with individual differences data. *Intelligence*, 80. doi: 10.1016/j.intell.2020.101449

The Dunning-Kruger Effect (DKE) refers to the simple idea that people's ignorance is often invisible to them. This meta-ignorance (i.e., ignorance of ignorance) that emerges as the result of a lack of expertise and relies on misbeliefs and background knowledge that are not appropriate for providing correct answers to a variety of questions. Studying DKE requires the analysis of self-assessed intelligence and comparing it with objective metrics. The raw data set revisited is a CSV file with 929 rows and 56 variables. The self-assessed intelligence (SAIQ) is on a discrete scale ranging from 1 to 25, while the intelligence objective metric (IQ) resulted from applying the Advanced Progressive Matrices.

### Variables in data set

```{r,  echo = FALSE}
# columns' names and structure
data$sex <- as.factor(data$sex)
data$study <- as.factor(data$study)
print(str(data))
```

### Balanced distribution of participants in the design

In the data, we observed the presence of grouping variables that might lead to possible distributional changes and variability of response by, e.g., sex and type of study. Furthermore, we noticed a slight imbalance in the distribution of individuals.

```{r,  echo = FALSE, fig.width = 10, fig.height=6}
p.1 <- ggplot(data = data) +
geom_mosaic(aes(x = product(study, sex), fill=study))+ scale_fill_manual(values = c("#00AFBB", "#E7B800", "#FC4E07"))+theme_light()+ guides(fill=FALSE)

p.2 <- data %>% 
  count(sex, study) %>%
  group_by(study) %>%
  mutate(yes_no_pct = round(n / sum(n) * 100)) %>%
  ungroup() %>%
  ggplot(mapping = aes(x = sex, y = study))+geom_tile(mapping = aes(fill = n))+  scale_fill_gradientn(colours = c("#FC4E07","#E7B800","#00AFBB"))+
  geom_text(colour = "white", aes(label = paste0('Percent of Total: ', round((n / sum(n) * 100)), '%'))) + 
  geom_text(colour = "white", nudge_y = .2, aes(label = paste0('Count: ', n))) + 
  geom_text(colour = "white", nudge_y = -.2, aes(label = paste0('Percent of Group: ', yes_no_pct, '%'))) + guides(fill=FALSE)+ theme_light()

png("F1.png", width = 15, height = 7, units = 'in', res = 300)
grid.arrange(p.1, p.2, ncol=2,   widths = c(1.5, 1.5), top="balancing size by age and study groups")
dev.off()
```

### Plots for each combination of variable types of interest

We now select some variables of interest for building the regression model. We now see the scatterplots and distributions to determine if we have met the assumption of linearity. By checking the scatterplot, we get insights into outliers, extreme values, variability in the data set and how such variability moves together with other variables for the regression model.

The proper use of linear regression requires the correct marginal distribution for inferential purposes (Hypothesis tests, confidence intervals, and prediction). We can see that the `Raven_IQ` (mean-centred), `sex`, `age`, and `study` covariates plotted against the dependent variable of interest `SAI_winsorized_IQ` has different distributional behaviour. We can observe that the covariate `Raven_IQ` shows a bimodality mainly at the junction of the covariate `study==motivation` and asymmetry in `study == couples`, which can reflect the lack of balancing of the data and can impact distributional hypotheses in the case of normality of the errors in traditional linear regression (i.e., normal symmetrical distribution with unimodality).

```{r,  echo = FALSE, fig.width = 14, fig.height=8}
# select covariates for exploration and modelling
db <- data %>%  dplyr::select(ID, study,sex,age,Raven_IQ,SAI_winsorized_IQ) 

dbggpair <- db %>% dplyr::select(!ID)
# ggpair: plot for each combination of interest variable types
p <- ggpairs(dbggpair, 
             mapping = ggplot2::aes(color = study, alpha=0.7,  shape=sex),
        lower = list(combo = wrap(ggally_facethist, bins = 10, alpha=0.7)), 
        diag = list(continuous = wrap("densityDiag"), 
                    mapping = ggplot2::aes(fill=sex, alpha=0.7)))
for(i in 1:p$nrow) {
  for(j in 1:p$ncol){
    p[i,j] <- p[i,j] + 
        # scale_fill_manual(values=c("red", "cornflowerblue", "purple")) +
      scale_fill_manual(values=c("#00AFBB", "#E7B800", "#FC4E07")) +
      scale_color_manual(values = c("#00AFBB", "#E7B800", "#FC4E07"))
        # scale_color_manual(values=c("red", "cornflowerblue", "purple"))  
  }
}

png("F2.png", width = 15, height = 9, units = 'in', res = 300)
p+ theme_light()
dev.off()
```

### Linearity of the Data

To perform a more thorough evaluation of linear regression supposition, we proceed with loess smoother for each variable in the model. Loess smoother selectively weights points in the regression. We can also overlay the straight regression line as the original paper of Gignac and Zajenkowski (2020) to compare how the loess smoothing (detection of nonlinearities or local behaviour in data) compares to linear regression. In the plot, the grey shading corresponds to the standard errors for the loess smoother fit, while the red shading corresponds to the standard errors for the linear fit. The blue line is a loess smoothed line, and the red one is a linear regression line. As long as the loess smoother roughly approximates the linear tendency, the assumption of linearity has been apparently met. Nonetheless, a more careful, simple visual inspection suggests that the linearity hypothesis seems not to be fully satisfied. The confidence bands are of short length, and loess shows local variations following the trend line, which may indicate nonlinearity in the data.

```{r,  echo = FALSE, fig.width = 10, fig.height=5}
png("F3.png", width = 10, height = 5, units = 'in', res = 300)
ggplot(data, aes(x=Raven_IQ,y=SAI_winsorized_IQ)) + stat_smooth(method="loess", alpha=.2) + stat_smooth(method="lm", formula=y ~ x, color="red", fill="red", alpha=.2) + geom_point( aes(colour=study, shape=sex), size=2.3)+ scale_color_manual(values = c("#00AFBB", "#E7B800", "#FC4E07"))+theme_linedraw()
dev.off()
```

Linearity assumption can also be checked by plotting the Residuals versus the Fitted observations. Ideally, the residual plot will show no fitted pattern. That is, the blue line should be approximately horizontal at zero. The presence of a pattern may indicate a problem with some aspects of the linear model. In our case, we note that linearity is a reasonable assumption. We can use `crPlots()` function of `car` package in `R`.

```{r,  echo = TRUE}
# The model - OLS estimation
lm.temp <- lm(SAI_winsorized_IQ~Raven_IQ, data=db)
```

```{r,  echo = FALSE, fig.width = 10, fig.height=5}
db$resid <- resid(lm.temp)
db$fitted <- fitted(lm.temp)  

png("F4.png", width = 15, height = 7, units = 'in', res = 300)
ggplot(db, aes(x=fitted,y=resid)) + stat_smooth(method="loess", alpha=.2) + 
  geom_point( aes(colour=study, shape=sex), size=2.3)+ scale_color_manual(values = c("#00AFBB", "#E7B800", "#FC4E07"))+geom_hline(yintercept=0, col="red", linetype="dashed")+xlab("Fitted values")+ylab("Residuals")+theme_linedraw()
dev.off()

# ggplot(lm.temp, aes(.fitted, .resid))+stat_smooth(method="loess", alpha=.2)+geom_hline(yintercept=0, col="red", linetype="dashed")+xlab("Fitted values")+ylab("Residuals")+ geom_point(shape=1)+  theme_linedraw()
```

The function `crPlots()` from the `car` package in `R` elaborates upon the partial-residual. The red line represents the line of best fit. The assumption is ok if the green line seems to be similarly linear as the red line. If the blue line appears curved relative to the red line, we likely have a linearity problem.

```{r,  echo = FALSE, fig.width = 10, fig.height=5}
db$pred.orig <- predict(lm.temp, type="terms")
db$partial.resid <- db$pred.orig+db$resid  

png("F5.png", width = 15, height = 7, units = 'in', res = 300)
ggplot(db, aes(x=Raven_IQ,y=partial.resid)) + stat_smooth(method="loess", alpha=.2) + 
  geom_point( aes(colour=study, shape=sex), size=2.3)+ scale_color_manual(values = c("#00AFBB", "#E7B800", "#FC4E07"))+
stat_smooth(method="lm", formula=y ~ x, color="red", fill="red", alpha=.2)+
  ylab(paste("Component+Residual(", responseName(lm.temp),")", sep=""))+theme_linedraw()# xlab("Raven_IQ")+ylab("Residuals")
dev.off()
```

### Checking homoscedasticity assumption

This assumption can be checked by examining the Scale-location plot, also known as the spread-location plot. Ideally, in the plot, we are interested in detecting a pattern in the residuals. If the blue line is flat and horizontal with equally and randomly spread data points, then this might be evidence of homoscedasticity. If the red line has a positive slope or the data points are not randomly spread out, the homoscedasticity assumption might not hold. This plot shows that residuals are not distributed equally along with the ranges of predictors. In our example, this is not entirely the case; rigorously, we need continued investigation of the homogeneity of variance.

```{r,  echo = FALSE, fig.width = 10, fig.height=5}
db$stdresid <- rstandard(lm.temp)
db$rstandard <-  sqrt(abs(db$stdresid ))

png("F6.png", width = 15, height = 7, units = 'in', res = 300)
ggplot(db, aes(x=fitted,y=rstandard)) + stat_smooth(method="loess", alpha=.2) + 
  geom_point( aes(colour=study, shape=sex), size=2.3)+ scale_color_manual(values = c("#00AFBB", "#E7B800", "#FC4E07"))+
  xlab("Fitted Value")+ylab(expression(sqrt("|Standardized residuals|")))+ggtitle("Scale-Location")+  theme_linedraw()
dev.off()
   # p3<-ggplot(lm.temp, aes(.fitted, sqrt(abs(.stdresid))))
   #  p3<-p3+stat_smooth(method="loess", na.rm = TRUE, alpha=.2)+xlab("Fitted Value")
   #  p3<-p3+ylab(expression(sqrt("|Standardized residuals|")))
   #  p3<-p3+ggtitle("Scale-Location")+ geom_point(shape=1)+  theme_linedraw()
   #  p3
```

We are looking for any lawful curves or skewness in the data that suggest that the regression model is better or worse at predicting for specific levels of predictors. In this way, absolute studentized residuals refer to the absolute values (ignoring over or underfitting) of the quotient resulting from the division of a residual by an estimate of its standard deviation.

```{r,  echo = FALSE, fig.width = 10, fig.height=5}
# abs of studentized residuals
db$resid2 <-  abs(rstudent(lm.temp))
db$fitval <- db$fitted

# non.pos <- db$fitval <= 0
# 	if (any(non.pos)){
# 		db$fitval <- db$fitval[!non.pos]
# 		db$resid2<- db$resid2[!non.pos]
# 		n.non.pos <- sum(non.pos)
# 		warning("\n", n.non.pos, " negative", if(n.non.pos > 1) " fitted values" else " fitted value", " removed")
# 	}

# spreadLevelPlot(lm.temp, lwd=3, col.lines=c("red","blue"))

png("F7.png", width = 15, height = 7, units = 'in', res = 300)
ggplot(db, aes(x=fitval,y=resid2)) + stat_smooth(method="loess", alpha=.2) + 
  geom_point( aes(colour=study, shape=sex), size=2.3)+ scale_color_manual(values = c("#00AFBB", "#E7B800", "#FC4E07"))+ scale_y_continuous(trans = log2_trans())+geom_smooth(method="rlm", color="red", fill="red", alpha=.2)+
  xlab("Fitted Value")+ylab("Absolute Studentized Residuals")+ggtitle("Spread-Level Plot")+  theme_linedraw()
dev.off()
mod = rlm(log(resid2) ~ log(fitval), data=db)
p <- 1 - (coefficients(mod))[2]
cat('\nSuggested power transformation: ', p,'\n')
```

In this case, even though there appear to be some small curves in the smoother, the linear fit seems fair straight across the scale. This is about the border of homoscedastic and heteroscedastic. However, we can observe patterns in the plot, and more investigations are necessary. Also, the function calculates spread-stabilizing power transformation of variation; in this case, the power transformation $y^* = response^\lambda$ suggested is for $\lambda = 2.20>1$, which is clear *evidence of non-constant variance* in the error term.

### Testing the homoscedasticity assumption.

We can conduct tests to evaluate non-constant error variance. The functions `ncvTest()` and `bptest()` compute a score test for the hypothesis of constant error variance against the alternative that the error variance changes with the level of the response (fitted values) or with a linear combination of predictors. The default function `bptest()` uses the studentized Breusch-Pagan proposed by R. Koenker (1981) [A Note on Studentizing a Test for Heteroscedasticity. Journal of Econometrics 17, 107--112.] In contrast, the function `ncvTest()` performs the original version of the Breusch-Pagan test [T.S. Breusch & A.R. Pagan (1979), A Simple Test for Heteroscedasticity and Random Coefficient Variation. Econometrica 47, 1287--1294].

```{r,  echo = TRUE}
# Breusch & Pagan test
bptest(lm.temp)

# Koenker version of Breusch & Pagan test  
ncvTest(lm.temp)
```

Using the usual nominal significance values of tests (for example, 5%), we **reject** the hypothesis of homoscedasticity of errors using the two tests.

### Cheking Independence Assumption

We expect the absence of patterns in the plot that depict the standardized residuals and the index of observations on ideal conditions. The presence of patterns pinpoints a lack of independence of observations or some endogeneity not controlled. By visual inspection of this plot, we observe groups formations for the original model proposed by Gignac and Zajenkowski 2020's paper.

```{r,  echo = FALSE, fig.width = 10, fig.height=5}
png("F8.png", width = 15, height = 7, units = 'in', res = 300)
  p2<-ggplot(db, aes(x=ID, y=stdresid))+geom_point( aes(colour=study, shape=sex), size=2.3)+ scale_color_manual(values = c("#00AFBB", "#E7B800", "#FC4E07"))+theme_linedraw()
  p2+labs(x='Observation ID', 
       y='Standardized Residuals', 
       title='Standardized Residuals vs ID of observations')
  dev.off()
```

We carried a Durbin Watson to examine the autocorrelation of errors. The null hypothesis in this regard states that these errors are not autocorrelated; that is, it is helpful to verify that we have not violated the independence assumption (presence of lack-of-fit).

```{r,  echo = FALSE, fig.width = 14, fig.height=8}
durbinWatsonTest(lm.temp)
```

Looking at the results, we see that the $p$-value is 0, so the errors are autocorrelated. We have violated the independence assumption. Such a violation would suggest, for example, the need to include a grouping variable, in our case, sex or study variables, to improve rigorousness.

### Checking the assumption of normality of residuals and homoscedasticity

Generalised additive model for location scale and shape (GAMLSS) is a semi-parametric framework to regression models [RIGBY, R. A.; STASINOPOULOS, D. M. Generalized additive models for location, scale and shape. Journal of the Royal Statistical Society: Series C (Applied Statistics), Wiley-Blackwell, v. 54, n. 3, p. 507--554, 2005.]. The gamlss package together with the gamlss.dist package provides extensions to the lm() and glm() functions from the stats package and the gam() function from the gam package for generalised additive models. In the GAMLSS approach, regression-based relationships are fitted without strong assumptions between the response and the covariates. Here, we can use this potential to find distributional evidence of model errors and lack-of-fit. Based on normality assumptions on errors and homoscedasticity, the linear model parametric fit based on OLS and likelihood are equivalent. Also, the hypothesis of normality of residuals is valid for inferential purposes. Therefore, we can fit the model proposed by Gignac and Zajenkowski (2020) under the GAMLSS perspective.

```{r,  echo = FALSE}
# The model - GAMLSS aproach
temp <- gamlss(SAI_winsorized_IQ ~ Raven_IQ, data=db, family=NO, trace=FALSE)
```

By a simple visual inspection, we notice no strong evidence of normality deviation. Nonetheless, we see that some density curves in the plots show different variability, which can be additional evidence of non-constant dispersion.

```{r,  echo = FALSE, fig.width = 10, fig.height=6}
y.lim<-range(db$SAI_winsorized_IQ)
png("F9.png", width = 15, height = 7, units = 'in', res = 300)
plotSimpleGamlss( x= Raven_IQ, 
                  y=SAI_winsorized_IQ,
                  model=temp, data=db, pch=1,
                   x.val=seq(65,160,6), val=200, N=500, 
                  ylim=c(80, 170),
                  xlim=c(58, 135),
                  cols=viridis_pal(begin = 0.3, end = 0.95, alpha=0.4, option ="D", direction=1)(100))
dev.off()
```

Worm-plots (Buuren, 2001) can be used to identify some data characteristics that the fitted model does not adequately capture. The yellow dots on the graph indicate how far the residuals depart from their expected null value. A high (low) variance pattern is related to positive (negative) sloped dots. When the dots are U-shaped (inverted U-shaped), there is evidence of positive (negative) skewness. There is evidence of large (small) kurtosis when the dots are S-shaped with left bent up (down); there is evidence of large (small) kurtosis.} When the model's specification is the right one, all the dots (or nearly all of them) should be inside the two dashed semicircles. A commonly used rule is to conclude that the model is correctly specified when at least $95\%$ of the dots fall between these limits.

```{r,  echo = FALSE, fig.width = 10, fig.height=6}
png("F10.png", width = 15, height = 7, units = 'in', res = 300)
wp(temp)
dev.off()
```

In our case, we observe a U-shape and several points outside the confidence band, indicating that the hypothesis of normality and constant variance does not hold. The residuals QQ plot is helpful to check the normality assumption visually. The normal probability plot of residuals should approximately follow a straight line. In our example, all the points fall approximately along this reference line to assume normality. The histogram with kernel density of residuals is not necessary under the normality assumption (bimodality).

```{r,  echo = FALSE, fig.width = 10, fig.height=5}
db$sresid <- studres(lm.temp)
temp <- qqnorm(db$sresid, plot.it = FALSE)
db$xqq <- temp$x
db$yqq <- temp$y

p.1 <-ggplot(db, aes(x=xqq, y=yqq))+geom_point( aes(colour=study, shape=sex), size=2.3)+ scale_color_manual(values = c("#00AFBB", "#E7B800", "#FC4E07"))+ geom_abline()+xlab("Theoretical Quantiles")+ylab("Standardized Residuals")
p.1<-p.1+ggtitle("Normal Q-Q")
p.1 = p.1+theme_linedraw()

# p.1<-ggplot(lm.temp)+stat_qq(aes(sample = .stdresid)) + geom_abline()+xlab("Theoretical Quantiles")+ylab("Standardized Residuals")+
# p.1<-p.1+ggtitle("Normal Q-Q")
# p.1 = p.1+theme_linedraw()

p.2<-ggplot(db, aes(x=sresid ))+geom_histogram(aes(y=..density..), colour="black", fill="white")+ geom_density(alpha=.2, fill="#E69F00")+xlab("Theoretical Quantiles")+ylab("Density")+ggtitle("Distribution of Studentized Residuals")+theme_linedraw() 

png("F11.png", width = 15, height = 7, units = 'in', res = 300)
grid.arrange(p.1, p.2, ncol=2,   widths = c(2, 1.5), top="")
dev.off()
```

We then go on to perform the robust and nonparametric Shapiro-Wilk Normality Test.

```{r,  echo = FALSE}
# distribution of studentized residuals
shapiro.test(db$sresid)
```

As the $p$-value\< 0.5, it is clear that the residuals distribution deviates from a Normal distribution significantly.

### More on heteroskedasticity evidence

With the idea of a regressogram to evaluate heteroscedasticity evidence, we propose constructing piecewise functions in the intervals $\{B_k:=[t_k,t_{k+1}):t_k=t_0+hk,k∈Z\}$ (uniform breaks) or $\{B_k:=[t_k,t_{k+1}):t_k-{\rm quantile \ empiric},k∈Z\}$ (quantile breaks) and calculate for the $y$ (ouput) the Coefficient of dispersion based on the MAD $$CV_{\text{MAD}}=\widehat{\theta}_{\rm MAD} = \frac{\text{MAD}}{\text{Mdn}}$$ proposed by [Ospina, R., & Marmolejo-Ramos, F. (2019). Performance of some estimators of relative variability. Frontiers in Applied Mathematics and Statistics, 5, 43.] across the midpoints of $x$ (input) in $B_k$ where \begin{equation}
\label{eq:MAD}
\text{MAD}=1.4826\,\cdot\text{Median}\{|x-\widehat{Q}_{2}|\},
\end{equation} and $\text{Mdn}=\widehat{Q}_{2}=\widehat{Q}(0.5)=F_n(0.5)^{-1}$ is an estimator of location that is robust as it has a high breakdown (The breakdown value is the smallest fraction of contamination that can cause the estimator to take on values far from its value on the uncontamined data). Here, we choose the bandwidth $h$ by the Sturges rule for histograms and we chose the percentils for quantile breaks. The plots are produced in terms of the behaviour of the covariates against the response and in terms of the fitted values against the residuals.

```{r,  echo = FALSE, message=FALSE,warning=FALSE}
## mad/median:
madmedian <- function(x){
  mdn <- median(x)
  MAD <- 1.4826 * median(abs(x - mdn))
  MAD/mdn
}

# old coefficient of variations
oldcv <- function(x) sd(x)/mean(x)

regres.madmedian<- function(x,y, x.lab="X",y.lab="Y",main="TITLE"){
  xy <- data.frame(x=x,y=y)
  xy <- xy[order(xy$x),]
  nbins <- nclass.Sturges(y)
  temp <- mean(y)
  z <- cut(xy$x,breaks=seq(min(xy$x),max(xy$x),length=nbins+1),
           labels=1:nbins,include.lowest=TRUE)
  xyz <- data.frame(xy,z=z)
  MEANS <- c(by(xyz$y,xyz$z,FUN=madmedian))
  x.seq <- seq(min(x), max(x), length = nbins + 1)
  midpts <- (x.seq[-1] + x.seq[-(nbins + 1)])/2
  d2 <- data.frame(midpts = midpts, MEANS = MEANS)
  p <- ggplot(xyz, aes(x, y)) + geom_point() + ggtitle(main) + 
    xlab(x.lab) + ylab(y.lab) + theme(text = element_text(size = 20))

    p <- p + geom_vline(xintercept = x.seq[-c(1, nbins + 
                                                1)], linetype = "dashed", color = "blue")+theme_linedraw()
    
    p2 <-  ggplot()+ geom_point(data = d2, aes(x = midpts, y = MEANS), 
                        color = "red", shape = 18, size = 3)+geom_vline(xintercept = x.seq[-c(1, nbins + 
                                                1)], linetype = "dashed", color = "blue")
  
    p2 <- p2 + geom_line(data = d2, aes(x = midpts, y = MEANS), 
                       color = "red", linetype = "dashed")+xlab("midpoints")+ylab("Mad/median coef.")+theme(text = element_text(size = 20))+theme_linedraw()
    
#     g <- ggplotGrob(p)
# g2 <- ggplotGrob(p2)
# g <- rbind(g, g2, size = "first")
# g$widths <- unit.pmax(g$widths, g2$widths)
# grid.newpage()
# grid.draw(g)
#  #    
#  #    
 output <- plot_grid(p, p2, ncol=1)
 output
}

regres.madmedian(x = db$Raven_IQ, y = db$SAI_winsorized_IQ,  
                 x.lab = "Raven_IQ", 
                 y.lab = "SAI_winsorized_IQ", main = "Regressogram MadMedian - - Uniform breaks")

regres.madmedian(x = db$fitted, y = db$resid,  
                 x.lab = "fitted", 
                 y.lab = "Residuals", main = "Regressogram MadMedian - Uniform breaks")

################################################

regresquan.madmedian<- function(x,y, x.lab="X",y.lab="Y",main="TITLE"){
  xy <- data.frame(x=x,y=y)
  xy <- xy[order(xy$x),]
  brks <- as.numeric(with(xy, quantile(x, probs = seq(0,1, .1))))
  nbins <- length(brks) # nclass.Sturges(y)
  temp <- mean(y)
  z <- cut(xy$x,breaks= brks,
           # seq(min(xy$x),max(xy$x),length=nbins+1),
  labels=1:(nbins-1),include.lowest=TRUE)

  # z <- cut(xy$x,breaks=seq(min(xy$x),max(xy$x),length=nbins+1),
  #          labels=1:nbins,include.lowest=TRUE)
xyz <- data.frame(xy,z=z)
  MEANS <- c(by(xyz$y,xyz$z,FUN=madmedian))
  x.seq <- seq(min(x), max(x), length = nbins)
  midpts <- (x.seq[-1] + x.seq[-(nbins)])/2
  d2 <- data.frame(midpts = midpts, MEANS = MEANS)
  p <- ggplot(xyz, aes(x, y)) + geom_point() + ggtitle(main) +
    xlab(x.lab) + ylab(y.lab) + theme(text = element_text(size = 20))

    p <- p + geom_vline(xintercept = x.seq[-c(1, nbins)], linetype = "dashed", color = "blue")+theme_linedraw()

    p2 <-  ggplot()+ geom_point(data = d2, aes(x = midpts, y = MEANS),
                        color = "red", shape = 18, size = 3)+geom_vline(xintercept = x.seq[-c(1, nbins )], linetype = "dashed", color = "blue")

    p2 <- p2 + geom_line(data = d2, aes(x = midpts, y = MEANS),
                       color = "red", linetype = "dashed")+xlab("midpoints")+ylab("Mad/median coef.")+theme(text = element_text(size = 20))+theme_linedraw()

#     g <- ggplotGrob(p)
# g2 <- ggplotGrob(p2)
# g <- rbind(g, g2, size = "first")
# g$widths <- unit.pmax(g$widths, g2$widths)
# grid.newpage()
# grid.draw(g)
#  #
#  #
 output2 <- plot_grid(p, p2, ncol=1)
 output2
}
# 
regresquan.madmedian(x = db$Raven_IQ, y = db$SAI_winsorized_IQ,
                 x.lab = "Raven_IQ",
                 y.lab = "SAI_winsorized_IQ", main = "Regressogram MadMedian - Quantile breaks")

regresquan.madmedian(x = db$fitted, y = db$resid,  
                 x.lab = "fitted", 
                 y.lab = "Residuals", main = "Regressogram MadMedian - Quantile breaks")

```

We can observe in plots that regardless of the way the breaks are constructed, the $CV_{\text{MAD}}$ is not constant, giving further evidence to our findings on heteroscedasticity.

### Global Test of Model Assumptions

Finally, we can use a global test of model assumptions using the `gvlma` package in `R`. The package allows the user to check out (almost) all the ideas discussed so far. This package relies on the methodology proposed by Pena and Slate (2006). [Global validation of linear model assumptions, Journal of the American Statistical Association, 101(473):341-354.]

```{r,  echo = FALSE}
# model proposed by Gignac and Zajenkowski 2020
gvmodel.gobble <- gvlma(lm.temp) 
summary(gvmodel.gobble)
```

We observe in the output that Skewness assumptions (symmetry by supposition in the NULL hypothesis) do not hold. Rejection of the null hypothesis ($p$-value \< .05) indicates the need to transform the data or employ robust statistical techniques. Note that the `gvlma` package is a recent contribution to the `R` user community and might not be fully known and used by psychologists. In this sense, the presented results help the user to understand the relevance of proceeding with careful attention to the data behaviour.

```{r,  echo = TRUE}
sessionInfo()
```