fig03.Rmd

---
title: 'A guide to robust statistical methods in neuroscience: Figure 3'
author: "Rand R. Wilcox & Guillaume A. Rousselet"
date: "`r Sys.Date()`"
output:
  pdf_document:
    fig_caption: no
    number_sections: no
    toc: yes
    toc_depth: 2
---

Illustrate homoscedastidicity (variance homogeneity) in panel A, and heteroscedasticity (variance heterogeneity) in panel B.

# Dependencies
```{r}
library(ggplot2)
library(cowplot)
library(tibble)
source("./code/Rallfun-v40.txt")
source("./code/theme_gar.txt")
```

# generate data
```{r}
set.seed(44)

n <- 100 # sample size
xvals <- seq(20,80,10)
nx <- length(xvals)
he.vals <- seq(from = 0.5, to = 3, length.out = nx)
x <- rep(xvals, each = n)
y <- matrix(rnorm(n*nx), nrow = n)
y.ho <- y
y.he <- y
for(C in 1:nx){
  # homoscedastic
  y.ho[,C] <- 2 + y[,C] / sd(y[,C]) + xvals[C] / 7
  # heteroscedastic
  y.he[,C] <- 2 + he.vals[C] * y[,C] / sd(y[,C]) + xvals[C] / 7
}
```

# base R figure
```{r}
par(mfrow=c(1,2))
plot(x,y.ho,xlab='Age',ylab='DV', ylim=c(0, 25))
abline(lm(as.vector(y.ho)~x))
plot(x,y.he,xlab='Age',ylab='DV', ylim=c(0, 25))
abline(lm(as.vector(y.he)~x))
par(mfrow=c(1,1))
```

# `ggplot2` figure

## panel A: homoscedasticity
```{r}
# plot parameters
axis_size <- 14
axis_title_size <- 16

dfA <- tibble(x = x, y = as.vector(y.ho)) # make data frame

pA <- ggplot(dfA, aes(x, y)) + theme_gar +
  geom_smooth(method=lm,   # Add linear regression line
              se=FALSE, # Don't add shaded confidence region
              colour = "grey50", size = 1.5) +  
  geom_point(size = 2, alpha = 0.1) +
  scale_x_continuous(breaks = x) +
  scale_y_continuous(breaks = seq(0,25,5), limits = c(0, 25)) +
  theme(axis.text = element_text(size = axis_size), 
        axis.title = element_text(size = axis_title_size),
        panel.grid.minor = element_blank()) +
  xlab("Age") +
  ylab("Dependent Variable")
pA
```

## panel B: heteroscedasticity
```{r}
dfB <- tibble(x = x, y = as.vector(y.he)) # make data frame

pB <- ggplot(dfB, aes(x, y)) + theme_gar +
  geom_smooth(method=lm,   # Add linear regression line
    se=FALSE, # Don't add shaded confidence region
    colour = "grey50", size = 1.5) +  
  geom_point(size = 2, alpha = 0.1) +
  scale_x_continuous(breaks = x) +
  scale_y_continuous(breaks = seq(0,25,5), limits = c(0, 25)) +
  theme(axis.text = element_text(size = axis_size), 
    axis.title = element_text(size = axis_title_size),
    panel.grid.minor = element_blank()) +
  xlab("Age") +
  ylab("Dependent Variable")
pB
```

## combine panels into one figure
```{r, eval=FALSE}
cowplot::plot_grid(pA, pB,
  labels=c("A", "B"), 
  ncol = 1, 
  nrow = 2,
  rel_widths = c(1, 1), 
  label_size = 20, 
  hjust = -0.5, 
  scale=.95,
  align = "h")

# save figure
ggsave(filename='./figures/figure3.pdf',width=5,height=7) 
```