-
Notifications
You must be signed in to change notification settings - Fork 0
/
GAMLSS Tutorial for Psychological Data.Rmd
449 lines (330 loc) · 27.4 KB
/
GAMLSS Tutorial for Psychological Data.Rmd
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
295
296
297
298
299
300
301
302
303
304
305
306
307
308
309
310
311
312
313
314
315
316
317
318
319
320
321
322
323
324
325
326
327
328
329
330
331
332
333
334
335
336
337
338
339
340
341
342
343
344
345
346
347
348
349
350
351
352
353
354
355
356
357
358
359
360
361
362
363
364
365
366
367
368
369
370
371
372
373
374
375
376
377
378
379
380
381
382
383
384
385
386
387
388
389
390
391
392
393
394
395
396
397
398
399
400
401
402
403
404
405
406
407
408
409
410
411
412
413
414
415
416
417
418
419
420
421
422
423
424
425
426
427
428
429
430
431
432
433
434
435
436
437
438
439
440
441
442
443
444
445
446
447
448
449
---
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()
grid.arrange(p.1, p.2, ncol=2, widths = c(1.5, 1.5), top="balancing size by age and study groups")
```
### 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"))
}
}
p+ theme_light()
```
### 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}
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()
```
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)
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()
# 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
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")
```
### 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 ))
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()
# 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"))
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()
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}
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')
```
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)
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))
```
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}
wp(temp)
```
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()
grid.arrange(p.1, p.2, ncol=2, widths = c(2, 1.5), top="")
```
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.
This tutorial was done in a computer with the following characteristics.
```{r, echo = TRUE}
sessionInfo()
```