# Statistics for Linguists using R - Quantitative and qualitative approaches

Jalal Al-Tamimi

<figure><img src="image_1.png" /></figure>

## To cite this version:

Jalal Al-Tamimi. Statistics for Linguists using R - Quantitative and qualitative approaches. , 2025, Jalal Al-Tamimi, 10.5281/zenodo.17321851. hal-05487648

HAL Id: hal-05487648

https://hal.science/hal-05487648v1

Submitted on 2 Feb 2026

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_2.png" /></figure>

Distributed under a Creative Commons CC BY 4.0 - Attribution - International License
Statistics for Linguists using R - Quantitative
and qualitative approaches

Jalal Al-Tamimi

2025-12-14
<table><tr><td>3.3</td><td>Linear Mixed-effects Models. Why random effects matter</td><td>181</td></tr><tr><td>3.4</td><td>Generalised Linear Mixed-effects Models</td><td>217</td></tr><tr><td>3.5</td><td>Cumulative Logit Link Mixed-effects Models</td><td>424</td></tr><tr><td>3.6</td><td>Generalised Additive Mixed-effects Models (GAMMs)</td><td>443</td></tr><tr><td>3.7</td><td>Other distributions</td><td>462</td></tr><tr><td>3.8</td><td>session info</td><td>463</td></tr></table>

III Qualitative data analysis and corpora 465

4 Creating a corpus and basic manipulations 467

<table><tr><td>4.1</td><td>Loading packages</td><td>467</td></tr><tr><td>4.2</td><td>janeaustenr</td><td>469</td></tr><tr><td>4.3</td><td>proustr</td><td>541</td></tr><tr><td>4.4</td><td>Inaugural Corpus USA</td><td>761</td></tr><tr><td>4.5</td><td>Universal Declaration of Human Rights</td><td>792</td></tr><tr><td>4.6</td><td>Twitter data</td><td>903</td></tr><tr><td>4.7</td><td>Single web page</td><td>1376</td></tr><tr><td>4.8</td><td>Multiple webpages</td><td>1393</td></tr><tr><td>4.9</td><td>session info</td><td>1461</td></tr></table>

5 Advanced manipulations of corpora 1465

<table><tr><td>5.1</td><td>Loading packages</td><td>1465</td></tr><tr><td>5.2</td><td>janeaustenr</td><td>1466</td></tr><tr><td>5.3</td><td>proustr</td><td>1481</td></tr><tr><td>5.4</td><td>Inaugural Corpus USA</td><td>1493</td></tr><tr><td>5.5</td><td>Universal Declaration of Human Rights</td><td>1505</td></tr><tr><td>5.6</td><td>Twitter data</td><td>1543</td></tr><tr><td>5.7</td><td>Single web page</td><td>1561</td></tr><tr><td>5.8</td><td>Multiple webpages</td><td>1572</td></tr><tr><td>5.9</td><td>session info</td><td>1587</td></tr></table>

6 Corpus analysis and statistics 1591

<table><tr><td>6.1</td><td>Loading packages</td><td>1591</td></tr><tr><td>6.2</td><td>janeaustenr</td><td>1592</td></tr><tr><td>6.3</td><td>proustr</td><td>1737</td></tr><tr><td>6.4</td><td>Inaugural Corpus USA</td><td>2111</td></tr><tr><td>6.5</td><td>Universal Declaration of Human Rights</td><td>2160</td></tr><tr><td>6.6</td><td>Twitter data</td><td>2339</td></tr><tr><td>6.7</td><td>Single web page</td><td>2627</td></tr><tr><td>6.8</td><td>Multiple webpages</td><td>2646</td></tr><tr><td>6.9</td><td>session info</td><td>2689</td></tr></table>
IV Data reduction techniques and machine learning 2693

7 Reducing Dimensionality: Correlations, Principal Component Analysis, cluster analysis and Multidimensional scaling 2695

<table><tr><td>7.1 Loading packages</td><td>2696</td></tr><tr><td>7.2 Read dataset</td><td>2696</td></tr><tr><td>7.3 Correlation tests</td><td>2697</td></tr><tr><td>7.4 Principal Component Analyses - PCA</td><td>2712</td></tr><tr><td>7.5 Kmeans Clustering</td><td>2769</td></tr><tr><td>7.6 Multidimensional scaling</td><td>2801</td></tr><tr><td>7.7 session info</td><td>2803</td></tr></table>

8 Decision trees and Random Forests 2807

<table><tr><td>8.1 What is predictive modelling</td><td>2807</td></tr><tr><td>8.2 Loading packages</td><td>2808</td></tr><tr><td>8.3 Read dataset</td><td>2808</td></tr><tr><td>8.4 GLM as a classification tool</td><td>2809</td></tr><tr><td>8.5 Issues with GLM (and regression analyses in general)</td><td>2813</td></tr><tr><td>8.6 Decision Trees</td><td>2822</td></tr><tr><td>8.7 Random Forests</td><td>2831</td></tr><tr><td>8.8 session info</td><td>2855</td></tr></table>
# Introduction

This book covers the basics of statistics and R programming for linguists. It is based on various courses I have taught at both Newcastle University, Newcastle, UK during various workshops and mainly at the Université Paris Cité, Paris, France during my teaching at undergraduate and postgraduate levels.

You can use this reference for citation:

Al-Tamimi, J. (2025). *Statistics for Linguists using R - Quantitative and qualitative approaches*. Bookdown. https://jalalal-tamimi.github.io/R-Stats-Linguistics-Book/. Doi: https://doi.org/10.5281/zenodo.17249197.

The book is designed to be accessible to those with little or no prior experience in statistics or programming. The book applies these concepts to real-world examples from linguistics, including phonetics, syntax, and semantics. The book also includes exercises and solutions to help readers practice and reinforce their understanding of the material. The book is divided into four main parts.

The first part covered in **chapter 1** introduces the basics of the R programming language and the tidyverse universe. It covers data manipulation, data cleaning and an introduction to the basics of descriptive statistics. We end the part with data visualisation.

The second part focuses on statistical analyses using R covering inferential statistics, and hypothesis testing. This part aims to consolidate the reader's understanding of fundamental statistical concepts and how they can be applied using R and is composed of two chapters. In **chapter 2**, we cover type of outcomes and the best statistical models adapted to each type of outcome, and cover type of errors in statistical designs, including the classical Type I and Type II errors, and Type M and S errors. We also cover statistical power, and effect sizes. This part then covers correlations and three main family distributions, Gaussian (via linear models), Binomial (via GLM) and Cumulative (via Ordinal or likert scale data). We introduce the Signal Detection Theory concepts in a way that allows the reader to understand the various notions linked with STD and GLM models.

This part continues with the **chapter 3**, which introduces Random effects structures via mixed effects modelling. We cover the notions of random effects, ran-
dom intercepts and random slopes, and modelling strategies. We look at LMMs, GLMM, CLMM and GAMMs (for additive models).

The third part covers qualitative research methods using R. It is composed of three chapters and includes examples of where to find textual corpora and creating a corpus (in **chapter 4**), how to manipulate a corpus, including part of speech tagging (in **chapter 5**) and how to conduct statistical analyses of textual data including poisson regressions, network analyses and word clouds (in **chapter 6**).

The last part of the book is composed of two chapters and introduces the reader in **chapter 7** with various techniques used for dimensionality reduction, including correlations thresholds, Principal Component Analysis (PCA), cluster analysis and Multi Dimensional Scaling (MDS). **Chapter 8** covers the basics of machine learning, exploring issues related with regression analyses, and introducing solutions in terms of decision trees and Random Forests. We look at two different implementations of Random Forests, via Conditional Inference Trees and via the basic implementation using the **ranger** package, for a faster computation. We end the chapter by introducing the **tidymodels** philosophy.

The structure is organised with each chapter being dedicated to a specific topic and can normally be covered in 1 or 2 sessions.

It is hoped that this book allows students to specialise in the field of statistical analyses applied to linguistic data and to be able to use R for their own research. The book is designed to be a practical guide that can be used in the classroom or for self-study. It is hoped that this book will help students to develop the skills they need to conduct their own research and to understand the research of others.
# Part I

## Overview and basics of the R programming language
<figure><img src="image_3.png" /></figure>
## -- Conflicts ------------------------------------ tidyverse_conflicts() --

## x dplyr::filter() masks stats::filter()

## x dplyr::lag() masks stats::lag()

## i Use the conflicted package (<http://conflicted.r-lib.org/>) to force all conflict

## Loading required package: languageR

##

## Loading required package: phonR

##

## Loading required package: summarytools

##

##

## Attaching package: 'summarytools'

##

##

## The following object is masked from 'package:tibble':

##

## view

##

##

## Loading required package: ggdist

##

## Loading required package: gtsummary

##

## Loading required package: cardx

## 1.2 Intro R Markdown

### 1.2.1 General

This is an R Markdown Notebook. When you execute code within the notebook, the results appear beneath the code.

Try executing this chunk by clicking the Run button within the chunk or by placing your cursor inside it and pressing `Ctrl+Shift+Enter`.

```r
plot(cars)
```
<figure><img src="image_4.png" /></figure>

Add a new chunk by clicking the *Insert Chunk* button on the toolbar or by pressing `Ctrl+Alt+I`.

When you save the notebook, an HTML file containing the code and output will be saved alongside it (click the *Preview* button or press `Ctrl+Shift+K` to preview the HTML file).

The preview shows you a rendered HTML copy of the contents of the editor. Consequently, unlike *Knit*, *Preview* does not run any R code chunks. Instead, the output of the chunk when it was last run in the editor is displayed.

### 1.2.2 Knitting to other formats

You can knit notebook into PDF or Word. Make sure to install the `tinytex` package using the following code: `install.packages("tinytex")`. If you already use `LATEX` with another distribution, then you do not need to install `tinytex` (though this works perfectly out of the bag). Then from the menu above (Knit), then choose *Knit to PDF* or *Knit to Word*. This is an excellent way to get your work fully written within RStudio: You will write a text (like this), have a specific structure and also have all your results in one place. This allows for transparency and replicability.

You may get an error message related to the default CRAN repository not being set up correctly. In this case, you will need to do the following:

1. Run the following code line: `usethis::edit_r_profile()`

2. Add the following to your profile: `options(repos = "https://cloud.r-project.org/")`

3. Restart your R distribution and then click on knit on the top
1.3 R and R Studio

1.3.1 R

R is the most influential statistical software that is widely used in data science. The R Project for Statistical Computing. R allows the user to take control of their analyses and being open about how the data were analysed, etc. R encourages transparency and reproducible research.

###Downloading base R

If you are a windows user, download the latest version here R version 4.3.1. If you are a MacX user, download the latest version here R version 4.3.1. Other Linux versions available here.

Using up-to-date versions of R is important as this allows you to use the latest developments of the software. You can have a look at what is new in this latest release here.

###Upgrading your current R installation

You can download the latest version from above and update. Or you can use the package installr and upgrade to the latest available version. If the package is not installed, use this: install.packages("installr") and then run with library/installr) then type installr in the console (what? what's a console?). We'll come to this later on!

1.3.2 R Studio

R Studio is one of the mostly used free and open-source integrated development
environment for R. It allows the user to have access to various information at
the same time, e.g., the Source, the Console, the Environment and the Files,
etc. When you open R studio, and if you have installed R appropriately, then
R Studio will “talk” to R by sending it messages to execute commands.

You can set up the layout to suit your needs. I always find the following layout
best for my needs:

1. The Source pane: the file where you write your code

2. The Console where actual code is run

3. The Environment pane, which shows you all variables/datasets, with the history of executed code, etc.

4. The Files/Viewer pane, which shows you the files in the current folder, the plots, the installed packages, and help files, etc.

If you click on Tools and Global options, then Pane Layout, you can change the
order, and add/remove options from the two panes below. You will see that I
use a specific formatting as this suits me best and also have a special colour-
coding used (Theme Modern et Editor theme = Tomorrow Night Bright). Use
the theme that works best for you!!
### 1.3.3 Other options?

#### Text Editors

I use Sublime Text to run Python, Praat and write in LATEX. I use R Markdown in R to publish my code and write notebooks. I am in the process of writing my first article with R Markdown for a fully reproducible research.

There are many development environments that can be used to “talk” to R: TinnR, Visual Studio, etc...

#### R GUIs

GUIs (for Graphical User Interface) for R are available. I list below a few. However, after trying some, I found it much easier to get to code directly in R. I don't remember all codes! I use these to my advantage, by saving my code into a script and using it later on in other scripts.

Some of the GUIs are meant to make R like excel or SPSS, while others are more specialised. Here is a list of some of these GUIs...

1. RCommander is the first GUI I used (and hated!). It is the one used in Discovering Statistics using R by Andy Field. There are compatibility issues between RCommander and RStudio... Install RCommander using `install.packages("Rcmdr")`. Then using R base, run RCommander from using `library(Rcmdr)`.

2. rattle is more of use for data mining and advanced statistics (use `library(rattle)` then `rattle()` to run)

3. Deducer. For basic and advanced statistics (run with `library(Deducer)` after installation)

4. RKWard. For basic and advanced statistics. Not available on CRAN and should be downloaded and installed.

5. Etc..

You can always start by using any of the above to familiarise yourself with the code, and then move to using R fully via code. My recommendation is to start coding first thing and search for help on how to write the specific code you are after.

## 1.4 Am I ready to use R now?

Well almost. There is one thing we need to consider: telling R where is our working directory. By default R saves this to your documents (or somewhere else). Here, this is generally OK, though when working on your own data, things get more complicated.

There are two schools of thought here. 1. Create R scripts that run the analyses and saves the output(s) directly to your working directory. Does not save the `.RData` image at the end 2. Create a project: a self-contained folder, where all
your scripts, figures, etc. will be automatically saved. Saves the `.RData` at the end

I subscribe to the second, as some of the computations I run take ages to finish.

### 1.4.1 Setting working directory

Click the menu `Session -> Set Workdir Directory -> Choose Directory` or use `setwd("path/to/directory")` (choose the location where you want to save the results)

You can also use `getwd()` to know where is your current working directory.

### 1.4.2 Creating a project

Look at the top-right hand where you can see **Projects (none)**. You can create a new project in a new path or based on a specific folder.

## 1.5 How to use packages?

Base R comes with many packages already installed. Look at **packages** to see which ones are already installed. There are currently 23067 packages on Cran: see here (to check all packages and their status) (repository for all packages). No one uses all packages so do not try to install all of them. Simply install what you need!! RMarkdown will let you know if you are running a specific code that lacks a package and asks you to download it.

### 1.5.1 Installation

The best option is to use the menu above (under Tools) and click **Install packages**, or type in `install.packages("package.name")`. Make sure to always have `install dependencies` ticked (using the first option).

### 1.5.2 Loading

Use the following to load a package: `library/package.name)`. Once the package is loaded, you can use any of its functions directly into your code. Sometimes you may need to specify to use a particular function from within a particular package, in this case use: `package.name::function`. We will most probably not use this today, but this is something you need to know about otherwise undesirable results may occur (or even errors!).

### 1.5.3 Finding packages and help

Under the Files pane (right bottom), click on the menu **Packages** and you will have access to all **installed** packages. Click on a package and you can see the as-
sociated help files. You can also type the following to find help: ?package.name.
??function e.g.,

?stats

## starting httpd help server ... done

??MASS

Or try clicking on the function name to find details of what to specify: e.g.,
scroll on lmer (assuming lme4 is installed). Do a Ctrl/Cmd + left mouse click
on a function to display options.

lme4::lmer

## function (formula, data = NULL, REML = TRUE, control = lmerControl(),

## start = NULL, verbose = 0L, subset, weights, na.action, offset,

## contrasts = NULL, devFunOnly = FALSE)

## {

## mc <- mcout <- match.call()

## missCtrl <- missing(control)

## if (!missCtrl && !inherits(control, "lmerControl")) {

## if (!is.list(control))

## stop("'control' is not a list; use lmerControl()")

## warning("passing control as list is deprecated: please use lmerControl() instead",

## immediate. = TRUE)

## control <- do.call(lmerControl, control)

## }

## mc$control <- control

## mc[[1]] <- quote(lme4::lFormula)

## lmod <- eval(mc, parent.frame(1L))

## mcout$formula <- lmod$formula

## lmod$formula <- NULL

## if (is.matrix(y <- model.response(lmod$fr)) && ncol(y) >

## 1) {

## stop("can't handle matrix-valued responses: consider using refit()")

## }

## devfun <- do.call(mkLmerDevfun, c(lmod, list(start = start,

## verbose = verbose, control = control)))

## if (devFunOnly)

## return(devfun)

## if (identical(control$optimizer, "none"))

## stop("deprecated use of optimizer==’none’; use NULL instead")

## opt <- if (length(control$optimizer) == 0) {

## s <- getStart(start, environment(devfun)$pp)

## list(par = s, fval = devfun(s), conv = 1000, message = "no optimization")

## }

## else {
18CHAPTER 1. INTRODUCTION TO R, THE TIDYVERSE - VISUALISATION

```R
## optimizeLmer(devfun, optimizer = control$optimizer, restart_edge = control$
## boundary tol = control$boundary tol, control = control$optCtrl,
## verbose = verbose, start = start, calc.derivs = control$calc.derivs,
## use.last.params = control$use.last.params)
## }
## cc <- checkConv(attr opt, "derivs"), opt$par, ctrl = control$checkConv,
## lbound = environment(devfun)$lower)
## mkMerMod environment(devfun), opt, lmod$reTrms, fr = lmod$fr,
## mc = mcout, lme4conv = cc)
## }
## <bytecode: 0x000002c961f2e0b0>
## <environment: namespace:lme4>
```

1.6 Let's get started with R

1.6.1 R as a calculator

1.6.1.1 Simple calculations

R can be used as a calculator. Try some of the following below:

1 + 2

## [1] 3

1+2*3

## [1] 7

Well wait a second! were you all expecting the result to be 7? how many
expected the result to be 9?

Check the following:

(1+2)*3

## [1] 9

1+(2*3)

## [1] 7

So parenthesis are important! Always use these to tell R (and any other software)
the order of operations. This is the order (remember PEMDAS):

1. Parentheses

2. Exponents

3. Multiplication and Division (from left to right)

4. Addition and Subtraction (from left to right)
### 1.6.1.2 Functions

There are many built-in functions in R to do some complicated mathematical calculations.

#### 1.6.1.2.1 Basic functions Run some of the following.

```
sqrt(3)
```

```
## [1] 1.732051
```

```
3^2
```

```
## [1] 9
```

```
log(3)
```

```
## [1] 1.098612
```

```
exp(3)
```

```
## [1] 20.08554
```

#### 1.6.1.2.2 Creating variables We can also create variables (aka temporary place holders).

```
x <- 2
```

```
y <- 5
```

```
b <- x*y
```

```
x
```

```
## [1] 2
```

```
y
```

```
## [1] 5
```

```
b
```

```
## [1] 10
```

```
b+log(y)*x^2
```

```
## [1] 16.43775
```

When you create a variable and assign to it a number (or characters), you can use it later on.

#### 1.6.1.2.3 Sequences We can also create sequences of numbers

```
seq(1, 10, 2)
```

```
## [1] 1 3 5 7 9
```
?seq

z <- 1:10

And we can do the following.. Can you explain what we have done here?

z2 <- z+1
z*z2

## [1] 2 6 12 20 30 42 56 72 90 110

Up to you... Write some more complex maths here just for fun!

## Add below

1.6.2 Objects

1.6.2.1 Basic objects

Objects are related to variables (we created above), but can also be dataframes,
and other things we create in R. All of these are stored in memory and are shown
below (under environment). You can check the type of the "object" below in
the list (look at "Type") or by using class().

Let's look at the variables we created so far.. We will create another one as
well...

class(b)

## [1] "numeric"

class(x)

## [1] "numeric"

class(y)

## [1] "numeric"

class(z)

## [1] "integer"

class(z2)

## [1] "numeric"

a <- "test"

class(a)

## [1] "character"
When we do calculations in R, we need to make sure we use numeric/integer variables only. Try some of the below. Uncomment the following `##x + two` and run the code.

x+y

## [1] 7

two <- "2"

##x + two

Can you explain the error?

We have tried to add a number to a (character) string which is clearly impossible. To do the maths, we need to change the class using any of the following commands: `as.character`, `as_integer`, `as_numeric`, `as_factor`, e.g.:

two <- as.numeric(two)

x + two

## [1] 4

### 1.6.2.2 Other functions and objects

#### 1.6.2.2.1 Some more calculations
We can create a vector of objects to do various things on.. We use the function `c()` and do various things on:

numbers <- c(1,4,5,12,55,13,45,38,77,836,543)

class(numbers)

## [1] "numeric"

mean(numbers)

## [1] 148.0909

sd(numbers)

## [1] 276.6375

median(numbers)

## [1] 38

min(numbers)

## [1] 1

max(numbers)

## [1] 836

range(numbers)

## [1] 1 836
sum(numbers)

## [1] 1629

1.6.2.2.2 Referring to a specific position Sometimes we may want to refer to a specific position in the list of numbers we just created... Use the following:

numbers[2]

## [1] 4

numbers[3:5]

## [1] 5 12 55

numbers[-4]

## [1] 1 4 5 55 13 45 38 77 836 543

numbers+numbers[6]

## [1] 14 17 18 25 68 26 58 51 90 849 556

Can you explain what we have done in the last operation?

1.7 Matrices and dataframes

1.7.1 Matrix

1.7.1.1 General

x <- 1:4

x <- as.matrix(x)

x

## [,1]

## [1,] 1

## [2,] 2

## [3,] 3

## [4,] 4

dim(x)

## [1] 4 1

dim(x) <- c(2,2)

dim(x)

## [1] 2 2
x

## [,1] [,2]
## [1,]   1    3
## [2,]   2    4

1.7.1.2 Referring to specific location

x[1,]

## [1] 1 3
x[,1]

## [1] 1 2
x[1,2]

## [1] 3
x[,] ## = x

## [,1] [,2]
## [1,]   1    3
## [2,]   2    4

1.7.2 Dataframes

A dataframe is the most important object we will be using over and over again...
It is an object that contains information in both rows and columns.

1.7.2.1 Creating a dataframe from scratch

In this exercise, we will create a 4*9 dataframe. The code below creates four
variables, and combines them together to make a dataframe. As you can see,
variables can also be characters. To create the dataframe, we use the functions
as.data.frame and cbind.

word <- c("a", "the", "lamp", "not", "jump", "it", "coffee", "walk", "on")
freq <- c(500, 600, 7, 200, 30, 450, 130, 33, 300) ## note this is completely made up!!
functionword <- c("y", "y", "n", "y", "n", "y", "n", "n", "y")
length <- c(1, 3, 4, 3, 4, 2, 6, 4, 2)
df <- as.data.frame(cbind(word,freq,functionword,length))

1.7.2.2 Deleting variables from the Environment

If you have created various variables you do not need any more, you can use rm
to remove these
rm(word,freq,functionword,length)

BUT wait, did I remove these from my dataframe? Well no.. We have removed
objects from within the R environment and not from the actual dataframe. Let's
check this up

df

| ## | word | freq | functionword | length |
|:---|:---|:---|:---|:---|
| ## 1 | a | 500 | y | 1 |
| ## 2 | the | 600 | y | 3 |
| ## 3 | lamp | 7 | n | 4 |
| ## 4 | not | 200 | y | 3 |
| ## 5 | jump | 30 | n | 4 |
| ## 6 | it | 450 | y | 2 |
| ## 7 | coffee | 130 | n | 6 |
| ## 8 | walk | 33 | n | 4 |
| ## 9 | on | 300 | y | 2 |

1.7.2.3 Saving and reading the dataframe

1.7.2.4 Reading and Saving in .csv

The code below allows you to save the dataframe and read it again. The exten-
sion .csv is for "comma delimited files". This is the best format to use as it is
simply a text file with no additional formatting.

write.csv(df, paste0("outputs/df.csv"))
dfNew <- read.csv(paste0("outputs/df.csv"))
df

| ## | word | freq | functionword | length |
|:---|:---|:---|:---|:---|
| ## 1 | a | 500 | y | 1 |
| ## 2 | the | 600 | y | 3 |
| ## 3 | lamp | 7 | n | 4 |
| ## 4 | not | 200 | y | 3 |
| ## 5 | jump | 30 | n | 4 |
| ## 6 | it | 450 | y | 2 |
| ## 7 | coffee | 130 | n | 6 |
| ## 8 | walk | 33 | n | 4 |
| ## 9 | on | 300 | y | 2 |

dfNew

| ## | X | word | freq | functionword | length |
|:---|:---|:---|:---|:---|:---|
| ## 1 | 1 | a | 500 | y | 1 |
| ## 2 | 2 | the | 600 | y | 3 |
| ## 3 | 3 | lamp | 7 | n | 4 |
| ## 4 | 4 | not | 200 | y | 3 |
<table><tr><td>## 5 5</td><td>jump</td><td>30</td><td>n</td><td>4</td></tr><tr><td>## 6 6</td><td>it</td><td>450</td><td>y</td><td>2</td></tr><tr><td>## 7 7</td><td>coffee</td><td>130</td><td>n</td><td>6</td></tr><tr><td>## 8 8</td><td>walk</td><td>33</td><td>n</td><td>4</td></tr><tr><td>## 9 9</td><td>on</td><td>300</td><td>y</td><td>2</td></tr></table>

The newly created object contains 5 columns rather than the 4 we initially created. This is normal. By default, R add a column that reflects the order of the list *before* it was saved. You can simply delete the column or keep as is (but be careful as this means you need to adjust any references to columns that we will use later on).

### 1.7.2.5 Reading and saving other formats

R allows us to read data in any format. If you have a `.txt`, `.sav`, `.xls`, `.xlsx`, etc., then there are packages specific to do that (e.g., package `xlsex` to read/save `.xlsex` files, or the function `haven` from the package `Tidyverse` to read/save `.sav` files).

You can use the built-in plugin in RStudio to **import** your dataset. See *Import Dataset* within the *Environment*.

In general, any specific formatting is kept, but sometimes variable names associated with numbers (as in `.sav` files) will be lost. Hence, it is always preferable to do minimal formatting on the data.. Start with a `.csv` file, import it to R and do the magic!

### 1.7.2.6 Checking the structure

The first thing we will do is to check the structure of our created dataset. We will use the originally created one (i.e., `df` and not the imported one (i.e., `dfNew`).

str(df)

```r
## 'data.frame': 9 obs. of 4 variables:
## $ word       : chr  "a" "the" "lamp" "not" ...
## $ freq       : chr  "500" "600" "7" "200" ...
## $ functionword: chr  "y" "y" "n" "y" ...
## $ length     : chr  "1" "3" "4" "3" ...
```

The function `str` gives us the following information:

1. How many observations (i.e., rows) and variables (i.e., columns)

2. The name of each variable (look at `$` and what comes after it)

3. Within each variable, we have the class with number of levels

### 1.7.2.7 Changing the class of a variable

As we can see, the four created variables were added to the dataframe as `factors`. We need to change the `class` of the `numeric` variables: `freq` and `length`. Let's
do that:

```
df$freq <- as.numeric(df$freq)
df$length <- as.numeric(df$length)
str(df)
```

```
## 'data.frame': 9 obs. of 4 variables:
## $ word      : chr  "a" "the" "lamp" "not" ...
## $ freq      : num  500 600 7 200 30 450 130 33 300
## $ functionword: chr  "y" "y" "n" "y" ...
## $ length    : num  1 3 4 3 4 2 6 4 2
```

### 1.7.2.8 Referring to particular variables, observations

As you can see from the above, we can refer to a particular variable in the dataframe by its name and adding $. There are additional options to do that. Let's see what we can do. Can you tell what each of the below does? chat to your neighbour....

```
df[1]
##      word
## 1      a
## 2      the
## 3      lamp
## 4      not
## 5      jump
## 6      it
## 7      coffee
## 8      walk
## 9      on
```

```
df[,1]
```

```
## [1] "a" "the" "lamp" "not" "jump" "it" "coffee" "walk"
## [9] "on"
```

```
df[1,]
```

```
##      word freq functionword length
## 1      a   500          y      1
```

```
df[1,1]
```

```
## [1] "a"
```

Here are the answers:

1. Refers to the full column 1

2. Refers to first variable

3. Refers to first row
4. Refers to first observation in first column

Practice a bit and use other specifications to obtain specific observations, columns or rows...

1.7.3 Descriptive statistics

1.7.3.1 Basic summaries, tables

We can use the function summary to do some basic summaries

summary(df)

<table><thead><tr><th>##</th><th>word</th><th>freq</th><th>functionword</th><th>length</th></tr></thead><tbody><tr><td>##</td><td>Length:9</td><td>Min. : 7</td><td>Length:9</td><td>Min. :1.000</td></tr><tr><td>##</td><td>Class :character</td><td>1st Qu.: 33</td><td>Class :character</td><td>1st Qu.:2.000</td></tr><tr><td>##</td><td>Mode :character</td><td>Median :200</td><td>Mode :character</td><td>Median :3.000</td></tr><tr><td>##</td><td></td><td>Mean :250</td><td></td><td>Mean :3.222</td></tr><tr><td>##</td><td></td><td>3rd Qu.:450</td><td></td><td>3rd Qu.:4.000</td></tr><tr><td>##</td><td></td><td>Max. :600</td><td></td><td>Max. :6.000</td></tr></tbody></table>

We can create a table with the function table

table(df$functionword, df$freq)

<table><tr><td>##</td><td>7</td><td>30</td><td>33</td><td>130</td><td>200</td><td>300</td><td>450</td><td>500</td><td>600</td></tr><tr><td>##</td><td>n</td><td>1</td><td>1</td><td>1</td><td>1</td><td>0</td><td>0</td><td>0</td><td>0</td></tr><tr><td>##</td><td>y</td><td>0</td><td>0</td><td>0</td><td>0</td><td>1</td><td>1</td><td>1</td><td>1</td></tr></table>

1.7.3.2 Basic manipulations

1.7.3.2.1 Creating variables We sometimes need to create and/or delete new variables.. Do you know how to do that?

Let's look at the structure again:

str(df)

<table><tr><td>'data.frame'</td><td>9 obs. of 4 variables:</td></tr><tr><td>$ word</td><td>chr "a" "the" "lamp" "not" ...</td></tr><tr><td>$ freq</td><td>num 500 600 7 200 30 450 130 33 300</td></tr><tr><td>$ functionword</td><td>chr "y" "y" "n" "y" ...</td></tr><tr><td>$ length</td><td>num 1 3 4 3 4 2 6 4 2</td></tr></table>

We said earlier that we can refer to a specific variable by using `$` + the name of the variable. Let's use this again and add a new name of variable not in the list of variables above

df newVariable

## NULL
What does NULL mean? The variable does not exist! Let's do something else

```
df$newVariable <- NA
```

Ah no error messages! Let's check the structure

```
str(df)
```

```
## 'data.frame':  9 obs. of 5 variables:
## $ word         : chr  "a" "the" "lamp" "not" ...
## $ freq         : num  500 600 7 200 30 450 130 33 300
## $ functionword: chr  "y" "y" "n" "y" ...
## $ length       : num  1 3 4 3 4 2 6 4 2
## $ newVariable  : logi  NA NA NA NA NA NA ...
```

So we now have five variables and the last one is named "newVariable" and assigned "NA". "NA" is used in R to refer to missing data or is a place holder. We can replace these with any calculations, or anything else. Let's do that:

```
df$newVariable <- log(df$freq)
str(df)
```

```
## 'data.frame':  9 obs. of 5 variables:
## $ word         : chr  "a" "the" "lamp" "not" ...
## $ freq         : num  500 600 7 200 30 450 130 33 300
## $ functionword: chr  "y" "y" "n" "y" ...
## $ length       : num  1 3 4 3 4 2 6 4 2
## $ newVariable  : num  6.21 6.4 1.95 5.3 3.4 ...
```

We replaced "NA" with the log of the frequencies. Let's check that this is correct only for one observation. Can you dissect the code below? what did I use to ask R to compute the log of the frequency (freq)? Remember rows and columns

```
log(df[1,2])
```

```
## [1] 6.214608
```

```
df[1,5]
```

```
## [1] 6.214608
```

So they are the same values.

### 1.7.3.2.2 Changing column names
Now we need to change the name of the variable to reflect the computations. "newVariable" is meaningless as a name, but "logFreq" is informative.

```
colnames(df)[5] <- "logFreq"
str(df)
```

```
## 'data.frame':  9 obs. of 5 variables:
## $ word         : chr  "a" "the" "lamp" "not" ...
```
```r
## $ freq      : num  500 600 7 200 30 450 130 33 300
## $ functionword: chr  "y" "y" "n" "y" ...
## $ length     : num  1 3 4 3 4 2 6 4 2
## $ logFreq    : num  6.21 6.4 1.95 5.3 3.4 ...
```

As can be seen from the above, using the command `colnames(df)[5] <- "logFreq"` allows us to change the column name in position 5 of the dataframe. If we were to change all of the columns names, we could use `colnames(df) <- c("col1","col2",...)`.

### 1.7.3.2.3 Deleting variables
Let us now create a new compound variable that we later delete. This new compound variable will the multiplication of two numeric variables. The result is meaningless of course, but will be used for this exercise.

```r
df$madeUpVariable <- df$freq*df$length
str(df)

## 'data.frame': 9 obs. of 6 variables:
## $ word      : chr  "a" "the" "lamp" "not"
## $ freq      : num  500 600 7 200 30 450 130 33 300
## $ functionword: chr  "y" "y" "n" "y" ...
## $ length     : num  1 3 4 3 4 2 6 4 2
## $ logFreq    : num  6.21 6.4 1.95 5.3 3.4 ...
## $ madeUpVariable: num  500 1800 28 600 120 900 780 132 600
```

Let us now delete this variable given that we are not interested in. Do you know how to do that? Think about how we referred to a variable before? We use `df[colNumber]`. What if we use `df[-colNumebr]`, what would be the result?

```
df[-6]
```

```
## word freq functionword length logFreq
## 1    a   500          y          1 6.214608
## 2    the  600          y          3 6.396930
## 3    lamp  7            n          4 1.945910
## 4    not  200          y          3 5.298317
## 5    jump  30           n          4 3.401197
## 6    it   450          y          2 6.109248
## 7    coffee 130         n          6 4.867534
## 8    walk  33           n          4 3.496508
## 9    on   300          y          2 5.703782
```

This shows all columns *minus* the one we are not interested in. If we rewrite the variable `df` and assign to it the newly created dataframe we just used above (with the minus sign), then the column we are not interested in will be deleted.

```
df <- df[-6]
str(df)
```
30CHAPTER 1. INTRODUCTION TO R, THE TIDYVERSE - VISUALISATION

```r
## 'data.frame': 9 obs. of 5 variables:
## $ word : chr "a" "the" "lamp" "not" ...
## $ freq : num 500 600 7 200 30 450 130 33 300
## $ functionword: chr "y" "y" "n" "y" ...
## $ length : num 1 3 4 3 4 2 6 4 2
## $ logFreq : num 6.21 6.4 1.95 5.3 3.4 ...
```

### 1.7.3.2.4 Changing names of observations
Let's say that we want to change the names of our observations. For instance, the variable "functionword" has the levels "y" and "n". Let us change the names to become "yes" and "no". We first need to change the factor level variable into character and then change the observations. Then we need to transform back to a factor

```r
df$functionword <- as.character(df$functionword)
df$functionword[df$functionword == "y"] <- "yes"
df$functionword[df$functionword == "n"] <- "no"
df$functionword <- as.factor(df$functionword)
str(df)
```

```r
## 'data.frame': 9 obs. of 5 variables:
## $ word : chr "a" "the" "lamp" "not" ...
## $ freq : num 500 600 7 200 30 450 130 33 300
## $ functionword: Factor w/ 2 levels "no","yes": 2 2 1 2 1 2 1 1 2
## $ length : num 1 3 4 3 4 2 6 4 2
## $ logFreq : num 6.21 6.4 1.95 5.3 3.4 ...
```

### 1.7.3.2.5 Checking levels of factors
We can also check the levels of factor and change the reference value. This is useful when doing any type of statistics or when plotting the data. We use `levels`, `relevel` and `ref`

```r
levels(df$functionword)

## [1] "no" "yes"

df$functionword <-relevel(df$functionword, ref = "yes")
levels(df$functionword)
```

```r
## [1] "yes" "no"
```

We can also use the following code to change the order of the levels of a multilevel factor

```r
levels(df$word)
```

```r
## NULL

df$word <- factor(df$word, levels = c("a", "coffee", "jump", "lamp", "not", "it", "on", "walk"))
levels(df$word)
```

```r
## [1] "a" "coffee" "jump" "lamp" "not" "it" "on" "walk"
```
## [9] "the"

1.7.3.2.6 Subsetting the dataframe We may sometimes need to subset
the dataframe and use parts of it. We use the function subset or which.

df_Yes1 <- df[which(df$functionword == 'yes'),]

##or

df_Yes2 <- subset(df, functionword="yes")

str(df_Yes1)

## 'data.frame': 5 obs. of 5 variables:

## $ word : Factor w/ 9 levels "a", "coffee", "jump", ...: 1 9 5 6 7

## $ freq : num 500 600 200 450 300

## $ functionword: Factor w/ 2 levels "yes", "no": 1 1 1 1 1

## $ length : num 1 3 3 2 2

## $ logFreq : num 6.21 6.4 5.3 6.11 5.7

str(df_Yes2)

## 'data.frame': 5 obs. of 5 variables:

## $ word : Factor w/ 9 levels "a", "coffee", "jump", ...: 1 9 5 6 7

## $ freq : num 500 600 200 450 300

## $ functionword: Factor w/ 2 levels "yes", "no": 1 1 1 1 1

## $ length : num 1 3 3 2 2

## $ logFreq : num 6.21 6.4 5.3 6.11 5.7

When we subset the data, the levels of a factor are kept as they are.

levels(df_Yes1$functionword)

## [1] "yes" "no"

levels(df_Yes2$functionword)

## [1] "yes" "no"

But we only have one level of our factor..

df_Yes1$functionword

## [1] yes yes yes yes yes

## Levels: yes no

df_Yes2$functionword

## [1] yes yes yes yes yes

## Levels: yes no

By default, R keeps the levels of the factor as they are unless we change it by
using the following:
```
df_Yes1$functionword <- factor(df_Yes1$functionword)
df_Yes2$functionword <- factor(df_Yes2$functionword)
df_Yes1$functionword
```

```
## [1] yes yes yes yes yes
## Levels: yes
df_Yes2$functionword
```

```
## [1] yes yes yes yes yes
## Levels: yes
```

## 1.8 The Tidyverse

### 1.8.1 Introduction

The Tidyverse is a family of packages used to speed up the use of R.

#### The tidyverse

##### Components

<figure><img src="image_5.png" /></figure>

The tidyverse is a collection of R packages that share common philosophies and are designed to work together. This site is a work-in-progress guide to the tidyverse and its packages.

You need to first install it (if you haven't already done so) and then load it. To install, use `Tools > Install packages` or `install.packages()` then add `tidyverse` to your `library()`.
Look at how many packages are installed within the *Tidyverse*. The messages you see are telling you which packages are loaded and which functions are in conflict (i.e., these are functions from other packages that are found within the *Tidyverse*). If you want to use the original function, simply add `package_name::function`.

### 1.8.2 Using piping

The difference between base R and the Tidyverse's way of doing things is that base R can sometimes be more complex, while tidyverse is more straightforward and allows you to "see" within a dataframe easily. You need to learn how to use the "pipe" in `magrittr` that is part of the *Tidyverse*.

<figure><img src="image_6.png" /></figure>

Pipes are written in R as `>%` (note you must use a percentage sign before and after the pipe). To demonstrate what pipes do, have a look at the following pseudo code. You can use a shortcut in your keyboard, type `Ctrl+Shift+m` to add a `pipe` (for mac users, it is `Cmd+Shift+m`).

```
data_frame <- function(data_frame, another_argument)

data_frame %>% function(another_argument)
```

Since R version 4.1.0, there is a native `pipe` `|>`. It seems to be doing almost
the same thing as the `%>%`. We will still use `%>%` as this is integrated within the `Tidyverse`.

### 1.8.3 Demo subseting

Below are two code lines for how to subset the dataframe using base R and piping from the `magrittr` package.

With base R, we always need to refer to the dataset twice: once at the beginning and then to look into the dataset to select a variable.

```r
df_Yes1 <- df[which(df$functionword == 'yes'),]
df_Yes1
```

<table><thead><tr><th>##</th><th>word</th><th>freq</th><th>functionword</th><th>length</th><th>logFreq</th></tr></thead><tbody><tr><td>## 1</td><td>a</td><td>500</td><td>yes</td><td>1</td><td>6.214608</td></tr><tr><td>## 2</td><td>the</td><td>600</td><td>yes</td><td>3</td><td>6.396930</td></tr><tr><td>## 4</td><td>not</td><td>200</td><td>yes</td><td>3</td><td>5.298317</td></tr><tr><td>## 6</td><td>it</td><td>450</td><td>yes</td><td>2</td><td>6.109248</td></tr><tr><td>## 9</td><td>on</td><td>300</td><td>yes</td><td>2</td><td>5.703782</td></tr></tbody></table>

With the pipe, you only need to specify the dataset once: By adding the pipe, you can already look into the dataset and select the variable you need.

```r
df %>% filter(functionword == 'yes')
```

<table><thead><tr><th>##</th><th>word</th><th>freq</th><th>functionword</th><th>length</th><th>logFreq</th></tr></thead><tbody><tr><td>## 1</td><td>a</td><td>500</td><td>yes</td><td>1</td><td>6.214608</td></tr><tr><td>## 2</td><td>the</td><td>600</td><td>yes</td><td>3</td><td>6.396930</td></tr><tr><td>## 3</td><td>not</td><td>200</td><td>yes</td><td>3</td><td>5.298317</td></tr><tr><td>## 4</td><td>it</td><td>450</td><td>yes</td><td>2</td><td>6.109248</td></tr><tr><td>## 5</td><td>on</td><td>300</td><td>yes</td><td>2</td><td>5.703782</td></tr></tbody></table>

And this is with the base R pipe (combined with code from the `Tidyverse` family)

```r
df |> filter(functionword == 'yes')
```

<table><thead><tr><th>##</th><th>word</th><th>freq</th><th>functionword</th><th>length</th><th>logFreq</th></tr></thead><tbody><tr><td>## 1</td><td>a</td><td>500</td><td>yes</td><td>1</td><td>6.214608</td></tr><tr><td>## 2</td><td>the</td><td>600</td><td>yes</td><td>3</td><td>6.396930</td></tr><tr><td>## 3</td><td>not</td><td>200</td><td>yes</td><td>3</td><td>5.298317</td></tr><tr><td>## 4</td><td>it</td><td>450</td><td>yes</td><td>2</td><td>6.109248</td></tr><tr><td>## 5</td><td>on</td><td>300</td><td>yes</td><td>2</td><td>5.703782</td></tr></tbody></table>

As you can see, using the pipe (either within the `Tidyverse` or with base R) is a quick and easy way to do various operations.

Out of convenience and because we will use other packages integrated within the `Tidyverse`, we will use its pipe.
ReCap:

* `>%` is called a “pipe”

* It passes the previous line into the `data` argument of the next line

* It does not save any changes after output

* If you want to save the output of a particular manipulation, simply save it with `xx < -`

### 1.8.4 Basic manipulations

We will use the pipe with the `Tidyverse` to obtain summaries. We will use an R built-in dataset. Type `data()` to see the full list of datasets installed by default in R. You can use `data/package = .packages(all.available = TRUE))` to see all datasets installed within all packages.

#### 1.8.4.1 First steps

Here is a list of all available datasets

```
data()

data/package = .packages(all.available = TRUE))
```

#### 1.8.4.2 Loading dataset

We will use the dataset `english` from the package `languageR`. This is a package that contains many linguistically-oriented datasets. See details of the dataset here. Or by typing `?languageR::english` (or simply `?english` if the package is already loaded) in the console.

You can load the dataset after loading the package. Simply refer to it by its name.

`?english`

#### 1.8.4.3 View

To see the dataset, run the code below to visualise it.

```
english %>%
  View()
```

#### 1.8.4.4 Structure

We can use `str()` to look at the structure of the dataset. Here we have a relatively large dataset with 4568 observations (`=rows`) and 36 variables (`=columns`).
english %>

str()

## 'data.frame': 4568 obs. of 36 variables:

<table><tr><td>$ RTlexdec</td><td>: num 6.54 6.4 6.3 6.42 6.45 ...</td></tr><tr><td>$ RTnaming</td><td>: num 6.15 6.25 6.14 6.13 6.2 ...</td></tr><tr><td>$ Familiarity</td><td>: num 2.37 4.43 5.6 3.87 3.93 3.27 3.73 5.67 3.1</td></tr><tr><td>$ Word</td><td>: Factor w/ 2197 levels "ace", "act", "add", ...: 467</td></tr><tr><td>$ AgeSubject</td><td>: Factor w/ 2 levels "old", "young": 2 2 2 2 2 2 2</td></tr><tr><td>$ WordCategory</td><td>: Factor w/ 2 levels "N", "V": 1 1 1 1 1 1 1 1 1 1</td></tr><tr><td>$ WrittenFrequency</td><td>: num 3.91 4.52 6.51 5.02 4.89 ...</td></tr><tr><td>$ WrittenSpokenFrequencyRatio</td><td>: num 1.022 0.35 2.089 -0.526 -1.045 ...</td></tr><tr><td>$ FamilySize</td><td>: num 1.39 1.39 1.61 1.95 2.2 ...</td></tr><tr><td>$ DerivationalEntropy</td><td>: num 0.141 0.427 0.062 0.43 0.359 ...</td></tr><tr><td>$ InflectionalEntropy</td><td>: num 0.0211 0.942 1.4434 0 1.7539 ...</td></tr><tr><td>$ NumberSimplexSynsets</td><td>: num 0.693 1.099 2.485 1.099 2.485 ...</td></tr><tr><td>$ NumberComplexSynsets</td><td>: num 0 0 1.95 2.64 2.48 ...</td></tr><tr><td>$ LengthInLetters</td><td>: int 3 5 6 4 4 4 3 3 5 ...</td></tr><tr><td>$ Ncount</td><td>: int 8 5 0 8 3 9 6 13 3 3 ...</td></tr><tr><td>$ MeanBigramFrequency</td><td>: num 7.04 9.54 9.88 8.31 7.94 ...</td></tr><tr><td>$ FrequencyInitialDiphone</td><td>: num 12 12.6 13.3 12.1 11.9 ...</td></tr><tr><td>$ ConspelV</td><td>: int 10 20 10 5 17 19 10 13 1 7 ...</td></tr><tr><td>$ ConspelN</td><td>: num 3.74 7.87 6.69 6.68 4.76 ...</td></tr><tr><td>$ ConphonV</td><td>: int 41 38 13 6 17 21 13 7 11 14 ...</td></tr><tr><td>$ ConphonN</td><td>: num 8.84 9.78 7.04 3.83 4.76 ...</td></tr><tr><td>$ ConfriendsV</td><td>: int 8 20 10 4 17 19 10 6 0 7 ...</td></tr><tr><td>$ ConfriendsN</td><td>: num 3.3 7.87 6.69 3.53 4.76 ...</td></tr><tr><td>$ ConffV</td><td>: num 0.693 0 0 0.693 0 ...</td></tr><tr><td>$ ConffN</td><td>: num 2.71 0 0 6.63 0 ...</td></tr><tr><td>$ ConfbV</td><td>: num 3.5 2.94 1.39 1.1 0 ...</td></tr><tr><td>$ ConfbN</td><td>: num 8.83 9.61 5.82 2.56 0 ...</td></tr><tr><td>$ NounFrequency</td><td>: int 49 142 565 150 170 125 582 2061 144 522 ...</td></tr><tr><td>$ VerbFrequency</td><td>: int 0 0 473 0 120 280 110 76 4 86 ...</td></tr><tr><td>$ CV</td><td>: Factor w/ 2 levels "C", "V": 1 1 1 1 1 1 1 1 1 1</td></tr><tr><td>$ Obstruent</td><td>: Factor w/ 2 levels "cont", "obst": 2 2 2 2 2 2 2</td></tr><tr><td>$ Frication</td><td>: Factor w/ 4 levels "burst", "frication", ...: 1 2 2</td></tr><tr><td>$ Voice</td><td>: Factor w/ 2 levels "voiced", "voiceless": 1 2 2</td></tr><tr><td>$ FrequencyInitialDiphoneWord</td><td>: num 10.13 9.05 12.42 10.05 11.8 ...</td></tr><tr><td>$ FrequencyInitialDiphoneSyllable</td><td>: num 10.41 9.15 13.13 11 12.16 ...</td></tr><tr><td>$ CorrectLexdec</td><td>: int 27 30 30 30 26 28 30 28 25 29 ...</td></tr></table>

1.8.4.5 See first 6 rows

english %>

head()
<table><thead><tr><th></th><th>RTlexdec</th><th>RTnaming</th><th>Familiarity</th><th>Word</th><th>Age</th><th>Subject</th><th>WordCategory</th><th>WrittenFrequency</th></tr></thead><tbody><tr><td>## 1</td><td>6.543754</td><td>6.145044</td><td>2.37</td><td>doe</td><td>young</td><td>N</td><td></td><td>3.912023</td></tr><tr><td>## 2</td><td>6.397596</td><td>6.246882</td><td>4.43</td><td>whore</td><td>young</td><td>N</td><td></td><td>4.521789</td></tr><tr><td>## 3</td><td>6.304942</td><td>6.143756</td><td>5.60</td><td>stress</td><td>young</td><td>N</td><td></td><td>6.505784</td></tr><tr><td>## 4</td><td>6.424221</td><td>6.131878</td><td>3.87</td><td>pork</td><td>young</td><td>N</td><td></td><td>5.017280</td></tr><tr><td>## 5</td><td>6.450597</td><td>6.198479</td><td>3.93</td><td>plug</td><td>young</td><td>N</td><td></td><td>4.890349</td></tr><tr><td>## 6</td><td>6.531970</td><td>6.167726</td><td>3.27</td><td>prop</td><td>young</td><td>N</td><td></td><td>4.770685</td></tr></tbody></table>

<table><thead><tr><th></th><th>WrittenSpokenFrequencyRatio</th><th>FamilySize</th><th>DerivationalEntropy</th></tr></thead><tbody><tr><td>## 1</td><td>1.0216512</td><td>1.386294</td><td>0.14144</td></tr><tr><td>## 2</td><td>0.3504830</td><td>1.386294</td><td>0.42706</td></tr><tr><td>## 3</td><td>2.0893560</td><td>1.609438</td><td>0.06197</td></tr><tr><td>## 4</td><td>-0.5263339</td><td>1.945910</td><td>0.43035</td></tr><tr><td>## 5</td><td>-1.0445451</td><td>2.197225</td><td>0.35920</td></tr><tr><td>## 6</td><td>0.9248014</td><td>1.386294</td><td>0.06268</td></tr></tbody></table>

<table><thead><tr><th></th><th>InflectionalEntropy</th><th>NumberSimplexSynsets</th><th>NumberComplexSynsets</th><th>LengthInLetters</th></tr></thead><tbody><tr><td>## 1</td><td>0.02114</td><td>0.6931472</td><td>0.000000</td><td>3</td></tr><tr><td>## 2</td><td>0.94198</td><td>1.0986123</td><td>0.000000</td><td>5</td></tr><tr><td>## 3</td><td>1.44339</td><td>2.4849066</td><td>1.945910</td><td>6</td></tr><tr><td>## 4</td><td>0.00000</td><td>1.0986123</td><td>2.639057</td><td>4</td></tr><tr><td>## 5</td><td>1.75393</td><td>2.4849066</td><td>2.484907</td><td>4</td></tr><tr><td>## 6</td><td>1.74730</td><td>1.6094379</td><td>1.386294</td><td>4</td></tr></tbody></table>

<table><thead><tr><th></th><th>Ncount</th><th>MeanBigramFrequency</th><th>FrequencyInitialDiphone</th><th>ConspelV</th><th>ConspelN</th><th>ConphonV</th></tr></thead><tbody><tr><td>## 1</td><td>8</td><td>7.036333</td><td>12.02268</td><td>10</td><td>3.737670</td><td>41</td></tr><tr><td>## 2</td><td>5</td><td>9.537878</td><td>12.59780</td><td>20</td><td>7.870930</td><td>38</td></tr><tr><td>## 3</td><td>0</td><td>9.883931</td><td>13.30069</td><td>10</td><td>6.693324</td><td>13</td></tr><tr><td>## 4</td><td>8</td><td>8.309180</td><td>12.07807</td><td>5</td><td>6.677083</td><td>6</td></tr><tr><td>## 5</td><td>3</td><td>7.943717</td><td>11.92678</td><td>17</td><td>4.762174</td><td>17</td></tr><tr><td>## 6</td><td>9</td><td>8.349620</td><td>12.19724</td><td>19</td><td>6.234411</td><td>21</td></tr></tbody></table>

<table><thead><tr><th></th><th>ConphonN</th><th>ConfriendsV</th><th>ConfriendsN</th><th>ConffV</th><th>ConffN</th><th>ConfbV</th><th>ConfbN</th></tr></thead><tbody><tr><td>## 1</td><td>8.837826</td><td>8</td><td>3.295837</td><td>0.6931472</td><td>2.708050</td><td>3.496508</td><td>8.833900</td></tr><tr><td>## 2</td><td>9.775825</td><td>20</td><td>7.870930</td><td>0.0000000</td><td>0.000000</td><td>2.944439</td><td>9.614738</td></tr><tr><td>## 3</td><td>7.040536</td><td>10</td><td>6.693324</td><td>0.0000000</td><td>0.000000</td><td>1.386294</td><td>5.817111</td></tr><tr><td>## 4</td><td>3.828641</td><td>4</td><td>3.526361</td><td>0.6931472</td><td>6.634633</td><td>1.098612</td><td>2.564949</td></tr><tr><td>## 5</td><td>4.762174</td><td>17</td><td>4.762174</td><td>0.0000000</td><td>0.000000</td><td>0.000000</td><td>0.000000</td></tr><tr><td>## 6</td><td>6.249975</td><td>19</td><td>6.234411</td><td>0.0000000</td><td>0.000000</td><td>1.098612</td><td>2.197225</td></tr></tbody></table>

<table><thead><tr><th></th><th>NounFrequency</th><th>VerbFrequency</th><th>CV</th><th>Obstruent</th><th>Frication</th><th>Voice</th></tr></thead><tbody><tr><td>## 1</td><td>49</td><td></td><td>0</td><td>C</td><td>obst</td><td>burst</td><td>voiced</td></tr><tr><td>## 2</td><td>142</td><td></td><td>0</td><td>C</td><td>obst</td><td>frication</td><td>voiceless</td></tr><tr><td>## 3</td><td>565</td><td></td><td>473</td><td>C</td><td>obst</td><td>frication</td><td>voiceless</td></tr><tr><td>## 4</td><td>150</td><td></td><td>0</td><td>C</td><td>obst</td><td>burst</td><td>voiceless</td></tr><tr><td>## 5</td><td>170</td><td></td><td>120</td><td>C</td><td>obst</td><td>burst</td><td>voiceless</td></tr><tr><td>## 6</td><td>125</td><td></td><td>280</td><td>C</td><td>obst</td><td>burst</td><td>voiceless</td></tr></tbody></table>

<table><thead><tr><th></th><th>FrequencyInitialDiphone</th><th>FrequencyInitialDiphone</th><th>Syllable</th><th>CorrectLexdec</th></tr></thead><tbody><tr><td>## 1</td><td>10.129308</td><td></td><td>10.409763</td><td>27</td></tr><tr><td>## 2</td><td>9.054388</td><td></td><td>9.148252</td><td>30</td></tr><tr><td>## 3</td><td>12.422026</td><td></td><td>13.127395</td><td>30</td></tr></tbody></table>
<table><tr><td>## 4</td><td>10.048151</td><td>11.003649</td><td>30</td></tr><tr><td>## 5</td><td>11.796336</td><td>12.163092</td><td>26</td></tr><tr><td>## 6</td><td>11.991567</td><td>12.436772</td><td>28</td></tr></table>

1.8.4.6 See last 6 rows

english %>%
tail()

### RTlexdec RTnaming Familiarity Word AgeSubject WordCategory

<table><tr><td>## 4563</td><td>6.608770</td><td>6.503839</td><td>3.70</td><td>spy</td><td>old</td><td>V</td></tr><tr><td>## 4564</td><td>6.753998</td><td>6.446513</td><td>2.40</td><td>jag</td><td>old</td><td>V</td></tr><tr><td>## 4565</td><td>6.711022</td><td>6.506979</td><td>3.17</td><td>hash</td><td>old</td><td>V</td></tr><tr><td>## 4566</td><td>6.592332</td><td>6.386879</td><td>3.87</td><td>dash</td><td>old</td><td>V</td></tr><tr><td>## 4567</td><td>6.565561</td><td>6.519884</td><td>4.97</td><td>flirt</td><td>old</td><td>V</td></tr><tr><td>## 4568</td><td>6.667300</td><td>6.496624</td><td>3.03</td><td>hawk</td><td>old</td><td>V</td></tr></table>

### WrittenFrequency WrittenSpokenFrequencyRatio FamilySize

<table><tr><td>## 4563</td><td>5.023881</td><td>0.9703580</td><td>1.609438</td></tr><tr><td>## 4564</td><td>2.079442</td><td>-1.6863990</td><td>1.386294</td></tr><tr><td>## 4565</td><td>3.663562</td><td>0.4367177</td><td>1.609438</td></tr><tr><td>## 4566</td><td>5.043425</td><td>0.5043947</td><td>1.945910</td></tr><tr><td>## 4567</td><td>3.135494</td><td>0.0628009</td><td>1.945910</td></tr><tr><td>## 4568</td><td>4.276666</td><td>1.0498221</td><td>1.945910</td></tr></table>

### DerivationalEntropy InflectionalEntropy NumberSimplexSynsets

<table><tr><td>## 4563</td><td>0.08753</td><td>1.64317</td><td>1.609438</td></tr><tr><td>## 4564</td><td>0.30954</td><td>1.85123</td><td>1.098612</td></tr><tr><td>## 4565</td><td>0.15110</td><td>0.77890</td><td>1.386294</td></tr><tr><td>## 4566</td><td>0.63316</td><td>1.65739</td><td>2.564949</td></tr><tr><td>## 4567</td><td>0.99953</td><td>1.75885</td><td>1.609438</td></tr><tr><td>## 4568</td><td>0.95422</td><td>1.81367</td><td>1.945910</td></tr></table>

### NumberComplexSynsets LengthInLetters Ncount MeanBigramFrequency

<table><tr><td>## 4563</td><td>0.6931472</td><td>3</td><td>5</td><td>6.838235</td></tr><tr><td>## 4564</td><td>0.0000000</td><td>3</td><td>18</td><td>6.229554</td></tr><tr><td>## 4565</td><td>1.7917595</td><td>4</td><td>11</td><td>8.825582</td></tr><tr><td>## 4566</td><td>1.6094379</td><td>4</td><td>10</td><td>8.356139</td></tr><tr><td>## 4567</td><td>0.6931472</td><td>5</td><td>1</td><td>8.751224</td></tr><tr><td>## 4568</td><td>3.0910425</td><td>4</td><td>4</td><td>7.426055</td></tr></table>

### FrequencyInitialDiphone ConspelV ConspelN ConphonV ConphonN ConfriendsV

<table><tr><td>## 4563</td><td>11.50982</td><td>18</td><td>8.917981</td><td>42</td><td>9.516132</td><td>17</td></tr><tr><td>## 4564</td><td>8.49433</td><td>20</td><td>4.744932</td><td>20</td><td>4.744932</td><td>19</td></tr><tr><td>## 4565</td><td>13.49254</td><td>25</td><td>5.141664</td><td>23</td><td>4.890349</td><td>21</td></tr><tr><td>## 4566</td><td>11.32815</td><td>25</td><td>5.141664</td><td>23</td><td>4.890349</td><td>21</td></tr><tr><td>## 4567</td><td>10.59918</td><td>7</td><td>4.624973</td><td>14</td><td>5.164786</td><td>7</td></tr><tr><td>## 4568</td><td>13.49254</td><td>3</td><td>2.772589</td><td>11</td><td>5.609472</td><td>3</td></tr></table>

### ConfriendsN ConffV ConffN ConfbV ConfbN NounFrequency

<table><tr><td>## 4563</td><td>8.916774</td><td>0.000000</td><td>0.000000</td><td>3.218876</td><td>8.7181729</td><td>219</td></tr></table>
<table><thead><tr><th></th><th></th><th></th><th></th><th></th><th></th><th></th></tr></thead><tbody><tr><td>##</td><td>4564</td><td>4.736198</td><td>0.000000</td><td>0.000000</td><td>0.000000</td><td>0.000000</td></tr><tr><td>##</td><td>4565</td><td>4.882802</td><td>1.609438</td><td>3.688879</td><td>1.098612</td><td>0.6931472</td></tr><tr><td>##</td><td>4566</td><td>4.882802</td><td>1.609438</td><td>3.688879</td><td>1.098612</td><td>0.6931472</td></tr><tr><td>##</td><td>4567</td><td>4.624973</td><td>0.000000</td><td>0.000000</td><td>2.079442</td><td>4.3040651</td></tr><tr><td>##</td><td>4568</td><td>2.772589</td><td>0.000000</td><td>0.000000</td><td>2.197225</td><td>5.5529596</td></tr><tr><td>##</td><td></td><td>VerbFrequency</td><td>CV</td><td>Obstruent</td><td>Frication</td><td>Voice FrequencyInitialDiphoneWord</td></tr><tr><td>##</td><td>4563</td><td>88</td><td>C</td><td>obst</td><td>frication</td><td>voiceless</td></tr><tr><td>##</td><td>4564</td><td>7</td><td>C</td><td>obst</td><td>frication</td><td>voiced</td></tr><tr><td>##</td><td>4565</td><td>7</td><td>C</td><td>obst</td><td>frication</td><td>voiceless</td></tr><tr><td>##</td><td>4566</td><td>231</td><td>C</td><td>obst</td><td>burst</td><td>voiced</td></tr><tr><td>##</td><td>4567</td><td>66</td><td>C</td><td>obst</td><td>frication</td><td>voiceless</td></tr><tr><td>##</td><td>4568</td><td>47</td><td>C</td><td>obst</td><td>frication</td><td>voiceless</td></tr><tr><td>##</td><td></td><td>FrequencyInitialDiphoneSyllable</td><td>CorrectLexdec</td><td></td><td></td><td></td></tr><tr><td>##</td><td>4563</td><td></td><td></td><td>12.492844</td><td></td><td>30</td></tr><tr><td>##</td><td>4564</td><td></td><td></td><td>8.390041</td><td></td><td>29</td></tr><tr><td>##</td><td>4565</td><td></td><td></td><td>12.665546</td><td></td><td>29</td></tr><tr><td>##</td><td>4566</td><td></td><td></td><td>9.287764</td><td></td><td>29</td></tr><tr><td>##</td><td>4567</td><td></td><td></td><td>10.932142</td><td></td><td>29</td></tr><tr><td>##</td><td>4568</td><td></td><td></td><td>9.148252</td><td></td><td>30</td></tr></tbody></table>

### 1.8.4.7 Selecting variables

Here, we select a few variables to use. For variables or columns, use the function `select`

```
english %>%
  select(RTlexdec, RTnaming, Familiarity) %>%
  head(10)
```

<table><thead><tr><th></th><th></th><th></th><th></th></tr></thead><tbody><tr><td>##</td><td>RTlexdec</td><td>RTnaming</td><td>Familiarity</td></tr><tr><td>##</td><td>1</td><td>6.543754</td><td>6.145044</td><td>2.37</td></tr><tr><td>##</td><td>2</td><td>6.397596</td><td>6.246882</td><td>4.43</td></tr><tr><td>##</td><td>3</td><td>6.304942</td><td>6.143756</td><td>5.60</td></tr><tr><td>##</td><td>4</td><td>6.424221</td><td>6.131878</td><td>3.87</td></tr><tr><td>##</td><td>5</td><td>6.450597</td><td>6.198479</td><td>3.93</td></tr><tr><td>##</td><td>6</td><td>6.531970</td><td>6.167726</td><td>3.27</td></tr><tr><td>##</td><td>7</td><td>6.370586</td><td>6.123808</td><td>3.73</td></tr><tr><td>##</td><td>8</td><td>6.266859</td><td>6.096050</td><td>5.67</td></tr><tr><td>##</td><td>9</td><td>6.608648</td><td>6.117657</td><td>3.10</td></tr><tr><td>##</td><td>10</td><td>6.284843</td><td>6.179188</td><td>4.43</td></tr></tbody></table>

### 1.8.4.8 Selecting observations

If we want to select observations, we use the function `filter`. We will use `select` to select particular variables and then use `filter` to select specific observations. This example shows how the pipe chain works, by combining multiple functions and using pipes
english %>
select(RTlexdec, RTnaming, Familiarity, AgeSubject) %>
filter(AgeSubject == "old") %>
head(10)

<table><thead><tr><th>RTlexdec</th><th>RTnaming</th><th>Familiarity</th><th>AgeSubject</th></tr></thead><tbody><tr><td>1453</td><td>6.664894</td><td>6.422597</td><td>2.37</td><td>old</td></tr><tr><td>1454</td><td>6.677209</td><td>6.636603</td><td>4.43</td><td>old</td></tr><tr><td>1455</td><td>6.538617</td><td>6.487075</td><td>5.60</td><td>old</td></tr><tr><td>1456</td><td>6.546943</td><td>6.404402</td><td>3.87</td><td>old</td></tr><tr><td>1457</td><td>6.637428</td><td>6.423409</td><td>3.93</td><td>old</td></tr><tr><td>1458</td><td>6.757444</td><td>6.529273</td><td>3.27</td><td>old</td></tr><tr><td>1459</td><td>6.598073</td><td>6.471728</td><td>3.73</td><td>old</td></tr><tr><td>1460</td><td>6.572464</td><td>6.424058</td><td>5.67</td><td>old</td></tr><tr><td>1461</td><td>6.817349</td><td>6.470645</td><td>3.10</td><td>old</td></tr><tr><td>1462</td><td>6.662877</td><td>6.549937</td><td>4.43</td><td>old</td></tr></tbody></table>

### 1.8.4.9 Changing order of levels

Use some of the code above to manipulate the dataframe but now using code from the *Tidyverse*. As you will see, once you know how to manipulate a dataset with base R, you can easily apply the same techniques with the *Tidyverse*. The *Tidyverse* provides additional ways to manipulate a dataframe.

For example, if I want to check levels of a variable and change the reference level, I will use the following code:

levels(english$AgeSubject)

## [1] "old" "young"

To change levels of `AgeSubject`, we need to save a new dataset (do not override the original dataset!!). The `mutate` function means we are manipulating an object.

english2<- english %>
mutate(AgeSubject = factor(AgeSubject, levels = c("young", "old")))
levels(english2$AgeSubject)

## [1] "young" "old"

### 1.8.4.10 Changing reference value

You can change the reference value by using `fct_relevel`. This is useful if you have many levels in one of the factors you are working with and you simply need to change the reference.
english2<- english %>%
  mutate(AgeSubject = fct_relevel(AgeSubject, "old"))
levels(english2$AgeSubject)

## [1] "old" "young"

## 1.8.5 String manipulation

We use the `str_func` variants from the `stringr` package within the Tidyverse.

topics <- c("Phonetics", "Phonology", "Morphology", "Syntax", "Semantics", "Pragmatics", "Psychol"

### 1.8.5.1 str_sub

Extracting substrings (2nd to 4th)

str_sub(topics, 2, 4)

## [1] "hon" "hon" "orp" "ynt" "ema" "rag" "syc"

### 1.8.5.2 str_detect

Detecting a particular pattern

str_detect(topics, "p")

## [1] FALSE FALSE TRUE FALSE FALSE FALSE
str_detect(topics, "n")

## [1] TRUE TRUE FALSE TRUE TRUE FALSE TRUE
str_detect(topics, "pho")

## [1] FALSE FALSE TRUE FALSE FALSE FALSE
str_detect(topics, "ling")

## [1] FALSE FALSE FALSE FALSE FALSE TRUE

### 1.8.5.3 str locate

Locate specific characters

str locate(topics, "n")

## start end

## [1,] 4 4

## [2,] 4 4

## [3,] NA NA

## [4,] 3 3
## [5,]      5      5
## [6,]     NA    NA
## [7,]      9      9

### 1.8.5.4 str locate all

Locate all instances of characters

str locate all ( topics , "n")

## [[1]]
##       start end
## [1,]    4      4
##
## [[2]]
##       start end
## [1,]    4      4
##
## [[3]]
##       start end
##
## [[4]]
##       start end
## [1,]    3      3
##
## [[5]]
##       start end
## [1,]    5      5
##
## [[6]]
##       start end
##
## [[7]]
##       start end
## [1,]    9      9

### 1.8.5.5 str_replace

replaces a single instance

str_replace ( topics , "o", "0")

## [1] "PhOnetics"      "PhOnology"      "MOrphology"
## [4] "Syntax"         "Semantics"       "Pragmatics"
## [7] "PsychOlinguistics"
```R
str_detect(topics, "[y]$")

## [1] FALSE TRUE TRUE FALSE FALSE FALSE FALSE

str_detect(topics, "[cs]$")

## [1] TRUE FALSE FALSE FALSE TRUE TRUE TRUE
```

### 1.8.6.2 Other characters

The `.` is used as a place holder asking for any character in the sequence

```R
str_detect(topics, "Ph.n")
```

```R
## [1] TRUE TRUE FALSE FALSE FALSE FALSE

[a-z] will detect all characters between "a" and "z"

str_detect(topics, "[a-z]")
```

```R
## [1] TRUE TRUE TRUE TRUE TRUE TRUE TRUE

[M-P] will detect all characters between "M" and "P"

str_detect(topics, "[M-P]")
```

```R
## [1] TRUE TRUE TRUE FALSE FALSE TRUE TRUE

[:lower:] will detect all characters in lower case

str_detect(topics, "[:lower:]")
```

```R
## [1] TRUE TRUE TRUE TRUE TRUE TRUE TRUE

str_detect(topics, "[:upper:]")
```

```R
## [1] TRUE TRUE TRUE TRUE TRUE TRUE TRUE
```

We can also use `[0-9]`, `[:digit:]`, `[:alpha:]`, `[:num:]`, `[:punct:]`, `[:graph:]`, `[:blank:]`, `[:space:]`, `[:print:]`, etc. See PDF of cheat sheet!

## 1.8.7 Advanced manipulations

Sometimes, you may have a dataset that comes in a wide format (i.e., columns contain data from participants) and you want to change to long format (i.e., each row contains one observation with minimal number of columns). Let's look at the functions `pivot_longer` and `pivot_wider`

### 1.8.7.1 Columns to rows

Let's use the `english` dataset to transform it from wide to long and see the first 10 rows.
english %>%
  select (Word, RTlexdec, RTnaming, Familiarity) %>%
  pivot_longer (cols = c (RTlexdec, RTnaming, Familiarity),
               names_to = "variable",
               values_to = "values") %>%
  head (10)

## # A tibble: 10 x 3
##   Word      variable      values
##   <fct>    <chr>         <db1>
##  1 doe     RTlexdec       6.54
##  2 doe     RTnaming       6.15
##  3 doe     Familiarity   2.37
##  4 whore   RTlexdec       6.40
##  5 whore   RTnaming       6.25
##  6 whore   Familiarity   4.43
##  7 stress  RTlexdec       6.30
##  8 stress  RTnaming       6.14
##  9 stress  Familiarity   5.6
## 10 pork    RTlexdec       6.42

### 1.8.7.2 Rows to columns

Let's use the same code above and change the code from long format, back to wide format. `Pivot_wider` allows you to go back to the original dataset. You will need to use `unnest` to get all rows in the correct place. Try without it to see the result.

english %>%
  select (Word, RTlexdec, RTnaming, Familiarity) %>%
  pivot_longer (cols = c (RTlexdec, RTnaming, Familiarity),
               names_to = "variable",
               values_to = "values") %>%
  pivot_wider (names_from = "variable",
               values_from = "values") %>%
  head (10)

## Warning: Values from `values` are not uniquely identified; output will contain
## list-cols.

## * Use `values_fn = list` to suppress this warning.

## * Use `values_fn = {summary Fun}` to summarise duplicates.

## * Use the following dplyr code to identify duplicates.

##  {data} |>
##  dplyr::summarise(n = dplyr::n(), .by = c(Word, variable)) |>
##  dplyr::filter(n > 1L)

## # A tibble: 10 x 4
46CHAPTER 1. INTRODUCTION TO R, THE TIDYVERSE - VISUALISATION

## Word RTlexdec RTnaming Familiarity
## <fct> <list> <list> <list>
## 1 doe <dbl [2]> <dbl [2]> <dbl [2]>
## 2 whore <dbl [2]> <dbl [2]> <dbl [2]>
## 3 stress <dbl [2]> <dbl [2]> <dbl [2]>
## 4 pork <dbl [2]> <dbl [2]> <dbl [2]>
## 5 plug <dbl [2]> <dbl [2]> <dbl [2]>
## 6 prop <dbl [2]> <dbl [2]> <dbl [2]>
## 7 dawn <dbl [2]> <dbl [2]> <dbl [2]>
## 8 dog <dbl [2]> <dbl [2]> <dbl [2]>
## 9 arc <dbl [2]> <dbl [2]> <dbl [2]>
## 10 skirt <dbl [2]> <dbl [2]> <dbl [2]>

But wait, where are the results? They are added in lists. We need to use the function unnest() to obtain the full results.

english %>
select (Word, RTlexdec, RTnaming, Familiarity) %>
pivot_longer (cols = c (RTlexdec, RTnaming, Familiarity),
              names_to = "variable",
              values_to = "values") %>
pivot_wider (names_from = "variable",
             values_from = "values") %>
unnest () %>
head (10)

## Warning: Values from `variable` are not uniquely identified; output will contain
list-cols.

## * Use `values_fn = list` to suppress this warning.

## * Use `values_fn = {summary Fun}` to summarise duplicates.

## * Use the following dplyr code to identify duplicates.

## {data} |
## dplyr::summarise(n = dplyr::n(), .by = c(Word, variable)) |
## dplyr::filter(n > 1L)

## Warning: `cols` is now required when using `unnest()`.

## i Please use `cols = c(RTlexdec, RTnaming, Familiarity)`.

## # A tibble: 10 x 4

<table><thead><tr><th>Word</th><th>RTlexdec</th><th>RTnaming</th><th>Familiarity</th></tr></thead><tbody><tr><td>&lt;fct&gt;</td><td>&lt;dbl&gt;</td><td>&lt;dbl&gt;</td><td>&lt;dbl&gt;</td></tr><tr><td>1 doe</td><td>6.54</td><td>6.15</td><td>2.37</td></tr><tr><td>2 doe</td><td>6.66</td><td>6.42</td><td>2.37</td></tr><tr><td>3 whore</td><td>6.40</td><td>6.25</td><td>4.43</td></tr><tr><td>4 whore</td><td>6.68</td><td>6.64</td><td>4.43</td></tr><tr><td>5 stress</td><td>6.30</td><td>6.14</td><td>5.6</td></tr><tr><td>6 stress</td><td>6.54</td><td>6.49</td><td>5.6</td></tr><tr><td>7 pork</td><td>6.42</td><td>6.13</td><td>3.87</td></tr></tbody></table>
<table><tr><td>##</td><td>8</td><td>pork</td><td>6.55</td><td>6.40</td><td>3.87</td></tr><tr><td>##</td><td>9</td><td>plug</td><td>6.45</td><td>6.20</td><td>3.93</td></tr><tr><td>##</td><td>10</td><td>plug</td><td>6.64</td><td>6.42</td><td>3.93</td></tr></table>

Ah that is better. But we get warnings. What does the warnings tell us? These are simple warnings and not errors. You can use the suggestions the *Tidyverse* makes. By default, we are told that the results are shown as lists of columns (what we are after). The second warning tells you to use a specific specification with unnest().

## 1.8.8 Basic descriptive statistics

### 1.8.8.1 Basic summaries

We can use `summary()` to obtain basic summaries of the dataset. For numeric variables, this will give you the minimum, maximum, mean, median, 1st and 3rd quartiles; for factors/characters, this will be the count. If there are missing values, you will get number of NAs. Look at the summaries of the dataset below.

```
english %>
summary()
```

<table><thead><tr><th>RTlexdec</th><th>RTnaming</th><th>Familiarity</th><th>Word</th><th>AgeSubject</th></tr></thead><tbody><tr><td>Min. :6.205</td><td>Min. :6.022</td><td>Min. :1.100</td><td>arm : 4</td><td>old :2284</td></tr><tr><td>1st Qu.:6.426</td><td>1st Qu.:6.149</td><td>1st Qu.:3.000</td><td>barge : 4</td><td>young:2284</td></tr><tr><td>Median :6.550</td><td>Median :6.342</td><td>Median :3.700</td><td>bark : 4</td><td></td></tr><tr><td>Mean :6.550</td><td>Mean :6.323</td><td>Mean :3.796</td><td>bear : 4</td><td></td></tr><tr><td>3rd Qu.:6.653</td><td>3rd Qu.:6.490</td><td>3rd Qu.:4.570</td><td>beef : 4</td><td></td></tr><tr><td>Max. :7.188</td><td>Max. :6.696</td><td>Max. :6.970</td><td>bind : 4</td><td></td></tr><tr><td></td><td></td><td></td><td>(Other):4544</td><td></td></tr><tr><th colspan="2">WordCategory</th><th>WrittenFrequency</th><th>WrittenSpokenFrequencyRatio</th><th>FamilySize</th></tr><tr><td>N:2904</td><td>Min. : 0.000</td><td>Min. :-6.55393</td><td>Min. :0.6931</td><td></td></tr><tr><td>V:1664</td><td>1st Qu.: 3.761</td><td>1st Qu.:-0.07402</td><td>1st Qu.:1.0986</td><td></td></tr><tr><td></td><td>Median : 4.832</td><td>Median : 0.68118</td><td>Median :1.7918</td><td></td></tr><tr><td></td><td>Mean : 5.021</td><td>Mean : 0.67763</td><td>Mean :1.8213</td><td></td></tr><tr><td></td><td>3rd Qu.: 6.247</td><td>3rd Qu.: 1.44146</td><td>3rd Qu.:2.3026</td><td></td></tr><tr><td></td><td>Max. :11.357</td><td>Max. : 5.63071</td><td>Max. :5.5175</td><td></td></tr></tbody></table>

<table><thead><tr><th colspan="2">DerivationalEntropy</th><th colspan="2">InflectionalEntropy</th><th colspan="2">NumberSimplexSynsets</th></tr></thead><tbody><tr><td>Min. :0.00000</td><td>Min. :0.0000</td><td>Min. :0.0000</td><td>Min. :0.0000</td><td>Min. :0.0000</td><td>Min. :0.0000</td></tr><tr><td>1st Qu.:0.03932</td><td>1st Qu.:0.7442</td><td>1st Qu.:0.7442</td><td>1st Qu.:1.099</td><td>1st Qu.:1.099</td><td>1st Qu.:1.099</td></tr><tr><td>Median :0.41097</td><td>Median :1.0982</td><td>Median :1.0982</td><td>Median :1.609</td><td>Median :1.609</td><td>Median :1.609</td></tr><tr><td>Mean :0.54089</td><td>Mean :1.1186</td><td>Mean :1.1186</td><td>Mean :1.708</td><td>Mean :1.708</td><td>Mean :1.708</td></tr><tr><td>3rd Qu.:0.89323</td><td>3rd Qu.:1.6325</td><td>3rd Qu.:1.6325</td><td>3rd Qu.:2.197</td><td>3rd Qu.:2.197</td><td>3rd Qu.:2.197</td></tr><tr><td>Max. :5.20728</td><td>Max. :2.4514</td><td>Max. :2.4514</td><td>Max. :4.357</td><td>Max. :4.357</td><td>Max. :4.357</td></tr></tbody></table>

<table><thead><tr><th colspan="2">NumberComplexSynsets</th><th>LengthInLetters</th><th>Ncount</th><th>MeanBigramFrequency</th></tr></thead><tbody></tbody></table>
##

### 1.8.8.2 Summary for a specific variable

```
english %>
summarise(count = n(),
           range_RTlexdec = range(RTlexdec),
           mean_RTlexdec = mean(RTlexdec),
           sd_RTlexdec = sd(RTlexdec),
           var_RTlexdec = var(RTlexdec),
           min_RTlexdec = min(RTlexdec),
           max_RTlexdec = max(RTlexdec),
           quart1_RTlexdec = quantile(RTlexdec, 0.25),
           quart1_RTlexdec = quantile(RTlexdec, 0.75),
           median_RTlexdec = median(RTlexdec))
```

```
## Warning: Returning more (or less) than 1 row per `summarise()` group was deprecated in
## dplyr 1.1.0.
```

```
## i Please use `reframe()` instead.
```

```
## i When switching from `summarise()` to `reframe()`, remember that `reframe()`
## always returns an ungrouped data frame and adjust accordingly.
```

```
## Call `lifecycle::last_lifecycle_warnings()` to see where this warning was
## generated.
```

```
## count range_RTlexdec mean_RTlexdec sd_RTlexdec var_RTlexdec min_RTlexdec
## 1  4568      6.205325    6.550097  0.1569188  0.02462351   6.205325
## 2  4568      7.187808    6.550097  0.1569188  0.02462351   6.205325
## max_RTlexdec quart1_RTlexdec median_RTlexdec
## 1      7.187808      6.653211      6.550466
## 2      7.187808      6.653211      6.550466
```

As you can see, we can add use `summarise` to obtain summaries of the dataset. We asked here for the mean, sd, variance, minimum and maximum values, etc. In the dataset `english`, we have many numeric variables, and if we want to obtain summaries for all of numeric variables, we can use `summarise_all`.

### 1.8.8.3 Summarise_all

If you want to add another level of summaries, e.g., for length, you can either add them as another level (with a new variable name) or use `summarise_all` to do that for you. We need to select only numeric variables to do that. This is the function to only select numeric variables `where(is.numeric)`. If you do not use it, you will get an error message

```
english %>
select(where(is.numeric)) %>
summarise_all(funs(mean = mean, sd = sd, var = var, min = min, max = max,
```
50CHAPTER 1. INTRODUCTION TO R, THE TIDYVERSE - VISUALISATION

range = range, median = median, Q1 = quantile(., prob s = 0.25), Q

## Warning: `funs()` was deprecated in dplyr 0.8.0.

## i Please use a list of either functions or lambdas:

##

## # Simple named list: list(mean = mean, median = median)

##

## # Auto named with `tibble::lst()`: tibble::lst(mean, median)

##

## # Using lambdas list(~ mean(. , trim = .2), ~ median(. , na .rm = TRUE))

## Call `lifecycle::last_lifecycle_warnings()` to see where this warning was

## generated.

## RTlexdec_mean RTnaming_mean Familiarity_mean WrittenFrequency_mean

## 1 6.550097 6.322505 3.795582 5.021145

## 2 6.550097 6.322505 3.795582 5.021145

## WrittenSpokenFrequencyRatio_mean FamilySize_mean DerivationalEntropy_mean

## 1 0.6776251 1.821324 0.5408901

## 2 0.6776251 1.821324 0.5408901

## InflectionalEntropy_mean NumberSimplexSynsets_mean NumberComplexSynsets_mean

## 1 1.118614 1.707797 1.567777

## 2 1.118614 1.707797 1.567777

## LengthInLetters_mean Ncount_mean MeanBigramFrequency_mean

## 1 4.342382 6.265762 8.489792

## 2 4.342382 6.265762 8.489792

## FrequencyInitialDiphone_mean ConspelV_mean ConspelN_mean ConphonV_mean

## 1 11.96292 11.71235 5.605324 18.25963

## 2 11.96292 11.71235 5.605324 18.25963

## ConphonN_mean ConfriendsV_mean ConfriendsN_mean ConffV_mean ConffN_mean

## 1 6.317727 10.41813 5.265387 0.410887 1.308101

## 2 6.317727 10.41813 5.265387 0.410887 1.308101

## ConfbV_mean ConfbN_mean NounFrequency_mean VerbFrequency_mean

## 1 1.556996 3.889859 600.1883 881.0201

## 2 1.556996 3.889859 600.1883 881.0201

## FrequencyInitialDiphoneWord_mean FrequencyInitialDiphoneSyllable_mean

## 1 10.35905 10.7892

## 2 10.35905 10.7892

## CorrectLexdec_mean RTlexdec_sd RTnaming_sd Familiarity_sd WrittenFrequency_sd

## 1 27.05166 0.1569188 0.1784815 1.149326 1.843559

## 2 27.05166 0.1569188 0.1784815 1.149326 1.843559

## WrittenSpokenFrequencyRatio_sd FamilySize_sd DerivationalEntropy_sd

## 1 1.165004 0.8205113 0.5573451

## 2 1.165004 0.8205113 0.5573451

## InflectionalEntropy_sd NumberSimplexSynsets_sd NumberComplexSynsets_sd

## 1 0.5660556 0.6693279 1.308999

## 2 0.5660556 0.6693279 1.308999
<table>
  <tbody>
    <tr>
      <td>##</td>
      <td>LengthInLetters_sd</td>
      <td>Ncount_sd</td>
      <td>MeanBigramFrequency_sd</td>
      <td></td>
      <td></td>
      <td></td>
    </tr>
    <tr>
      <td>## 1</td>
      <td></td>
      <td>0.8418279</td>
      <td>4.893778</td>
      <td></td>
      <td>0.6982397</td>
      <td></td>
    </tr>
    <tr>
      <td>## 2</td>
      <td></td>
      <td>0.8418279</td>
      <td>4.893778</td>
      <td></td>
      <td>0.6982397</td>
      <td></td>
    </tr>
    <tr>
      <td>##</td>
      <td>FrequencyInitialDiphone_sd</td>
      <td>ConspelV_sd</td>
      <td>ConspelN_sd</td>
      <td>ConphonV_sd</td>
      <td>ConphonN_sd</td>
      <td></td>
    </tr>
    <tr>
      <td>## 1</td>
      <td></td>
      <td>1.11962</td>
      <td>7.090395</td>
      <td>1.99571</td>
      <td>11.83507</td>
      <td>1.980334</td>
    </tr>
    <tr>
      <td>## 2</td>
      <td></td>
      <td>1.11962</td>
      <td>7.090395</td>
      <td>1.99571</td>
      <td>11.83507</td>
      <td>1.980334</td>
    </tr>
    <tr>
      <td>##</td>
      <td>ConfriendsV_sd</td>
      <td>ConfriendsN_sd</td>
      <td>ConffV_sd</td>
      <td>ConffN_sd</td>
      <td>ConfbV_sd</td>
      <td>ConfbN_sd</td>
    </tr>
    <tr>
      <td>## 1</td>
      <td></td>
      <td>6.923388</td>
      <td>2.076256</td>
      <td>0.6881427</td>
      <td>2.413261</td>
      <td>1.125484</td>
      <td>3.131148</td>
    </tr>
    <tr>
      <td>## 2</td>
      <td></td>
      <td>6.923388</td>
      <td>2.076256</td>
      <td>0.6881427</td>
      <td>2.413261</td>
      <td>1.125484</td>
      <td>3.131148</td>
    </tr>
    <tr>
      <td>##</td>
      <td>NounFrequency_sd</td>
      <td>VerbFrequency_sd</td>
      <td>FrequencyInitialDiphoneWord_sd</td>
      <td></td>
      <td></td>
      <td></td>
    </tr>
    <tr>
      <td>## 1</td>
      <td></td>
      <td>1858.115</td>
      <td>6852.356</td>
      <td></td>
      <td>1.583577</td>
      <td></td>
    </tr>
    <tr>
      <td>## 2</td>
      <td></td>
      <td>1858.115</td>
      <td>6852.356</td>
      <td></td>
      <td>1.583577</td>
      <td></td>
    </tr>
    <tr>
      <td>##</td>
      <td>FrequencyInitialDiphoneSyllable_sd</td>
      <td>CorrectLexdec_sd</td>
      <td>RTlexdec_var</td>
      <td>RTnaming_var</td>
      <td></td>
      <td></td>
    </tr>
    <tr>
      <td>## 1</td>
      <td></td>
      <td>1.598599</td>
      <td>4.302697</td>
      <td>0.02462351</td>
      <td>0.03185564</td>
      <td></td>
    </tr>
    <tr>
      <td>## 2</td>
      <td></td>
      <td>1.598599</td>
      <td>4.302697</td>
      <td>0.02462351</td>
      <td>0.03185564</td>
      <td></td>
    </tr>
    <tr>
      <td>##</td>
      <td>Familiarity_var</td>
      <td>WrittenFrequency_var</td>
      <td>WrittenSpokenFrequencyRatio_var</td>
      <td></td>
      <td></td>
      <td></td>
    </tr>
    <tr>
      <td>## 1</td>
      <td></td>
      <td>1.32095</td>
      <td>3.39871</td>
      <td></td>
      <td>1.357234</td>
      <td></td>
    </tr>
    <tr>
      <td>## 2</td>
      <td></td>
      <td>1.32095</td>
      <td>3.39871</td>
      <td></td>
      <td>1.357234</td>
      <td></td>
    </tr>
    <tr>
      <td>##</td>
      <td>FamilySize_var</td>
      <td>DerivationalEntropy_var</td>
      <td>InflectionalEntropy_var</td>
      <td></td>
      <td></td>
      <td></td>
    </tr>
    <tr>
      <td>## 1</td>
      <td></td>
      <td>0.6732388</td>
      <td>0.3106336</td>
      <td></td>
      <td>0.320419</td>
      <td></td>
    </tr>
    <tr>
      <td>## 2</td>
      <td></td>
      <td>0.6732388</td>
      <td>0.3106336</td>
      <td></td>
      <td>0.320419</td>
      <td></td>
    </tr>
    <tr>
      <td>##</td>
      <td>NumberSimplexSynsets_var</td>
      <td>NumberComplexSynsets_var</td>
      <td>LengthInLetters_var</td>
      <td></td>
      <td></td>
      <td></td>
    </tr>
    <tr>
      <td>## 1</td>
      <td></td>
      <td>0.4479998</td>
      <td>1.713479</td>
      <td></td>
      <td>0.7086742</td>
      <td></td>
    </tr>
    <tr>
      <td>## 2</td>
      <td></td>
      <td>0.4479998</td>
      <td>1.713479</td>
      <td></td>
      <td>0.7086742</td>
      <td></td>
    </tr>
    <tr>
      <td>##</td>
      <td>Ncount_var</td>
      <td>MeanBigramFrequency_var</td>
      <td>FrequencyInitialDiphone_var</td>
      <td>ConspelV_var</td>
      <td></td>
      <td></td>
    </tr>
    <tr>
      <td>## 1</td>
      <td></td>
      <td>23.94906</td>
      <td>0.4875387</td>
      <td></td>
      <td>1.25355</td>
      <td>50.27371</td>
    </tr>
    <tr>
      <td>## 2</td>
      <td></td>
      <td>23.94906</td>
      <td>0.4875387</td>
      <td></td>
      <td>1.25355</td>
      <td>50.27371</td>
    </tr>
    <tr>
      <td>##</td>
      <td>ConspelN_var</td>
      <td>ConphonV_var</td>
      <td>ConphonN_var</td>
      <td>ConfriendsV_var</td>
      <td>ConfriendsN_var</td>
      <td></td>
    </tr>
    <tr>
      <td>## 1</td>
      <td></td>
      <td>3.982859</td>
      <td>140.0688</td>
      <td>3.921722</td>
      <td>47.9333</td>
      <td>4.31084</td>
    </tr>
    <tr>
      <td>## 2</td>
      <td></td>
      <td>3.982859</td>
      <td>140.0688</td>
      <td>3.921722</td>
      <td>47.9333</td>
      <td>4.31084</td>
    </tr>
    <tr>
      <td>##</td>
      <td>ConffV_var</td>
      <td>ConffN_var</td>
      <td>ConfbV_var</td>
      <td>ConfbN_var</td>
      <td>NounFrequency_var</td>
      <td></td>
    </tr>
    <tr>
      <td>## 1</td>
      <td></td>
      <td>0.4735403</td>
      <td>5.823826</td>
      <td>1.266714</td>
      <td>9.804086</td>
      <td>3452593</td>
    </tr>
    <tr>
      <td>## 2</td>
      <td></td>
      <td>0.4735403</td>
      <td>5.823826</td>
      <td>1.266714</td>
      <td>9.804086</td>
      <td>3452593</td>
    </tr>
    <tr>
      <td>##</td>
      <td>VerbFrequency_var</td>
      <td>FrequencyInitialDiphoneWord_var</td>
      <td></td>
      <td></td>
      <td></td>
      <td></td>
    </tr>
    <tr>
      <td>## 1</td>
      <td></td>
      <td>46954789</td>
      <td></td>
      <td>2.507715</td>
      <td></td>
      <td></td>
    </tr>
    <tr>
      <td>## 2</td>
      <td></td>
      <td>46954789</td>
      <td></td>
      <td>2.507715</td>
      <td></td>
      <td></td>
    </tr>
    <tr>
      <td>##</td>
      <td>FrequencyInitialDiphoneSyllable_var</td>
      <td>CorrectLexdec_var</td>
      <td>RTlexdec_min</td>
      <td></td>
      <td></td>
      <td></td>
    </tr>
    <tr>
      <td>## 1</td>
      <td></td>
      <td>2.555518</td>
      <td>18.51321</td>
      <td></td>
      <td>6.205325</td>
      <td></td>
    </tr>
    <tr>
      <td>## 2</td>
      <td></td>
      <td>2.555518</td>
      <td>18.51321</td>
      <td></td>
      <td>6.205325</td>
      <td></td>
    </tr>
    <tr>
      <td>##</td>
      <td>RTnaming_min</td>
      <td>Familiarity_min</td>
      <td>WrittenFrequency_min</td>
      <td></td>
      <td></td>
      <td></td>
    </tr>
    <tr>
      <td>## 1</td>
      <td></td>
      <td>6.021751</td>
      <td>1.1</td>
      <td></td>
      <td>0</td>
      <td></td>
    </tr>
    <tr>
      <td>## 2</td>
      <td></td>
      <td>6.021751</td>
      <td>1.1</td>
      <td></td>
      <td>0</td>
      <td></td>
    </tr>
    <tr>
      <td>##</td>
      <td>WrittenSpokenFrequencyRatio_min</td>
      <td>FamilySize_min</td>
      <td>DerivationalEntropy_min</td>
      <td></td>
      <td></td>
      <td></td>
    </tr>
    <tr>
      <td>## 1</td>
      <td></td>
      <td>-6.553933</td>
      <td>0.6931472</td>
      <td></td>
      <td>0</td>
      <td></td>
    </tr>
    <tr>
      <td>## 2</td>
      <td></td>
      <td>-6.553933</td>
      <td>0.6931472</td>
      <td></td>
      <td>0</td>
      <td></td>
    </tr>
    <tr>
      <td>##</td>
      <td>InflectionalEntropy_min</td>
      <td>NumberSimplexSynsets_min</td>
      <td>NumberComplexSynsets_min</td>
      <td></td>
      <td></td>
      <td></td>
    </tr>
  </tbody>
</table>
## 1 0 0 0

## 2 0 0 0

## LengthInLetters_min Ncount_min MeanBigramFrequency_min

## 1 2 0 5.390053

## 2 2 0 5.390053

## FrequencyInitialDiphone_min ConspelV_min ConspelN_min ConphonV_min

## 1 4.14313 0 0 0

## 2 4.14313 0 0 0

## ConphonN_min ConfriendsV_min ConfriendsN_min ConffV_min ConffN_min ConfbV_min

## 1 0 0 0 0 0 0

## 2 0 0 0 0 0 0

## ConfbN_min NounFrequency_min VerbFrequency_min

## 1 0 0 0

## 2 0 0 0

## FrequencyInitialDiphoneWord_min FrequencyInitialDiphoneSyllable_min

## 1 3.091042 3.367296

## 2 3.091042 3.367296

## CorrectLexdec_min RTlexdec_max RTnaming_max Familiarity_max

## 1 1 7.187808 6.695675 6.97

## 2 1 7.187808 6.695675 6.97

## WrittenFrequency_max WrittenSpokenFrequencyRatio_max FamilySize_max

## 1 11.35666 5.630714 5.517453

## 2 11.35666 5.630714 5.517453

## DerivationalEntropy_max InflectionalEntropy_max NumberSimplexSynsets_max

## 1 5.20728 2.45141 4.356709

## 2 5.20728 2.45141 4.356709

## NumberComplexSynsets_max LengthInLetters_max Ncount_max

## 1 6.111467 7 22

## 2 6.111467 7 22

## MeanBigramFrequency_max FrequencyInitialDiphone_max ConspelV_max ConspelN_max

## 1 10.28277 14.65437 32 10.49202

## 2 10.28277 14.65437 32 10.49202

## ConphonV_max ConphonN_max ConfriendsV_max ConfriendsN_max ConffV_max

## 1 66 10.59975 31 10.30304 3.332205

## 2 66 10.59975 31 10.30304 3.332205

## ConffN_max ConfbV_max ConfbN_max NounFrequency_max VerbFrequency_max

## 1 10.34744 4.189655 10.59975 35351 242066

## 2 10.34744 4.189655 10.59975 35351 242066

## FrequencyInitialDiphoneWord_max FrequencyInitialDiphoneSyllable_max

## 1 13.9249 13.92962

## 2 13.9249 13.92962

## CorrectLexdec_max RTlexdec_range RTnaming_range Familiarity_range

## 1 30 6.205325 6.021751 1.10

## 2 30 7.187808 6.695675 6.97

## WrittenFrequency_range WrittenSpokenFrequencyRatio_range FamilySize_range

## 1 0.00000 -6.553933 0.6931472
<table>
  <tbody>
    <tr>
      <td>## 2</td>
      <td>11.35666</td>
      <td>5.630714</td>
      <td>5.5174529</td>
    </tr>
    <tr>
      <td>##</td>
      <td>DerivationalEntropy_range InflectionalEntropy_range</td>
      <td></td>
      <td></td>
    </tr>
    <tr>
      <td>## 1</td>
      <td>0.00000</td>
      <td>0.00000</td>
      <td></td>
    </tr>
    <tr>
      <td>## 2</td>
      <td>5.20728</td>
      <td>2.45141</td>
      <td></td>
    </tr>
    <tr>
      <td>##</td>
      <td>NumberSimplexSynsets_range NumberComplexSynsets_range LengthInLetters_range</td>
      <td></td>
      <td></td>
    </tr>
    <tr>
      <td>## 1</td>
      <td>0.000000</td>
      <td>0.000000</td>
      <td>2</td>
    </tr>
    <tr>
      <td>## 2</td>
      <td>4.356709</td>
      <td>6.111467</td>
      <td>7</td>
    </tr>
    <tr>
      <td>##</td>
      <td>Ncount_range MeanBigramFrequency_range FrequencyInitialDiphone_range</td>
      <td></td>
      <td></td>
    </tr>
    <tr>
      <td>## 1</td>
      <td>0</td>
      <td>5.390053</td>
      <td>4.14313</td>
    </tr>
    <tr>
      <td>## 2</td>
      <td>22</td>
      <td>10.282767</td>
      <td>14.65437</td>
    </tr>
    <tr>
      <td>##</td>
      <td>ConspelV_range ConspelN_range ConphonV_range ConphonN_range ConfriendsV_range</td>
      <td></td>
      <td></td>
    </tr>
    <tr>
      <td>## 1</td>
      <td>0</td>
      <td>0.00000</td>
      <td>0.00000</td>
    </tr>
    <tr>
      <td>## 2</td>
      <td>32</td>
      <td>10.49202</td>
      <td>66</td>
    </tr>
    <tr>
      <td>##</td>
      <td>ConfriendsN_range ConffV_range ConffN_range ConfbV_range ConfbN_range</td>
      <td></td>
      <td></td>
    </tr>
    <tr>
      <td>## 1</td>
      <td>0.00000</td>
      <td>0.000000</td>
      <td>0.00000</td>
    </tr>
    <tr>
      <td>## 2</td>
      <td>10.30304</td>
      <td>3.332205</td>
      <td>4.189655</td>
    </tr>
    <tr>
      <td>##</td>
      <td>NounFrequency_range VerbFrequency_range FrequencyInitialDiphoneWord_range</td>
      <td></td>
      <td></td>
    </tr>
    <tr>
      <td>## 1</td>
      <td>0</td>
      <td>0</td>
      <td>3.091042</td>
    </tr>
    <tr>
      <td>## 2</td>
      <td>35351</td>
      <td>242066</td>
      <td>13.924902</td>
    </tr>
    <tr>
      <td>##</td>
      <td>FrequencyInitialDiphoneSyllable_range CorrectLexdec_range RTlexdecarf</td>
      <td></td>
      <td></td>
    </tr>
    <tr>
      <td>## 1</td>
      <td>3.367296</td>
      <td>1</td>
      <td>6.550466</td>
    </tr>
    <tr>
      <td>## 2</td>
      <td>13.929620</td>
      <td>30</td>
      <td>6.550466</td>
    </tr>
    <tr>
      <td>##</td>
      <td>RTnamingarf Familiarityarf WrittenFrequencyarf</td>
      <td></td>
      <td></td>
    </tr>
    <tr>
      <td>## 1</td>
      <td>6.342023</td>
      <td>3.7</td>
      <td>4.832298</td>
    </tr>
    <tr>
      <td>## 2</td>
      <td>6.342023</td>
      <td>3.7</td>
      <td>4.832298</td>
    </tr>
    <tr>
      <td>##</td>
      <td>WrittenSpokenFrequencyRatioarf FamilySizearf</td>
      <td></td>
      <td></td>
    </tr>
    <tr>
      <td>## 1</td>
      <td>0.681184</td>
      <td>1.791759</td>
      <td></td>
    </tr>
    <tr>
      <td>## 2</td>
      <td>0.681184</td>
      <td>1.791759</td>
      <td></td>
    </tr>
    <tr>
      <td>##</td>
      <td>DerivationalEntropyarf InflectionalEntropy_range</td>
      <td></td>
      <td></td>
    </tr>
    <tr>
      <td>## 1</td>
      <td>0.410975</td>
      <td>1.09821</td>
      <td></td>
    </tr>
    <tr>
      <td>## 2</td>
      <td>0.410975</td>
      <td>1.09821</td>
      <td></td>
    </tr>
    <tr>
      <td>##</td>
      <td>NumberSimplexSynsetsarf NumberComplexSynsetsarf</td>
      <td></td>
      <td></td>
    </tr>
    <tr>
      <td>## 1</td>
      <td>1.609438</td>
      <td>1.386294</td>
      <td></td>
    </tr>
    <tr>
      <td>## 2</td>
      <td>1.609438</td>
      <td>1.386294</td>
      <td></td>
    </tr>
    <tr>
      <td>##</td>
      <td>LengthInLettersarf Ncountarf MeanBigramFrequencyarf</td>
      <td></td>
      <td></td>
    </tr>
    <tr>
      <td>## 1</td>
      <td>4</td>
      <td>5</td>
      <td>8.55855</td>
    </tr>
    <tr>
      <td>## 2</td>
      <td>4</td>
      <td>5</td>
      <td>8.55855</td>
    </tr>
    <tr>
      <td>##</td>
      <td>FrequencyInitialDiphonearf ConspelVarf ConspelNarf</td>
      <td></td>
      <td></td>
    </tr>
    <tr>
      <td>## 1</td>
      <td>12.02268</td>
      <td>11</td>
      <td>5.710427</td>
    </tr>
    <tr>
      <td>## 2</td>
      <td>12.02268</td>
      <td>11</td>
      <td>5.710427</td>
    </tr>
    <tr>
      <td>##</td>
      <td>ConphonVarf ConphonNarf ConfriendsVarf ConfriendsNarf</td>
      <td></td>
      <td></td>
    </tr>
    <tr>
      <td>## 1</td>
      <td>16</td>
      <td>6.340359</td>
      <td>10</td>
    </tr>
    <tr>
      <td>## 2</td>
      <td>16</td>
      <td>6.340359</td>
      <td>10</td>
    </tr>
    <tr>
      <td>##</td>
      <td>ConffVarf ConffNarf ConfbVarf ConfbNarf NounFrequencyarf</td>
      <td></td>
      <td></td>
    </tr>
    <tr>
      <td>## 1</td>
      <td>0</td>
      <td>0</td>
      <td>1.386294</td>
    </tr>
    <tr>
      <td>## 2</td>
      <td>0</td>
      <td>0</td>
      <td>1.386294</td>
    </tr>
  </tbody>
</table>
VerbFrequency_median FrequencyInitialDiphoneWord_median

1 30 10.51651

2 30 10.51651

FrequencyInitialDiphoneSyllable_median CorrectLexdec_median RTlexdec_Q1

1 10.97207 29 6.425525

2 10.97207 29 6.425525

RTnaming_Q1 Familiarity_Q1 WrittenFrequency_Q1 WrittenSpokenFrequencyRatio_Q1

1 6.148682 3 3.7612 -0.07401695

2 6.148682 3 3.7612 -0.07401695

FamilySize_Q1 DerivationalEntropy_Q1 InflectionalEntropy_Q1

1 1.098612 0.039325 0.7441725

2 1.098612 0.039325 0.7441725

NumberSimplexSynsets_Q1 NumberComplexSynsets_Q1 LengthInLetters_Q1 Ncount_Q1

1 1.098612 0 4 2

2 1.098612 0 4 2

MeanBigramFrequency_Q1 FrequencyInitialDiphone_Q1 ConspelV_Q1 ConspelN_Q1

1 8.099924 11.27694 6 4.519056

2 8.099924 11.27694 6 4.519056

ConphonV_Q1 ConphonN_Q1 ConfriendsV_Q1 ConfriendsN_Q1 ConffV_Q1 ConffN_Q1

1 10 5.267858 4 4.158883 0 0

2 10 5.267858 4 4.158883 0 0

ConfbV_Q1 ConfbN_Q1 NounFrequency_Q1 VerbFrequency_Q1

1 0.6931472 0 28.75 0

2 0.6931472 0 28.75 0

FrequencyInitialDiphoneWord_Q1 FrequencyInitialDiphoneSyllable_Q1

1 9.556808 10.00048

2 9.556808 10.00048

CorrectLexdec_Q1 RTlexdec_Q3 RTnaming_Q3 Familiarity_Q3 WrittenFrequency_Q3

1 27 6.653211 6.489699 4.57 6.247074

2 27 6.653211 6.489699 4.57 6.247074

WrittenSpokenFrequencyRatio_Q3 FamilySize_Q3 DerivationalEntropy_Q3

1 1.44146 2.302585 0.893225

2 1.44146 2.302585 0.893225

InflectionalEntropy_Q3 NumberSimplexSynsets_Q3 NumberComplexSynsets_Q3

1 1.632455 2.197225 2.564949

2 1.632455 2.197225 2.564949

LengthInLetters_Q3 Ncount_Q3 MeanBigramFrequency_Q3

1 5 9 8.972658

2 5 9 8.972658

FrequencyInitialDiphone_Q3 ConspelV_Q3 ConspelN_Q3 ConphonV_Q3 ConphonN_Q3

1 12.69734 17 6.997358 24 7.491413

2 12.69734 17 6.997358 24 7.491413

ConfriendsV_Q3 ConfriendsN_Q3 ConffV_Q3 ConffN_Q3 ConfbV_Q3 ConfbN_Q3

1 15 6.642077 0.6931472 1.386294 2.564949 6.241734

2 15 6.642077 0.6931472 1.386294 2.564949 6.241734

NounFrequency_Q3 VerbFrequency_Q3 FrequencyInitialDiphoneWord_Q3
<table><tr><td>## 1</td><td>424.75</td><td>164.25</td><td>11.31995</td></tr><tr><td>## 2</td><td>424.75</td><td>164.25</td><td>11.31995</td></tr><tr><td colspan="4">## FrequencyInitialDiphoneSyllable_Q3 CorrectLexdec_Q3</td></tr><tr><td>## 1</td><td></td><td>11.70264</td><td>30</td></tr><tr><td>## 2</td><td></td><td>11.70264</td><td>30</td></tr></table>

As you can see, in this example, we see the chains of commands in the Tidyverse. We can continue to add commands each time we want to investigate something in particular. Keep adding pipes and commands. The most important point is that the dataset `english` did not change at all. If you want to create a new dataset with the results, simply use the assignment function `<-` at the beginning or `->` at the end and give a name to the new dataset.

### 1.8.8.4 Group_by

### 1.8.8.5 One variable

What if you want to obtain all results summarised by a specific grouping? Let's obtain the results grouped by the levels of `AgeSubject`.

```
english %>%
  group_by(AgeSubject) %>%
  summarise(count = n(),
            range_RTlexdec = range(RTlexdec),
            mean_RTlexdec = mean(RTlexdec),
            sd_RTlexdec = sd(RTlexdec),
            var_RTlexdec = var(RTlexdec),
            min_RTlexdec = min(RTlexdec),
            max_RTlexdec = max(RTlexdec),
            quart1_RTlexdec = quantile(RTlexdec, 0.25),
            quart1_RTlexdec = quantile(RTlexdec, 0.75),
            median_RTlexdec = median(RTlexdec))
```

```
## Warning: Returning more (or less) than 1 row per `summarise()` group was deprecated in
## dplyr 1.1.0.
```

```
## i Please use `reframe()` instead.
```

```
## i When switching from `summarise()` to `reframe()`, remember that `reframe()`
## always returns an ungrouped data frame and adjust accordingly.
```

```
## Call `lifecycle::last_lifecycle_warnings()` to see where this warning was
## generated.
```

```
## `summarise()` has grouped output by 'AgeSubject'. You can override using the
## `.groups` argument.
```

```
## # A tibble: 4 x 10
```

```
## # Groups: AgeSubject [2]
```

```
## AgeSubject count range_RTlexdec mean_RTlexdec sd_RTlexdec var_RTlexdec
## <fct> <int> <db1> <db1> <db1> <db1>
## 1 old 2284 6.40 6.66 0.116 0.0134
```
56CHAPTER 1. INTRODUCTION TO R, THE TIDYVERSE - VISUALISATION

## 2 old 2284 7.19 6.66 0.116 0.0134
## 3 young 2284 6.21 6.44 0.106 0.0113
## 4 young 2284 6.88 6.44 0.106 0.0113

## # i 4 more variables: min_RTlexdec <db1>, max_RTlexdec <db1>,
## # quart1_RTlexdec <db1>, median_RTlexdec <db1>

1.8.8.6 Multiple variables

What if you want to obtain all results summarised by multiple groupings? Let's
obtain the results grouped by the levels of AgeSubject, WordCategory and
Voice and we want to save the output.

english %>%
group_by(AgeSubject, WordCategory, Voice) %>%
  summarise(count = n(),
    range_RTlexdec = range(RTlexdec),
    mean_RTlexdec = mean(RTlexdec),
    sd_RTlexdec = sd(RTlexdec),
    var_RTlexdec = var(RTlexdec),
    min_RTlexdec = min(RTlexdec),
    max_RTlexdec = max(RTlexdec),
    quart1_RTlexdec = quantile(RTlexdec, 0.25),
    quart1_RTlexdec = quantile(RTlexdec, 0.75),
    median_RTlexdec = median(RTlexdec)) -> dfMeans

## Warning: Returning more (or less) than 1 row per `summarise()` group was deprecated
## dplyr 1.1.0.

## i Please use `reframe()` instead.

## i When switching from `summarise()` to `reframe()`, remember that `reframe()`

## always returns an ungrouped data frame and adjust accordingly.

## Call `lifecycle::last_lifecycle_warnings()` to see where this warning was
generated.

## `summarise()` has grouped output by 'AgeSubject', 'WordCategory', 'Voice'. You
## can override using the `.groups` argument.

dfMeans

## # A tibble: 16 x 12
## # Groups: AgeSubject, WordCategory, Voice [8]

## AgeSubject WordCategory Voice count range_RTlexdec mean_RTlexdec sd_RTlexdec
## <fct> <fct> <fct> <int> <db1> <db1> <db1>
## 1 old N voiced 676 6.45 6.67 0.121
## 2 old N voiced 676 7.15 6.67 0.121
## 3 old N voice~ 776 6.45 6.66 0.120
## 4 old N voice~ 776 7.19 6.66 0.120
## 5 old V voiced 354 6.40 6.66 0.112
## 6 old V voiced 354 7.14 6.66 0.112
<table><tr><td>## 7 old</td><td>V</td><td>voice~</td><td>478</td><td>6.46</td><td>6.65</td><td>0.101</td></tr><tr><td>## 8 old</td><td>V</td><td>voice~</td><td>478</td><td>7.01</td><td>6.65</td><td>0.101</td></tr><tr><td>## 9 young</td><td>N</td><td>voiced</td><td>676</td><td>6.21</td><td>6.45</td><td>0.116</td></tr><tr><td>## 10 young</td><td>N</td><td>voiced</td><td>676</td><td>6.88</td><td>6.45</td><td>0.116</td></tr><tr><td>## 11 young</td><td>N</td><td>voice~</td><td>776</td><td>6.21</td><td>6.44</td><td>0.102</td></tr><tr><td>## 12 young</td><td>N</td><td>voice~</td><td>776</td><td>6.79</td><td>6.44</td><td>0.102</td></tr><tr><td>## 13 young</td><td>V</td><td>voiced</td><td>354</td><td>6.24</td><td>6.43</td><td>0.106</td></tr><tr><td>## 14 young</td><td>V</td><td>voiced</td><td>354</td><td>6.75</td><td>6.43</td><td>0.106</td></tr><tr><td>## 15 young</td><td>V</td><td>voice~</td><td>478</td><td>6.21</td><td>6.43</td><td>0.0982</td></tr><tr><td>## 16 young</td><td>V</td><td>voice~</td><td>478</td><td>6.83</td><td>6.43</td><td>0.0982</td></tr></table>

```r
## # i 5 more variables: var_RTlexdec <dbl>, min_RTlexdec <dbl>,
## #   max_RTlexdec <dbl>, quart1_RTlexdec <dbl>, median_RTlexdec <dbl>
```

## 1.8.9 Summary tables

### 1.8.9.1 Using summarytools

Using the package `summarytools` you can obtain nice table ready for publications!

```r
english %>%
  dfSummary(graph.col = TRUE, style = "grid") %>%
  stview(method = "render")
```

#### 1.8.9.1.1 Summary of full dataset with graphs

```r
english %>%
  descr() %>%
  stview(method = "render")
```

#### 1.8.9.1.2 Descriptive statistics for numeric variables

```r
## Non-numerical variable(s) ignored: Word, AgeSubject, WordCategory, CV, Obstruent, Frication, V
```

```r
english %>%
  freq() %>%
  stview(method = "render")
```

#### 1.8.9.1.3 Frequencies for non-numeric variables

```r
## Variable(s) ignored: RTlexdec, RTnaming, Familiarity, WrittenFrequency, WrittenSpokenFrequency
```

### 1.8.9.2 Using gtsummary

We use the `gtsummary` package to create nice tables for publication. This package is very useful for creating publication ready tables with descriptive statistics,
<table><thead><tr><th>Characteristic</th><th>N = 4,568<sup>1</sup></th></tr></thead><tbody><tr><td>AgeSubject</td><td></td></tr><tr><td>old</td><td>2,284 (50%)</td></tr><tr><td>young</td><td>2,284 (50%)</td></tr><tr><td>RTlexdec</td><td>6.55 (6.43, 6.65)</td></tr><tr><td>Familiarity</td><td>3.70 (3.00, 4.57)</td></tr><tr><td>WordCategory</td><td></td></tr><tr><td>N</td><td>2,904 (64%)</td></tr><tr><td>V</td><td>1,664 (36%)</td></tr><tr><td>Voice</td><td></td></tr><tr><td>voiced</td><td>2,060 (45%)</td></tr><tr><td>voiceless</td><td>2,508 (55%)</td></tr></tbody></table>

¹n (%); Median (Q1, Q3)

regression models, etc. See this link for more details on the various options possible with this package.

english %>%
tbl_summary(include = c(AgeSubject, RTlexdec, Familiarity, WordCategory, Voice))

### 1.8.9.2.1 Basic summaries

#### 1.8.9.2.2 With statistics (comparing means)
We use the two numeric outcomes and assess differences between the two groups.

english %>%
tbl_summary(include = c(RTlexdec, RTnaming),
by = AgeSubject, # split table by AgeSubject
missing = "no" # don't list missing data separately
) %>%
add_n() %>% # add column with total number of non-missing observations
add_p() %>% # test for a difference between groups
modify_header label = "**Variable**" %>% # update the column header
bold labels()

### 1.8.10 Visualisation

In the *tidyverse*, the package for making elegant plots is called `ggplot2`. It works a lot like how pipes work, but since it was originally designed as a separate package, it uses `+` instead of `%>-%`.
<table><thead><tr><th>Variable</th><th>N</th><th>old N = 2,284<sup>1</sup></th><th>young N = 2,284<sup>1</sup></th><th>p-value<sup>2</sup></th></tr></thead><tbody><tr><td>RTlexdec</td><td>4,568</td><td>6.64 (6.57, 6.73)</td><td>6.43 (6.36, 6.51)</td><td>&lt;0.001</td></tr><tr><td>RTnaming</td><td>4,568</td><td>6.49 (6.45, 6.53)</td><td>6.15 (6.12, 6.18)</td><td>&lt;0.001</td></tr></tbody></table>

¹Median (Q1, Q3)

²Wilcoxon rank sum test

### 1.8.10.1 First steps

1.8.10.1.1 **Empty plot area** Let's produce a basic plot with nothing drawn on it. This is the basic plotting area in R. We need to then add layers on top of it to show our plot

```r
english %>%
  ggplot() +
  themeBW()
```

1.8.10.1.2 **Adding x and y values** Let's add the x and y values from our dataset. X = subjective familiarity rating, y = RT in Visual Lexical Decision task

```r
english %>%
  ggplot(aes(x = Familiarity,
             y = RTlexdec)) +
```
theme()

<figure><img src="image_7.png" /></figure>

There are no differences between the two. We need to tell ggplot2 to add a
geometric function for plotting

1.8.10.1.3 Adding geom s
Geoms are integrated within ggplot2 to obtain
various types of plots.

english %>

```r
ggplot(aes(x = Familiarity,
           y = RTlexdec)) +
  themeBW() +
  geom_point()
```
<figure><img src="image_8.png" /></figure>

1.8.10.1.4 Adding line of best fit We will add a line of best fit. This is used to evaluate presence/absence of a relationship between two numeric variables

english %>

```
ggplot(aes(x = Familiarity,
           y = RTlexdec)) +
  themeBW() +
  geom_point() +
  geom.smooth(method = "lm") ## line of best fit based on the lm() method
```

```r
## `geom.smooth()` using formula = 'y ~ x'
```
<figure><img src="image_9.png" /></figure>

The result shows a nice negative correlation! RT lexical decision decreases when familiarity rating increases.

We can ask, are there differences related to the word category, i.e., verb vs noun?

1.8.10.1.5 By word category We change colour by levels of word category;

```r
ggplot(aes(x = Familiarity,
           y = RTlexdec,
           colour = WordCategory)) + ## add colour to the base aesthetics
  themeBW() +
  geom_point() +
  geom.smooth(method = "lm")
```

```r
## `geom.smooth()` using formula = 'y ~ x'
```
<figure><img src="image_10.png" /></figure>

1.8.10.1.6 Making final touches
Let's add a title and a subtitle, change
x and y labels, change size of overall plot, and colours of the categories.

english %>

```r
ggplot(aes(x = Familiarity,
           y = RTlexdec,
           colour = WordCategory)) + ## add colour to the base aesthetics
theme_bw() +
geom_point() +
geom.smooth(method = "lm") +
labs(x = "Familiarity rating", y = "RT Lexical Decision", title = "Familiarity rating vs RT in theme(text = element_text(size = 15)) + ## increase size of plot
theme(legendbottom, legend.title = elementقد)) + ## remove legend title and scale_color_manual labels = c("Nouns", "Verbs"), values = c("blue", "red")) ## change colours
```

## `geom.smooth()` using formula = 'y ~ x'
<figure><img src="image_11.png" /></figure>

To choose colours, use the addin `colour picker` from above. See this link for full list of colours available. Use colours that are colour-blind friendly here

### 1.8.10.2 Additional plots

We looked above at one example of plots (with points). We could use additional types of plots.

#### 1.8.10.2.1 A bar plot Will show barplots of the dataset

```r
e English %>

ggplot(aes(x = RTlexdec,
           colour = AgeSubject)) +
  themeBW() +
  geom_bar()
```
<figure><img src="image_12.png" /></figure>

And another view with error bars! This is a nice example that shows how you
can combine multiple chains with the pipe:

• Group by Age of subject

• Compute mean and SD

• use ggplot2 syntax to plot a barplot and error bars

english %>
group_by(AgeSubject) %>
summarise(
  sd = sd(RTlexdec),
  RTlexdecM = mean(RTlexdec)
) %>
ggplot(aes(x = AgeSubject,
           y = RTlexdecM)) +
themeBW(
  geom_col(fill = "lightgray", color = "black") +
  geom_errorbar(aes(ymin = RTlexdecM-sd, ymax = RTlexdecM+sd), width = 0.2)
<figure><img src="image_13.png" /></figure>

1.8.10.2.2 A histogram This looks at the distribution of the variable. We
look at a histogram

english %>

```
ggplot(aes(x = RTlexdec,
           colour = AgeSubject)) +
  themeBW() +
  geom_histogram(fill = "white") +
  scale_color_manual(values = c("red", "blue"))
```

## `stat_bin()` using `bins = 30`. Pick better value `binwidth`.
<figure><img src="image_14.png" /></figure>

**1.8.10.2.3 A density plot** This looks at the distribution of the variable. We see that the two variables have different means. We can superpose the density plot on top of the histogram or have the density plot on its own.

```r
ggplot(aes(x = RTlexdec,
           colour = AgeSubject)) +
themeBW() +
geom_histogram(aes(y = ..density..), fill = "white") +
scale_color_manual(values = c("red", "blue")) +
geom_density()
```

**1.8.10.2.3.1 Histogram and density plot**

```r
## Warning: The dot-dot notation (`..density..`) was deprecated in ggplot2 3.4.0.
## i Please use `after_stat(density)` instead.
## This warning is displayed once every 8 hours.
## Call `lifecycle::last_lifecycle_warnings()` to see where this warning was
## generated.
```

```r
## `stat_bin()` using `bins = 30`. Pick better value `binwidth`.
```
<figure><img src="image_15.png" /></figure>

```r
ggplot(aes(x = RTlexdec,
           colour = AgeSubject)) +
  themeBW() +
  geom_density()
```
<figure><img src="image_16.png" /><figcaption>1.8.10.2.3.2 Density plot only</figcaption></figure>

1.8.10.2.4 A boxplot This allows you to see various information, including the Median, SD, Quartiles (25% and 75%) and outliers. Looking at the medians, we see clear difference between the two distributions.

english %>
ggplot(aes(x = AgeSubject,
y = RTlexdec)) +
themeBW() +
geom_boxplot()
<figure><img src="image_17.png" /></figure>

1.8.10.2.5 **A Violin plot** This allows you to see various information, including the Median, SD, Quartiles (25% and 75%) and outliers. Looking at the medians, we see clear difference between the two distributions.

english %>

```r
ggplot(aes(x = AgeSubject,
           y = RTlexdec)) +
  themeBW() +
  geom_violin()
```
<figure><img src="image_18.png" /></figure>

1.8.10.3 Advanced plots

We will use the package ggdist see herewith allows us to obtain various plots
adding the distributions.

1.8.10.3.1 A dotplot This dotplot allow to display the distribution of the
dataframe showing their density

english %>
ggplot(aes(x = AgeSubject,
           y = RTlexdec)) +
themeBW() +
stat_dotsinterval()
<figure><img src="image_19.png" /></figure>

1.8.10.3.2 **Rain cloud plots**
Rain clouds are used to combine both a density plot with dotplots in the same figure. These can be informative as to the overall distribution of the dataset.

```r
ggplot(aes(x = AgeSubject,
           y = RTlexdec,
           fill = AgeSubject)) +
  themeBW() +
  stat_slab(aes(thickness = after_stat(pdf*n)), scale = 0.7) +
  stat_dotsInterval(side = "bottom", scale = 0.7, slab_linewidth = NA) +
  scale_fill_brewer palette = "Set2")
```
<figure><img src="image_20.png" /></figure>

1.8.10.3.3 Eye plots and half-eye plots We can also obtain density plots with eye plots and half-eye plots. These are similar to rain clouds but al- low you to see the distribution of the data in a different way. For these plots, using stat_slabinterval() with side = "both", "top", "bottom", "right", or "left" you can adjust positions of the eye plots. Note that with "top", or "bottom", the y value is a categorical predictor and the x is a nu- meric outcome.

1.8.10.3.3.1 Eye plot We can use statEye() or stat_slabinterval(side = "both") to obtain eye plots.

english %>

```r
ggplot(aes(x = AgeSubject,
           y = RTlexdec,
           fill = AgeSubject)) +
  theme_bw() +
  statEye()
```
74CHAPTER 1. INTRODUCTION TO R, THE TIDYVERSE - VISUALISATION

<figure><img src="image_21.png" /></figure>

```r
ggplot(aes(x = AgeSubject,
           y = RTlexdec,
           fill = AgeSubject)) +
themeBW() +
stat_slabinterval(side = "both")
```
<figure><img src="image_22.png" /></figure>

1.8.10.3.3.2 **Half-eye plot** We can use `stat EMME` or `stat_slabinterval()` to obtain half-eye plots. Use `stat_slabinterval(side = "left")` or `stat_slabinterval(side = "right")` to obtain half-eye plots on the left or right side of the plot, respectively.

```r
ggplot(aes(x = AgeSubject,
           y = RTlexdec,
           fill = AgeSubject)) +
  themeBW() +
  stat EMME
```
<figure><img src="image_23.png" /></figure>

```r
ggplot(aes(x = AgeSubject,
           y = RTlexdec,
           fill = AgeSubject)) +
themeBW() +
stat_slabinterval()
```
<figure><img src="image_24.png" /></figure>

1.8.10.4 Facet_grid

The plots we used so far allowed to plot data as a function of one categorical
variable, e.g., AgeSubject. What if we wanted to show the different patterns
emerging when combining AgeSubject (old vs young), WordCategory (Noun or
Verb), CV (Consonant or Vowel) and Voice (Voiced and Voiceless) ? What if
we also wanted to modify the labels and order of levels of variables?

We will start slowly below to show how we can combine two categorical variables
and extend them to additional ones

1.8.10.4.1 Two categorical variables

1.8.10.4.1.1 First steps Here we obtain a boxplot with two categorical
variables AgeSubject and WordCategory

english %>

```
ggplot(aes(x = AgeSubject,
           y = RTlexdec)) +
  themeBW() +
  geom_boxplot() +
  facet_grid(~ WordCategory)
```
<figure><img src="image_25.png" /></figure>

1.8.10.4.1.2 Changing order of levels within a variable and its labels

What would you do to change both order of levels within a variable and its
labels? We want to change order for AgeSubject to be Young vs Old (rather
than old vs young) and change labels of WordCategory from N vs V to Noun
vs Verb?

english %>%
  mutate(AgeSubject = factor(AgeSubject, levels = c("young", "old"), labels = c("Young" ,
  WordCategory = factor(WordCategory, labels = c("Noun", "Verb"))) %>%
  ggplot(aes(x = AgeSubject,
            y = RTlexdec)) +
  themeBW() +
  geom_boxplot() +
  facet_grid(~ WordCategory)
80CHAPTER 1. INTRODUCTION TO R, THE TIDYVERSE - VISUALISATION

<figure><img src="image_26.png" /></figure>

1.8.10.4.3 Comparing two numeric outcomes What if we want to compare performance in relation to reaction time for the lexical decision task (RTlexdec) and reaction time for naming (RTnaming). We want to see if there are differences related to the AgeSubject, WordCategory. We use pivot_longer here to do change the format of our table and then change names and use facet_grid.

english %>
select(RTlexdec, RTnaming, AgeSubject, WordCategory) %>
pivot_longer cols = c(RTlexdec, RTnaming),
  names_to = "variable",
  values_to = "values") %>
mutate(AgeSubject = factor(AgeSubject, levels = c("young", "old"), labels = c("Young" "Old")),
  WordCategory = factor(WordCategory, labels = c("Noun", "Verb"))) %>
ggplot(aes(x = variable,
  y = values)) +
themeBW () +
geom_boxplot () +
facet_grid(AgeSubject ~ WordCategory, margins = TRUE, scales = "free")
<figure><img src="image_27.png" /></figure>

1.8.10.4.4 Exporting images
When you use R markdown, your figures are already embedded within the generated output. If you are using an R script and/or want to add the figure in a different document, you can use the following code:

jpeg(filename = paste0("outputs/test.jpeg"), width = 15, height = 15, units = "cm", res = 300)

english %>

select(RTlexdec, RTnaming, AgeSubject, WordCategory) %>

pivot_longer cols = c(RTlexdec, RTnaming),
  names_to = "variable",
  values_to = "values") %>

mutate(AgeSubject = factor(AgeSubject, levels = c("young", "old"), labels = c("Young", "Old")),
  WordCategory = factor(WordCategory, labels = c("Noun", "Verb"))) %>

ggplot(aes(x = variable,
  y = values)) +

themeBW() +

geom_boxplot() +

facet_grid(AgeSubject ~ WordCategory, margins = TRUE, scales = "free")

dev.off()

## pdf

## 2

The image is automatically saved into your working directory and you can import it to your word (LATEX) document.
You can use any device to save the output. Jpeg, PNG, PDF, TIFF, etc.. From an R script, you can run the code and then the image will appear within the "Plots" area. Simply click on export and you will be able to save the image.

### 1.8.10.5 Conclusion

As you can see, visualisations in R using the *Tidyverse* provide you with many options and you can explore these further.

See here for a full list of geoms. This will help you in thinking about visualisation.

See extensions to ggplot2 here for additional plugins to enhance plots.

## 1.9 PhonR

### 1.9.1 Intro

These two packages allow to do various things with phoneticy data. PhonR former allows to normalise data and generate nice plots; the Vowels is a suite of normalisation procedures and plotting functions (it used to be an online platform, but this is not working anymore! see here or here). See help for Vowels.

We'll heavily use PhonR as it provides standard and many advanced options for normalisation and plotting, etc..

I am borrowing many examples from the PhonR link. You can adapt the codes below to your own data

```r
data(indoVowels) # provides dataframe 'indo'
str(indo)
```

```r
## 'data.frame': 1725 obs. of 5 variables:
## $ subj : chr "F02" "F02" "F02" "F02" ...
## $ gender: Factor w/ 2 levels "f","m": 1 1 1 1 1 1 1 1 1 1 ...
## $ vowel : Factor w/ 5 levels "a","e","i","o",..: 1 1 1 1 1 1 1 1 1 1 ...
## $ f1   : int 1052 904 606 652 1129 968 1124 924 602 607 ...
## $ f2   : int 1800 1911 1829 1901 1863 1699 1882 1876 1815 1950 ...
```

### 1.9.2 Normalisation

```r
#bark
with(indo, normBark(cbind(f1, f2))) %>
head(20)
```

```r
##         f1      f2
## [1,] 8.833918 12.30457
```
## [2,] 7.932374 12.70532
## [3,] 5.801590 12.41154
## [4,] 6.162236 12.67016
## [5,] 9.268799 12.53488
## [6,] 8.333415 11.91881
## [7,] 9.241219 12.60285
## [8,] 8.059612 12.58146
## [9,] 5.769617 12.36011
## [10,] 5.809568 12.84072
## [11,] 7.822817 11.93057
## [12,] 7.803347 11.68344
## [13,] 7.646099 11.80408
## [14,] 7.803347 11.64686
## [15,] 5.769617 12.31943
## [16,] 8.370116 11.80408
## [17,] 5.656923 12.76815
## [18,] 7.593087 11.85573
## [19,] 5.057221 11.67127
## [20,] 6.368064 10.73797

#erb

with(indo, normErb(cbind(f1, f2))) %>%
  head(20)

## f1 f2
## [1,] 16.00664 20.28137
## [2,] 14.86545 20.77643
## [3,] 12.02853 20.41328
## [4,] 12.52702 20.73290
## [5,] 16.54921 20.56558
## [6,] 15.37622 19.80679
## [7,] 16.51493 20.64962
## [8,] 15.02810 20.62316
## [9,] 11.98390 20.34983
## [10,] 12.03966 20.94429
## [11,] 14.72493 19.82124
## [12,] 14.69990 19.51791
## [13,] 14.49727 19.66592
## [14,] 14.69990 19.47305
## [15,] 11.98390 20.29968
## [16,] 15.42270 19.66592
## [17,] 11.82595 20.85428
## [18,] 14.42874 19.72933
## [19,] 10.96836 19.50298
## [20,] 12.80763 18.35952
#mel

with(indo, normMel(cbind(f1, f2))) %>%
  head(20)

<table><thead><tr><th>##</th><th>f1</th><th>f2</th></tr></thead><tbody><tr><td>## [1,]</td><td>1033.9416</td><td>1434.625</td></tr><tr><td>## [2,]</td><td>934.4759</td><td>1483.584</td></tr><tr><td>## [3,]</td><td>702.8431</td><td>1447.623</td></tr><tr><td>## [4,]</td><td>741.8551</td><td>1479.260</td></tr><tr><td>## [5,]</td><td>1082.4152</td><td>1462.673</td></tr><tr><td>## [6,]</td><td>978.5693</td><td>1388.149</td></tr><tr><td>## [7,]</td><td>1079.3300</td><td>1470.997</td></tr><tr><td>## [8,]</td><td>948.4413</td><td>1468.375</td></tr><tr><td>## [9,]</td><td>699.3861</td><td>1441.367</td></tr><tr><td>## [10,]</td><td>703.7057</td><td>1500.294</td></tr><tr><td>## [11,]</td><td>922.4677</td><td>1389.558</td></tr><tr><td>## [12,]</td><td>920.3353</td><td>1360.086</td></tr><tr><td>## [13,]</td><td>903.1289</td><td>1374.443</td></tr><tr><td>## [14,]</td><td>920.3353</td><td>1355.743</td></tr><tr><td>## [15,]</td><td>699.3861</td><td>1436.427</td></tr><tr><td>## [16,]</td><td>982.6160</td><td>1374.443</td></tr><tr><td>## [17,]</td><td>687.2023</td><td>1491.327</td></tr><tr><td>## [18,]</td><td>897.3346</td><td>1380.607</td></tr><tr><td>## [19,]</td><td>622.3732</td><td>1358.640</td></tr><tr><td>## [20,]</td><td>764.1399</td><td>1249.344</td></tr></tbody></table>

#lobanov

with(indo, normLobanov(cbind(f1, f2), group=subj)) %>%
  head(20)

<table><thead><tr><th>##</th><th>f1</th><th>f2</th></tr></thead><tbody><tr><td>## 1</td><td>1.90935053</td><td>0.014251925</td></tr><tr><td>## 2</td><td>1.24911468</td><td>0.180502189</td></tr><tr><td>## 3</td><td>-0.08027913</td><td>0.057686678</td></tr><tr><td>## 4</td><td>0.12492931</td><td>0.165524688</td></tr><tr><td>## 5</td><td>2.25285162</td><td>0.108610183</td></tr><tr><td>## 6</td><td>1.53462208</td><td>-0.137020839</td></tr><tr><td>## 7</td><td>2.23054635</td><td>0.137067435</td></tr><tr><td>## 8</td><td>1.33833574</td><td>0.128080934</td></tr><tr><td>## 9</td><td>-0.09812334</td><td>0.036718177</td></tr><tr><td>## 10</td><td>0.09897984</td><td>0.409751192</td></tr><tr><td>## 11</td><td>1.73839969</td><td>0.024780619</td></tr><tr><td>## 12</td><td>1.72083448</td><td>-0.071462024</td></tr><tr><td>## 13</td><td>1.58031278</td><td>-0.024893003</td></tr><tr><td>## 14</td><td>1.72083448</td><td>-0.085432730</td></tr><tr><td>## 15</td><td>0.06970449</td><td>0.183115290</td></tr></tbody></table>
## 16 2.24779086 -0.024893003
## 17 -0.01226651 0.377152877
## 18 1.53347221 -0.004713094
## 19 0.51288465 0.309606914
## 20 1.88650426 -0.100056446

#lognormed

with(indo, normLog(cbind(f1, f2))) %>%
  head(20)

## f1 f2
## [1,] 3.022016 3.255273
## [2,] 2.956168 3.281261
## [3,] 2.782473 3.262214
## [4,] 2.814248 3.278982
## [5,] 3.052694 3.270213
## [6,] 2.985875 3.230193
## [7,] 3.050766 3.274620
## [8,] 2.965672 3.273233
## [9,] 2.779596 3.258877
## [10,] 2.783189 3.290035
## [11,] 2.947924 3.230960
## [12,] 2.946452 3.214844
## [13,] 2.934498 3.222716
## [14,] 2.946452 3.212454
## [15,] 2.779596 3.256237
## [16,] 2.988559 3.222716
## [17,] 2.769377 3.285332
## [18,] 2.930440 3.226084
## [19,] 2.712650 3.214049
## [20,] 2.831870 3.152594

# using the convenience wrapper:

#bark

with(indo, normVowels('bark', f1=f1, f2=f2)) %>%
  head(20)

## f1 f2
## [1,] 8.833918 12.30457
## [2,] 7.932374 12.70532
## [3,] 5.801590 12.41154
## [4,] 6.162236 12.67016
## [5,] 9.268799 12.53488
## [6,] 8.333415 11.91881
## [7,] 9.241219 12.60285
## [8,] 8.059612 12.58146
## [9,] 5.769617 12.36011
## [10,] 5.809568 12.84072
## [11,] 7.822817 11.93057
## [12,] 7.803347 11.68344
## [13,] 7.646099 11.80408
## [14,] 7.803347 11.64686
## [15,] 5.769617 12.31943
## [16,] 8.370116 11.80408
## [17,] 5.656923 12.76815
## [18,] 7.593087 11.85573
## [19,] 5.057221 11.67127
## [20,] 6.368064 10.73797

```
#wattfab
with(indo, normVowels('wattfabricius', f1=f1, f2=f2,
    vowel=vowel, group=subj)) %>
head(20)
```

```
## f1 f2
## [1,] 1.923869 0.9674299
## [2,] 1.653211 1.0270880
## [3,] 1.108236 0.9830162
## [4,] 1.192360 1.0217134
## [5,] 2.064684 1.0012899
## [6,] 1.770252 0.9131463
## [7,] 2.055541 1.0115017
## [8,] 1.689786 1.0082769
## [9,] 1.100921 0.9754918
## [10,] 1.122104 1.1089047
## [11,] 1.639714 0.9678748
## [12,] 1.634168 0.9326173
## [13,] 1.589801 0.9496774
## [14,] 1.634168 0.9274993
## [15,] 1.112861 1.0258790
## [16,] 1.800542 0.9496774
## [17,] 1.086980 1.0969627
## [18,] 1.575012 0.9570701
## [19,] 1.241886 1.0669316
## [20,] 1.634187 0.9261513
```

### 1.9.3 Calculating vowel space area

```
poly.area <- with(indo, vowelMeansPolygonArea(f1, f2, vowel, poly.order = c("i",
"e", "a", "o", "u"), group = subj))
hull.area <- with(indo, convexHullArea(f1, f2, group = subj))
rbind(poly.area, hull.area)
```
## F02 F04 F08 F09 M01 M02 M03
## poly.area 485051.4 337364.0 434816 302064.9 197746.1 229501.7 215713.3
## hull.area 1254575.0 866109.5 1020835 751327.0 517212.5 666246.0 477518.5
## M04
## poly.area 177131.1
## hull.area 568364.0

with(indo, plotVowels(f1, f2, var.sty.by = vowel, var.col.by = vowel, pretty = TRUE, legend.kwd =

1.9.4 The plotVowels function

<figure><img src="image_28.png" /></figure>

with(indo, plotVowels(f1, f2, var.sty.by = vowel, var.col.by = gender, pretty = TRUE, legend.kwd =
<figure><img src="image_29.png" /></figure>

1.9.5 Token means

with(indo, plotVowels(f1, f2, vowel, group = gender, plotAPP = FALSE, plot.means = TRUE, var.col.by = gender, var.sty.by = vowel, pretty = TRUE, xlim = c(2900, 600), ylim = c(1000, 250), legend.kwd = "bottomleft"))

<figure><img src="image_30.png" /></figure>
1.9.6 Ellipses and areas

<figure><img src="image_31.png" /></figure>

#areas

with(indo, plotVowels(f1, f2, vowel, group = subj, plotAPP = FALSE, plot.means = TRUE,
 pch.means = vowel, cex.means = 2, var.col.by = vowel, var.sty.by = subj, ellipse line = TRUE,
 ellipse.fill = TRUE, fill opacity = 0.1, pretty = TRUE, legend.kwd = "bottomright"))
<figure><img src="image_32.png" /><figcaption>with(indo, plotVowels(f1, f2, vowel, group = subj, var.col.by = subj, var.sty.by = subj, plot_tokens = FALSE, plot.means = TRUE, pch.means = vowel, cex.means = 2, pretty = TRUE, poly.line = TRUE, poly.order = c("i", "e", "a", "o", "u"), legend.kwd = "bottomleft", legend.args = list(segmentLen = 3, cex = 1.2, lwd = 2), col = c("gray70", "gray40", "gray10"), lty = c("solid", "solid", "dotdash", "dotdash", "dotted", "longdash", "dashed", "dotted")))</figcaption></figure>
<figure><img src="image_33.png" /></figure>

1.10 session info

sessionInfo()

## R version 4.5.2 (2025-10-31 ucrt)
## Platform: x86_64-w64-mingw32/x64
## Running under: Windows 11 x64 (build 26200)
##
## Matrix products: default
## LAPACK version 3.12.1
##
## locale:
## [1] LC_COLLATE=English_United Kingdom.utf8
## [2] LC_CTYPE=English_United Kingdom.utf8
## [3] LC_MONETARY=English_United Kingdom.utf8
## [4] LC_NUMERIC=C
## [5] LC_TIME=English_United Kingdom.utf8
##
## time zone: Europe/Paris
## tzcode source: internal
##
## attached base packages:
## [1] stats graphics grDevices utils datasets methods base
##
## other attached packages:
## [1] cardx_0.3.0 gtsummary_2.4.0 ggdist_3.3.3 summarytools_1.1.4
## [5] phonR_1.0-7 languageR_1.6 lubridate_1.9.4 forcats_1.0.1
## [9] stringr_1.6.0 dplyr_1.1.4 purrr_1.2.0 readr_2.1.6
## [13] tidyr_1.3.1 tibble_3.3.0 ggplot2_4.0.1 tidyverse_2.0.0

##

## loaded via a namespace (and not attached):

## [1] tidyselect_1.2.1 farver_2.1.2 S7_0.2.1
## [4] fastmap_1.2.0 digest_0.6.39 timechange_0.3.0
## [7] lifecycle_1.0.4 magrittr_2.0.4 compiler_4.5.2
## [10] rlang_1.1.6 tools_4.5.2 utf8_1.2.6
## [13] yaml_2.3.10 gt_1.1.0 knitr_1.50
## [16] labeling_0.4.3 sp_2.2-0 plyr_1.8.9
## [19] xml2_1.5.0 RColorBrewer_1.1-3 withr_3.0.2
## [22] grid_4.5.2 scales_1.4.0 MASS_7.3-65
## [25] tinytex_0.58 cli_3.6.5 r markdown_2.30
## [28] crayon_1.5.3 reformulas_0.4.2 generics_0.1.4
## [31] rstudioapi_0.17.1 reshape2_1.4.5 tzdb_0.5.0
## [34] commonmark_2.0.0 minqa_1.2.8 pander_0.6.6
## [37] splines_4.5.2 matrixStats_1.5.0 base64enc_0.1-3
## [40] vctrs_0.6.5 boot_1.3-32 Matrix_1.7-4
## [43] bookdown_0.45 litedown_0.8 hms_1.1.4
## [46] rapportools_1.2 magick_2.9.0 splancs_2.01-45
## [49] glue_1.8.0 nloptr_2.2.1 codetools_0.2-20
## [52] distributional_0.5.0 stringi_1.8.7 gtable_0.3.6
## [55] quadprog_1.5-8 lme4_1.1-37 pillar_1.11.1
tcltk_4.5.2
## [58] htmltools_0.5.8.1 R6_2.6.1 lattice_0.22-7
## [61] Rdpack_2.6.4 rbibutils_2.4 backports_1.5.0
## [64] markdown_2.0 broom_1.0.10 pryr_0.1.6
## [67] cards_0.7.0 nlme_3.1-168 checkmate_2.3.3
## [70] Rcpp_1.1.0 xfun_0.54 fs_1.6.6
## [73] mgcv_1.9-4
## [76] pkgconfig_2.0.3
# Part II

## Statistical analyses using R
<figure><img src="image_34.png" /></figure>
# Chapter 2

## Correlation, LM, GLM, STD and CLM

This chapter introduces you to the basics of inferential statistics. We cover correlation tests to identify specific patterns in the data. We then use linear models to evaluate group differences, and generalised linear models to evaluate binary outcomes. We also introduce the basics of signal detection theory (SDT) to evaluate the accuracy of our models. At the end of this chapter, we introduce the basics of cumulative link models (CLM) to evaluate ordinal (likert) outcomes on two datasets.

### 2.1 Loading packages

```r
## Use the code below to check if you have all required packages installed. If some are not installed, install them.
requiredPackages = c('tidyverse', 'broom', 'knitr', 'Hmisc', 'corrplot', 'emmeans', 'ggsignif', 'library(p, character_only = TRUE)
}
```

### 2.2 Inferential statistics

#### 2.2.1 Introduction

What do we mean by inferential statistics?

We want to evaluate if there are differences observed between two groups of datapoints and if these differences are statistically significant.
For this, we need to pay attention to the format of our outcome and predictors. An outcome is the response variable (or dependent variable); a predictor is our independent variable(s). This is fed into our equation

$$y = x\beta + \varepsilon$$

where:

* $y \to$ outcome (DV) $\Rightarrow$ known

* $x \to$ fixed effect (IV) $\Rightarrow$ known

* $\beta \to$ coefficient of fixed effect $\Rightarrow$ unknown

* $\varepsilon \to$ random error term $\Rightarrow$ unknown

### 2.2.2 Hypothesis testing

In the approach that this book follows, we use a frequentist approach to hypothesis testing. By frequentist we mean when modelling the relationship between an outcome and a predictor, we are interested in testing if the predictor has a clear and statistically significant influence on the outcome. When evaluating this, we try to reject the null hypothesis ($H_0$) in favour of the alternative hypothesis ($H_1$).

The null hypothesis ($H_0$) is the hypothesis we make that there is no statistical evidence in favour of an effect nor of a difference between the groups. The alternative hypothesis ($H_1$) is the hypothesis we make to reject the null hypothesis and there is indeed a statistical evidence in favour of an effect or a difference between the groups.

### 2.2.3 Type of outcomes

It is important to know the class of the outcome before doing any pre-data analyses or inferential statistics. Outcome classes can be one of:

1. **Numeric**: As an example, we have length/width of leaf; height of mountain; fundamental frequency of the voice; etc. These are generally modelled with a **Gaussian** distribution. These are **true** numbers and we can use summaries, t-tests, linear models, etc. Integers are a family of numeric variables and can still be considered as a normal numeric variable. Count data are numbers related to a category. In general, these are related to the second type of outcomes, and should be analysed using a poisson logistic regression

2. **Categorical** (Unordered): Observations for two or more categories. As an example, we can have gender of a speaker (male or female); responses to a True vs False perception tests; Colour (unordered) categorisation, e.g., red, blue, yellow, orange, etc.. For these we can use a Generalised Linear Model (binomial or multinomial) or a simple chi-square test.
3. **Categorical** (Ordered): When you run a rating experiment, where the outcome is either **numeric** (i.e., 1, 2, 3, 4, 5) or **categories** (i.e., disagree, neutral, agree). The **numeric** option is NOT a true number as for the participant, these are categories. Cumulative Logit models (or Generalised Linear Model with a cumulative function) are used. The mean is meaningless here, and the median is a preferred descriptive statistic.

### 2.2.4 Types of Errors?

In the literature, we discuss various types of errors in statistical models. There are four different types, and the excellent Sonderegger (2023) book (Sonderegger, M. (2023). *Regression Modeling for Linguistic Data*. The MIT Press.) summarises these nicely. Here are the four types:

1. **Type I** (or “false positive” or error $\alpha$) $\Rightarrow$ falsely concluding there is an effect when none exists. (generally due to inaccurate modelling strategies)

2. **Type II** (or “false negative” or error $\beta$)) $\Rightarrow$ falsely concluding there is no effect when one in fact exists (generally due to inaccurate modelling strategies)

<table><thead><tr><th colspan="4">Statistical analysis result (sample)</th></tr><tr><th></th><th>Reject H<sub>0</sub></th><th>Don't reject H<sub>0</sub></th><th></th></tr></thead><tbody><tr><th>Reality<br>(population)</th><td>H<sub>0</sub> is true<br>H<sub>0</sub> is false</td><td>Type I error (α)<br>Correct decision (significant)</td><td>Correct decision (null result)<br>Type II error (β)</td></tr></tbody></table>

3. **Type S** $\Rightarrow$ Inaccurate sign (generally due to *hidden* multicollinearity and low power = 1- $\beta$)

4. **Type M** $\Rightarrow$ Inaccurate magnitude (generally due to *hidden* multicollinearity and low power = 1- $\beta$)

It is important to note that the **Type I** and **Type II** errors are the most commonly discussed in the literature and are related to each other. Modelling is an art, and depending on model specifications, we can easily increase **Type I** and/or **Type II** errors. **Type I** is usually increased in the case of using inaccurately specified model, e.g., not using a random effects structure when the data comes from multiple speakers and/or multiple items. **Type II** error is usually increased when there is an over specification of the model and/or not accounting for dependencies in the data. **Type S** and **Type M** are generally caused due to collinearity in the data, which causes instability in the model and leads to cases of suppression; a point discussed in **chapter 8**.
## 2.3 Effect size and power analysis

### 2.3.1 Effect size

#### 2.3.1.1 What is the purpose of Effect Size?

‘Effect size’ is simply a way of quantifying the size of the difference between two groups.

There are two family groups of Effect Size:

1. Standardised mean differences (also called group difference indices; e.g., Cohen’s *d*), The mean difference Effect Sizes are computed when the focus of the statistical analysis uses mean outcome scores to compare potential group differences.

2. Strength of association (also referred to as relationship variance-accounted-for or variance explained Effect Sizes; e.g., *R*², η², η²*p*, ω², etc. Strength of association Effect Size indices are calculated across a wide variety of research designs (e.g., ex post facto or causal comparative, experimental, quasi-experimental, prediction [regression and correlational]) using the general linear model (GLM) for the statistical analyses. Strength of association Effect Sizes indicate how much of the variance associated with the dependent variable (e.g., students’ English test scores) can be accounted for (or explained) by the independent variable(s) (e.g., gender).

Both Effect Size families can be further categorised into biased (uncorrected) and unbiased (corrected) measures:

* Biased measures relates to **"sampling error"** and how the Effect Size derived from statistical procedures computed with smaller groups of participants (samples) drawn from a larger population will be higher than would be found either in the original population or in participant samples studied at a later date.

For example, assume a school counselling researcher randomly selects a sample of 15 girls and 15 boys from a 10th-grade physical education class in one high school. This 30-person sample then is assumed to represent the entire population of male and female students attending all Grade 10 physical education classes across a district's five high schools. Because there is always sampling error (i.e., problems with selecting the students to participate in the research study), uncorrected Effect Sizes (e.g., Cohen's *d*, *R*², η²) will be higher (i.e., positively biased or overestimated) for the 30 students than if the researcher actually could have included in the study every 10th-grade physical education student across all five high schools.

* Unbiased ones (e.g., adjusted *R*², Hays' ω², etc.) statistically correct for this positive bias by calculating the strength of association indices so that they better reflect the true Effect Sizes for the entire population and those calculated in future samples.
The issue is that because these formulas must adjust for sampling error occurring both in the present study and in subsequent investigations, Effect Size estimates for future research samples result in more shrinkage (i.e., are smaller or more conservative) than those Effect Size estimates based on an entire population. Therefore, interpretation of Effect Size s is often less than straightforward and must be exercised with caution.

### 2.3.1.2 What is the relationship between Effect Size and Significance?

Effect size quantifies the size of the difference between two groups, and may therefore be said to be a true measure of the significance of the difference.

The statistical significance is usually calculated as a ‘p-value’, the probability that a difference of at least the same size would have arisen by chance, even if there really were no difference between the two populations. For differences between the means of two groups, this p-value would normally be calculated from a ‘t-test’. By convention, if $p < 0.05$ (i.e. below 5%), the difference is taken to be large enough to be ‘significant’; if not, then it is ‘not significant’.

There are a number of problems with using ‘significance tests’ in this way. The main one is that the p-value depends essentially on two things:

a) The size of the effect

b) The size of the sample.

One would get a ‘significant’ result either if the effect were very big (despite having only a small sample) or if the sample were very big (even if the actual effect size were tiny). It is important to know the statistical significance of a result, since without it there is a danger of drawing firm conclusions from studies where the sample is too small to justify such confidence. However, statistical significance does not tell you the most important thing: **the size of the effect**. One way to overcome this confusion is to report the effect size, together with an estimate of its likely ‘margin for error’ or ‘confidence interval’.

### 2.3.1.3 What is Confidence Interval?

The confidence interval is an alternative way to indicate the variability in estimates from small samples. The default calculation here is a ‘95% confidence interval’. If multiple samples of two groups of the same size as these, taken from a population in which the true difference was the value in column J, there would be variation in the differences found. However, for every 100 samples taken, for 95 of them (on average) the difference would be between the lower and upper confidence limits. The confidence interval is usually interpreted as a ‘margin of uncertainty’ around the estimate of the difference between experimental and control groups.
### 2.3.1.4 Which Effect Size measure to use?

We highlighted above various effect size measures, e.g., Cohen's d, R², η², η²p, ω². Which one to choose from? Well, all depends on the statistical approach you are using and whether you want to report on standardised or unstandardised measures.

However, as a rule of thumb, you can do the following:

1. In general, if we want to talk about the sample we have, we can use any biased method (like the eta-squared) because it estimates only for the sample and not the population.

2. If we want to generalise our results, we have to use an unbiased method.

Literature recommends using an Effect Size measure accompanied with the Confidence Interval, because telling only that a result is statistically significant doesn't inform us of the contribution of the independent variable to explain the variation in the dependent variable.

### 2.3.1.5 How to compute these in R

We use the package `effectsize`

#### 2.3.1.5.1 Standardised effect sizes (Cohen's d)
It can be computed using the following formula:

$$d = \frac{Mean_{factor1} - Mean_{factor2}}{SD_{pooled}}$$

With the pooled standard deviation computed as:

$$SD_{pooled} = \sqrt{\frac{SD_{factor1}^2 * (Number_{factor1} - 1) + SD_{factor2}^2 * (Number_{factor2} - 1)}{Number_{factor1} + Number_{factor2} - 2}}$$

Fortunately, instead of calculating it manually, we can use the function `cohens_d()` from the package `effectsize`.

```
options(es.useSymbols = TRUE) # get nice symbols when printing! (On Windows, requires .R symbols)
```

```
english %>%
  lm(RTlexdec ~ AgeSubject, data = .) %>%
  summary()
```

#### 2.3.1.5.1.1 Linear model

```
##
```

```
## Call:
```

```
## lm(formula = RTlexdec ~ AgeSubject, data = .)
```

```
##
```

```
## Residuals:
```
```
## Min      1Q   Median      3Q   Max
## -0.25776 -0.08339 -0.01669  0.06921  0.52685
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept)   6.660958  0.002324 2866.44 <2e-16 ***
## AgeSubjectyoung -0.221721  0.003286 -67.47 <2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.1111 on 4566 degrees of freedom
## Multiple R-squared: 0.4992, Adjusted R-squared: 0.4991
## F-statistic: 4552 on 1 and 4566 DF, p-value: < 2.2e-16
```

What the model above tells us is that there is a very highly statistically significant impact of AgeSubject on RTlexdec (p < 0.0001).

### 2.3.1.5.1.2 Cohen's d

However, it does not tell us how large this effect is. That is where we can use the `cohens_d()` function.

```
# Cohen's d
english %>%
  cohens_d(RTlexdec ~ AgeSubject, data = .) %>%
  interpret_cohens_d()
```

```
## Cohen's d |      95% CI | Interpretation
## --------------------------
## 2.00 | [1.93, 2.07] |      large
##
## - Estimated using pooled SD.
## - Interpretation rule: cohen1988
```

One nice option in the package `effectsize` is to use an *Automated Interpretation of Indices of Effect Size* vignette. We can use the original publication to interpret the effect size or use an alternative. See the link above for more details.

We can also standardise our model's coefficients, which allows us to compare the effect sizes of different predictors in the same model and between models, variables and studies.

```
model <- english %>%
  lm(RTlexdec ~ AgeSubject, data = .)

datawizard::standardize(model)
```

### 2.3.1.5.1.3 Standardised coefficients
##
## Call:
## lm(formula = RTlexdec ~ AgeSubject, data = data_std)
##
## Coefficients:
## (Intercept) AgeSubjectyoung
## 0.7065 -1.4130

parameters::standardize_parameters(model)

## # Standardization method: refit
##
## Parameter | Std. Coef. | 95% CI
##
## (Intercept) | 0.71 | [ 0.68, 0.74]
## AgeSubject [young] | -1.41 | [-1.45, -1.37]

Standardised coefficients are interpreted as the change in standard deviations of the outcome variable for a one standard deviation change in the predictor variable. We can initially standardise our outcomes using the function `scale(center = TRUE, scale = TRUE)` and then run our model, OR use the function `standardize()` from the package `datawizard` or `standardize_parameters()` from the package `parameters` to do that.

### 2.3.1.5.2 Strength of association effect sizes
Here, we saw already that there are biased and unbiased measures. We will document the classical measures used.

#### 2.3.1.5.2.1 $R^2$ $R^2$ (R squared):
It is the proportion of variance explained by the model. This is a biased effect size measure. It is computed as:

$$R^2 = \frac{SS_{model}}{SS_{total}} = 1 - \frac{SS_{residual}}{SS_{total}}$$

Where:

* $SS_{model}$ is the sum of squares of the model

* $SS_{total}$ is the total sum of squares

* $SS_{residual}$ is the sum of squares of the residuals.

We can get the $R^2$ and the adjusted $R^2$ from the summary of the model:

`summary(model)$r.squared`

`## [1] 0.49923`

`summary(model)$adj.r.squared`

`## [1] 0.4991203`
2.3.1.5.2.2 η² (eta squared) The eta-squared effect size measure is a biased measure and describes the ratio of variance explained in the dependent variable by a predictor (or factor) while controlling for other predictors (or in a simple way: it is the proportion of total variation attributable to the factor). The range is from 0 to 1 and it sums up to 1. Eta squared is a biased estimator of the variance explained by the model in the population (it only estimates effect size in the sample). On average it overestimates the variance explained in the population. It is computed as:

$$ \eta^2 = \frac{SS_{model}}{SS_{total}} $$

* $SS_{model}$ is the sum of squares of the model

* $SS_{total}$ is the total sum of squares

* $SS_{residual}$ is the sum of squares of the residuals.

In a sense, both $R^2$ and $η^2$ are computed in the same way.

To obtain it directly from the model, we can use the function `eta_squared()` from the package `effectsize` and specify `partial = FALSE`.

```r
# eta squared
english %>%
  aov(RTlexdec ~ AgeSubject, data = .) %>%
  eta_squared(. , partial = FALSE) %>%
  interpret_eta_squared()
```

```r
## # Effect Size for ANOVA (Type I)
```

```r
##
```

```r
## Parameter | ² | 95% CI | Interpretation
```

```r
## ----------------------------------------------------------------------
```

```r
## AgeSubject | 0.50 | [0.48, 1.00] | large
```

```r
##
```

```r
## - One-sided CIs: upper bound fixed at [1.00].
```

```r
## - Interpretation rule: field2013
```

2.3.1.5.2.3 η²p (partial eta squared) η²p (partial eta squared): “proportion of total variation attributable to the factor, partiallying out (excluding) other factors from the total non error variation”. Partial eta squared is normally higher than eta squared (except in simple one-factor models). The range is from 0 to 1, and it sums up to more than 1 (except in one-way ANOVA). It is computed as:

$$ \eta^2 p = \frac{SS_{factor}}{SS_{factor} + SS_{error}} $$

Where:

* $SS_{factor}$ is the sum of squares of the factor

* $SS_{error}$ is the sum of squares of the error.
To obtain it directly from the model, we can use the function `eta()`
from the package `effectsize` and specify `partial = TRUE`.

```r
# partial eta squared
english %>%
  aov(RTlexdec ~ AgeSubject, data = .) %>%
  etaüstred(., partial = TRUE) %>%
  interpret_etaüstred()
```

## For one-way between subjects designs, partial eta squared is equivalent
to eta squared. Returning eta squared.

## # Effect Size for ANOVA
##
## Parameter | % | 95% CI | Interpretation
##
## AgeSubject | 0.50 | [0.48, 1.00] | large
##
## - One-sided CIs: upper bound fixed at [1.00].
## - Interpretation rule: field2013

2.3.1.5.2.4 $\omega^2$ (omega squared) $\omega^2$ (omega squared): Provides a relatively unbiased estimate of the variance explained in the population by a predictor variable. It takes random error into account more so than eta squared, which is incredibly biased to be too large. The range is from 0 to 1 and it sums up to 1.

It is computed as:

$$ \omega^2 = \frac{SS_{model} - (df_{model} * MS_{error})}{SS_{total} + MS_{error}} $$

Where:

* $SS_{model}$ is the sum of squares of the model

* $SS_{total}$ is the total sum of squares

* $df_{model}$ is the degrees of freedom of the model

* $MS_{error}$ is the mean square of the error.

To obtain it directly from the model, we can use the function `omega()`
from the package `effectsize` and specify `partial = FALSE`.

```r
# omega squared
english %>%
  aov(RTlexdec ~ AgeSubject, data = .) %>%
  omega_squared(., partial = FALSE) %>%
  interpret_omega_squared()
```

## # Effect Size for ANOVA (Type I)
##
Parameter | ² | 95% CI | Interpretation
--- | --- | --- | ---
AgeSubject | 0.50 | [0.48, 1.00] | large

One-sided CIs: upper bound fixed at [1.00].
Interpretation rule: field2013

### 2.3.1.5.2.5 Partial ω² (partial omega squared)
Partial ω² (omega squared): “proportion of total variation attributable to the factor, partialling out (excluding) other factors from the total non error variation”. The range is from 0 to 1 and it sums up to more than 1 (except in one-way ANOVA). It is the equivalent of partial eta-squared but as an unbiased measure. It is computed as:

$$ω^2 p = \frac{df_{factor} * MS_{factor} - MS_{error}}{df_{factor} * MS_{factor} + df_{total} - df_{factor} * MS_{error}}$$

Where:

* $df_{factor}$ is the degrees of freedom of the factor

* $MS_{factor}$ is the mean square of the factor

* $df_{total}$ is the total degrees of freedom

* $MS_{error}$ is the mean square of the error.

To obtain it directly from the model, we can use the function `omega_squared()` from the package `effectsize` and specify `partial = TRUE`.

```r
# omega squared partial
english %>%
  aov(RTlexdec ~ AgeSubject, data = .) %>%
  omega_squared(., partial = TRUE) %>%
  interpret_omega_squared()
```

For one-way between subjects designs, partial omega squared is equivalent to omega squared. Returning omega squared.

### # Effect Size for ANOVA

Parameter | ² | 95% CI | Interpretation
--- | --- | --- | ---
AgeSubject | 0.50 | [0.48, 1.00] | large

One-sided CIs: upper bound fixed at [1.00].
Interpretation rule: field2013

### 2.3.1.5.2.6 Odds ratio
We will see in the sections 2.6 and 3.4 that we use odds ratio as an exponentiation of the logodds (logit) coefficients. Odds ratios can be interpreted as the odds for an event to occur. In comparison to
the logodds (logit) scale, which spans -6 and 6 (well $-\infty$ - $\infty$), odds ratio are very (very high) and can reach 999 for a logodds of 6.9. Of course, in some cases, odds ratio are not easy to interpret, and it is often “easier” to transform the logodds into proportions. See this excellent link where you can understand the links better.

## 2.3.2 Power analysis and sample size

Usually, when we design an experiment, we always ask the question: what is the optimal sample size to use so that I have confidence in the results I will obtain? To answer this question, we usually say: well it all depends on your discipline and what is feasible for a particular study. However, this can easily lead to misconceptions and issues in designing our studies and more importantly, in the interpretation of the results!

Hence, we demonstrate in the following power analysis and the sample size. These two are crucial as they are interlinked. We saw in **section 2.2.4** that there are four types of errors, with **Type M** and **Type S** depending heavily on the power and on the sample size.

We use the package *pwr*, which allows us to compute both power and sample size (but not at the same time!!). Looking at the vignette of the *pwr* package, they cite the various functions available, such as:

* *pwr.p.test*: one-sample proportion test

* *pwr.2p.test*: two-sample proportion test

* *pwr.2p2n.test*: two-sample proportion test (unequal sample sizes)

* *pwr.t.test*: two-sample, one-sample and paired t-tests

* *pwr.t2n.test*: two-sample t-tests (unequal sample sizes)

* *pwr.anova.test*: one-way balanced ANOVA

* *pwr.r.test*: correlation test

* *pwr.chisq.test*: chi-squared test (goodness of fit and association)

* *pwr.f2.test*: test for the general linear model

### 2.3.2.1 Cohen's *d* interpretation rules and *t.test*

Based on J. Cohen (1988) rule of thumb, we can interpret Cohen's *d* effect size measure as:

* *d* < 0.2 - Very small

* 0.2 <= *d* < 0.5 - Small

* 0.5 <= *d* < 0.8 - Medium

* *d* >= 0.8 - Large

We can use the function *pwr.t.test()* to compute the power of a *t.test*.
### 2.3.2.1.1 pwr.t.test - one-sample and two-sample t tests for means

Based on the example below, we know that the variance between the two age groups with respect to RTlexdec is significantly different; meaning that there is no equality of variance in this case. This is important to know to be able to compute the t.test appropriately.

Here are the results of the t.test between the two groups:

```
english %>
    t.test(RTlexdec ~ AgeSubject, data = ., var equal = FALSE)
```

##

## Welch Two Sample t-test

##

## data: RTlexdec by AgeSubject

## t = 67.468, df = 4534.6, p-value < 2.2e-16

## alternative hypothesis: true difference in means between group old and group young is not equal

## 95 percent confidence interval:

## 0.2152787 0.2281642

## sample estimates:

## mean in group old mean in group young

## 6.660958 6.439237

```
english %>
    cohens_d(RTlexdec ~ AgeSubject, data = .) %>
    interpret_cohens_d()
```

### 2.3.2.1.2 Cohen's d effect size

## Cohen's d | 95% CI | Interpretation

## ------------------------------------------

## 2.00 | [1.93, 2.07] | large

##

## - Estimated using pooled SD.

## - Interpretation rule: cohen1988

## cohen.ES(test = "t", size = "large")

##

## Conventional effect size from Cohen (1982)

##

## test = t

## size = large

## effect.size = 0.8

## pwr.t.test(n = 2284, d = 2, sig.level = 2.2e-16)

##
## Two-sample t test power calculation
##
## n = 2284
## d = 2
## sig.level = 2.2e-16
## power = 1
## alternative = two.sided
##
## NOTE: n is number in *each* group

The results here show that with a sample size of 2284, a large effect size of 2 and a significance level of 2.2e-16, we have a power of 1 (or 100%). This is a massive power and effect size, but we know that the results are inflated de to not using a mixed effects model.

Now, let's imagine we use another sample.

We're interested to know if there is a difference in the mean reaction time in a lexical decision task between young and old participants. Let's say we randomly ask 30 young and 30 old participants to perform the test and calculate the mean reaction time for each group. We'll test for a difference in means using a two-sample t-test.

How powerful is this experiment if we want to detect a "medium" effect in either direction with a significance level of 0.001?

Unfortunately, we don't have any idea what the number of participants in the *english* dataset. It is composed of 2284 monomorphemic English nouns and verbs that are averaged across old and young subjects.

cohen.ES(test = "t", size = "large")

##
## Conventional effect size from Cohen (1982)
##
## test = t
## size = large
## effect.size = 0.8

pwr.t.test(n = 30, d = 0.8, sig.level = 0.01)

##
## Two-sample t test power calculation
##
## n = 30
## d = 0.8
## sig.level = 0.01
## power = 0.6677563
## alternative = two.sided
###

### NOTE: n is number in *each* group

The power is only at 66.78%! Not very powerful. How many students should we observe for a test with 80% power?

pwr.t.test(d = 0.8, power = 0.80, sig.level = 0.01)

###

### Two-sample t test power calculation

###

### n = 38.18831

### d = 0.8

### sig.level = 0.01

### power = 0.8

### alternative = two.sided

###

### NOTE: n is number in *each* group

We need at least 38 students per group. If we lower our p value, the sample size will increase drastically

### 2.3.2.1.3 Plots

The `pwr` package comes with a handy plotting function that we can use to evaluate the impact of the increase in power, effect size and sample size for a specific p value

#### 2.3.2.1.3.1 p < 0.05

What if we use a specific p value < 0.05 and a power of 0.80? How does the change in effect size impact the sample size? We simulate variable effect size measures at 0.2, 0.5, 0.8 and 2 based on the interpretation above.

Based on J. Cohen (1988) rule of thumb, we can interpret Cohen's *d* effect size measure as:

* d < 0.2 - Very small

* 0.2 <= d < 0.5 - Small

* 0.5 <= d < 0.8 - Medium

* d >= 0.8 - Large

p.out <- pwr.t.test(d = 0.2, power = 0.80, sig.level = 0.05)

p.out

###

### Two-sample t test power calculation

###

### n = 393.4057

### d = 0.2
## sig.level = 0.05
## power = 0.8
## alternative = two.sided
##
## NOTE: n is number in *each* group
plot(p.out)

## Warning: Using `size` aesthetic for lines was deprecated in ggplot2 3.4.0.
## i Please use `linewidth` instead.
## i The deprecated feature was likely used in the pwr package.
## Please report the issue to the authors.
## This warning is displayed once every 8 hours.
## Call `lifecycle::last_lifecycle_warnings()` to see where this warning was
## generated.

Two-sample t test power calculation

<figure><img src="image_35.png" /></figure>

p.out <- pwr.t.test(d = 0.5, power = 0.80, sig.level = 0.05)
p.out

##
## Two-sample t test power calculation
##
## n = 63.76561
## d = 0.5
## sig.level = 0.05
## power = 0.8
## alternative = two.sided
Two-sample t test power calculation

<figure><img src="image_36.png" /></figure>

p.out <- pwr.t.test(d = 2, power = 0.80, sig.level = 0.05)

p.out

##
## Two-sample t test power calculation
##
## n = 5.089995
## d = 2
## sig.level = 0.05
## power = 0.8
## alternative = two.sided
##
## NOTE: n is number in *each* group
plot(p.out)
<figure><img src="image_37.png" /></figure>

### 2.3.2.1.3.2 p < 0.01
We do the same as above but we change the p value < 0.01 and a power of 0.80? How does the change in effect size impact the sample size? We simulate variable effect size measures at 0.2, 0.5, 0.8 and 2 based on the interpretation above.

Based on J. Cohen (1988) rule of thumb, we can interpret Cohen's *d* effect size measure as:

* d < 0.2 - Very small

* 0.2 <= d < 0.5 - Small

* 0.5 <= d < 0.8 - Medium

* d >= 0.8 - Large

```R
p.out <- pwr.t.test(d = 0.2, power = 0.80, sig.level = 0.01)
p.out
```

##

### Two-sample t test power calculation

##

### n = 585.6093

### d = 0.2

### sig.level = 0.01

### power = 0.8

### alternative = two.sided

##
## NOTE: n is number in *each* group

plot(p.out)

<figure><img src="image_38.png" /></figure>

p.out <- pwr.t.test(d = 0.5, power = 0.80, sig.level = 0.01)

p.out

##
## Two-sample t test power calculation
##
## n = 95.10364
## d = 0.5
## sig.level = 0.01
## power = 0.8
## alternative = two.sided

##
## NOTE: n is number in *each* group

plot(p.out)
### Two-sample t test power calculation

<figure><img src="image_39.png" /></figure>

```
p.out <- pwr.t.test(d = 0.8, power = 0.80, sig.level = 0.01)
```

```
p.out
```

```
##
## Two-sample t test power calculation
##
```

```
##
```

```
##           n = 38.18831
```

```
##           d = 0.8
```

```
##       sig.level = 0.01
```

```
##           power = 0.8
```

```
##       alternative = two.sided
```

```
##
```

```
## NOTE: n is number in *each* group
```

```
plot(p.out)
```
Two-sample t test power calculation

<figure><img src="image_40.png" /></figure>

p.out <- pwr.t.test(d = 2, power = 0.80, sig.level = 0.01)

p.out

##

Two-sample t test power calculation

##

n = 7.66057

d = 2

sig.level = 0.01

power = 0.8

alternative = two.sided

##

## NOTE: n is number in *each* group

plot(p.out)
<figure><img src="image_41.png" /></figure>

The plot above shows the sample size required for the specific power, p value and effect size. Play around with these values to have an idea of the impact of these parameters on the sample size.

### 2.3.2.2 Linear model

#### 2.3.2.2.1 General example
We borrow from the vignette of the `pwr` package an example of a linear model.

##### 2.3.2.2.1.1 Power of the test
We want to compute the power we have with a sample size of 40 words, an effect size of 0.3 and a significance level of 0.001.

```
pwr.f2.test(u = 1, v = 40 - 1 - 1, f2 = 0.3/(1 - 0.3), sig.level = 0.001)
```

```
##
## Multiple regression power calculation
##
##
## u = 1
## v = 38
## f2 = 0.4285714
## sig.level = 0.001
## power = 0.7104043
```

Where:

* u = number of predictors
* v = degrees of freedom of the residuals (n - u - 1)

* f2 = effect size (Cohen's f2)

* sig.level = significance level (alpha)

### 2.3.2.2.1.2 Sample size
We want to compute the sample size needed to have a power of 0.8, with a significance level of 0.001 and an effect size of 0.3 (medium effect size according to Cohen's conventions).

pwr.f2.test(u = 1, f2 = 0.3/(1 - 0.3), sig.level = 0.001, power = 0.8)

##
## Multiple regression power calculation
##
## u = 1
## v = 43.35263
## f2 = 0.4285714
## sig.level = 0.001
## power = 0.8

Where:

* u = number of predictors

* f2 = effect size (Cohen's f2)

* sig.level = significance level (alpha)

* power = desired power (1 - beta)

Recall $n = v + u + 1$. Therefore, we need $43.35263 + 1 + 1 = 45.35263 \sim 45$ words.

### 2.3.2.2.2 Using our LM model
In our case, we use the dataset english from the languageR package. We used AgeSubject as our predictor to predict RTlexdec. What if we want to compute the power of our model that we ran previously? Here is the summary of our model

summary(model)

##
## Call:
## lm(formula = RTlexdec ~ AgeSubject, data = .)
##
## Residuals:
##       Min     1Q   Median     3Q      Max
## -0.25776 -0.08339 -0.01669  0.06921  0.52685
##
## Coefficients:

Estimate Std. Error t value Pr(>\|t\|)
## (Intercept) 6.660958 0.002324 2866.44 <2e-16 ***
## AgeSubjectyoung -0.221721 0.003286 -67.47 <2e-16 ***
## ---
## Signif. codes: 0 '**' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.1111 on 4566 degrees of freedom
## Multiple R-squared: 0.4992, Adjusted R-squared: 0.4991
## F-statistic: 4552 on 1 and 4566 DF, p-value: < 2.2e-16

### 2.3.2.2.2.1 Power analysis from our model
Recall that the F test has numerator and denominator degrees of freedom. The numerator degrees of freedom, u, is the number of coefficients you'll have in your model (minus the intercept).

The effect size, f2, is R²/(1 - R²), where R² is the coefficient of determination, aka the “proportion of variance explained”. To determine effect size you hypothesize the proportion of variance your model explains, or the R².

Once we computed these three values, we can compute the power of our model based on its summary.

u <- summary(model)$fstatistic[2]
v <- summary(model)$fstatistic[3]
f2 <- summary(model)$r.squared/(1 - summary(model)$r.squared)
# p.value <- summary(model)$coefficients[2,4]
p.value <- 0.0001
pwrModel <- pwr.f2.test(u = u, v = v, f2 = f2, sig.level = p.value)
pwrModel

##
## Multiple regression power calculation
##
## u = 1
## v = 4566
## f2 = 0.9969247
## sig.level = 1e-04
## power = 1

### 2.3.2.2.2.2 Sample size
For the optimal sample size, we use the same values for u, f2 and p.value as above, but we specify the power we want to have. Here, we use the power we obtained from our model above.

pwrSample <- pwr.f2.test(u = u, f2 = f2, sig.level = p.value, power = pwrModel$power)
pwrSample

##
## Multiple regression power calculation
##
## u = 1
## v = 1e+09
## f2 = 0.9969247
## sig.level = 1e-04
## power = 1

Based on the output, we can compute the optimal sample size needed for our model as

n <- pwrSample$v + u + 1
n

## numdf
## 1e+09

The optimal sample size to reach this level of power is **1,000,000,000 words!**
BUT we do not have that many words in our dataset. This is kind of telling us that there are serious issues in our model. In fact, we see that we have a power of 1, which is very high power. However, we know very well that our linear model is inaccurate. We did not tell our model that the data comes from multiple words, and that we really need to compute a mixed effects regression (see **Section 3.3** from **chapter 3**).

## 2.4 Correlation tests

### 2.4.1 Correlation does not imply causation!

There is a general misconception that when one conducts a correlation test (see below), we can explain causation. An excellent example is borrowed from wikipedia. It is not because two variables are correlated (positively or negatively) that one explains the reasons behind the second!

If you want to evaluate the cause of an effect then inferential statistics that we cover after the correlation tests are the ones to go for.
<figure><img src="image_42.png" /></figure>

2.4.2 Basic correlations

Let us start with a basic correlation test. We want to evaluate if two numeric
variables are correlated with each other.

We use the function cor to obtain the pearson correlation and cor.test to run
a basic correlation test on our data with significance testing

cor(english$RTlexdec, english$RTnaming, method = "pearson")

[1] 0.7587033
cor.test(english$RTlexdec, english$RTnaming)

##
## Pearson's product-moment correlation
##
## data: english$RTlexdec and english$RTnaming
## t = 78.699, df = 4566, p-value < 2.2e-16
## alternative hypothesis: true correlation is not equal to 0
## 95 percent confidence interval:
## 0.7461195 0.7707453
## sample estimates:
## cor
## 0.7587033

What these results are telling us? There is a positive correlation between RTlexdec and RTnaming. The correlation coefficient (R²) is 0.76 (limits between -1 and 1). This correlation is statistically significant with a t value of 78.699, degrees of freedom of 4566 and a p-value < 2.2e-16.

What are the degrees of freedom? These relate to number of total observations - number of comparisons. Here we have 4568 observations in the dataset, and two comparisons, hence 4568 - 2 = 4566.

For the p value, there is a threshold we usually use. This threshold is p = 0.05. This threshold means we have a minimum to consider any difference as significant or not. 0.05 means that we have a probability to find a significant difference that is at 5% or lower. IN our case, the p value is lower that 2.2e-16. How to interpret this number? this tells us to add 15 0s before the 2!! i.e., 0.0000000000000002. This probability is very (very!!) low. So we conclude that there is a statistically significant correlation between the two variables.

The formula to calculate the t value is below.

$$t = \frac{\bar{x} - \mu_0}{s / \sqrt{n}}$$

$\bar{x}$ = sample mean $0$ = population mean $s$ = sample standard deviation $n$ = sample size

The p value is influenced by various factors, number of observations, strength of the difference, mean values, etc.. You should always be careful with interpreting p values taking everything else into account.
### 2.4.3 Using the package corrplot

Above, we did a correlation test on two predictors. What if we want to obtain a nice plot of all numeric predictors and add significance levels?

#### 2.4.3.1 Correlation plots

```r
corr <-
  english %>
  select(where(is.numeric)) %>
  cor()
corrplot(corr, method = 'ellipse', type = 'upper')
```

<figure><img src="image_43.png" /></figure>
### 2.4.3.2 More advanced

Let's first compute the correlations between all numeric variables and plot these
with the p values

```r
## correlation using "corrplot"
## based on the function 'rcorr' from the 'Hmisc' package
## Need to change dataframe into a matrix
corr <-
  english %>
  select(where(is.numeric)) %>
  data.matrix(english) %>
  rcorr(type = "pearson")

# use corrplot to obtain a nice correlation plot!
corrplot(corr$r$, p.mat = corr$P,
  addCoef.col = "black", diag = FALSE, type = "upper", tl.srt = 55)
```

<figure><img src="image_44.png" /></figure>

```r
english %>
  group_by(AgeSubject) %>
  summarise(mean = mean(RTlexdec),
    sd = sd(RTlexdec))
```

```r
## # A tibble: 2 x 3
## AgeSubject mean sd
```
```
## <fct> <db1> <db1>
## 1 old 6.66 0.116
## 2 young 6.44 0.106
```

## 2.5 Linear Models

Up to now, we have looked at descriptive statistics, and evaluated summaries, correlations in the data (with p values).

We are now interested in looking at group differences. But before moving into linear models, let's assume we simply wanted to use a student's t.test to evaluate the differences between two groups. In our case, we have two age groups within the variable AgeSubject (young vs old). Before running a t.test, there are some assumptions one needs to take into account:

### 2.5.1 Normality of distribution

#### 2.5.1.1 Shapiro test

To check normality of distribution, we use the shapiro.test on the numeric outcome. Given that our predictor now has two levels.

```
shapiro.test(english$RTlexdec)
```

```
##
## Shapiro-Wilk normality test
##
## data: english$RTlexdec
## W = 0.98595, p-value < 2.2e-16
```

How to interpret this statistically significant result? This tells us that the distribution of the data is statistically different from a normal distribution.

#### 2.5.1.2 Density plot

We can also use a density plot to evaluate normality. The results show that both levels have bell shaped distributions but they are both skewed to the right (positively skewed data). See this website for a nice example of left and right skewed data distributions.

```
english %>
 ggplot(aes(x = RTlexdec))+
   geom_density()+
   facetWrap(~ AgeSubject)
```
<figure><img src="image_45.png" /></figure>

## 2.5.2 Homogeneity of variance

If our data were normally distributed, we can use the `bartlett` test. If our data were non-normally distributed, we would use the Levene Test (either with `var.test` from base-R or `leveneTest` from the car package We can check both.

### 2.5.2.1 Bartlett test

english %>%
bartlett.test(RTlexdec ~ AgeSubject, data = .)

##

## Bartlett test of homogeneity of variances

##

## data: RTlexdec by AgeSubject

## Bartlett's K-squared = 15.883, df = 1, p-value = 6.737e-05

### 2.5.2.2 Levene test

english %>%
  var.test(RTlexdec ~ AgeSubject, data = .)

#### 2.5.2.2.1 Base R

##
### F test to compare two variances

##

## data: RTlexdec by AgeSubject

## F = 1.1817, num df = 2283, denom df = 2283, p-value = 6.737e-05

## alternative hypothesis: true ratio of variances is not equal to 1

## 95 percent confidence interval:

## 1.088581 1.282737

## sample estimates:

## ratio of variances

## 1.181678

english %>%

car::leveneTest(RTlexdec ~ AgeSubject, data = .)

### 2.5.2.2.2 From the car package

## Levene's Test for Homogeneity of Variance (center = median)

## Df F value Pr(>F)

## group 1 5.1615 0.02314 *

## 4566

## ---

## Signif. codes: 0 '**' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

In all cases, the statistically significant result indicates that there is evidence that the variance of two levels of the factor `AgeSubject` is statistically significant; i.e., the variances are not equal. This is important as we will use this in our t-test later on.

### 2.5.3 First steps

Up to now, we have looked at descriptive statistics, and evaluated summaries, correlations in the data (with p values), and checked the normality of distribution of our data.

We are now interested in looking at group differences.

Let us start with a simple t-test

#### 2.5.3.1 T-test

We then run a t-test on the data. We specify the formula as `y ~ x` and add `var equal = FALSE` (based on the normality tests above)

english %>%

t.test(RTlexdec ~ AgeSubject, data = ., var equal = FALSE)

##

## Welch Two Sample t-test
##
## data: RTlexdec by AgeSubject
## t = 67.468, df = 4534.6, p-value < 2.2e-16
## alternative hypothesis: true difference in means between group old and group young
## 95 percent confidence interval:
## 0.2152787 0.2281642
## sample estimates:
## mean in group old mean in group young
## 6.660958 6.439237

To interpret the t-test, we say that there is evidence for a statistically significant
difference between the two groups: young and old: t(4534.6) = 67.468, p
< 2.2e-16. The mean of RTlexdec for the young group is significantly lower
than that of the old group.

2.5.3.2 Linear Model

Let us run a linear model on the same data

english %>%
  lm(RTlexdec ~ AgeSubject, data = .) %>%
  summary()

##
## Call:
## lm(formula = RTlexdec ~ AgeSubject, data = .)
##
## Residuals:
## Min 1Q Median 3Q Max
## -0.25776 -0.08339 -0.01669 0.06921 0.52685
##
## Coefficients:
## Estimate Std. Error t value Pr(>\|t\|)
## (Intercept)   6.660958  0.002324  2866.44 <2e-16 ***
## AgeSubjectyoung -0.221721  0.003286  -67.47 <2e-16 ***
## ---
## Signif. codes: 0 '**' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.1111 on 4566 degrees of freedom
## Multiple R-squared: 0.4992, Adjusted R-squared: 0.4991
## F-statistic: 4552 on 1 and 4566 DF, p-value: < 2.2e-16

Any comments? discuss with your neighbour.

The results of the linear model are exactly the same, albeit with a difference
in the sign of the difference. This indicates that the β coefficient for young is
significantly lower than that of old by -0.221721: t = -67.47, p < <2e-16.
The dataset `english` contains two levels for `AgeSubject`. We will run an ANOVA and a linear model on the data.

### 2.5.3.3 Basic ANOVA

We can use the function `aov` to run an Analysis of Variance on the full dataset

```
english %>%
  aov(RTlexdec ~ AgeSubject, data = .) %>%
  summary()
```

```
## Df Sum Sq Mean Sq F value Pr(>F)
## AgeSubject      1  56.14   56.14    4552 <2e-16 ***
## Residuals     4566  56.31    0.01
## ---
## Signif. codes: 0 '**' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
```

### 2.5.3.4 Linear model

We can use the function `lm` to run a linear model

```
english %>%
  lm(RTlexdec ~ AgeSubject, data = .) %>%
  summary()
```

```
##
## Call:
## lm(formula = RTlexdec ~ AgeSubject, data = .)
##
## Residuals:
##       Min     1Q   Median     3Q     Max
## -0.25776 -0.08339 -0.01669  0.06921  0.52685
##
## Coefficients:
##                 Estimate Std. Error t value Pr(>|t|)
## (Intercept)    6.660958  0.002324 2866.44 <2e-16 ***
## AgeSubjectyoung -0.221721  0.003286 -67.47 <2e-16 ***
## ---
## Signif. codes: 0 '**' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.1111 on 4566 degrees of freedom
## Multiple R-squared: 0.4992, Adjusted R-squared: 0.4991
## F-statistic: 4552 on 1 and 4566 DF, p-value: < 2.2e-16
```

But wait... How is the linear model comparable to the analysis of variance we ran above? This linear model derives the analysis of variance we saw above, use `anova` on your linear model.

Here are the results of the initial Analysis of variance:
english %>
aov(RTlexdec ~ AgeSubject, data = .) %>
summary()

## Df Sum Sq Mean Sq F value Pr(F>
## AgeSubject 1 56.14 56.14 4552 <2e-16
## Residuals 4566 56.31 0.01
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

And here are the results of the linear model with the anova function

english %>
lm(RTlexdec ~ AgeSubject, data = .) %>
anova()

## Analysis of Variance Table
##
## Response: RTlexdec
## Df Sum Sq Mean Sq F value Pr(F>
## AgeSubject 1 56.141 56.141 4552 < 2.2e-16
## Residuals 4566 56.314 0.012
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

They are exactly the same... The underlying of an Analysis of variance is a linear model. We will continue with a linear model to understand it better

### 2.5.4 Introduction to Linear models

The basic assumption of a Linear model is to create a regression analysis on the data. We have an outcome (or dependent variable) and a predictor (or an independent variable). The formula of a linear model is as follows outcome ~ predictor that can be read as “outcome as a function of the predictor”. We can add “1” to specify an intercept, but this is by default added to the model

#### 2.5.4.1 Model estimation

english2 <- english %>
mutate(AgeSubject = factor(AgeSubject, levels = c("young", "old")))
mdl.lm <- english2 %>
lm(RTlexdec ~ AgeSubject, data = .)
#lm(RTlexdec ~ AgeSubject, data = english)
mdl.lm #also print(mdl.lm)

##
## Call:
```r
## lm(formula = RTlexdec ~ AgeSubject, data = .)
##
## Coefficients:
## (Intercept) AgeSubjectold
##       6.4392       0.2217

summary(mdl.lm)
```

```r
##
## Call:
## lm(formula = RTlexdec ~ AgeSubject, data = .)
##
## Residuals:
##       Min     1Q   Median     3Q     Max
## -0.25776 -0.08339 -0.01669  0.06921  0.52685
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept)   6.439237  0.002324  2771.03 <2e-16 ***
## AgeSubjectold 0.221721  0.003286   67.47 <2e-16 ***
## ---
## Signif. codes: 0 '**' 0.001 '**' 0.01 '**' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.1111 on 4566 degrees of freedom
## Multiple R-squared: 0.4992, Adjusted R-squared: 0.4991
## F-statistic: 4552 on 1 and 4566 DF, p-value: < 2.2e-16
```

### 2.5.5 Role of coding schemes

#### 2.5.5.1 Intro

There are other coding schemes that can be used. See Schad et al. (2020): "How to capitalize on a priori contrasts in linear (mixed) models: A tutorial". *Journal of Memory and Language*, vol. 110, 104038.

We have:

1. The treatment coding

2. The contrast or sum coding

3. The polynomial coding

4. The Repeated coding

5. The helmert coding

The first two are the most commonly used coding schemes, but see the paper for why one can use the others.
### 2.5.5.2 Treatment coding

By default, R uses a treatment coding scheme. By this we mean that the intercept has the outcome of a specific level in the dataset (usually, the first in alphabetical order, unless you have changed that!). In our example above, we changed the reference level of the variable `AgeSubject` to be "young" vs "old". What we see as an result in our LM is the intercept (`=Young`) vs "AgeSubject-told". This will mean we are looking at the difference between the reference level and all other levels.

The code below allows you to see what the coding scheme is.

```r
english2$AgeSubject2 <- english2$AgeSubject

contrasts(english2$AgeSubject2) <- contr.treatment(2)
contrasts(english2$AgeSubject2)
```

```r
## 2
## young 0
## old 1
```

```r
mdl.lm.T <- english2 %>%
  lm(RTlexdec ~ AgeSubject2, data = .)

#lm(RTlexdec ~ AgeSubject, data = english)
mdl.lm.T #also print(mdl.lm)
```

```r
##
## Call:
## lm(formula = RTlexdec ~ AgeSubject2, data = .)
##
## Coefficients:
## (Intercept) AgeSubject22
##       6.4392      0.2217
```

```r
summary(mdl.lm.T)
```

```r
##
## Call:
## lm(formula = RTlexdec ~ AgeSubject2, data = .)
##
## Residuals:
##       Min     1Q   Median     3Q     Max
## -0.25776 -0.08339 -0.01669  0.06921  0.52685
```

```r
##
## Coefficients:
## Estimate Std. Error t value Pr(>\|t\|)
## (Intercept)  6.439237  0.002324  2771.03  <2e-16 ***
## AgeSubject22  0.221721  0.003286   67.47  <2e-16 ***
## ---
```
```
## Signif. codes: 0 '**' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.1111 on 4566 degrees of freedom
## Multiple R-squared: 0.4992, Adjusted R-squared: 0.4991
## F-statistic: 4552 on 1 and 4566 DF, p-value: < 2.2e-16
```

### 2.5.5.3 Contrast (or sum) coding

#### 2.5.5.3.1 Default
Let's change the coding scheme and see the difference in the output.

```
english2 <- english2 %>%
  mutate(AgeSubject2 = factor(AgeSubject2, levels = c("old", "young")))

contrasts(english2$AgeSubject2) <- contr.sum(2)
contrasts(english2$AgeSubject2)

[,1]
old      1
young    -1

mdl.lm.C <- english2 %>%
  lm(RTlexdec ~ AgeSubject2, data = .)

#lm(RTlexdec ~ AgeSubject, data = english)
mdl.lm.C #also print(mdl.lm)

##
## Call:
## lm(formula = RTlexdec ~ AgeSubject2, data = .)
##
## Coefficients:
## (Intercept) AgeSubject21
##       6.5501      0.1109

summary(mdl.lm.C)

##
## Call:
## lm(formula = RTlexdec ~ AgeSubject2, data = .)
##
## Residuals:
##       Min     1Q   Median     3Q     Max
## -0.25776 -0.08339 -0.01669  0.06921  0.52685
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept)  6.550097  0.001643  3986.29 <2e-16 ***
## AgeSubject21  0.110861  0.001643   67.47 <2e-16 ***
```
## ---
## Signif. codes: 0 '**' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.1111 on 4566 degrees of freedom
## Multiple R-squared: 0.4992, Adjusted R-squared: 0.4991
## F-statistic: 4552 on 1 and 4566 DF, p-value: < 2.2e-16

With contrast (or sum) coding, the intercept is almost the same as before (treatment = 6.439237 vs contrast = 6.550097). The intercept is now the average of all of the data points

english2 %>
mutate(meanRTlexdec = mean(RTlexdec)) %>
select(meanRTlexdec) %>
head(10)

<table><thead><tr><th></th><th>meanRTlexdec</th></tr></thead><tbody><tr><td>## 1</td><td>6.550097</td></tr><tr><td>## 2</td><td>6.550097</td></tr><tr><td>## 3</td><td>6.550097</td></tr><tr><td>## 4</td><td>6.550097</td></tr><tr><td>## 5</td><td>6.550097</td></tr><tr><td>## 6</td><td>6.550097</td></tr><tr><td>## 7</td><td>6.550097</td></tr><tr><td>## 8</td><td>6.550097</td></tr><tr><td>## 9</td><td>6.550097</td></tr><tr><td>## 10</td><td>6.550097</td></tr></tbody></table>

The coefficient for old is different. It is now nearly half of that in the treatment coding (0.221721 vs 0.110861). Why is this the case? In treatment coding, the distance between old and young was of 1 (1-0 = 1), in contrast coding, it is of 2 (1 - -1 = 2). The coefficient of the intercept is exactly the same; for the second, it is half of the one above (0.221721 / 1 = 0.221721; 0.221721 / 2 = 0.110861)!

How to interpret the coefficient for the level “old”? It is the distance from the average!

Let's change slightly the coding scheme

english2 <- english2 %>
mutate(AgeSubject2 = factor(AgeSubject, levels = c("young", "old")))

contrasts(english2$AgeSubject2) <- c(-0.5, 0.5)
contrasts(english2$AgeSubject2)

2.5.5.3.2 Modified

## [,1]
```r
## young -0.5
## old    0.5
```

```r
mdl.lmC.2 <- english2 %>
lm(RTlexdec ~ AgeSubject2, data = .)
#lm(RTlexdec ~ AgeSubject, data = english)
mdl.lmC.2 #also print(mdl.lmC.2)
```

```r
##
## Call:
## lm(formula = RTlexdec ~ AgeSubject2, data = .)
##
## Coefficients:
## (Intercept) AgeSubject21
##       6.5501       0.2217
summary(mdl.lmC.2)
```

```r
##
## Call:
## lm(formula = RTlexdec ~ AgeSubject2, data = .)
##
## Residuals:
##       Min     1Q   Median     3Q     Max
## -0.25776 -0.08339 -0.01669  0.06921  0.52685
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept)  6.550097  0.001643 3986.29 <2e-16 ***
## AgeSubject21  0.221721  0.003286  67.47 <2e-16 ***
## ---
## Signif. codes: 0 '**' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.1111 on 4566 degrees of freedom
## Multiple R-squared: 0.4992, Adjusted R-squared: 0.4991
## F-statistic: 4552 on 1 and 4566 DF, p-value: < 2.2e-16
```

Our intercept is still the average of all datapoints of the dependent variable; that of “old” is now at 0.221721 as expected. The distance between young and old is of 1 (0.5 - -0.5 = 1).

When and why do we use the different coding schemes? We use it to make the intercept more interpretable. This is especially the case when having more than two categorical predictors, or interactions between a numeric and a categorical predictor. Read the paper and the Bodo Winter's book for more details.
2.5.6 Further steps

2.5.6.1 Tidying the output

We use our original model with treatment coding.

```r
# from library(broom)
tidy(mdl.lm) %>%
  select(term, estimate) %>%
  mutate(estimate = round(estimate, 3))
```

```r
## # A tibble: 2 x 2
```

<table><thead><tr><th>term</th><th>estimate</th></tr></thead><tbody><tr><td>&lt;chr&gt;</td><td>&lt;dbl&gt;</td></tr><tr><td>1 (Intercept)</td><td>6.44</td></tr><tr><td>2 AgeSubjectold</td><td>0.222</td></tr></tbody></table>

```r
mycoefE <- tidy(mdl.lm) %>% pull(estimate)
```

Obtaining mean values from our model

```r
old
mycoefE[1]
```

```r
[1] 6.439237
```

```r
young
mycoefE[1] + mycoefE[2]
```

```r
[1] 6.660958
```

2.5.6.2 Nice table of our model summary

We can also obtain a nice table of our model summary. We can use the package `knitr` or `xtable`

```r
kable(summary(mdl.lm)$coef, digits = 3)
```

2.5.6.2.1 Directly from model summary

<table><thead><tr><th></th><th>Estimate</th><th>Std. Error</th><th>t value</th></tr></thead><tbody><tr><td>(Intercept)</td><td>6.439</td><td>0.002</td><td>2771.027</td></tr><tr><td>AgeSubjectold</td><td>0.222</td><td>0.003</td><td>67.468</td></tr></tbody></table>

```r
mdl.lmT <- tidy(mdl.lm)
kable(mdl.lmT, digits = 3)
```

2.5.6.2.2 From the tidy output

<table><thead><tr><th>term</th><th>estimate</th><th>std.error</th><th>statistic</th><th>p.value</th></tr></thead><tbody><tr><td>(Intercept)</td><td>6.439</td><td>0.002</td><td>2771.027</td><td>0</td></tr><tr><td>AgeSubjectold</td><td>0.222</td><td>0.003</td><td>67.468</td><td>0</td></tr></tbody></table>
```R
print tab_model (mdl .lm, file = paste0("outputs/mdl_lms.html")))
htmltools::includeHTML("outputs/mdl_lms.html")
```

### 2.5.6.2.3 Model's fit

#### 2.5.6.3 Dissecting the model

Let us dissect the model. If you use `str`, you will be able to see what is available under our linear model. To access some info from the model

```
str(mdl.lms)
```

##### 2.5.6.3.1 `str` and `coef`

```R
## List of 13
## $ coefficients: Named num [1:2] 6.439 0.222
## ... attr(*, "names")= chr [1:2] "(Intercept)" "AgeSubjectold"
## $ residuals : Named num [1:4568] 0.1045 -0.0416 -0.1343 -0.015 0.0114 ...
## ... attr(*, "names")= chr [1:4568] "1" "2" "3" "4"
## $ effects : Named num [1:4568] -442.7013 7.4927 -0.1352 -0.0159 0.0105 ...
## ... attr(*, "names")= chr [1:4568] "(Intercept)" "AgeSubjectold" "" ""
## $ rank : int 2
## $ fitted.values: Named num [1:4568] 6.44 6.44 6.44 6.44 6.44 ...
## ... attr(*, "names")= chr [1:4568] "1" "2" "3" "4"
## $ assign : int [1:2] 0 1
## $ qr :List of 5
## ..$ qr : num [1:4568, 1:2] -67.587 0.0148 0.0148 0.0148 0.0148 ...
## ... attr(*, "dimnames")=List of 2
## ..$ : chr [1:4568] "1" "2" "3" "4"
## ..$ : chr [1:2] "(Intercept)" "AgeSubjectold"
## ... attr(*, "assign")= int [1:2] 0 1
## ... attr(*, "contrasts")=List of 1
## ..$ AgeSubject: chr "contr.treatment"
## ..$ qraux: num [1:2] 1.01 1.01
## ..$ pivot: int [1:2] 1 2
## ..$ tol : num 1e-07
## ..$ rank : int 2
## ... attr(*, "class")= chr "qr"
## $ df.residual : int 4566
## $ contrasts :List of 1
## ..$ AgeSubject: chr "contr.treatment"
## $ xlevels :List of 1
## ..$ AgeSubject: chr [1:2] "young" "old"
## $ call : language lm(formula = RTlexdec ~ AgeSubject, data = .)
```
## $ terms :Classes 'terms', 'formula' language RTlexdec ~ AgeSubject
## ...- attr(*, "variables")= language list(RTlexdec, AgeSubject)
## ...- attr(*, "factors")= int [1:2, 1] 0 1
## ... ...- attr(*, "dimnames")=List of 2
## ... ... ...$ : chr [1:2] "RTlexdec" "AgeSubject"
## ... ... ...$ : chr "AgeSubject"
## ...- attr(*, "term.labels")= chr "AgeSubject"
## ...- attr(*, "order")= int 1
## ...- attr(*, "intercept")= int 1
## ...- attr(*, "response")= int 1
## ...- attr(*, ".Environment")=<environment: 0x000002c964c27b98>
## ...- attr(*, "predvars")= language list(RTlexdec, AgeSubject)
## ...- attr(*, "dataClasses")= Named chr [1:2] "numeric" "factor"
## ... ...- attr(*, "names")= chr [1:2] "RTlexdec" "AgeSubject"
## $ model : 'data.frame': 4568 obs. of 2 variables:
## ..$ RTlexdec : num [1:4568] 6.54 6.4 6.3 6.42 6.45 ...
## ..$ AgeSubject: Factor w/ 2 levels "young","old": 1 1 1 1 1 1 1 1 1 1 ...
## ...- attr(*, "terms")=Classes 'terms', 'formula' language RTlexdec ~ AgeSubject
## ... ...- attr(*, "variables")= language list(RTlexdec, AgeSubject)
## ... ...- attr(*, "factors")= int [1:2, 1] 0 1
## ... ... ...- attr(*, "dimnames")=List of 2
## ... ... ... ...$ : chr [1:2] "RTlexdec" "AgeSubject"
## ... ... ... ...$ : chr "AgeSubject"
## ... ...- attr(*, "term.labels")= chr "AgeSubject"
## ... ...- attr(*, "order")= int 1
## ... ...- attr(*, "intercept")= int 1
## ... ...- attr(*, "response")= int 1
## ... ... ...- attr(*, ".Environment")=<environment: 0x000002c964c27b98>
## ... ...- attr(*, "predvars")= language list(RTlexdec, AgeSubject)
## ... ...- attr(*, "dataClasses")= Named chr [1:2] "numeric" "factor"
## ... ... ...- attr(*, "names")= chr [1:2] "RTlexdec" "AgeSubject"
## - attr(*, "class")= chr "lm"

coef(mdl.lm)

## (Intercept) AgeSubjectold
## 6.4392366 0.2217215

## same as
## mdl.lm$coefficients

2.5.6.3.2 “coef” and “coefficients” What if I want to obtain the “Intercept”? Or the coefficient for distance? What if I want the full row for distance?

coef(mdl.lm)[1] # same as mdl.lm$coefficients[1]

## (Intercept)
## 6.439237
`coef(mdl.lm)[2] # same as mdl.lm$coefficients[2]`

## AgeSubjectold
## 0.2217215

`summary(mdl.lm)$coefficients[2, ] # full row`

## Estimate Std. Error t value Pr(>\|t\|)
## 0.22172146 0.00328631 67.46820211 0.00000000

`summary(mdl.lm)$coefficients[2, 4] #for p value`

## [1] 0

### 2.5.6.3.3 Residuals
What about residuals (difference between the observed value and the estimated value of the quantity) and fitted values? This allows us to evaluate how normal our residuals are and how different they are from a normal distribution.

`hist(residuals(mdl.lm))`

<figure><img src="image_46.png" /><figcaption>Histogram of residuals(mdl.lm)</figcaption></figure>

#### 2.5.6.3.3.1 Histogram

residuals(mdl.lm)

`qqnorm(residuals(mdl.lm)); qqline(residuals(mdl.lm))`
Normal Q-Q Plot

<figure><img src="image_47.png" /></figure>

Theoretical Quantiles

2.5.6.3.3.2 qqplots

plot(fitted mdl lm), residualsponds, cex = 4)

<figure><img src="image_48.png" /></figure>

fitted(mdl.lm)

2.5.6.3.3.3 Residuals vs fitted values

sim_residuals <- simulateResiduals(mdl.lm)
plot(sim_residuals)
### 2.5.6.3.3.4 Generating and plotting residual plot with DHARMa
DHARMa residual

<figure><img src="image_49.png" /></figure>

AIC(mdl.lm) # Akaike's Information Criterion, lower values are better

### 2.5.6.3.4 Goodness of fit?

## [1] -7110.962

BIC(mdl.lm) # Bayesian AIC

## [1] -7091.682

logLik(mdl.lm) # log likelihood

## 'log Lik.' 3558.481 (df=3)

Or use the following from broom

glance(mdl.lm)

## # A tibble: 1 x 12

## r.squared adj.r.squared sigma statistic p.value df logLik AIC BIC
## <db1> <db1> <db1> <db1> <db1> <db1> <db1> <db1> <db1>
## 1 0.499 0.499 0.111 4552. 0 1 3558. -7111. -7092.

## # i 3 more variables: deviance <db1>, df.residual <int>, nobs <int>

### 2.5.6.3.5 Significance testing
Are the above informative? of course not directly. If we want to test for overall significance of model. We run a null model (aka intercept only) and compare models.
```
mdl.lm.Null <- english %>%
  lm(RTlexdec ~ 1, data = .)
mdl.comp <- anova(mdl.lm.Null, mdl.lm)
mdl.comp
```

```
## Analysis of Variance Table
##
## Model 1: RTlexdec ~ 1
## Model 2: RTlexdec ~ AgeSubject
## Res.Df      RSS Df Sum of Sq      F      Pr(>F)
## 1          4567  112.456
## 2          4566   56.314   1      56.141   4552 < 2.2e-16 ***
## ---
## Signif. codes: 0 '**' 0.001 '**' 0.01 '**' 0.05 '.' 0.1 ' ' 1
```

The results show that adding the variable "AgeSubject" improves the model fit. We can write this as follows: Model comparison showed that the addition of AgeSubject improved the model fit when compared with an intercept only model ($F(1) = 4551.96, p < 0$) ($F(1) = 4552, p < 2.2e-16$)

## 2.5.7 Plotting fitted values

### 2.5.7.1 Trend line

Let's plot our fitted values but only for the trend line

```
english %>%
  ggplot(aes(x = AgeSubject, y = RTlexdec))+
  geom_boxplot()+
  themeBW() + theme(text = element_text(size = 15))+
  geom.smooth(aes(x = as.numeric(AgeSubject), y = predict(mdl.lm)), method = "lm", col
  labs(x = "Age", y = "RTLexDec", title = "Boxplot and predicted trend line", subtitle
```
<figure><img src="image_50.png" /></figure>

This allows us to plot the fitted values from our model with the predicted linear trend. This is exactly the same as our original data.

### 2.5.7.2 Predicted means and the trend line

We can also plot the predicted means and linear trend

```r
eyes

english %>

 ggplot(aes(x = AgeSubject, y = predict(mdl.lm))) +
   geom_boxplot(color = "blue") +
   themeBW() + theme(text = element_text(size = 15)) +
   geom.smooth(aes(x = as.numeric(AgeSubject), y = predict(mdl.lm)),
               method = "lm", color = "blue")) +
   labs(x = "Age", y = "RTLexDec", title = "Predicted means and trend line",
       subtitle = "with ggplot2")
```
<figure><img src="image_51.png" /></figure>

### 2.5.7.3 Raw data, predicted means and the trend line

We can also plot the actual data, the predicted means and linear trend

```r
e English %>
  ggplot(aes(x = AgeSubject, y = RTlexdec)) +
  geom_boxplot() +
  geom_boxplot(aes(x = AgeSubject, y = predict(mdl.lm)), color = "blue") +
  themeBW() + theme(text = element_text(size = 15)) +
  geom.smooth(aes(x = as.numeric(AgeSubject), y = predict(mdl.lm)), method = "lm", color = "blue")) +
  labs(x = "Species", y = "Length", title = "Boxplot raw data, predicted means (in blue)")
```
<figure><img src="image_52.png" /><figcaption>Boxplot raw data, predicted means (in blue) and trend with ggplot2</figcaption></figure>

### 2.5.7.4 Add significance levels and trend line on a plot?

We can use the p values generated from either our linear model to add significance levels on a plot. We use the code from above and add the significance level. We also add a trend line

```r
e English %>
  ggplot(aes(x = AgeSubject, y = RTlexdec)) +
  geom_boxplot() +
  geom_boxplot(aes(x = AgeSubject, y = predict(mdl.lm)), color = "blue") +
  themeBW() + theme(text = element_text(size = 15)) +
  geom.smooth(aes(x = as.numeric(AgeSubject), y = predict(mdl.lm)), method = "lm", color = "blue")) +
  labs(x = "Species", y = "Length", title = "Boxplot raw data, predicted means (in blue) and trend line") +
  geom_signif(comparison = list(c("old", "young")),
    map_signif_level = TRUE, test = function(a, b) {
      list(p.value = summary(mdl.lm)$coefficients[2, 4]))}
```
Boxplot raw data, predicted means (in blue) and tren
with significance testing

<figure><img src="image_53.png" /></figure>

2.5.8 What about pairwise comparison?

When having three of more levels in our predictor, we can use pairwise compar-
isons, with corrections to evaluate differences between each level.

summary(mdl.lm)

##

## Call:

## lm(formula = RTlexdec ~ AgeSubject, data = .)

##

## Residuals:

## Min 1Q Median 3Q Max
## -0.25776 -0.08339 -0.01669 0.06921 0.52685

##

## Coefficients:

## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 6.439237 0.002324 2771.03 <2e-16 ***
## AgeSubjectold 0.221721 0.003286 67.47 <2e-16 ***

## ---

## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

##

## Residual standard error: 0.1111 on 4566 degrees of freedom

## Multiple R-squared: 0.4992, Adjusted R-squared: 0.4991

## F-statistic: 4552 on 1 and 4566 DF, p-value: < 2.2e-16
```
mdl.lm %>% emmeans(pairwise ~ AgeSubject, adjust = "fdr") -> mdl.emmeans
```

```
## $emmeans
## AgeSubject emmean      SE   df lower.CL upper.CL
## young          6.439 0.00232 4566   6.435    6.444
## old           6.661 0.00232 4566   6.656    6.666
##
## Confidence level used: 0.95
##
## $contrasts
## contrast   estimate      SE   df t.ratio p.value
## young - old  -0.222 0.00329 4566 -67.468 <.0001
```

How to interpret the output? Discuss with your neighbour and share with the group.

Hint... Look at the emmeans values for each level of our factor "AgeSubject" and the contrasts.

### 2.5.9 Multiple predictors?

Linear models require a numeric outcome, but the predictor can be either numeric or a factor. We can have more than one predictor. The only issue is that this complicates the interpretation of results

```
english %>%
  lm(RTlexdec ~ AgeSubject * WordCategory, data = .) %>%
  summary()

##
## Call:
## lm(formula = RTlexdec ~ AgeSubject * WordCategory, data = .)
##
## Residuals:
##       Min    1Q   Median     3Q      Max
## -0.25079 -0.08273 -0.01516  0.06940  0.52285
##
## Coefficients:
##                 Estimate   Std. Error   t value   Pr(>\|t\|)
## (Intercept)      6.664955  0.002911  2289.950  <2e-16 ***
## AgeSubjectyoung  -0.220395  0.004116  -53.545   <2e-16 ***
## WordCategoryV     -0.010972  0.004822   -2.275   0.0229 *
## AgeSubjectyoung:WordCategoryV -0.003642  0.006820   -0.534   0.5934
## ---
## Signif. codes: 0 '**' 0.001 '**' 0.01 '**' 0.05 '.' 0.1 '' 1
##
```
## Residual standard error: 0.1109 on 4564 degrees of freedom
## Multiple R-squared: 0.5008, Adjusted R-squared: 0.5005
## F-statistic: 1526 on 3 and 4564 DF, p-value: < 2.2e-16

And with an Anova

english %>%
  lm(RTlexdec ~ AgeSubject * WordCategory, data = .) %>%
  anova()

## Analysis of Variance Table

##

## Response: RTlexdec

| | Df | Sum Sq | Mean Sq | F value | Pr(>F) |
| :--- | :--- | :--- | :--- | :--- | :--- |
| **AgeSubject** | 1 | 56.141 | 56.141 | 4564.2810 | < 2.2e-16 *** |
| **WordCategory** | 1 | 0.173 | 0.173 | 14.0756 | 0.0001778 *** |
| **AgeSubject:WordCategory** | 1 | 0.004 | 0.004 | 0.2851 | 0.5933724 |
| **Residuals** | 4564 | 56.138 | 0.012 | | |
| **---** | | | | | |
| **Signif. codes:** | 0 | *** | 0.001 | ** | 0.01 | * | 0.05 | . | 0.1 | ' | 1 |

The results above tell us that all predictors used are significantly different.

## 2.6 Generalised Linear Models

Here we will look at an example when the outcome is binary. This simulated data is structured as follows. We asked one participant to listen to 165 sentences, and to judge whether these are “grammatical” or “ungrammatical”. There were 105 sentences that were “grammatical” and 60 “ungrammatical”. Responses varied, with 110 “yes” responses and 55 “no” responses, spread across the grammatical.

This fictitious example can apply in any other situation. Let’s think Geography: 165 lands: 105 “flat” and 60 “non-flat”, etc. This applies to any case where you need to “categorise” the outcome into two groups.

### 2.6.1 Load and summaries

Let’s create this dataset from scratch (you can load the “grammatical.csv file as well) and do some basic summaries

grammatical <- as.data.frame(
  cbind("grammaticality" = c("grammatical" = rep("grammatical", 105),
                            "ungrammatical" = rep("ungrammatical", 60)),
  "response" = c("yes" = rep("yes", 100),
                "no" = rep("no", 5),
                "yes" = rep("yes", 10),
                "no" = rep("no", 50))),
  row.names = FALSE)
grammatical %>%
  head(20)

<table><thead><tr><th></th><th></th><th></th></tr></thead><tbody><tr><td>##</td><td>grammaticality</td><td>response</td></tr><tr><td>## 1</td><td>grammatical</td><td>yes</td></tr><tr><td>## 2</td><td>grammatical</td><td>yes</td></tr><tr><td>## 3</td><td>grammatical</td><td>yes</td></tr><tr><td>## 4</td><td>grammatical</td><td>yes</td></tr><tr><td>## 5</td><td>grammatical</td><td>yes</td></tr><tr><td>## 6</td><td>grammatical</td><td>yes</td></tr><tr><td>## 7</td><td>grammatical</td><td>yes</td></tr><tr><td>## 8</td><td>grammatical</td><td>yes</td></tr><tr><td>## 9</td><td>grammatical</td><td>yes</td></tr><tr><td>## 10</td><td>grammatical</td><td>yes</td></tr><tr><td>## 11</td><td>grammatical</td><td>yes</td></tr><tr><td>## 12</td><td>grammatical</td><td>yes</td></tr><tr><td>## 13</td><td>grammatical</td><td>yes</td></tr><tr><td>## 14</td><td>grammatical</td><td>yes</td></tr><tr><td>## 15</td><td>grammatical</td><td>yes</td></tr><tr><td>## 16</td><td>grammatical</td><td>yes</td></tr><tr><td>## 17</td><td>grammatical</td><td>yes</td></tr><tr><td>## 18</td><td>grammatical</td><td>yes</td></tr><tr><td>## 19</td><td>grammatical</td><td>yes</td></tr><tr><td>## 20</td><td>grammatical</td><td>yes</td></tr></tbody></table>

str(grammatical)

## 'data.frame': 165 obs. of 2 variables:
## $ grammaticality: chr "grammatical" "grammatical" "grammatical" "grammatical" ...
## $ response : chr "yes" "yes" "yes" "yes" ...

head(grammatical)

<table><thead><tr><th></th><th></th><th></th></tr></thead><tbody><tr><td>##</td><td>grammaticality</td><td>response</td></tr><tr><td>## 1</td><td>grammatical</td><td>yes</td></tr><tr><td>## 2</td><td>grammatical</td><td>yes</td></tr><tr><td>## 3</td><td>grammatical</td><td>yes</td></tr><tr><td>## 4</td><td>grammatical</td><td>yes</td></tr><tr><td>## 5</td><td>grammatical</td><td>yes</td></tr><tr><td>## 6</td><td>grammatical</td><td>yes</td></tr></tbody></table>

table(grammatical$response, grammatical$grammaticality)

<table><thead><tr><th></th><th></th><th></th></tr></thead><tbody><tr><td>##</td><td>grammatical</td><td>ungrammatical</td></tr><tr><td>## no</td><td>5</td><td>50</td></tr><tr><td>## yes</td><td>100</td><td>10</td></tr></tbody></table>
## 2.6.2 GLM - Categorical predictors

Let's run a first GLM (Generalised Linear Model). A GLM uses a special family "binomial" as it assumes the outcome has a binomial distribution. In general, results from a Logistic Regression are close to what we get from SDT (see above).

To run the results, we will change the reference level for both response and grammaticality. The basic assumption about GLM is that we start with our reference level being the "no" responses to the "ungrammatical" category. Any changes to this reference will be seen in the coefficients as "yes" responses to the "grammatical" category.

### 2.6.2.1 Model estimation and results

The results below show the logodds for our model.

```
grammatical <- grammatical %>
mutate(response = factor(response, levels = c("no", "yes")),
        grammaticality = factor(grammaticality, levels = c("ungrammatical", "grammatical")))
```

```
grammatical %>
group_by(grammaticality, response) %>
table()
```

```
## response
```

```
## grammaticality no yes
```

```
## ungrammatical 50 10
```

```
## grammatical 5 100
```

```
mdl.glm <- grammatical %>
glm(response ~ grammaticality, data = ., family = binomial)
summary(mdl.glm)
```

```
##
## Call:
##   glm(formula = response ~ grammaticality, family = binomial, data = .)
```

```
##
## Coefficients:
##                 Estimate Std. Error z value Pr(>\|z\|)
## (Intercept)      -1.6094   0.3464  -4.646 3.38e-06 ***
## grammaticality    grammatical    4.6052   0.5744   8.017 1.08e-15 ***
## ---
## Signif. codes: 0 '**' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
```

```
##
## (Dispersion parameter for binomial family taken to be 1)
```

```
##
## Null deviance: 210.050 on 164 degrees of freedom
## Residual deviance: 94.271 on 163 degrees of freedom
```

```
##
```