03-exploring-data.qmd

# Exploring and Transforming Data {#sec-chap03}


```{r}
#| label: setup
#| results: hold
#| include: false

base::source(file = "R/helper.R")
ggplot2::theme_set(ggplot2::theme_bw())
options(show.signif.stars = FALSE)
options(rgl.useNULL = TRUE) # Suppress separate window opened by {rgl}.

```

## Chapter section list {.unnumbered}

::::: {#obj-chap03}
:::: {.my-objectives}
::: {.my-objectives-header}
Table of content (TOC)
:::

::: {.my-objectives-container}

- Examing distribution (@sec-chap03-1)
    - Histograms (@sec-chap03-1-1)
    - Density estimation (@sec-chap03-1-2)
    - Quantile-comparison plots (@sec-chap03-1-3)
    - Boxplots (@sec-chap03-1-4)
- Examing relationships (@sec-chap03-2)
    - Scatterplots (@sec-chap03-2-1)
    - Parallel boxplots (@sec-chap03-2-2)
    - More on the plot function (@sec-chap03-2-3)
- Examining multivariate data (@sec-chap03-3)
    - Three-dimensional plots (@sec-chap03-3-1)
    - Scatterplot Matrices (@sec-chap03-3-2)
- Transforming data (@sec-chap03-4)
    - Logarithms (@sec-chap03-4-1)
    - Power transformations (@sec-chap03-4-2)
    - Transformations and Exploarative Data Analysis (EDA) {@sec-chap03-4-3}
    - Transforming resticted-range variables (@sec-chap03-4-4)
    - Other transformations (@sec-chap03-4-5)
- Point labelling and identification (@sec-chap03-5)
    - The identify() function (@sec-chap03-5-1)
    - Automatic point labelling with car::showLabels() (@sec-chap03-5-2)
    



:::
::::
:::::

## Examing distributions {#sec-chap03-1}

### Histograms {#sec-chap03-1-1}

> The shape of the histogram is determined in part by the number of bins and the location of their boundaries. With too few bins, the histogram can hide interesting features of the data, while with too many bins, the histogram is very rough, displaying spurious features of the data.

The help page for `ggplot2::geom_histogram()` recommends: "You should always override this [default] value, exploring multiple widths to find the best to illustrate the stories in your data."

There are several algorithm to calculate an optimal number of bins depending of the sample size and distribution. Fox/Weisberg mention the rule by Freedman and Diaconis [-@freedman1981a]:

$$
\frac{n^\frac{1}{3}(max-min)}{2(Q_{3}-Q_{1})}
$$ {#eq-chap03-FD}

:::::{.my-example}
:::{.my-example-header}
:::::: {#exm-chap03-histograms}
: Histograms
::::::
:::
::::{.my-example-container}

::: {.panel-tabset}

###### base R

:::::{.my-r-code}
:::{.my-r-code-header}
:::::: {#cnj-chap03-default-base-r-histogram}
: Default base R histogram
::::::
:::
::::{.my-r-code-container}
::: {#lst-chap03-default-base-r-histogram}
```{r}
#| label: default-base-r-histogram

Prestige <- carData::Prestige

base::with(Prestige, graphics::hist(income))
```
Default	base R histogram of income in the Canadian occupational-prestige data
:::

::::
:::::

###### FD rule

:::::{.my-r-code}
:::{.my-r-code-header}
:::::: {#cnj-chap03-compute-bins-FD}
: Compute number of bins with the Freedman-Diaconis rule [-@freedman1981a]
::::::
:::
::::{.my-r-code-container}
::: {#lst-compute-bins-FD}
```{r}
#| label: compute-bins-FD


Income <- Prestige$income
    
base::ceiling(
    base::length(Income)^(1/3) * 
    (base::max(Income) - base::min(Income)) / 
    (2 * (stats::quantile(Income, 0.75) 
          - stats::quantile(Income, 0.25)))[[1]]
)


```

Compute number of bins of `income` with the Freedman-Diaconis rule [-@freedman1981a]
:::

::::
:::::


###### base R 2

:::::{.my-r-code}
:::{.my-r-code-header}
:::::: {#cnj-chap03-fd-base-r-histogram2}
: Base R histogram with `breaks = "FD"`
::::::
:::
::::{.my-r-code-container}
::: {#lst-chap03-fd-base-r-histogram2}
```{r}
#| label: fd-base-r-histogram2

base::with(Prestige, graphics::hist(income, breaks = "FD"))
```

Base R histogram: Number of bins computed after Freedman and Diaconis [-@freedman1981a]
:::

::::
:::::



###### ggplot2

:::::{.my-r-code}
:::{.my-r-code-header}
:::::: {#cnj-chap03-default-ggplot2-histogram}
: Histogram with {**ggplot2**} (30 bins default)
::::::
:::
::::{.my-r-code-container}
::: {#lst-chap03-default-ggplot2-histogram}    
```{r}
#| label: default-ggplot2-histogram

Prestige |> 
    ggplot2::ggplot(
        ggplot2::aes(x = income)
    ) +
    ggplot2::geom_histogram(
        fill = "grey",
        color = "black"
        )
```

Histogram with {**ggplot2**} with default number of bins (30)
:::

::::
:::::

###### ggplot2 2

:::::{.my-r-code}
:::{.my-r-code-header}
:::::: {#cnj-chap03-default-ggplot2-histogram2}
: Histogram with {**ggplot2**} with bin number computed after Freedman & Diaconis
::::::
:::
::::{.my-r-code-container}
::: {#lst-chap03-default-ggplot2-histogram2}    
```{r}
#| label: default-ggplot2-histogram2

Prestige |> 
    ggplot2::ggplot(
        ggplot2::aes(x = income)
    ) +
    ggplot2::geom_histogram(
        fill = "grey",
        color = "black",
        bins = grDevices::nclass.FD(Prestige$income)
        )
```

Histogram with {**ggplot2**}. Bin number computed after Freedman and Diaconis [-@freedman1981a]
:::

::::
:::::

:::

::::
:::::

A special type of histograms are *stem-and-leaf displays*, which are histograms that encode the numeric data directly in their bars. I do not find them very usefull, because they are different to interpret and do not have essential advantages. They may be helpful for small data sets as a kind of paper-and-pencil method for visualizing.

Stem-and-leaf graphs can be produced with the `graphics::stem()` function or using the `aplpack::stem.leaf()` function. 

:::::{.my-resource}
:::{.my-resource-header}
:::::: {#lem-chap03-aplpack}
: Another Plot Package {aplpack}
::::::
:::
::::{.my-resource-container}
##### Functions for drawing some special plots

Besides plotting stem-and-leaf graphs with `aplpack::stem.leaf()` from the {**aplpack**} package, I noticed that there are some other interesting graphic functions in this package, which I should explore.

For instance `aplpack::plotsummary()` is an interesting function that shows important characteristics of the variables of a data set. For each variable a plot is computed consisting of a `r glossary("Bar Charts", "barplot")`, an `r glossary("ecdf")`, a `r glossary("Density Plots", "density trace")` and a `r glossary("Boxplots", "boxplot")`. 

See @pak-aplpack.
::::
:::::









### Density Estimation {#sec-chap03-1-2}

*Nonparametric density estimation* often produces a more satisfactory representation of the distribution of a numeric variable than a traditional histogram. Unlike a histogram, a nonparametric density estimate is continuous and so it doesn’t dissect the range of a numeric variable into discrete bins.

`r glossary("Kernel density estimation")` (KDE) is the application of kernel smoothing for probability density estimation. The bandwith controls	the	degree of smoothness of	the density	estimate:

The bandwidth of a density estimate is the continuous analog of the bin width of a histogram: If the bandwidth is too large, then the density estimate is smooth but biased as an estimator of the true density, while if the bandwidth is too small, then bias is low but the estimate is too rough and the variance of the estimator is large.

The `adaptiveKernel()` function in the {**car**} package employs an algorithm that uses different bandwidths depending on the local observed density of the data, with smaller bandwidths in dense regions and larger bandwidths in sparse regions.




:::::{.my-example}
:::{.my-example-header}
:::::: {#exm-chap03-density-estimation}
: Density estimation
::::::
:::
::::{.my-example-container}

::: {.panel-tabset}

###### stats::density() 1

:::::{.my-r-code}
:::{.my-r-code-header}
:::::: {#cnj-chap03-density1}
: Nonparametric fixed-bandwidth and adaptive-bandwidth kernel density estimates `(adjust = 1)`
::::::
:::
::::{.my-r-code-container}
::: {#lst-chap03-density1}
```{r}
#| label: density1

with(Prestige, {
    hist(
        income,
        freq = FALSE,
        ylim = c(0, 1.5e-4),
        breaks = "FD",
        main = ""
    )
    lines(density(income, from = 0), lwd = 3, lty = 2)
    lines(car::adaptiveKernel(income, from = 0, adjust = 1),
          lwd = 2,
          lty = 1) # solid line
    rug(income)
    legend(
        "topright",
        c("Fixed bandwidth", "Adaptive bandwidth"),
        lty = 2:1, # dashed with proportion 2:1
        lwd = 2,
        inset = .02
    )
    box()
})
```
Nonparametric fixed-bandwidth and adaptive-bandwidth kernel density estimates (adjust = 1) for the distribution of income in the `Prestige` data set; a density histogram of income is also shown. The rug-plot at the bottom of the graph shows the location of the income values.
:::

***

The command to draw the graph in @lst-chap03-density1 is relatively complicated and thus requires some explanation: 

- The `base::with()` function is used as usual to pass the data frame `Prestige` to the second argument. Here the second argument is a compound expression consisting of all the commands between the initial { and the final }. 
- The call to `graphics::hist()` draws a histogram in density scale, so the areas of all the bars in the histogram sum to 1. 
- The argument `main=""` suppresses the title for the histogram. 
- The `ylim` argument sets the range for the y-axis to be large enough to accommodate the adaptive-kernel density estimate. The value 1.5e-4 is in scientific notation, 1.5 × 10−4 = 0.00015. 
- The fixed-bandwidth and adaptive-bandwidth kernel estimates are computed, respectively, by `stats::density()` and `car::adaptiveKernel()`. 
- In each case, the result returned by the function is then supplied as an argument to the `graphics::lines()` function to add the density estimate to the graph. 
- The argument `from = 0` to both `density()` and `adaptiveKernel()` ensures that the density estimates don’t go below income = 0. 
- The call to `graphics::rug()` draws a rug-plot (one-dimensional scatterplot) of the data at the bottom of the graph. 
- The remaining two commands add a legend and a frame around the graph.




::::
:::::

Both nonparametric density estimates and the histogram suggest a mode at around $5,000, and all three show that the distribution of income is rightskewed. The fixed-bandwidth kernel estimate has more wiggle at the right where data are sparse, and the histogram is rough in this region, while the adaptive-kernel estimator is able to smooth out the density estimate in the low-density region. 


###### stats::density() 2

:::::{.my-r-code}
:::{.my-r-code-header}
:::::: {#cnj-chap03-density2}
: Nonparametric fixed-bandwidth and adaptive-bandwidth kernel density estimates `(adjust = 0.5)`
::::::
:::
::::{.my-r-code-container}
::: {#lst-chap03-density2}    
```{r}
#| label: density2

with(Prestige, {
    hist(
        income,
        freq = FALSE,
        ylim = c(0, 1.5e-4),
        breaks = "FD",
        main = ""
    )
    lines(density(income, from = 0), lwd = 3, lty = 2)
    lines(car::adaptiveKernel(income, from = 0, adjust = 0.5),
          lwd = 2,
          lty = 1) # solid
    rug(income)
    legend(
        "topright",
        c("Fixed bandwidth", "Adaptive bandwidth"),
        lty = 2:1, # dashed with proportion 2:1
        lwd = 2,
        inset = .02
    )
    box()
})



```

Nonparametric fixed-bandwidth and adaptive-bandwidth kernel density estimates (adjust = 0.5) for the distribution of `income` in the `Prestige` data set; a density histogram of `income` is also shown. The rug-plot at the bottom of the graph shows the location of the `income` values.
:::

::::
:::::

###### ggplot2::geom_density()

:::::{.my-r-code}
:::{.my-r-code-header}
:::::: {#cnj-chap03-density-ggplot2}
: Histogram, rug and smoothed density estimate
::::::
:::
::::{.my-r-code-container}
::: {#lst-density-ggplot2}
```{r}
#| label: density-ggplot2

Prestige |> 
    ggplot2::ggplot(
        ggplot2::aes(x = income)
    ) +
    ggplot2::geom_histogram(
        ggplot2::aes(y = ggplot2::after_stat(density)),
        color = "black",
        fill = "grey",
        bins = grDevices::nclass.FD(Prestige$income)
    ) +
    ggplot2::geom_density(
        adjust = 1,
        kernel = "gaussian",
        color = "red",
        linewidth = 1
    ) +
    ggplot2::geom_rug()
```

Histogram with a kernel density estimate overlaid, which is a smoothed version of the histogram. The rug-plot at the bottom of the graph shows the location of the `income` values from the `Prestige` dataset.
:::

::::
:::::

::: {.callout-warning #wrn-chap03-adaptive-kernel}
##### Don't know how to inlcude adaptive kernel density estimation

The {**car**} package has with `car::densityPlot()` an additional function, that both computes and draws density estimates with either a fixed or adaptive kernel. I do not know how to include the adaptive kernel density estimation from the {**car**} package to get a full reproduction of book’s Figure 3.3 with {**ggplot2**}.
:::




###### my_hist_dnorm()

:::::{.my-r-code}
:::{.my-r-code-header}
:::::: {#cnj-chap03-my-hist-dnorm}
: Using my own function to create a histogram with density estimate and overlaid dnorm curve
::::::
:::
::::{.my-r-code-container}
::: {#lst-chap03-my-hist-dnorm}
```{r}
#| label: my-hist-dnorm

my_hist_dnorm(Prestige, Prestige$income)
```

Histogram with density estimation and overlaid dnorm curve
:::

***

I have my own function `my_hist_dnorm()` adapted by supplying the Freeman-Diaconis rule [-@freedman1981a] as default value for the number of bins.

Additionally I have written with `my_nbins()` another small functions which returns the number of bins according to the Freeman-Diaconis rule.

::::
:::::


:::

::::
:::::

### Quantile-comparison plots {#sec-chap03-1-3}

A *quantile-comparison plot*, or quantile-quantile plot (`r glossary("Q-Q-plot")`), provides an effective graphical means of making the comparison, plotting the ordered data on the vertical axis against the corresponding quantiles of the reference distribution on the horizontal axis. If the data conform to the reference distribution, then the points in the quantile-comparison plot should fall close to a straight line, within sampling error.

:::::{.my-example}
:::{.my-example-header}
:::::: {#exm-chap03-qq-plots}
: Quantile-quantile plots
::::::
:::
::::{.my-example-container}

::: {.panel-tabset}

###### stats::qqnorm()

:::::{.my-r-code}
:::{.my-r-code-header}
:::::: {#cnj-chap03-stats-qqnorm}
: Base R: Quantile-Quantile Plot
::::::
:::
::::{.my-r-code-container}
::: {#lst-chap03-stats-qqnorm}
```{r}
#| label: stats-qqnorm

stats::qqnorm(Prestige$income,
              ylab = "Income")
stats::qqline(Prestige$income)
```
Normal quantile-comparison plot	for	`income` from the `Prestige` dataset.
:::

***
Many points, especially at the right of the graph, are far from the line of the theoretical quantiles. We have therefore evidence that the distribution of `income` is not like a sample from a normal population.

::::
:::::


###### car::qqPlot()

:::::{.my-r-code}
:::{.my-r-code-header}
:::::: {#cnj-chap03-car-qqplot}
: car: Quantile-quantile plot
::::::
:::
::::{.my-r-code-container}
::: {#lst-chap03-car-qqplot}    
```{r}
#| label: car-qqplot

car::qqPlot( ~ income, data = Prestige, id = list(n = 3))
```

Normal quantile-comparison plot for `income.` 
:::
***
The broken lines give a pointwise 95% confidence envelope around the fitted solid line. Three points were labeled automatically. Because many points, especially at the right of the graph, are outside the confidence bounds, we have evidence that the distribution of `income` is not like a sample from a normal population.

::::
:::::

The function `car::qqPLot()` has several advantages:

1) It produces a pointwise 95% confidence envelope around the fitted solid line.
2) It labels the most extreme data points.
3) The `car::qqPlot()` function can be used to plot the data against *any* reference distribution.

###### my_qq_plot()

:::::{.my-r-code}
:::{.my-r-code-header}
:::::: {#cnj-chap03-my-qq-plot}
: Q-Q-plot using my own function applying {**ggplot2**}
::::::
:::
::::{.my-r-code-container}
::: {#lst-chap03-my-qq-plot}
```{r}
#| label: my-qq-plot

my_qq_plot(Prestige, Prestige$income)
```

Q-Q-plot: Comparing data points against a normal distribution 
:::

***

My own function lacks the confidence interval cannot label the most extreme points. I have to think if and how I could add these features to `my_qq_plot()`.

::::
:::::

###### chi-square

:::::{.my-r-code}
:::{.my-r-code-header}
:::::: {#cnj-chap03-chi-square-illustration}
: car::qqPlot: Chi-square illustration
::::::
:::
::::{.my-r-code-container}
::: {#lst-chap03-chi-square-illustration}
```{r}
#| label: chi-square-illustration

base::set.seed(124) # for reproducibility
car::qqPlot(
    rchisq(100, 3), 
    distribution = "chisq", df = 3,
    id = FALSE)
```

Quantile-comparison plot of a sample of size n = 100 from the $χ^2(3)$ distribution against the distribution from which the sample was drawn.
:::
***

The argument `df = 3` to `car::qqPlot()` is passed by it to the `stats::dchisq()` and `stats::qchisq()` functions. The points should, and do, closely match the straight line on the graph, with the fit a bit worse for the larger values in the sample. The confidence envelope suggests that these deviations for large values are to be expected, as they reflect the greater variability of sampled values in the long right tail of the $X^2(3)$ density function.

::::
:::::



:::

::::
:::::

The `car::qqPlot()` function can be used to plot the data against any reference distribution for which there are `r glossary("quantile")` and density function in R. You have simply to specify the root word for the distribution. For	example, the root	for	the	normal distribution is "norm", with	density	function	`stats::dnorm()` and quantile function `stats::qnorm()`. See also [chapter 8](https://rstudio.github.io/r-manuals/r-intro/Probability-distributions.html) of the the Quarto manual of [An Introduction to R](https://rstudio.github.io/r-manuals/r-intro/). [@rcoreteam2024]

![Probability functions in R](img/chap03-r-probability-functions.png){#fig-03-1 
fig-alt='List of probability functions in R as a two part table with five columns: It shows the name of the "Distribution" (column 1), such as "normal" or "chi-square", the "Design or Mass Function" (column 2), such as "dnorm(x, mean = 0, sd = 1)" or "qchisq(n, df)", the "Quantile Function" (column 3), such as "qnorm(p, mean = 0, sd = 1)" the "Distribution Function" column 4), such as "pnorm(q, mean = 0, sd = 1)" or "pchisq(q, df)" and the "Random Number Function" (column 5), such as "rnorm(n, mean = 0, sd = 1)" or "rchisq(n, df)"' fig-align="center" 
width="70%"}

An illustration is shown with @lst-chap03-chi-square-illustration: The `rchisq()` function is used to generate a random sample from the chi-square distribution with three df (`r glossary("degrees of freedom")`) and then plotted the sample against the distribution from which it was drawn.

### Boxplots {#sec-chap03-1-4}

Although boxplots are most commonly used to compare distributions	among groups, they can also be drawn to summarize the distribution of a numeric variable in a single sample, providing a quick check of symmetry and the presence of outliers.

:::::{.my-example}
:::{.my-example-header}
:::::: {#exm-chap03-boxplot}
: Boxplots
::::::
:::
::::{.my-example-container}

::: {.panel-tabset}

###### graphics::boxplot

:::::{.my-r-code}
:::{.my-r-code-header}
:::::: {#cnj-chap03-boxplot-base-r}
: graphics::boxplot(): Boxplot of income
::::::
:::
::::{.my-r-code-container}
::: {#lst-chap03-boxplot-base-r}
```{r}
#| label: boxplot-base-r

graphics::boxplot(Prestige$income)
```
Boxplot of income, plotted with base R `graphics::boxplot()`
:::

::::
:::::


###### car::Boxplot

:::::{.my-r-code}
:::{.my-r-code-header}
:::::: {#cnj-chap03-boxplot-car}
: car::Boxplot(): Boxplot of income
::::::
:::
::::{.my-r-code-container}
::: {#lst-chap03-boxplot-car}    
```{r}
#| label: boxplot-car

car::Boxplot(~ income, data = Prestige)
```

Boxplot of income, plotted with `car::Boxplot()`. Several outlying cases were labeled automatically.
:::

***

`car::Boxplot()` adds automatic identification of outlying values (by default, up to 10), points that are shown individually in the boxplot. Points identified as outliers are those beyond the *inner* fences, which are 1.5 times the interquartile range below the first quartile and above the third quartile.

The	names shown	in the output are the cases that are labeled on the graph and are drawn from the row names of the `Prestige` data set.

::::
:::::





:::

::::
:::::

The `ggplot2::geom_boxplot()` draws boxplots but has no options to label the outliers. To reproduce Figure 3.6 I need to add some code to the compute the outliers and label them with {**ggrepel**}. 

:::::{.my-experiment}
:::{.my-experiment-header}
:::::: {#def-chap03-boxplot-outlier-labelled}
: {ggplot2} boxplot with outlier labelled
::::::
:::
::::{.my-experiment-container}

::: {.panel-tabset}

###### default

:::::{.my-r-code}
:::{.my-r-code-header}
:::::: {#cnj-chap03-boxplot-ggplot2}
: Boxplot of income (default)
::::::
:::
::::{.my-r-code-container}
::: {#lst-chap03-boxplot-ggplot2}
```{r}
#| label: boxplot-ggplot2


Prestige <- carData::Prestige

Prestige |> 
    ggplot2::ggplot(
        ggplot2::aes(y = income)
    ) +
    ggplot2::geom_boxplot() 
```

Boxplot of income, plotted with `ggplot2::geom_boxplot()`
:::

::::
:::::


###### outliers

:::::{.my-r-code}
:::{.my-r-code-header}
:::::: {#cnj-chap03-compute-boxplot-outliers}
: Boxplot of income with outliers labelled
::::::
:::
::::{.my-r-code-container}
::: {#lst-chap03-compute-boxplot-outliers}  
```{r}
#| label: compute-boxplot-outliers

my_find_boxplot_outlier <- function(x) {
  return(x < quantile(x, .25) - 1.5*IQR(x) | x > quantile(x, .75) + 1.5*IQR(x))
}


prestige1 <- Prestige |> 
  dplyr::mutate(outlier = my_find_boxplot_outlier(income)) |> 
  tibble::rownames_to_column(var = "ID") |>
  dplyr::mutate(outlier =
        dplyr::case_when(outlier == TRUE ~ ID,
                         outlier == FALSE ~ "")
  )

set.seed(42)

prestige1 |> 
    ggplot2::ggplot(
        ggplot2::aes(y = income)
    ) +
    ggrepel::geom_text_repel(
      ggplot2::aes(
        x = 0,
        label = outlier
        ),
      nudge_x      = 0.05,
      direction    = "y",
      hjust        = 0,
    ) +
    ggplot2::geom_boxplot()

```

Boxplot of income with outliers labelled plotted with `ggplot2::geom_boxplot()` and using `ggrepel::geom_text_repel()` 
:::

::::
:::::

:::

::::
:::::

***


## Examing relationships {#sec-chap03-2}

### Scatterplots {#sec-chap03-2-1}

Understanding and using scatterplots is at the heart of regression analysis!

`r glossary("Scatterplot", "Scatterplots")` are useful for studying the mean and variance functions in the regression of the y-variable on the x-variable. In addition, scatterplots can help identify `r glossary("outliers")`, points that have values of the response far different from the expected value, and `r glossary("leverage")` points, cases with extremely large or small values on the x-axis.

:::::{.my-example}
:::{.my-example-header}
:::::: {#exm-chap03-scatterplot}
: Scatterplots
::::::
:::
::::{.my-example-container}

::: {.panel-tabset}

###### graphics::plot()

:::::{.my-r-code}
:::{.my-r-code-header}
:::::: {#cnj-chap03-scatterplot-base-r}
: Simple scatterplot using base R `graphics::plot()` function
::::::
:::
::::{.my-r-code-container}
::: {#lst-chap03-scatterplot-base-r}
```{r}
#| label: scatterplot-base-r

graphics::plot(prestige ~ income, data = Prestige)
```
Scatterplot of `prestige` versus `income` for the Canadian occupational-prestige data
:::
***

The scatterplot in @lst-chap03-scatterplot-base-r is a *summary graph* for the regression problem in which `prestige` is the response and `income` is the predictor. As our eye moves from left to right across the graph, we see how the distribution of `prestige` changes as `income` increases. In technical terms, we are visualizing the *conditional distributions* of `prestige` given the value of `income`.

::::
:::::

The overall story here is that as `income` increases, so does `prestige`, at least up to about $10,000, after which the value of `prestige` stays more or less fixed on average at about 80.


###### car::scatterplot

:::::{.my-r-code}
:::{.my-r-code-header}
:::::: {#cnj-chap03-scatterplot-car}
: Enhanced scatterplot with `car::scatterplot()`
::::::
:::
::::{.my-r-code-container}
::: {#lst-chap03-scatterplot-car}    
```{r}
#| label: scatterplot-car

car::scatterplot(
    prestige ~ income,
    data = Prestige,
    id = list(n = 4)
    )
```

Enhanced scatterplot of prestige versus income produced by the `car::scatterplot()` function. Four points were identified using the `id` argument.
:::

***

- **Points** correspond to the pairs of (`income`, `prestige`) values in the Prestige data frame, which is supplied in the data argument. 
- **The thick solid straight line** in the scatterplot is the simple linear regression of `prestige` on `income` fit by ordinary least squares (`r glossary("OLS")`). 
- **The dashed line** is fit using a nonparametric *scatterplot smoother*, and it provides an estimate of the mean function that does not depend on a parametric regression model, linear or otherwise. 
- **The two dash-dotted lines** provide a nonparametric estimate of the square root of the variance function (i.e., the conditional standard deviation of the response), based on separately smoothing the positive and negative residuals from the fitted *nonparametric* mean function. 
- Also shown on each axis are **marginal boxplots** of the plotted variables, summarizing the univariate distributions of x and y. 
- The only optional argument is `id=list (n=4)`, to identify by row name the four **most extreme cases**, where by default “extreme” means farthest from the point of means using Mahalanobis distances.

All these features, and a few others that are turned off by default, can be modified by the user, including the color, size, and symbol for points; the color, thickness, and type for lines; and inclusion or exclusion of the least-squares fit, the mean smooth, variance smooth, and marginal boxplots. See help ("scatterplot") for the available options.



::::
:::::

In @lst-chap03-scatterplot-car the least squares line cannot match the obvious curve in the mean function that is apparent in the smooth fit. So modeling the relationship of `prestige` to `income` by simple linear regression is likely to be inappropriate.

:::::{.my-resource}
:::{.my-resource-header}
:::::: {#lem-resource-text}
: Nonparametric regression as an appendix to the R Companion
::::::
:::
::::{.my-resource-container}
@lst-chap03-scatterplot-car uses as an plot enhancement in the {**car**} package a smoother, designed to help us extract information from a graph. 

*Nonparametric regression*, in which smoothers are substituted for more traditional regression models, is described in an [online appendix to the R Companion](https://www.john-fox.ca/Companion/appendices/Appendix-Nonparametric-Regression.pdf).
::::
:::::


###### ggplot2::geom_point()

:::::{.my-r-code}
:::{.my-r-code-header}
:::::: {#cnj-chap03-scatterplot-ggplot2}
: Scatterplot with `ggplot2::geom_point()`
::::::
:::
::::{.my-r-code-container}
::: {#lst-chap03-scatterplot-ggplot2}
```{r}
#| label: scatterplot-ggplot2
#| results: hold

gg_Prestige <- carData::Prestige |> 
    ggplot2::ggplot(
        ggplot2::aes(
            x = income,
            y = prestige
        )
    ) +
    ggplot2::geom_point(
        shape = "circle open",
        size = 2
    ) +
    ggplot2::stat_smooth(
        formula = 'y ~ x',
        method = "loess",
        span = .75,
        level = 0.95
    ) +
    ggplot2::stat_smooth(
        formula = 'y ~ x',
        method = "lm",
        se = FALSE
    ) +
    ggplot2::coord_cartesian(
      xlim = c(0, 26000),
      ylim = c(10, 90),
      expand = TRUE,
      default = FALSE,
      clip = "on"
    )

save_data_file("chap03", gg_Prestige, "gg_Prestige.rds")

gg_Prestige

```

Scatterplot of `prestige` versus `income` produced by different functions of the {**ggplot2**} package. 
:::
***

The default smoother employed by `ggplot2::geom_smooth()` depends on the size of the largest group. `stats::loess()` would have been used here because we have less than 1,000 observation. But I have applied the `loess() `function explicitly to prevent a warning message.

The span in loess is the fraction of the data used to determine the fitted value of y at each x. A larger span therefore produces a smoother result, and too small a span produces a fitted curve with too much wiggle. The trick is to select the smallest span that produces a sufficiently smooth regression mean function. The default span of 2/3 works well most of the time but not always




::::
:::::

Besides the labelled points and the marginal boxplots, I could not create smoother with the similar width as in the {**car**} example. I assume it has to do with "the *nonparametric estimate* of the square root of the variance function (i.e., the conditional standard deviation of the response), based on separately smoothing the positive and negative residuals from the fitted nonparametric mean function. "

###### categorical

:::::{.my-r-code}
:::{.my-r-code-header}
:::::: {#cnj-chap03-scatterplot-car-categorical}
: Conditioning on a categorical variable
::::::
:::
::::{.my-r-code-container}
::: {#lst-chap03-scatterplot-car-categorical}
```{r}
#| label: scatterplot-car-categorical

Prestige$type <- factor(Prestige$type, levels = c("bc", "wc", "prof"))
car::scatterplot(
    prestige ~ income | type,
    data = Prestige,
    legend = list(coords = "bottomright", inset = 0.1),
    smooth = list(span = 0.9),
    col = ggokabeito::palette_okabe_ito(c(1, 3, 5))
)


xtabs( ~ type, data = Prestige)

glue::glue(" ")
glue::glue("bc = blue color, wc = white color, prof = professional")
```

Scatterplot of `prestige` versus `income`, coded with {**car**} by `type` of occupation. The span of the loess smoother is set to 0.9.
:::
***

- **The variables** for the coded scatterplot are given in a formula as `y ~ x | g`, which is interpreted as plotting y on the vertical axis, x on the horizontal axis, and marking points according to the value of g (or “y versus x given g”).
- **The legend** for the graph, automatically generated by the scatterplot () function, is placed by default above the plot; we specify the legend argument to move the legend to the lower-right corner of the graph.
- We select a **larger span**, span=0.9, for the scatterplot smoothers than the default (span=2/3) because of the small numbers of cases in the occupational groups.
- As **color palette** I used the colorblind friendly Okabe Ito scale from {**ggokabeito**} instead of the default colors from {**car**} .

The default smoother employed by `car::scatterplot()` is the loess smoother. The span in loess is the fraction of the data used to determine the fitted value of y at each x. A larger span therefore produces a smoother result, and too small a span produces a fitted curve with too much wiggle. The trick is to select the smallest span that produces a sufficiently smooth regression mean function. The default span in {**car**} of 2/3 works well most of the time but not always. (In {**ggplot2**} the default span is 3/4.)

::::
:::::

The nonlinear relationship in @lst-chap03-scatterplot-car has disappeared, and we now have three reasonably linear regressions with different slopes. The slope of the relationship between `prestige` and `income` looks steepest for blue-collar occupations and least steep for professional and managerial occupations.

###### categorical2

:::::{.my-r-code}
:::{.my-r-code-header}
:::::: {#cnj-chap03-scatterplot-ggplot2-categorical}
: Numbered R Code Title
::::::
:::
::::{.my-r-code-container}
::: {#lst-chap03-scatterplot-ggplot2-categorical}
```{r}
#| label: scatterplot-ggplot2-categorical

Prestige |> 
    tidyr::drop_na(type) |> 
    ggplot2::ggplot(
        ggplot2::aes(
            x = income,
            y = prestige,
            group = type,
            color = type,
            shape = type
        ), 
    ) +
    ggplot2::geom_point(
        size = 2
    ) +
    ggplot2::stat_smooth(
        formula = 'y ~ x',
        method = "lm",
        se = FALSE
    ) +
    ggplot2::stat_smooth(
        linetype = "dashed",
        formula = 'y ~ x',
        method = "loess",
        se = FALSE, 
        span = .95,
    ) +
    ggokabeito::scale_color_okabe_ito(
        order = c(1, 3, 5, 2, 4, 6:9)
    ) +
    ggplot2::scale_shape_manual(values = 1:3)


```

Scatterplot of `prestige` versus `income`, coded with {**ggplot2**} by `type` of occupation. The span of the loess smoother is set to 0.95.
:::

:::::: {#tdo-chap03-scatterplot-ggplot2-categorical}
:::::{.my-checklist}
:::{.my-checklist-header}
TODO: Better legend inside graphic panel
:::
::::{.my-checklist-container}
- Add legend for linetype
- Create legend of occupation type without line, just the shape
- Put legend bottom right inside the graph panel
::::
:::::
How to improve the legend of @lst-chap03-scatterplot-ggplot2-categorical
:::


::::
:::::

###### Jitter (car)

:::::{.my-r-code}
:::{.my-r-code-header}
:::::: {#cnj-chap03-scatterplot-jitter-car}
: Scatterplots of `vocabulary` by `education` using {**car**}
::::::
:::
::::{.my-r-code-container}
::: {#lst-chap03-scatterplot-jitter-car}
```{r}
#| label: scatterplot-jitter-car
#| cache: true

Vocab <- carData::Vocab

car::scatterplot(vocabulary ~ education, data = Vocab, main = "(a)")

car::scatterplot(
  vocabulary ~ education,
  data = Vocab,
  jitter = list(x = 2, y = 2),
  cex = 0.01,
  col = "darkgray",
  smooth = list(
    span = 1 / 3,
    col.smooth = "black",
    col.spread = "black"
  ),
  regLine = list(col = "black"),
  main = "(b)"
)
```

Scatterplots of vocabulary by education with {**car**}: (a) using defaults, (b) jittering the points.
:::
***

- The argument `jitter=list (x=2, y=2)` adds small random numbers to the x and y coordinates of each plotted point. The values x=2 and y=2 specify the degree of jittering relative to the default amount, in this case twice as much jittering. The amounts of jitter used were determined by trial and error to find the choices that provide the most information. 
- The argument `cex=0.01` reduces the size of the circles for the plotted points to 1% of their normal size, and `col="darkgray"` sets their color to gray, more appropriate choices when plotting more than 30,000 points. As a consequence of jittering and using smaller and lighter plotting symbols, we see clearly that the density of points for education = 12, high school graduates, is higher than for other years, and that the data for education < 8 are very sparse. 
- We use the *smooth argument* to set the span for the default loess smoother to 1/3, half the default value of 2/3. Setting a smaller span uses less data to estimate the fitted curve at each value of education, allowing us to resolve greater detail in the regression function. Here we observe a dip in the regression function when education ≈ 11, so individuals who just missed graduating from high school perform somewhat worse than expected by a straight-line fit on the vocabulary test. Similarly, there is a small positive bulge in the regression function corresponding to college graduation, approximately 16 years of education. The specifications `col.smooth="black"` and `col.var="black"` set the color of the loess and variability lines, making them more visible in the graph; the default is the color of the points, now gray. 
- Finally, as before, the *main argument* sets the title of the graph.

::::
:::::

###### Jitter (ggplot2)

:::::{.my-r-code}
:::{.my-r-code-header}
:::::: {#cnj-chap03-scatterplot-jitter-ggplot2}
: Scatterplots of `vocabulary` by `education` using {**ggplot2**
::::::
:::
::::{.my-r-code-container}
::: {#lst-chap03-scatterplot-jitter-ggplot2}
```{r}
#| label: scatterplot-jitter-ggplot2
#| cache: true

Vocab <- carData::Vocab

Vocab |> 
    ggplot2::ggplot(
        ggplot2::aes(
            x = education,
            y = vocabulary
        )
    ) +
    ggplot2::geom_jitter(
        color = "grey",
        size = .01,
        alpha = 0.2
    ) +
    ggplot2::stat_smooth(
        formula = 'y ~ x',
        method = "loess",
        span = .7,
        level = 0.95
    ) +
    ggplot2::stat_smooth(
        formula = 'y ~ x',
        method = "lm"
    ) 
```

Scatterplots of vocabulary by education with {**car**}
:::

***
 
At first, using `span = 0.33` as in @lst-chap03-scatterplot-jitter-car I've got many warnings that I do not understand. (The same happened with the {**car**} approach, but these warnings did not appear in the outcome.) It seems that these warnings come when there are too many identical values, e.g. when the variable behaves more like a discrete variable than a continuous variable as it would be required for a proper smoothing (see: ).

Despite the advice in a StackOverflow posting [R warnings, simpleLoess, pseudoinverse etc.](https://forum.posit.co/t/r-warnings-simpleloess-pseudoinverse-etc/8651) I changed span to 0.7 and the warnings disappeared. 

> Increasing the `span` parameter has the effect of "squashing out", … the piles of repeated values where they occur. … I would definitely not increase span to achieve the squashing: it is a lot better to use a tiny amount of jitter for that purpose; span should be dictated by other considerations….

Changing the width and/or height of jitter didn't remove the warnings, but I got the advise in the warnings that I should enlarge `span`.

::::
:::::


:::

::::
:::::

The {**car**} packages has the ability to display "extreme" points, mostly defined by the `r glossary("Mahalanobis")` distance. Unlike simple Euclidean distances, which are inappropriate when the variables are scaled in different units, Mahalanobis distances take into account the variation and correlation of x and y.

> The Mahalanobis distance (MD) is the distance between two points in multivariate space. In a regular Euclidean space, variables (e.g. x, y, z) are represented by axes drawn at right angles to each other; The distance between any two points can be measured with a ruler. For uncorrelated variables, the Euclidean distance equals the MD. However, if two or more variables are correlated, the axes are no longer at right angles, and the measurements become impossible with a ruler. In addition, if you have more than three variables, you can’t plot them in regular 3D space at all. The MD solves this measurement problem, as it measures distances between points, even correlated points for multiple variables. [@glenn2017]

#### Scatterplots enhanced with marginal plots {#sec-chap03-2-1-1}

To reproduce the marginal boxplots without {**car**} and using the {**tidyverse**} approach I learned about the {**ggExtra**} package.
It now only allows scatterplots with marginal boxplots but also with histograms, density, densigram (a histogram overlaid with a density distribution) and violins.

:::::{.my-experiment}
:::{.my-experiment-header}
:::::: {#def-chap03-marginal-plots}
: Using {ggExtra} for adding marginal graphs to scatterplots
::::::
:::
::::{.my-experiment-container}

::: {.panel-tabset}

###### boxplot

:::::{.my-r-code}
:::{.my-r-code-header}
:::::: {#cnj-chap03-marginal-boxplot}
: Scatterplot of `prestige` versus `income` with marginal boxplots
::::::
:::
::::{.my-r-code-container}
::: {#lst-chap03-marginal-boxplot}
```{r}
#| label: marginal-boxplot

gg_Prestige <- base::readRDS("data/chap03/gg_Prestige.rds")


ggExtra::ggMarginal(gg_Prestige, type = "boxplot")
```
  
Scatterplot of `prestige` versus `income` with marginal boxplots
:::


::::
:::::

###### histogram

:::::{.my-r-code}
:::{.my-r-code-header}
:::::: {#cnj-chap03-marginal-histogram}
: Scatterplot of `prestige` versus `income` with marginal histograms
::::::
:::
::::{.my-r-code-container}
::: {#lst-chap03-marginal-histogram}
```{r}
#| label: marginal-histogram

gg_Prestige <- base::readRDS("data/chap03/gg_Prestige.rds")


ggExtra::ggMarginal(gg_Prestige, type = "histogram")
```
  
Scatterplot of `prestige` versus `income` with marginal histograms
:::


::::
:::::

###### density

:::::{.my-r-code}
:::{.my-r-code-header}
:::::: {#cnj-chap03-marginal-density}
: Scatterplot of `prestige` versus `income` with marginal density distribution
::::::
:::
::::{.my-r-code-container}
::: {#lst-chap03-marginal-density}
```{r}
#| label: marginal-density

gg_Prestige <- base::readRDS("data/chap03/gg_Prestige.rds")


ggExtra::ggMarginal(gg_Prestige, type = "density")
```
  
Scatterplot of `prestige` versus `income` with marginal density distribution
:::


::::
:::::

###### densigram

:::::{.my-r-code}
:::{.my-r-code-header}
:::::: {#cnj-chap03-marginal-densigram}
: Scatterplot of `prestige` versus `income` with marginal densigrams
::::::
:::
::::{.my-r-code-container}
::: {#lst-chap03-marginal-densigram}
```{r}
#| label: marginal-densigram

gg_Prestige <- base::readRDS("data/chap03/gg_Prestige.rds")


ggExtra::ggMarginal(gg_Prestige, type = "densigram")
```
  
Scatterplot of `prestige` versus `income` with marginal densigrams
:::

I reported that `..density..` was deprecated since {**ggplot2**} 3.4.0. (See [my comment](https://github.com/daattali/ggExtra/issues/175#issuecomment-2149784674))


::::
:::::

###### violin

:::::{.my-r-code}
:::{.my-r-code-header}
:::::: {#cnj-chap03-marginal-violin}
: Scatterplot of `prestige` versus `income` with marginal violin plots
::::::
:::
::::{.my-r-code-container}
::: {#lst-chap03-marginal-violin}
```{r}
#| label: marginal-violin

gg_Prestige <- base::readRDS("data/chap03/gg_Prestige.rds")


ggExtra::ggMarginal(gg_Prestige, type = "violin")
```
  
Scatterplot of `prestige` versus `income` with marginal violin plots
:::


::::
:::::



:::

::::
:::::

#### Scatterplot enhanced with point labelling {#sec-chap03-2-1-2}

The	`car::scatterplot()` function can	draw	scatterplots with	a	wide variety of enhancements and options. One enhancement is to include marginal plots as I have done with the use of {**ggExtra**} in the previous section @sec-chap03-2-1-1.

Another enhancement is to label "extreme" points with different measures and in different conditions. I have already developed this enhancement for 

- Outliers in boxplots using 1.5 times the `r glossary("IQR")`  subtracting to the lower quartile (.25) and adding to the upper quartile (.75) (see: @lst-chap03-compute-boxplot-outliers) and for
- `r glossary("Mahalanobis")` distances using p-values < .001 (see: @lst-chap04-davis-ggrepel).

Here now I will not use the `r glossary("p-value")` for the Mahalanobis distances but a fix amount of points to label. To go conform with the example I will at first apply a fixed `n = 4` but in the next challenge I will try to develop a general function that distinguishes between type (scatterplot, boxplot) and method (p-value, fixed number of points).

:::::{.my-experiment}
:::{.my-experiment-header}
:::::: {#def-chap03-replicating-books-figure-3-8}
: Label 4 points of farest Mahalanobis distances in scatterplots 
::::::
:::
::::{.my-experiment-container}

::: {.panel-tabset}



###### labelled-4 

:::::{.my-r-code}
:::{.my-r-code-header}
:::::: {#cnj-chap03-prestige-labelled-n-points}
: Scatterplot	of	`prestige`	versus	`income`	with four	points	identified
::::::
:::
::::{.my-r-code-container}
::: {#lst-chap03-prestige-labelled-n-points}
```{r}
#| label: prestige-labelled-4-points

set.seed(420)
df_temp <- carData::Prestige |> 
  dplyr::select(prestige, income) |> 
  tidyr::drop_na() 

df <- df_temp |>
  dplyr::mutate(mahal =
        {
          stats::mahalanobis(
                df_temp,
                base::colMeans(df_temp),
                stats::cov(df_temp)
                )
        }
  ) |>
  dplyr::mutate(p =
    pchisq(
      q = mahal,
      df = 1,
      lower.tail = FALSE
      )
    ) |>
  tibble::rownames_to_column(var = "ID") |>
  dplyr::mutate(label = "") |> 
  dplyr::arrange(desc(mahal)) |> 
  dplyr::mutate(label =
      dplyr::case_when(dplyr::row_number() <= 4 ~ ID)
)
  
gg_enhanced <- df |>
  ggplot2::ggplot(
        ggplot2::aes(
                x = income,
                y = prestige
        )
  ) +
  ggplot2::geom_point(
    shape = 1,
    size = 2
  ) +
  ggrepel::geom_text_repel(
    ggplot2::aes(label = label),
    na.rm = TRUE
  ) +
  ggplot2::stat_smooth(
    ggplot2::aes(linetype = "solid"),
    formula = y ~ x,
    method = lm,
    se = FALSE,
    color = "black"
  ) +
    ggplot2::geom_smooth(
      ggplot2::aes(
          linetype = "dashed"
        ),
      formula = y ~ x,
      method = loess,
      se = TRUE,
      color = "black"
  ) +
  ggplot2::scale_linetype_discrete(
    name = "",
    label = c(
      "Least Squares",
      "Unbiased Reporting   "
    )
  ) +
  ggplot2::guides(
    linetype = ggplot2::guide_legend(position = "inside")
  ) +
  ggplot2::theme(
    legend.title = ggplot2::element_blank(),
    legend.position.inside = c(.8, .15),
    legend.box.background = ggplot2::element_rect(
                              color = "black",
                              linewidth = 1
                              )
                 )

ggExtra::ggMarginal(p = gg_enhanced, type = "boxplot")
```

Replication of book’s Figure 3.8 with {**ggplot2**}, {**ggrepel**} and {**ggExtra**}.
:::
***

This is my attempt to replicate book’s Figure 3.8. Instead of using functions of tthe he {**car**} package I used packages from the {**tidyverse**} approach: {**gplot2**} with {**ggrepel**} and {**ggExtra**} (see @pak-ggrepel and @pak-ggExtra). 



::::
:::::

Using the `r glossary("p-value")` of .001 we would get 5 (instead of 4) points labelled.

###### function

**Still to do!**

:::::{.my-r-code}
:::{.my-r-code-header}
:::::: {#cnj-ID-text}
: Numbered R Code Title
::::::
:::
::::{.my-r-code-container}
::: {#lst-chaplisting-ID}
```{r}
#| label: function-missing
#| eval: false

my_find_boxplot_outlier <- function(df) {
  return(df < quantile(df, .25) - 1.5*IQR(df) | df > quantile(df, .75) + 1.5*IQR(df))
  
my_find_boxplot_outlier(carData::Prestige$income)

prestige1 <- Prestige |> 
  dplyr::mutate(outlier = my_find_boxplot_outlier(income)) |> 
  tibble::rownames_to_column(var = "ID") |>
  dplyr::mutate(outlier =
        dplyr::case_when(outlier == TRUE ~ ID,
                         outlier == FALSE ~ "")
  )

}

# label_points <-  function(
#     df, 
#     n = 2, 
#     type = c("scatterplot", "boxplot"),
#     output = c("dataframe", "labels", "both")
#     ) {
#     if (type == "scatterplot") {
#           df |>
#           dplyr::mutate(mahal =
#                 {
#                   stats::mahalanobis(
#                         df,
#                         base::colMeans(df),
#                         stats::cov(df)
#                         )
#                 }
#           )
#       } |>
#     if (n == 0) ( {
#         dplyr::mutate(p =
#           pchisq(
#             q = mahal,
#             df = 1,
#             lower.tail = FALSE
#             )
#           ) 
#       } ) |>
#       tibble::rownames_to_column(var = "ID") |>
#       dplyr::mutate(label = "") |> 
#       dplyr::arrange(desc(mahal)) |> 
#       if (n != 0) (
#         dplyr::mutate(
#           label = dplyr::case_when(dplyr::row_number() <= 4 ~ ID)
#           )
#       ) 
# }


set.seed(420)
df_temp <- carData::Prestige |> 
  dplyr::select(prestige, income) |> 
  tidyr::drop_na()

test <-  label_points(
  df = df_temp,
  n = 4,
  type = "scatterplot",
  output = "dataframe")
```

Listing title
:::

::::
:::::



:::

::::
:::::



### Parallel boxplots {#sec-chap03-2-2}

:::::{.my-example}
:::{.my-example-header}
:::::: {#exm-chap03-parallel-boxplots}
: Parallel boxplots
::::::
:::
::::{.my-example-container}

::: {.panel-tabset}

###### graphics::boxplot()

:::::{.my-r-code}
:::{.my-r-code-header}
:::::: {#cnj-chap03-parallel-boxplot-base-r}
: Parallel boxplots with `graphics::boxplot()`
::::::
:::
::::{.my-r-code-container}
::: {#lst-chap03-parallel-boxplot-base-r}
```{r}
#| label: parallel-boxplot-base-r

Vocab <- carData::Vocab

graphics::boxplot(vocabulary ~ education, data = Vocab)
```
Boxplots with `graphics::boxplot()` of `vocabulary` separately for each value of years of `education.`
:::

::::
:::::


###### car::Boxplot()

:::::{.my-r-code}
:::{.my-r-code-header}
:::::: {#cnj-chap03-parallel-boxplot-car}
: Parallel boxplots with `car::Boxplot()`
::::::
:::
::::{.my-r-code-container}
::: {#lst-chap03-parallel-boxplot-car}    
```{r}
#| label: parallel-boxplot-car

car::Boxplot(vocabulary ~ education, data = Vocab, id = FALSE)
```
Boxplots with `car::Boxplot()` of `vocabulary` separately for each value of years of `education.`
:::

***
This command draws a separate boxplot for vocabulary for all cases with the same value of education, so we condition on the value of education. The formula has the response variable on the left of the `~` and a discrete conditioning variable — typically, but not necessarily, a factor — on the right.

Setting `id = FALSE` prevents labelling up to ten outliers, which would be very distracting because of the many boxplots. Therefore we get using `car::Boxplot()` an identical graph as in @lst-chap03-parallel-boxplot-base-r with the base R version of `graphics::boxplot()`.

With `xtabs(~ education, data = Vocab)` we get the distribution of education:

```{r}
xtabs(~ education, data = Vocab)
```


::::
:::::

In this example, the conditioning predictor education is a discrete numeric variable. For the subsamples with larger sample sizes, 8 ≤ education ≤ 18, rather than a steady increase in vocabulary with education, there appear to be jumps every 2 or 3 years, at years 10, 12, 15, and 18.

###### Ornstein

:::::{.my-r-code}
:::{.my-r-code-header}
:::::: {#cnj-chap03-canadian-directorats}
: Interlocking directorates	among	248	major	Canadian	corporations
::::::
:::
::::{.my-r-code-container}
::: {#lst-chap03-canadian-directorats}
```{r}
#| label: canadian-directorats

car::brief(carData::Ornstein)
```

Brief summary of the Ornstein data set about "Interlocking Directorates Among Major Canadian Firms"
:::

***

The variables in the data set include the assets of each corporation, in millions of dollars; the corporation’s sector of operation, a factor with 10 levels; the factor nation, indicating the country in which the firm is controlled, with levels "CAN" (Canada), "OTH" (other), "UK", and "US"; and interlocks, the number of interlocking directorate and executive positions maintained between each company and others in the data set.

::::
:::::

###### boxplot

:::::{.my-r-code}
:::{.my-r-code-header}
:::::: {#cnj-chap03-ornstein-boxplot}
: Parallel boxplots of interlocks by nation of control
::::::
:::
::::{.my-r-code-container}
::: {#lst-chap03-ornstein-boxplot}
```{r}
#| label: ornstein-boxplot
#| results: hold


Ornstein <- carData::Ornstein

car::Boxplot(
  interlocks ~ nation, 
  data = Ornstein, 
  main = "(a)")



```

Parallel boxplots of interlocks by nation of control, for Ornstein’s interlocking-directorate data.
:::

***

Because the names of the companies are not given in the original data source, the points are labeled by case numbers. The firms are in descending order by assets, and thus the identified points are among the largest companies.

::::
:::::

###### errorbars

:::::{.my-r-code}
:::{.my-r-code-header}
:::::: {#cnj-chap03-plotrix-errorbars}
: A mean/standard deviation plot
::::::
:::
::::{.my-r-code-container}
::: {#lst-chap03-plotrix-errorbars}
```{r}
#| label: plotrix-errorbars

means <- car::Tapply(interlocks ~ nation, mean, data = Ornstein)
sds <- car::Tapply(interlocks ~ nation, sd, data = Ornstein)

plotrix::plotCI(
  1:4,
  means,
  sds,
  xaxt = "n",
  xlab = "Nation of Control",
  ylab = "interlocks",
  main = "(b)",
  ylim = range(Ornstein$interlocks)
)
lines(1:4, means)
axis(1, at = 1:4, labels = names(means))
```

A mean/standard deviation plot for Ornstein’s interlocking-directorate data
:::

***

- The `car::Tapply()` function adds a formula interface to the standard `base::tapply()` function. 
- The first call to `car::Tapply()` computes the mean of interlocks for each level of nation, and the second computes within-nation standard deviations. 
- The basic graph is drawn by `plotrix::plotCI()`. The first argument to this function specifies the coordinates on the horizontal axis, the second the coordinates on the vertical axis, and the third the vector of SDs. 
- The standard graphical argument `xaxt="n"` suppresses the x-axis tick marks and labels, and the ylim argument is used here to match the vertical axis of panel (b) to that of panel (a). 
- The `graphics::lines()` function joins the means, and the `graphics::axis()` function labels the horizontal axis with the names of the groups. 
- The first argument to `graphics::axis()` specifies the side of the graph where the axis is to be drawn: `side = 1` (as in the example) is below the graph, 2 at the left, 3 above, and 4 at the right.

::::
:::::

We discourage the use of bar charts for means because interpretation of the length of the bars, and therefore the visual metaphor of the graph, depends on whether or not a meaningful origin (zero) exists for the measured variable and whether or not the origin is included in the graph. 

The error bars can also lead to misinterpretation, because neither standard-error bars nor standard-deviation bars are the appropriate measure of variation for comparing means between groups: They make no allowance or correction for multiple testing, among other potential problems.

So is the mean/SD graph misleading because rather than showing the outliers, the graph inflates both the means and the SDs, particularly for Canada, and disguises the skewness that is obvious in the boxplots.



:::

::::
:::::

### More on the graphics::plot() command {#sec-chap03-2-3}

- There are different option to call the plot() function:
  - **Basic scatterplot**: `plot(y ~ x)` or `plot(x, y)` with x on the horizontal and y on the vertical axis. 
  - **Including reference to data**: `plot(y ~ x, data = D)`,	`plot(D$x, D$y)`,	or	`with(D, plot(x, y)`.
  - **Index plot**: `plot(x)` produces a scatterplot with x on the *vertical* axis and case numbers of the horizontal axis.
  - **Boxplot**: `plot(y ~ factor)` is the same as `boxplot(y ~ factor)` using the standard `graphics::boxplot()` function and not  `car::Boxplot()`.
  -**Model plots**: Depending on the type of object, for instance with a `lm`-object `plot(model-object)` draws several graphs that are commonly associated with linear regression models fit by least squares. In contrast, `plot(density(x))` draws a plot of the density estimate for the numeric variable x.
  
The `plot()` function can take many additional optional arguments that control the appearance of the graph, the labeling of its axes, the fonts used, and so on. We can set some of these options globally with the `par()` function or just for the current graph via arguments to the `plot()` function. But the parameter names are not very intuitive (e.g., `cex` (magnifying plotting text and symbols, or `crt` rotating of single characters etc.) 

Plots can be built up sequentially by first creating a basic graph and then adding to it. 

::: {.callout-note #nte-chap03-plot-vs.ggplot2}
Generally I prefer the {**ggplot2**} over the base R `plot()` function. But it is important to be able to understand graphs developed with `graphics::plot()` as they appear quite often in examples and textbooks.
:::

## Examining multivariate data {#sec-chap03-3}

Because we can only perceive objects in three spatial dimensions, examining multivariate data is intrinsically more difficult than examining univariate or bivariate data.

### Three-dimensional plots {#sec-chap03-3-1}

Perspective, motion, and illumination can convey a sense of depth, enabling us to examine data in three dimensions on a two-dimensional computer display. The most effective software of this kind allows the user to rotate the display, to mark points, and to plot surfaces such as regression mean functions. The {**rgl**} package (@rgl) links R to the [OpenGL](https://www.opengl.org/) threedimensional graphics library often used in animated films. The `car::scatter3d()` function {**rgl**} to provide a three-dimensional generalization of the `car::scatterplot()` function.

:::::{.my-r-code}
:::{.my-r-code-header}
:::::: {#cnj-chap03-3d-scatter}
: Three-dimensional	scatterplot
::::::
:::
::::{.my-r-code-container}
::: {#lst-chap03-3d-scatter}
```{r}
#| label: 3d-scatter
#| results: hold

Duncan <- carData::Duncan

car::scatter3d(prestige	~ income + education,	
               data = Duncan,	
               id = list(n = 3)
               )
rgl::rglwidget()
```

Three-dimensional scatterplot for Duncan’s occupational-prestige data, showing the least-squares regression plane. Three unusual points were labeled automatically.
:::

::: {.callout-warning #wrn-chap03-3d-xquartz}
The graph with the proposed code of the book is in my macOS installation shown with an extra XQuartz window only inside RStudio.
:::

:::::{.my-solution}
:::{.my-solution-header}
:::::: {#sol-chap03-3d-rgl}
: How to embed 3d graphics into Quarto
::::::
:::
::::{.my-solution-container}
I found a solution from a [StackOverflow post](https://stackoverflow.com/a/63597059/7322615) to embed the 3D-graph into the Quarto document. The post refers to a R Markdown document. Reading the comments I learned that there may be some glitches with Quarto. But luckily I did not experience these mentioned problems.

My solution:

- `options(rgl.useNULL = TRUE)` to prevent opening an extra window with Xquartz.
- using `rgl::rglwidget()` after the code for the 3D-plot.

(I did not try the option with calling `rgl::setupKnitr(autoprint = TRUE)` as the first line.)
::::
:::::


The graph shows the least-squares regression plane for the regression of the variable on the vertical or y-axis, `prestige`, on the two variables on the horizontal (or x- and z-) axes, `income` and `education`; three cases (minister, conductor, and railroad engineer) are identified as the most unusual based on their `r glossary("Mahalanobis")` distances from the centroid (i.e., the point of means) of the three variables. The three-dimensional scatterplot can be rotated by left-clicking and dragging with the mouse. Color is used by default, with perspective, sophisticated lighting, translucency, and fog-based depth cueing.

::::
:::::

The `car::scatter3d()` function 

- can also plot other regression surfaces (e.g., nonparametric regressions), 
- can identify points interactively and according to other criteria, 
- can plot concentration ellipsoids, and 
- can rotate the plot automatically. 

Because three-dimensional dynamic graphs depend on color, perspective, motion, and so on for their effectiveness, it is important to read the help file for `car::scatter3d()` and to experiment with the examples therein. 

:::::{.my-resource}
:::{.my-resource-header}
:::::: {#lem-chap03-3d-plots}
: 3-dimensional graphic packages in R
::::::
:::
::::{.my-resource-container}
Besides several functions for 3D graphics in {**car**} there are also facilities in R for drawing static three-dimensional graphs:

- `graphics::pers()` draws perspective plots of a surface over the x–y plane.
- {**lattice**} has generic functions to draw static 3d scatter plot and surfaces `cloud()` and `wireframe()`

{**rggobi**} is also mentionded and should link R to the [GGobi](http://ggobi.org/index.html) system for visualizing data in three and more dimensions. See also the [GGobi book](http://ggobi.org/book.html). But I believe that {**rggobi**} and the GGobi system is outdated, because {**rggobi**} is archived and not at CRAN anymore and links from the [GGobi page](http://ggobi.org/index.html) are rotten. But maybe [GGally](https://ggobi.github.io/ggally/) is the successor, because it is on the same GitHub account and the URL stats with "ggobi".
::::
:::::

### Scatterplot matrices {#sec-chap03-3-2}

`r glossary("ScatterplotMatrix", "Scatterplot matrices")` are graphical analogs of correlation matrices, displaying bivariate scatterplots of all pairs of numeric variables in a data set as a twodimensional graphical array. Because the panels of a scatterplot matrix are just two-dimensional scatterplots, each panel is the appropriate summary graph for the regression of the y-axis variable on the x-axis variable.

:::::{.my-example}
:::{.my-example-header}
:::::: {#exm-chap03-scatterplot-matrices}
: Scatterplot Matrices
::::::
:::
::::{.my-example-container}

::: {.panel-tabset}

###### graphics::pairs

:::::{.my-r-code}
:::{.my-r-code-header}
:::::: {#cnj-chap03-scatterplot-matrices-pairs}
: Scatterplot matrices with `graphics::pairs()`
::::::
:::
::::{.my-r-code-container}
::: {#lst-chap03-scatterplot-matrices-pairs}
```{r}
#| label: scatterplot-matrices-pairs

graphics::pairs(
  formula = ~ prestige + income + education + women,
  data = carData::Prestige
)
```
Scatterplot matrix for the Canadian occupational-prestige data plotted with `graphics::pairs()`
:::

::::
:::::


###### lattice::splom()

:::::{.my-r-code}
:::{.my-r-code-header}
:::::: {#cnj-chap03-scatterplot-matrices-lattice-splom}
: Scatterplot matrices with `lattice::splom()`
::::::
:::
::::{.my-r-code-container}
::: {#lst-chap03-scatterplot-matrices-lattice-splom}  
```{r}
#| label: scatterplot-matrices-lattice-splom
#| fig-width: 7
#| fig-height: 7

lattice::splom(
  x = ~ Prestige[1:4],
  data = carData::Prestige
)
```

Scatterplot matrix for the Canadian occupational-prestige data plotted with `lattice::splom()`
:::

::::
:::::

###### car::scatterplotMatrix()

:::::{.my-r-code}
:::{.my-r-code-header}
:::::: {#cnj-chap03-scatterplot-matrices-car-scatterplotMatrix}
: Scatterplot matrices with `car::scatterplotMatrix()`
::::::
:::
::::{.my-r-code-container}
::: {#lst-chap03-scatterplot-matrices-car-scatterplotMatrix}
```{r}
#| label: scatterplot-matrices-car-scatterplotMatrix

car::scatterplotMatrix( ~ prestige + income + education + women,
                   data = carData::Prestige)
```

Scatterplot matrix for the Canadian occupational-prestige data plotted with `car::scatterplotMatrix()` with density estimates on the diagonal
:::

As in `car::scatterplot()` (@lst-chap03-scatterplot-car), mean and variability smoothers and a least-squares regression line are added by default to each panel of a scatterplot matrix and are controlled respectively by optional `smooth` and `regLine` arguments. (The `car::scatterplot()` and `car::scatterplotMatrix()` functions have many other arguments in common.)

The `id` argument for marking points is also the same for both functions, except that interactive point-marking isn’t supported for scatterplot matrices. 

The diagonal argument to car::scatterplotMatrix() controls the contents of the diagonal panels. ()


::::
:::::

###### GGally::ggpairs()

:::::{.my-r-code}
:::{.my-r-code-header}
:::::: {#cnj-chap03-scatterplot-matrices-GGally-ggpairs}
: Scatterplot matrices with `GGally::ggpairs()
::::::
:::
::::{.my-r-code-container}
::: {#lst-chap03-scatterplot-matrices-GGally-ggpairs}
```{r}
#| label: scatterplot-matrices-GGally-ggpairs

carData::Prestige |> 
  dplyr::select(4, 2, 1, 3) |> 
  GGally::ggpairs()
```

Scatterplot matrix for the Canadian occupational-prestige data plotted with `GGally::ggpairs()` with density estimates on the diagonal and correlation coefficients on the upper half of the panel.
:::

In contrast to the `car::scatterplotMatrix()` function `GGally::ggpairs()` also shows the correlation coefficients instead of reverting x/y axis for the other half of the scatterplots.

::::
:::::

:::::{.my-resource}
:::{.my-resource-header}
:::::: {#lem-chap03-GGally}
: Scatterplot matrices with {**GGally**}
::::::
:::
::::{.my-resource-container}

- For theroetical background see "The Generalized Pairs Plot" [@emerson2013].
- For practical usage: [Extra - Scatterplot Matrices With GGally](https://unc-libraries-data.github.io/R-Open-Labs/Extras/ggally/ggally.html) 
- For the general function `GGally::ggmatrix()` for managing multiple plots in a matrix-like layout: [ggmatrix(): Plot matrix](https://ggobi.github.io/ggally/articles/ggmatrix.html)


::::
:::::



:::

::::
:::::

## Transforming data {#sec-chap03-4}

### Logarithms {#sec-chap03-4-1}

#### Introduction

Logarithmic scale corresponds to viewing variation through relative or percentage changes, rather than through absolute changes. In many instances humans perceive relative changes rather than absolute changes, and so logarithmic scales are often theoretically appropriate as well as useful in practice [@varshney2013].

- The log of a positive number x to the base b (where b is also a positive number), written $log_{b}x$, is the exponent to which b must be raised to produce x. That is, if $y = log_{b}x$, then $b_{y} = x$. 
- Thus, for example, 
  - $log_{10} 100 = 2$ because $10^2 = 100$; 
  - $log_{10} 0.01 = −2$ because $10^{−2} = 1/10^2 = 0.01$; 
  - $log_{2} 8 = 3$ because $2^3 = 8$; and 
  - $log_{2} 1/8 = −3$ because $2^{−3} = 1/8$. 
- So-called **natural logs** use the base $e \approx 2.71828$. 
- Thus, for example, $log_{e} e = 1$ because $e^1 = e$. 
- Logs to the bases 2 and 10 are often used in data analysis, because powers of 2 and 10 are familiar. Logs to the base 10 are called **common logs**.
- Regardless of the base b, $log_{b} 1 = 0$, because $b^0 = 1$. 
- Regardless of the base, the log of zero is undefined (or taken as log 0 = −∞), and the logs of negative numbers are undefined.

Summarizing we can say that there are different types of logarithms (abbreviated as "logs"):

- **Natural logs** have as base $e \approx 2.718$.
- **Common logs** have as base 10.
- **Others**, for instance with base 2.

The	R	functions	`base::log()`,	`base::log10()`,	and	`base::log2()`	compute	the	natural,	base-10,	and base-2	logarithms,	respectively. 

Some examples with 7

- Natural log = `log(7)` = `r log(7)`.
- Common log = `log10(7)` = `r log10(7)`.
- Base 2-log = `log2(7)` = `r log2(7)`.
- General form = `log(x, base = exp(1))` or `logb(x, base = exp(1))` for instance: 
  - `log10(7) = log(x = 7, base = 10) = log(7,10)` = `r log(x = 7, base = 10)`
- The **exponential function** `base::exp()` computes powers of $e$, and thus `exp(1)` = $e$ = `r exp(1)`. 
- Exponentiating is the inverse of taking logarithms, and so the equalities `x = exp(log(x))` = `log(exp(x))` hold for any positive number x and base b, for instance 
  - `exp(log(50))` = `r exp(log(50))`,
  - `log(exp(50))` = `r log(exp(50))`,
  - `10^(log10(50))` = `r 10^(log10(50))`
  - `log10(10^50)` = `r log10(10^50)`

The various log functions in R can be applied to numbers, numeric vectors, numeric matrices, and numeric data frames. They return 

- NA for missing values, 
- NaN (Not a Number) for negative values, and 
- -Inf (-$\infty$) for zeros.

::: {.callout-important #imp-chap03-log-rule}
##### A rule for using logarithms

**Although it isn’t necessarily true that linearizing transformations also stabilize variation and make marginal distributions more symmetric, in practice, this is often the case.**

For any strictly positive variable with no fixed upper bound whose values cover two or more orders of magnitude (that is, powers of 10), replacing the variable by its log is likely to be helpful. 

Conversely, if the range of a variable is considerably less than an order of magnitude, then transformation by logarithms, or indeed any simple transformation, is unlikely to make much of a difference. 

Variables that are intermediate in range may or may not benefit from a log or other transformation.

For many important aspects of statistical analyses, such as fitted values, tests, and model selection, **the base of logarithms is inconsequential**, in that identical conclusions will be reached regardless of the base selected. 

The interpretation of regression coefficients, however, may be simpler or more complicated depending on the base. For example, increasing the log base-2 of x by 1 implies doubling x, but increasing the natural log of x by 1 implies multiplying x by e.
:::

#### Practice

:::::{.my-example}
:::{.my-example-header}
:::::: {#exm-chap03-log-transformations}
: Log transformations
::::::
:::
::::{.my-example-container}

::: {.panel-tabset}

###### car (Ornstein)

:::::{.my-r-code}
:::{.my-r-code-header}
:::::: {#cnj-chap03-car-assets-ornstein}
: Distribution of assets in the Ornstein data set
::::::
:::
::::{.my-r-code-container}
::: {#lst-chap03-car-assets-ornstein}
```{r}
#| label: car-assets-ornstein
#| results: hold
#| fig-height: 5

graphics::par(mfrow = c(1, 2), mar = c(5, 4, 6, 2) + 0.1)

car::densityPlot( ~ assets,
             data = carData::Ornstein,
             xlab = "assets",
             main = "(a)"
             )
car::densityPlot(
          ~ base::log10(assets),
          data = carData::Ornstein,
          adjust = 0.65,
          xlab = base::expression(log[10] ~ "(assets)"),
          main = "(b)"
          )
car::basicPowerAxis(
      0,
      base = 10,
      side = "above",
      at = 10 ^ (2:5),
      axis.title = ""
      )
```
Distribution of assets (millions of dollars) in the Ornstein data set (a) before and (b) after log transformation.
:::

***

- The command `graphics::par(mfrow=c (1, 2))` sets the global graphics parameter `mfrow`, which divides the graphics window into an array of subplots with one row and two columns. 
- We also set the `mar` (margins) graphics parameter to increase space at the top of the plots and 
- call the `car::basicPowerAxis()` function to draw an additional horizontal axis at the top of panel (b) showing the original units of assets. 
- The `xlab` argument in the second command is set equal to an expression, allowing us to typeset 10 as the subscript of log.

::::
:::::

###### ggplot2 (Ornstein)

:::::{.my-r-code}
:::{.my-r-code-header}
:::::: {#cnj-chap03-ggplot2-assets-ornstein}
: Distribution of assets in the Ornstein data set
::::::
:::
::::{.my-r-code-container}
::: {#lst-chap03-ggplot2-assets-ornstein}
```{r}
#| label: ggplot2-assets-ornstein
#| fig-height: 3

gg1 <- carData::Ornstein |> 
  ggplot2::ggplot(
    ggplot2::aes(x = assets)
  ) +
  ggplot2::geom_density() +
  ggplot2::ylab("Density") +
  ggplot2::geom_rug() +
  ggplot2::ggtitle("(a)") +
  ggplot2::theme_bw()

gg2 <- carData::Ornstein |> 
  ggplot2::ggplot(
    ggplot2::aes(x = base::log10(assets))
  ) +
  ggplot2::geom_density(
    adjust = 0.65
  ) +
  ggplot2::labs(
    y = "Density",
    x = base::expression(log[10](assets))
  ) +
  ggplot2::geom_rug() +
  ggplot2::scale_x_continuous(
    sec.axis = ggplot2::sec_axis(
      transform = ~ 10^.,
      breaks = c(100, 1000, 1e+5),
      labels = \(x) ifelse(x <= 1000, x, scales::scientific(x))
      )
    ) +
  ggplot2::theme_bw()


gridExtra::grid.arrange(gg1, gg2, ncol = 2)

```

Distribution of assets (millions of dollars) in the Ornstein data set (a) before and (b) after log transformation.
:::

***
 @lst-chap03-ggplot2-assets-ornstein is my try to replicate Figure 3.15 with {**ggplot2**}. At first I didn't succeed because I used as formula for the transformation `.^10` instead of `10^.`. I have to confess that this was not a typo as the [answer in SO](https://stackoverflow.com/a/78620195/7322615) supposed, but I didn't understand what the formula effectively meant.


::::
:::::


###### car (UN)

:::::{.my-r-code}
:::{.my-r-code-header}
:::::: {#cnj-chap03-car-infant-mortality-un}
: Infant mortality rate and gross domestic product per capita, from the UN data set, produced with {car}
::::::
:::
::::{.my-r-code-container}
::: {#lst-chap03-car-infant-mortality-un}  
```{r}
#| label: car-infant-mortality-un

car::scatterplot(
  infantMortality ~ ppgdp,
  data = carData::UN,
  xlab = "GDP per Capita",
  ylab = "Infant Mortality Rate (per 1000 births)",
  main = "(a)"
)

car::scatterplot(
  infantMortality ~ ppgdp,
  data = carData::UN,
  xlab = "GDP per capita",
  ylab = "Infant Mortality Rate (per 1000 births)",
  main = "(b)",
  log = "xy",
  id = list(n = 3)
)


```

Infant mortality rate and gross domestic product per capita, from the UN data set: (a) untransformed data and (b) both variables log-transformed, using the {car} package.
:::

***
The UN data set in the {**carData**} package, for example, contains data obtained from the United Nations on several characteristics of 213 countries and other geographic areas recognized by the UN around 2010; the data set includes the infant mortality rate (`infantMortality`, infant deaths per 1,000 live births) and the per-capita gross domestic product (`ppgdp`, in U.S. dollars) of most of the countries. 

There are only 193 countries in the UN data set with both variables observed. The R commands used here automatically reduce the data to the 193 fully observed cases on these two variables.

***

- The argument `log="xy"`, which also works with `graphics::plot()` and `car::scatterplotMatrix()`, uses logs on both axes but labels the axes in the units of the original variables.
- The optional argument `id=list(n=3)` causes the three most “unusual” points according to the Mahalanobis distances of the data points to the means of the two variables to be labeled with their row names. These are three African countries that have relatively large infant mortality rates for their levels of percapita GDP. The names of the identified points and their row numbers in the data set are printed in the R console.

::::
:::::

###### ggplot2 (UN)

:::::{.my-r-code}
:::{.my-r-code-header}
:::::: {#cnj-chap03-ggplot2-infant-mortality-un}
: Infant mortality rate and gross domestic product per capita, from the UN data set, produced with {ggplot2}
::::::
:::
::::{.my-r-code-container}
::: {#lst-chap03-ggplot2-infant-mortality-un}
```{r}
#| label: ggplot2-infant-mortality-un
#| fig-height: 10

gg1_un <- carData::UN |> 
  ggplot2::ggplot(
    ggplot2::aes(
      x = ppgdp,
      y = infantMortality
    ),
  ) +
  ggplot2::geom_point(na.rm = TRUE) +
  ggplot2::geom_smooth(
    formula = y ~ x,
    na.rm = TRUE,
    se = TRUE,
    method = "lm"
  ) +
    ggplot2::geom_smooth(
    formula = y ~ x,
    na.rm = TRUE,
    se = TRUE,
    method = "loess"
  ) +
  ggplot2::labs(    title = "(a)") +
  ggplot2::xlab("GDP per Capita") +
  ggplot2::ylab("Infant Mortality Rate (per 1000 births)") +
  ggplot2::scale_y_continuous(
    limits = c(0, 120), 
    breaks = seq(0, 120, by = 20)
    ) +
  ggplot2::scale_x_continuous(
    limits = c(0, 1e+05), 
    breaks = seq(0, 1e+05, by = 2e+04)
    )

gg2_un <- carData::UN |> 
  ggplot2::ggplot(
    ggplot2::aes(
      x = log(ppgdp),
      y = log(infantMortality)
    ),
  ) +
  ggplot2::geom_point(na.rm = TRUE) +
  ggplot2::geom_smooth(
    formula = y ~ x,
    na.rm = TRUE,
    se = TRUE,
    method = "lm"
  ) +
    ggplot2::geom_smooth(
    formula = y ~ x,
    na.rm = TRUE,
    se = TRUE,
    method = "loess"
  ) +
  ggplot2::labs(title = "(b)") +
  ggplot2::xlab("GDP per Capita") +
  ggplot2::ylab("Infant Mortality Rate (per 1000 births)") +
  ggplot2::theme(
    axis.title.y.right = ggplot2::element_blank(),
    axis.ticks.y.right = ggplot2::element_blank()
  ) +
  ggplot2::scale_y_continuous(
    position = "right",
    breaks = base::log(c(2, 5, 10, 20, 50, 100)),
    labels = NULL,
    sec.axis = ggplot2::sec_axis(
      name = ggplot2::derive(),
      transform = ~ base::exp(1)^ .,
      breaks = c(2, 5, 10, 20, 50, 100)
    )
  )  +
  ggplot2::scale_x_continuous(
    position = "top",
    breaks = log(c(100, 500, 1000, 5000, 10000, 50000, 100000)),
    labels = NULL,
    sec.axis = ggplot2::sec_axis(
      name = ggplot2::derive(),
      transform = ~ exp(1)^ .,
      breaks = c(100, 500, 1000, 5000, 10000, 50000, 100000)
    )
  )



gg1 <- ggExtra::ggMarginal(gg1_un, type = "boxplot")
gg2 <- ggExtra::ggMarginal(gg2_un, type = "boxplot")

gridExtra::grid.arrange(gg1, gg2, nrow = 2)

```
Infant mortality rate and gross domestic product per capita, from the UN data set: (a) untransformed data and (b) both variables log-transformed, using {ggplot2}
:::

***

With the exception of showing the most interesting points is this my replication of Figure 3.16 or here in my notes of @lst-chap03-car-infant-mortality-un.

To replicate Figure 2.16 I had to learn new things. Most of them was maninly an exercise in adapting the axes. I learned several things:

- adding a second axis with transformed values
- adapting axis title and ticks
- setting breaks and labels
- positioning axes
- learning about the many detailed parameters of the `ggplot2::theme()` function

::::
:::::

###### fit model

:::::{.my-r-code}
:::{.my-r-code-header}
:::::: {#cnj-chap03-fit-model-log-data}
: Fitting the linear regression model to the log-transformed data
::::::
:::
::::{.my-r-code-container}
::: {#lst-chap03-fit-model-log-data}
```{r}
#| label: fit-model-log-data

car::brief(
  stats::lm(
    base::log(infantMortality) ~ base::log(ppgdp), 
    data = carData::UN)
  )

```

Summary	of the log-data	fitted model
:::

***

Larger values of `ppgdp` are associated with smaller values of `infantMortality`, as indicated by the negative sign of the coefficient for `log(ppgdp)`.



::::
:::::

Take for instance a hypothetical increase in `ppgdp` in a country of 1%, to $1.01 \times ppgdp$, the expected `infantMortality` would be (approximately):

$$
\alpha_{0}(1.01 \times ppgdp)^{\beta_{1}} = 1.01^{\beta_{1}} \times \alpha_{0}(ppgdp^{\beta_{1}})
$$

For the estimated slope $b^1 = −0.617$, we have $1.01−0.617 = 0.994%, and so the estimated infantMortality would be 0.6% smaller — a substantial amount.

:::

::::
:::::

### Power transformation {#sec-chap03-4-2}

#### Introduction

Transformations have four different goals:

1. **Approximating linearity** as demonstrated in @lst-chap03-car-infant-mortality-un and @lst-chap03-ggplot2-infant-mortality-un. Particularly in multiple linear regression, data with approximately linear mean functions for all 2D plots of the response and regressors, and among the regressors, are very useful.
2. **Smoothing variablity**, again seen in @lst-chap03-car-infant-mortality-un and @lst-chap03-ggplot2-infant-mortality-un.
3. **Inducing a simpler structure** to get a better understanding of the model. See later examples.
4. **Approaching normality** so that the transformed data are as close to normally distributed as possible. This can be useful for statistical methods such as factor analysis that depend on multivariate normality for their justification.

**Simple Power Transformation**

We define a simple power transformation as the replacement of a variable x by $x^λ$, where, in the current general context, x can be either a response variable or a numeric predictor variable, and where λ is a power to be selected based on the data, usually in the range [−1, 1] but occasionally in the range [−3, 3]. The simple power transformations require that all elements of x be strictly positive, x > 0.

This use of power transformations is distinct from **polynomial regression**. In polynomial regression, the highest-order power, say p, is a small positive integer, usually 2 or 3, and rather than replacing a predictor x by $x^p$, we also include all lower-order powers of x in the regression model. For example, for p = 3, we would include the regressors x, $x^2$, and $x^3$ . Polynomial regression in R can be accomplished with the `stats::poly()` function.

The power $λ = −1$ is the **reciprocal or inverse transformation**, 1/x. If, for example, x represents the time to completion of a task (e.g., in hours), then $x^{−1}$ represents the rate of completion (tasks completed per hour). The power λ = 1/3 can convert a volume measure (e.g., in cubic meters) to a linear measure (meters), which may be appropriate in some problems. Likewise, λ= 3 can convert a linear measure (meters) to a volume (cubic meters); λ= 1/2, an area (square meters) to a linear measure (meters); and λ = 2, a linear measure (meters) to an area (square meters).

**Scaled Power Transformations**



:::::{.my-theorem}
:::{.my-theorem-header}
:::::: {#thm-chap03-box-cox-transformation}
: Scaled Power Transformation by Box & Cox
::::::
:::
::::{.my-theorem-container}
$$
\begin{align*}
T_{BC}(x, \lambda) = x^{(\lambda)} =& \\
=& \frac{x^{\lambda}-1}{\lambda} \text{when }\lambda \ne 0 \\
=& log(x) \text{when } \lambda = 0
\end{align*}
$$ {#eq-chap03-box-cox-transformation}

***

BC stands for Box-Cox power transformation [@box1964].

- For $\lambda \ne 0$, the scaled power transformations are essentially $x^{\lambda}$, because the scaled-power family only subtracts 1 and divides by the constant $\lambda$. 
- One can show that as $\lambda$ approaches zero, $T_{BC}(x, \lambda)$ approaches log(x), and so the scaled-power family includes the log transformation when λ = 0, even though the literal zeroth power, $x^0 = 1$, is useless as a transformation. 
- Also, the scaled power $T_{BC}(x, \lambda)$ preserves the order of the x-values, while the basic powers preserve order only when λ is positive and reverse the order of x when λ is negative.
::::
:::::

The scaled power family is generally called the **Box-Cox power family**. Usual practice is to use the Box-Cox power family to estimate the λ parameter of the transformation, but then to use in statistical modeling the equivalent simple power xλ if the selected value of λ is far from zero, and log(x) for λ close to zero. This estimation method is similar to the method of maximum likelihood.

#### Practice

:::::{.my-example}
:::{.my-example-header}
:::::: {#exm-chap03-power-transformations}
: Power transformations
::::::
:::
::::{.my-example-container}

::: {.panel-tabset}

###### p1

:::::{.my-r-code}
:::{.my-r-code-header}
:::::: {#cnj-chap03-power-transformation-infant-mortality-un}
: Demonstration of a transformation of the Box-Cox power family
::::::
:::
::::{.my-r-code-container}
::: {#lst-chap03-power-transformation-infant-mortality-un}
```{r}
#| label: power-transformation-infant-mortality-un

p1 <- car::powerTransform(
  infantMortality ~ 1, 
  data = carData::UN, 
  family = "bcPower")
summary(p1)

```

Estimating a normalizing transformation of `infantMortality` in the UN data

:::

***

- The first argument to `car::powerTransform()` is a formula, with the variable to be transformed on the left of the ~ symbol and the variables, if any, for conditioning on the right. In this instance, there are no conditioning variables, indicated by 1 on the right, so we are simply transforming `infantMortality` unconditionally.
- The argument `family = "bcPower"`, to specify the Box-Cox power family, isn’t strictly necessary because "bcPower" is the default.
- Other transformations are: "bcnPower" for a two-parameter modification of the Box-Cox family that allows negative responses [@hawkins2017], and the "yjPower" family [@yeo2000], another modifiation of the Box-Cox family that allows a few negative values. 
- All three families are documented at [Box-Cox, Box-Cox with Negatives Allowed, Yeo-Johnson and Basic Power Transformations](https://search.r-project.org/CRAN/refmans/car/html/bcPower.html).

::::
:::::

The first part of the output gives information about the estimated power. 

- **Est Power**: The value under "Est Power" is the point estimate found using the procedure suggested by Box and Cox. 
- **Rounded Pwr**: The second value, labeled "Rounded Pwr", is a rounded version of the estimate. 
- **Wald Lwr Bnd / Wald Upr Bnd**: The remaining two values are the bounds of the 95% Wald confidence interval for λ. The confidence interval includes $\lambda = 0$, which in turn suggests using a log-transform because of our preference for simple transformations. In the example, zero is included in the confidence interval, suggesting the use in practice of the log transformation.
- **Tests**: The remainder of the output provides tests concerning λ, first testing `H0: λ= 0` for the log transformation and then testing `H0: λ= 1`, which is equivalent to no transformation. The test for the log transformation has a very large p-value, indicating that the log transformation is consistent with the data, while the tiny p-value for λ = 1 indicates that leaving `infantMortality` untransformed is inconsistent with the goal of making the variable normally distributed.


###### λ = 0.5

:::::{.my-r-code}
:::{.my-r-code-header}
:::::: {#cnj-chap03-test-transform-un}
: Test UN data for trasnformation $\lambda = .5$
::::::
:::
::::{.my-r-code-container}
::: {#lst-chap03-test-transform-un}  
```{r}
#| label: test-transform-un 

car::testTransform(p1, lambda = 1/2)
```

Test for another value of $\lambda$, in this case $\lambda = .5$
:::

::::
:::::

The very small `r glossary("p-value")` for the test of λ = 1/2 suggests that we shouldn’t use this transformation.

###### add variable

:::::{.my-r-code}
:::{.my-r-code-header}
:::::: {#cnj-chap03-add-transformed-variable}
: Add transformed variable to data set
::::::
:::
::::{.my-r-code-container}
::: {#lst-chap03-add-transformed-variable}
```{r}
#| label: add-transformed-variable

UN <- transform(
  carData::UN, 
  infantMortality.tran = infantMortality^p1$lambda
  )

skimr::skim(UN)
```

Add transformed variable to data set
:::

***

- The updated UN data set contains all the variables in the original data set plus an additional variable called `infantMortality.tran`. 
- The value of `p1$lambda` in `infantMortality^p1$lambda` is the estimated power transformation parameter.

::::
:::::

###### p2

:::::{.my-r-code}
:::{.my-r-code-header}
:::::: {#cnj-chap03-multivariate-transformation-prestige-data}
: Mutivariate transformation
::::::
:::
::::{.my-r-code-container}
::: {#lst-chap03-multivariate-transformation-prestige-data}
```{r}
#| label: multivariate-transformation-prestige-data

base::summary(p2 <- car::powerTransform(
  base::cbind(income, education) ~ 1,
  data = carData::Prestige,
  family = "bcPower"
  )
)
```

Multivariate transformation of `income` and `education` toward bivariate normality
:::

***

Let $λ = (λ1, λ2)$ represent the vector of power transformation parameters. Looking at the likelihood-ratio tests, untransformed variables, corresponding to $λ = (1, 1)′$, and all logarithms, corresponding to $λ = (0, 0)′$, both have small p-values, implying that neither of these choices is appropriate. 

The suggested procedure is to use the rounded values of the power parameters (1/3, 1) in further work.

::::
:::::

The multivariate extension of Box-Cox power transformations was provided by Velilla [-@velilla1993] and is also implemented by the `car::powerTransform()` function.

A useful early step in any regression analysis is considering transformations toward linearity and normality of all the quantitative predictors, without at this point including the response.

###### λ = 1/3 and (1, 0)

:::::{.my-r-code}
:::{.my-r-code-header}
:::::: {#cnj-chap03-test-transform-prestige}
: Test UN data for transformation $\lambda = .33$ and log (0, 1)
::::::
:::
::::{.my-r-code-container}
::: {#lst-chap03-test-transform-prestige}
```{r}
#| label: test-transform-prestige
#| results: hold

glue::glue(" ")
glue::glue("############################################################")
glue::glue("test transformation with lambda = 0.33")
car::testTransform(p2, lambda = c(1/3, 1))

glue::glue(" ")
glue::glue("############################################################")
glue::glue("test log transformation with lambda = (0, 1)")

car::testTransform(p2, lambda = c(0, 1))
```

Test UN data for transformation $\lambda = .33$ and log (0, 1)
:::

***

The nominal p-value is about 0.3, but interpreting this result is complicated by the use of the observed data to select the null hypothesis. As with a one dimensional transformation, the value of estimation and of the rounded estimates are stored respectively in `p2$lambda` (p2[["lambda"]][["income"]], p2[["lambda"]][["education"]] = 0.262, 0.424) and `p2$roundlam` (p2[["roundlam"]][["income"]], p2[["roundlam"]][["education"]] = 0.33, 1).

Because we generally prefer logarithms we have also tried $\lambda = (0, 1)$.

Interpretation of the p-values is again difficult because the data were used to generate the hypothesis.



::::
:::::

The relatively small p-value suggests that the choice of $λ = (1/3, 1)$ is to be preferred to $λ = (0, 1)$.

###### log

:::::{.my-r-code}
:::{.my-r-code-header}
:::::: {#cnj-chap03-log-transformation-scatterplot}
: Log transformation of `income`
::::::
:::
::::{.my-r-code-container}
::: {#lst-chap03-log-transformation-scatterplot}
```{r}
#| label: log-transformation-scatterplot

car::scatterplotMatrix(
  ~ prestige + log(income) + education + women,
  smooth = list(span = 0.7),
  data = carData::Prestige
)

```

Scatterplot matrix for variables in the Prestige data set, with log-transformed `income.`
:::

***
The panels in this scatterplot matrix show little nonlinearity, although in the plot of `prestige` versus `log(income)` the two occupations with the lowest incomes appear not to fit the pattern in the rest of the data.

::::
:::::

###### cube

:::::{.my-r-code}
:::{.my-r-code-header}
:::::: {#cnj-chap03-cube-transformation-scatterplot}
: Cube transformation of `income`
::::::
:::
::::{.my-r-code-container}
::: {#lst-chap03-cube-transformation-scatterplot}
```{r}
#| label: cube-transformation-scatterplot

car::scatterplotMatrix(
  ~ prestige + income^3 + education + women,
  smooth = list(span = 0.7),
  data = carData::Prestige
)

```

Scatterplot matrix for variables in the Prestige data set, with cube-transformed `income.`
:::


::::
:::::

###### p3

The variable type in the Prestige data set is a factor that divides occupations into blue-collar ("bc"), white-collar ("wc"), and professional/managerial ("prof") categories. Examination of the data suggests that the distributions of the numeric variables in the data set vary by occupational type. We can then seek transformations of the numeric variables that make their within-group distributions as close as possible to multivariate normal with common covariance structure (i.e., similar variances and covariances within types of occupations).

:::::{.my-r-code}
:::{.my-r-code-header}
:::::: {#cnj-chap03-prestige-types}
: Conditioning on groups
::::::
:::
::::{.my-r-code-container}
::: {#lst-chap03-prestige-types}
```{r}
#| label: prestige-types

summary(p3 <- car::powerTransform(cbind(income, education) ~ type,
                             data = carData::Prestige))

car::testTransform(p3, c(0, 1))
```

Listing title
:::

::::
:::::

###### cond

:::::{.my-r-code}
:::{.my-r-code-header}
:::::: {#cnj-chap03-scatterplot-cond-types}
: Revised scatterplot matrix
::::::
:::
::::{.my-r-code-container}
::: {#lst-chap03-scatterplot-cond-types}
```{r}
#| label: scatterplot-cond-types

car::scatterplotMatrix(
  ~ prestige + log(income) + education + women | type,
  data = carData::Prestige,
  smooth = FALSE,
  ellipse = list(levels = 0.5),
  legend = list(coords = "center")
)
```

Scatterplot matrix with log-transformed income and conditioning on type of occupation, for the Prestige data set. The ellipse in each panel is a 50% concentration ellipse assuming bivariate normality.
:::

***

- We suppress the loess smooths with the argument `smooth = FALSE` to minimize clutter in the graph. 
- We move the legend from its default position in the top-right corner of the first diagonal panel to the center of the panel. 
- We set the argument ellipse=list (levels=0.5) to get separate 50% concentration ellipses for the groups in the various off-diagonal panels of the plot (see, e.g., [@friendly2013]). --- If the data in a panel are bivariately normally distributed, then the ellipse encloses approximately 50% of the points.




::::
:::::

- Approximate within-group linearity is apparent in most of the panels. 
- The size of the ellipse in the vertical and horizontal directions reflects the standard deviations of the two variables, and its tilt reflects their correlation. The within-group variability of prestige, log (income), education, and women is therefore reasonably similar across levels of type; the tilts of the ellipses in each panel are somewhat more variable, suggesting some differences in within-group correlations. 
-The least-squares lines in the panels also show some differences in within-group slopes. 
- The plotted points are the same in @lst-chap03-log-transformation-scatterplot resp. @lst-chap03-cube-transformation-scatterplot and @lst-chap03-scatterplot-cond-types. Only the symbols used for the points, and the fitted lines, are different. In particular, the two points with low income no longer appear exceptional when we control for type. 
- A small difference between the two graphs is that type is missing for four occupations, and thus @lst-chap03-scatterplot-cond-types has only 98 of the 102 points shown in @lst-chap03-log-transformation-scatterplot and @lst-chap03-cube-transformation-scatterplot.

###### p4

We may wish to apply the Box and Cox transformation methodology to a variable with a small number of zero or even negative values. For example, 28 of the 240 companies in the Ornstein data set maintained zero interlocking directors and officers with other firms; the numbers of interlocks for the remaining companies were positive, ranging between 1 and 107.

The Box-Cox family with negatives [@hawkins2017] is defined to be the Box-Cox power family applied to z(x,γ) rather than to x: $T_{BCN} (x, λ,γ) = T_{BC} (z(x,γ),λ)$. If x is strictly positive and γ = 0, then z(x, 0) = x and the transformation reduces to the regular Box-Cox power family. The Box-Cox family with negatives is defined for any γ > 0 when x has nonpositive entries. 

The important improvement of this method over just adding a start is that for any γ and large enough x, the transformed values with power λ can be interpreted as they would be if the Box-Cox family had been used and no nonpositive data were present.

:::::{.my-r-code}
:::{.my-r-code-header}
:::::: {#cnj-chap03-special-transformation}
: Transformations With a Few Negative or Zero Values
::::::
:::
::::{.my-r-code-container}
::: {#lst-chap03-special-transformation}
```{r}
#| label: special-transformation

p4 <-
  car::powerTransform(
    interlocks ~ 1, 
    data = carData::Ornstein, 
    family = "bcnPower")
summary(p4)
```

Using car::powerTransform() with the "bcnPower" family, invoking the `car:::bcnPower()` function
:::

***

Employing the "bcnPower" family estimates both the power λ and the shift γ transformation parameters using Box and Cox’s maximum-likehood-like method. 

As is often the case with the "bcnPower" family, the estimate of the shift γ is very close to the boundary of zero, in which case `car::powerTransform()` fixes γ = 0.1 and then estimates the power λ. In the example, it is close to the cube root. The tests shown in the output are for the fixed value of γ = 0.1; had γ been estimated, the tests would have been based on the log-likelihood profile averaging over γ. 

::::
:::::



:::

::::
:::::


### Transformations and EDA {#sec-chap03-4-3}

#### Introduction

Exploratory Data Analysis (`r glossary("ExDA")`) proposed by Tukey and Mosteller [@tukey1977; @mosteller2019], is a collection of method that can be used to understand data, mostly without requiring models.

The {**car**} package implements some methods for selecting a transformation. These methods are of some continuing utility, and they can also help us to understand how, why, and when transformations work.

Mosteller and Tukey characterize the family of power transformations $x^λ$ as a *ladder of powers and roots*, with no transformation, λ = 1, as the central rung of the ladder. Transformations like λ = 1/2 (square root), λ = 0 (taken, as in the BoxCox family, as the log transformation), and λ = −1 (the inverse transformation) entail moving successively down the ladder of powers, and transformations like λ = 2 (squaring) and λ = 3 (cubing) entail moving up the ladder.

#### Transforming for symmetry

Transformation down the ladder of powers serves to spread out the small values of a variable relative to the large ones and consequently can serve to correct a positive skew. Negatively skewed data are much less common but can be made more symmetric by transformation up the ladder of powers. The `symbox()` function in the {**car**} package implements a suggestion of Tukey and Mosteller to induce approximate symmetry by trial and error.

See: @lst-chap03-symmetry.

#### Transforming toward linearity

If a relationship is nonlinear but monotone (i.e., strictly increasing or strictly decreasing) and simple (in the sense that the curvature of the regression function is relatively constant), then Mosteller and Tukey’s bulging rule, illustrated in Figure 3.20, can guide the selection of linearizing power transformations of the two variables. 

When the bulge of the regression function in a scatterplot points *down* and to the *right*, for example — as suggested by the line in the lower-right quadrant of Figure 3.20 we move y *down* the ladder of powers, x *up* the ladder of powers, or both, to straighten the relationship between the two variables. The bulging rule generalizes to the other quadrants of Figure 3.20.

![Mosteller and Tukey’s bulging rule for finding linearizing transformations](img/chap03-bulging-rule-min.png
){#fig-bulgin-rule
fig-alt="Circle shows the bulging rule: When the bulge points down, transform y down the ladder of powers; when the bulge points up, transform y up; when the bulge points left transform x down; when the bulge points right, transform x up."
fig-align="center"}

To illustrate, let's inspect @lst-chap03-car-infant-mortality-un once more:

- The nonlinear relationship between `infantMortality` and `ppgdp` @lst-chap03-car-infant-mortality-un (a) is simple and monotone.
- The bulge in the scatterplot points down and to the left, suggesting the transformation of `infantMortality` or `ppgdp`, or both, down the ladder of powers. 
- In this case, the log transformation of both variables produces a nearly linear relationship, as shown in  @lst-chap03-car-infant-mortality-un (b).

#### Transforming to equalize spread

In @lst-chap03-parallel-boxplot-car, we examined the relationship between number of interlocks and nation of control among the 248 large Canadian corporations in Ornstein’s interlocking-directorate data. 

That graph revealed an association between *level* and *spread*: Nations with a relatively high level of interlocks (Canada, other) evidence more variation than nations with fewer interlocks on average (United States, United Kingdom). The term **level** is meant to imply an average or typical value, such as the mean, and **spread** is meant to imply variability, as measured, for example, by the SD. 

These terms are used in exploratory data analysis in place of more technical terms such as mean and SD because level and spread may be measured using other statistics that are resistant to outlying values. 

A spread-level plot [@tukey1977] is a scatterplot of the logarithm of the interquartile range, which measures spread, versus the logarithm of the within-group median, which measures level. Interquartile ranges and medians are insensitive to a few outliers, and so the spread-level plot provides a robust representation of the dependence of spread on level. 

- If there is an association between spread and level in the spread-level plot, Tukey suggested fitting a straight line to the plot, with estimated slope $b$, and then transforming the response variable $y$ by the power transformation $y^{1−b}$.
- If $b \approx 1$, then a log transformation is used. 

See: @lst-chap03-equalize-spread.

### Transforming restricted-range variables {#sec-chap03-4-4}

Variables with ranges bounded both below and above may require different transformations because of `r glossary("ceiling")` or `r glossary("floor")` effects, in which values tend to accumulate near the boundaries of the variable. 

The most common kinds of restricted-range variables are percentages and proportions, but other variables, such as the number of questions correct on a test of fixed length, are also “disguised proportions.” 

Because small changes in a proportion near zero or 1 are potentially much more important than small changes near 0.5, promising transformations spread out both large and small proportions relative to proportions near 0.5. If a proportion has a more limited observed range, say between 0.25 and 0.75, transformations are unlikely to be helpful. The idea is to transform proportions so the transformed values behave more like a sample from a normal distribution. Usual methods of statistical inference for normally distributed data can then be used with the transformed proportions as the response. 

Statistical models specifically designed for proportion data — for example, logistic-regression and loglinear models (discussed in @sec-chap06) — are the preferred methods for working with proportions, and so transformation for normality is now rarely used with proportion response variables if the counts on which the proportions are based are available. 

Transforming predictor variables that are expressed as proportions or percents may prove useful. Suppose that $u_{i}$ is a sample count of number of “successes” for the $i$th of $n$ cases and that $N_{i}$ is the corresponding number of trials. Each $u_{i}$ can have its own probability of success $θ_{i}$. From the binomial distribution, the standard deviation of the sample proportion $x_{i} = u_{i} / N_{i}$ is ????, and therefore varies depending on the value of $θ_{i}$. If we transform $x_{i} = u_{i} / N_{i}$ using the arcsine square root, then the SD of the transformed variable is approximately for any value of $θ_{i}$, and hence, $T_{asinsqrt}(x)$ is a variance-stabilizing transformation. 

**Logistic regression** uses a different transformation of a proportion called the logit or log-odds transformation. The **logit transformation** is not defined for sample proportions exactly equal to zero or 1, however[^03-exploring-data-2]. 

[^03-exploring-data-2]: This restriction does not affect logistic regression, where it is the probability of successθ, and not the observed sample proportion of successes x, that is transformed.

We can get around this limitation if we want to use sample logits in a graph or statistical analysis by remapping proportions from the interval [0, 1] to [0.025, 0.975], for example, taking the logit of $0.025 + 0.95 \times x$ rather than the logit of x. The `logit()` function in the {**car**} package takes care of remapping proportions or percentages when there are zeros or ones—or 0% or 100% for percentage data — printing a warning if remapping is required.

See: @lst-chap03-logit-transformation.

This restriction does not affect logistic regression, where it is the *probability* of success θ, and not the observed sample *proportion* of successes x, that is transformed. Even better, if we have access to the original data from which the proportions are calculated, we can avoid proportions of zero or 1 by computing *empirical logits*, $log[(u + 1/2)/(N + 1)]$, where, as before, $u$ is the number of successes in $N$ trials. 

To illustrate, we apply the logit and arcsine-square-root transformations to the distribution of the gender composition of occupations (percentage women) in the Prestige data set: See: @lst-chap03-transformation-restricted-range-variables.



:::::{.my-example}
:::{.my-example-header}
:::::: {#exm-eda-transformation}
: Transformation and Exploratory Data Analysis
::::::
:::
::::{.my-example-container}

::: {.panel-tabset}

###### symmetry

:::::{.my-r-code}
:::{.my-r-code-header}
:::::: {#cnj-chap03-symmetry}
: Transforming for symmetry
::::::
:::
::::{.my-r-code-container}
::: {#lst-chap03-symmetry}
```{r}
#| label: symmetry

car::symbox( ~ ppgdp, data = carData::UN)

```
Boxplots of various power transformations of per-capita GDP in the United Nations data.
:::

***

By default, `car::symbox()` uses the function `car:::bcPower()` and displays boxplots of the transformed variable for several transformations down the ladder of powers. 

::::
:::::

The graphic shows that the ubiquitous log transformation of `ppgdp` does the best job of making the distribution of the variable symmetric.

###### equalize spread

:::::{.my-r-code}
:::{.my-r-code-header}
:::::: {#cnj-chap03-equalize-spread}
: Spread-level plot
::::::
:::
::::{.my-r-code-container}
::: {#lst-chap03-equalize-spread}
```{r}
#| label: equalize-spread

car::spreadLevelPlot(
  interlocks + 1 ~ nation, 
  data = carData::Ornstein)
```

Spread-level plot for the relationship between number of interlocks and nation of control in Ornstein’s interlocking-directorate data
:::

***

- Because Tukey’s method requires a strictly positive response variable, we add a start of 1 to interlocks to avoid zero values. 
- In addition to drawing a graph, the `car::spreadLevelPlot()` function prints a table showing the first quartile (Tukey’s lower hinge), median, third quartile (upper hinge), and interquartile range (hinge spread ) in each group, along with the suggested power transformation of interlocks + 1. 

::::
:::::

In the example, the suggested power is p = 0.15, close to the log transformation.

###### parallel boxplots

:::::{.my-r-code}
:::{.my-r-code-header}
:::::: {#cnj-chap03-parallel-boxplots-ornstein}
: Parallel boxplots for the log-transformed Ornstein data
::::::
:::
::::{.my-r-code-container}
::: {#lst-chap03-parallel-boxplots-ornstein}
```{r}
#| label: parallel-boxplots-ornstein

par(mar = c(5.1, 4.1, 4.1, 4.1))
car::Boxplot(
  log10(interlocks + 1) ~ nation, 
  data = carData::Ornstein
  )
car::basicPowerAxis(
  power = 0,
  base = 10,
  at = c(1, 3, 6, 11, 21, 51, 101),
  start = 1,
  axis.title = "Interlocks"
)
```

Parallel boxplots for $log_{10}(interlocks + 1)$ by nation of control
:::

***

- Many global graphics parameters in R are set or queried by the base::par() function. 
- The "mar" setting is for the plot margins as a numberical vector of the form  `c(bottom, left, top, right)` which gives the number of lines of margin to be specified on the four sides of the plot. The default is `c(5, 4, 4, 2) + 0.1`.
- We specify two arguments in the call to `car::Boxplot()`: 
  - a formula with the base-10 logarithm of interlocks + 1 on the left-hand side and the factor nation on the right-hand side, and
  - the argument data = carData::Ornstein, supplying the data frame in which these variables reside. 
- The remaining commands make the graph more elaborate, first by increasing the right-side margin of the plot and then by adding a second axis on the right, labeled in the original units of interlocks, produced by the `car::basicPowerAxis()` function[^03-exploring-data-1]. 
  - The argument `power = 0` to `car::basicPowerAxis()` specifies the log transformation; 
  - `base = 10` gives the base used for the logs; 
  - `at = c(1, 3, 6, 11, 21, 51, 101)` indicates where the tick marks are to appear on the interlocks + 1 scale; and 
  - `start = 1` gives the start that was used, so that the tick labels can be adjusted to the original interlocks-scale.
  
[^03-exploring-data-1]: The functions `car:::bcPowerAxis()`, `car:::bcnPowerAxis()`, and `car:::probabilityAxis()` may be used similarly to produce axes on the untransformed scale corresponding to Box-Cox, Box-Cox with negatives, and logit transformations.

::::
:::::

The spreads in the transformed data for the four groups are much more similar than the spreads in the untransformed data (shown in @lst-chap03-ornstein-boxplot). 

###### logit

:::::{.my-r-code}
:::{.my-r-code-header}
:::::: {#cnj-chap03-logit-transformation}
: Logit transformation
::::::
:::
::::{.my-r-code-container}
::: {#lst-chap03-logit-transformation}
```{r}
#| label: logit-transformation
#| results: hold

glue::glue(" ")
glue::glue("############################################################")
glue::glue("Logit transformation for the range 0.1 to 0.9")
glue::glue("############################################################")
car::logit(base::seq(0.1, 0.9, 0.1))
glue::glue(" ")
glue::glue("############################################################")
glue::glue("Logit transformation for the range 0 to 1")
glue::glue("############################################################")
car::logit(base::seq(0, 1, 0.1))
```

Logit transformation with remapping proportions or percentages when there are zeros or ones
:::

::::
:::::


###### restricted range



:::::{.my-r-code}
:::{.my-r-code-header}
:::::: {#cnj-chap03-transformation-restricted-range-variables}
: Transformation of restricted-range Variables
::::::
:::
::::{.my-r-code-container}
::: {#lst-chap03-transformation-restricted-range-variables}
```{r}
#| label: transformation-restricted-range-variables
#| results: hold

graphics::par(mfrow = c(1, 3))

car::densityPlot(
  ~ women,
  data = carData::Prestige,
  from = 0,
  to = 100,
  main = "(a) Untransformed"
  )

car::densityPlot( ~ car::logit(women),
             data = carData::Prestige,
             adjust = 0.7,
             main = "(b) Logit"
             )

car::densityPlot( ~ base::asin(base::sqrt(women / 100)),
             data = carData::Prestige,
             adjust = 0.7,
             main = "(c) Arcsine square-root"
             )

```

Distribution of women in the Canadian occupational-prestige data: (a) untransformed, (b) logit transformed, and (c) arcsine square root transformed.
:::

::::
:::::


:::

::::
:::::

### Other transformations {#sec-chap03-4-5}

The discussed transformation methods are my no means a complete catalogue. There are many other possible transformationsd, such as:

- People’s height in cm and weight in kg are often combined in the body mass index, $BMI = weight / (height^2)$, which is intended to approximate body composition in units of $kg/cm^2$. 
- A variable like month.number can be replaced in a regression model by two regressors, `sin(month.number/12)` and `cos(month.number/12)`, to model seasonal time trends. 
- In a study of highway accident rates, a variable giving the number of traffic signals in a highway segment is converted to the number of signals per mile by dividing by the length of the segment. 

::: {.callout-important #imp-cahp03-transformation-when-why}
##### Lesson learned about transformations

One should think carefully about how variables are expressed in the substantive context in which they are used in research.

As a general rule, measures of total size often need to be converted to size per unit, for example converting GDP to a per-capita basis by dividing by population.

Transformations are very important and there are many different options. The {**car**} package provides some tools to learn not only about the optimal method of transformations but also about parameters for it’s application.
:::

## Point labelling and identification {#sec-chap03-5}

Identifying extreme points can be especially valuable in graphs used for model building and diagnostics.
    
### The graphics::identify() function {#sec-chap03-5-1}

The `graphics::identify()` function in base R has a (clumsy) interface for interactive marking of interesting points.

:::::{.my-procedure}
:::{.my-procedure-header}
:::::: {#prp-chap03-graphics-identify}
: Displaying interesting points with graphics::identify()
::::::
:::
::::{.my-procedure-container}
`graphics::identify()` is	interactive, so you can’t use	it	conveniently in R	Markdown or Quarto documents.

1. Enter the plot command to create the graphics where you want to label interesting points. You have to do this at the console! For instance: `with(carData::Freedman, plot(density, crime))`. This command produces the graphic in the RStudio "Plots" tab. 
2. Enter at the console the command for `graphics::identify()`. It needs at least three arguments: The x/y coordinates for the points and a vector with the labels. If the label vector is missing `identify()` uses the row number. The RStudio "Plots" tab is now interactive. You can see this as the message "Locator active (Esc to finish)" appears at the upper left of the RStudio "Plots" tab and a button "Finish" at the upper right. Additionally you see two icons at the upper right of the console window. When you hover over the green circle with the $\infty$ sign a tool tip appears: "Session suspend timeout paused: Waiting for event: locator_ completed R is executing". The other icon is a stop sign.
3. Click at (or near) the points you want to label. 
4. A “pin” icon flashes momentarily near the point, but point labels aren’t shown until you exit from point identification mode.
5. After you have clicked on all points you are interested in, press "Esc" or the "Finish" button.
6. The labels appear near the clicked positions and the console prints the row numbers of the points.
7. To include the plot with the labelled points in your report you have either export the graph with the options "Save as image…", "Save as PDF… or "Copy to clipboard…" or to create  the graph programmatically (for instance using {**ggplot2**} and {**ggrepel**}).

![Scatterplots	of	crime	by	population	density	for	Freedman’s	data, with a few high-density cities identified manually by the graphics::indentify() function](img/chap03-rplot-identify-min.png){#fig-rplot-identify
fig-alt="Scatterplots of crime by population density for Freedman’s data, with a few high-density cities identified manually by the graphics::indentify() function"
fig-align="center"}

::::
:::::


### Automatic point labelling with car::showLabels() {#sec-chap03-5-2}

The graphics functions in the {**car**} package employ a common general mechanism for point identification, using the `car::showLabels()` function to identify potentially noteworthy points.

:::{.my-bulletbox}
:::: {.my-bulletbox-header}
::::: {.my-bulletbox-icon}
:::::
:::::: {#bul-chap03-arguments-showLabels}
::::::
: Arguments of `car::showLabels()` function
::::
:::: {.my-bulletbox-body}

- **Point identification** is controlled by the `id` argument, which takes the values TRUE, FALSE, or a detailed list of specifications. In	most	cases,	the	user	only	needs	to	set	`id = TRUE`	or to	specify	the	number	of	points	to	be	identified, for	example,	`id = list(n = 5)`.
- **labels**: By default, points are labeled with a data frame’s row labels or by row numbers if there are no row labels. You can use whatever point labels you like by setting labels to a character vector with the same number of elements as there are data points. 
- **n**: the number of points to label. 
- **cex**: the relative size of the labels; cex=0.5, for example, produces labels half the default size. 
- **col**: the color of the labels; if not set, the color is determined automatically. 
-**location**: where the labels are drawn. The default is `location = "lr"` to draw labels to the left of points in the right half of the graph and to the right of points in the left half. Another option is `location="ab"`, to draw labels above points below the middle of the graph and below points above the middle. Finally, `location="avoid"` tries to avoid overplotting labels. 
- **method**: Setting `method = "identify"` (the default) enables interactive point identification, as with `graphics::identify()`. Several automatic methods are available for determining how interesting points are to be identified. For example, in `car::scatterplot()`, the default is `method = "mahal"`, in which the `r glossary("Mahalanobis")` distance of each point to the centriod (point of averages) is computed, and the identified points are those with the largest Mahalanobis distances. For `car::residualPlot()`, the default is `method = "r"`, selecting noteworthy points according to the absolute values of their vertical coordinates (i.e., the residuals) and labeling the points with the largest values. There are many other options for the method argument; see help ("showLabels") for their description and the help pages for particular car graphics functions for how they use `car::showLabels()`.


::::
:::



::: {.callout-important #imp-chap03-show-labels-with-ggplot}
##### Developing my own show_labels() command

I am trying to replicate the labelling of interesting points with {**ggplot2**}. The main idea for this additional work is to use the {**tidyverse**} approach and to stick with the more exhaustive functionality of the grammar of graphics.

This is not only programmatically a challenge but I have to learn and understand what type of graphics shows what type of "interesting" points (outliers, extreme values, high-leverage points) and how to compoute these points.

So far I have code proposals for

- all `r glossary("outliers")` in `r glossary("boxplots")` (@lst-chap03-compute-boxplot-outliers) and
- `r glossary("Mahalanobis")` distances with `r glossary("p-value")` of .001 in a linear model (@lst-chap04-davis-ggrepel)

I am still missing the option to display a certain number of extreme points.


:::