Alongside type, objects in R also have a class, which can be viewed with the class() function. For atomic types, their class is the same as their type, but for more complicated types of variables, a class essentially describes a special instance of a type. Objects within a specific class typically have their own set of functions or methods that can only be applied for that specific class. A common example of a class is factor. Factors are a special class of integer vectors designed to represent categorical and ordinal data. Variables that are of the factor class have levels, which differentiates factor variables from ordinary integer variables. For example, a factor could represent group membership in a randomized experiment indicating inclusion in the control group or the treatment group for each individual. The levels of this factor could be called "control" and "treatment" for example. Supposing that 0 indicates that an individual belongs in the control group and 1 indicates that an individual belongs in the treatment group, we could create such a factor by writing

```
group <- c(0, 0, 1, 0, 1, 0, 1, 1)
factor(group, labels = c("control", "treatment"))
```
[1] control control treatment control treatment control treatment treatment Levels: control treatment

Factors are important when fitting statistical models or when performing statistical tests, because if we would simply use the corresponding integer values, they may be erroneously interpreted as continuous. The levels of the factor also often help to provide more informative output. In the next section we will discuss more complicated data types including data frames, which can contain data of various types simultaneously.

#### **7 Basic R Objects**

In this section we will discuss the concepts of data frame and tibble. Let's assume that your data can be stored in a two-dimensional array where the columns represent variables and each row represents a case of measurement. In R, a data frame is a concept for storing such a two-dimensional array of data. Let's first study how data frames work in R and we will then move on to see how another concept called tibble extends the capabilities of a data frame.

A data frame which has been loaded into R under the name grades and printed in the R console will look as follows.

grades group grade 1 2 4.67 2 2 4.90 3 3 2.63 4 4 3.39 5 4 6.89

In the above printout we can see that R prints the data frame just as we would expect the data to look like. One issue with a standard data frame is that if there are very many columns or rows, then the printout may be difficult to read. In RStudio, a neater (and more flexible) way of inspecting the data is by using a command called View():

#### View(grades)

To go further with the data, we need tools for data manipulation. We begin with the very basic tools which are commonplace in any data analysis workflow. More comprehensive knowledge about data transforming, cleaning etc. can be found in Chap. 2.

A data frame can be constructed directly in the R by using the function called data.frame() and then listing variable names and their values. The grades data above has two columns, group and grade, and each row represents a student. The variable group is indicates the number of a group in which each person has studied. The variable grade stores the final grade that the student has received from a course.

#### **8 Working with Dataframes**

To extract a column from a data frame, one needs to start with the name of the data frame object and connect the column name with the name of the data frame object by using the dollar symbol (\$). Thus, extracting the column group from grades data can be accomplished by writing

grades\$group

[1] 22344

Next, let's use a couple of basic functions which are often needed when developing R code. First of all, to get a summary of every variable of a data frame, we can call

summary(grades)

The result show us the minimum, maximum, median, mean and quartiles of both variables. You may note that as group was intended to be a categorical variable, so computing its mean value in the data does not make sense. We can change that behavior by converting the group column into a factor.

```
grades$group <- as.factor(grades$group)
summary(grades)
```
group grade 2:2 Min. :2.630 3:1 1st Qu.:3.390 4:2 Median :4.670 Mean :4.496 3rd Qu.:4.900 Max. :6.890

On the other hand, we might want to calculate the mean or the sample standard deviation of the variable grade. This can be done with functions called mean() and sd(), respectively.

```
mean(grades$grade)
```
[1] 4.496

sd(grades\$grade)

[1] 1.630178

The functions above can only take numeric vectors as input. If we tried using another type of argument, we would encounter an error message.

sd(grades)

Error in is.data.frame(x): 'list' object cannot be coerced to type 'double'

The above message can actually teach us a couple of things. First of all, the error comes from the function is.data.frame(), which is called somewhere in the definition of sd(). Second, the actual error message tells us that the object which we gave to the function as its argument is a list object. List is an object type of R on top of which data.frame objects have been built. We will describe lists in greater detail later. Further, the message tells us that R has tried to coerce the argument object into the double type. This means that the object we supplied to the function is not of the right type and cannot easily be converted into the proper format.

#### *8.1 tibble*

A tibble is an expansion of data.frame objects which is used in the tidyverse programming paradigm [10]. To use tidyverse, we advise to load the tidyverse [11] (meta)package which loads all key tidyverse packages.

```
# load tidyverse to use as_tibble
library("tidyverse")
# convert a data frame as tibble
grades2 <- as_tibble(grades)
```
Next, let's see what a tibble looks like when printed in the console

```
grades2
# A tibble: 5 x 2
 group grade
 <fct> <dbl>
1 2 4.67
2 2 4.9
3 3 2.63
4 4 3.39
5 4 6.89
```
Tibbles behave similarly compared to data frames when printed, but they also describe the dimensions of the data and the types of the columns right under the column names. For instance, <fct> refers to a factor column and <dbl> refers to a double column. In order to discover the column types when using data frames, one would need to apply the class() or typeof() function to the columns, or write str(grades) to see the types of the columns and the structure of the data. Another useful property of tibble tables is that if a tibble has a large number of observations or variables, then only the rows or the columns which can fit on to the screen are printed.

Tidyverse and tibbles also support so called lazy evaluation, which is useful when your data is stored in a database, for instance. With lazy evaluation, the commands that you use on your data would be evaluated directly in the database (if possible). Without lazy evaluation, the entire data would be downloaded onto your computer only after which the commands would be evaluated. Lazy evaluation can perform many tasks faster and it can also alleviate memory usage of the computer.

#### **9 Pipes**

Piping is a fairly recent concept in R and was very rarely used in R code just a few years ago. The concept of a pipe originates from a package called magrittr, and pipes are commonly used under the tidyverse programming paradigm, but the pipe was also later added to base R. The notations for a tidyverse pipe and a native R pipe are %>% and |>, respectively. The idea of a pipe is that you can connect multiple function calls sequentially while keeping the code more readable. Pipes also serve to unnest standard R code which often involves using many nested parentheses, and can quickly become hard to read as one has to read the code based on the order of the operations instead of reading it linearly. For example, consider the following code where we apply a sequence of operations on a numeric vector x.

```
x <- 1:10
round(mean(diff(log(x))), digits = 2)
```
[1] 0.26

This code computes the rounded mean differences of the logarithms of the vector x, however this description does not match the order of operations, where the logarithm is computed first. To accomplish the same result using pipes we would write

```
x <- 1:10
x |> log() |> diff() |> mean() |> round(digits = 2)
```
[1] 0.26

Here, the order of operations can be easily read from left to right. Next, we will discuss the use of pipes in more detail.

#### *9.1 magrittr pipe %>%*

Let's have a look at an example, where we call the summary() function for the grades2 data using the magrittr pipe %>%, and we also define that results should be printed with two digits.

```
grades2 %>%
 summary(digits = 2)
group grade
2:2 Min. :2.6
3:1 1st Qu.:3.4
4:2 Median :4.7
      Mean :4.5
      3rd Qu.:4.9
      Max. :6.9
```
In the code above, the object grades2 is taken by the pipe operator %>% and forwarded to the first argument of the summary() function. The summary() function also has a second argument, which is defined by digits = 2. Thus, the pipe only takes the object mentioned before the pipe operator and forwards it to the function after the pipe as the first free argument. It is very common and recommended to structure R code so that there is only one pipe per row and that a new line is started after each pipe.

Although the above example is easy to understand as we already know the summary() function, there is also a more general way to compute summarized information following the tidyverse style. The function summarise() can be used to compute arbitrary statistics from the data, for example the number of observations (via the function n()) and the mean and the sample standard deviation of the variable grade.

```
grades2 %>%
  summarise(
    n = n(),
    mean = mean(grade),
    sd = sd(grade)
  )
# A tibble: 1 x 3
      n mean sd
  <int> <dbl> <dbl>
1 5 4.50 1.63
```
The summarise() function also produces a tibble enabling further operations via piping, if desired.

#### *9.2 Native pipe |>*

In R version 4.1, the native pipe |> was introduced to the R language, which does not require any external packages to use. In most scenarios, it does not matter whether the native or the magrittr pipe is used. However, there are two technical differences between the magrittr pipe and the native pipe. First, the magrittr pipe is actually a function

```
class(`%>%`)
```
[1] "function"

while the native pipe is not, and simply converts the written code into a non-piped version, i.e., into a form that one would write without using the pipe:

```
x <- 1:5
quote(x |> sum())
```
sum(x)

We see that providing x to the sum function via the native pipe is exactly the same as writing sum(x) directly. What this means in practice is that the native pipe may have better performance for example when passing a large dataset through a large number of pipes. The reason for this is that the magrittr pipe incurs an additional computational function call overhead each time it is called. The second difference between the pipes is that parentheses have to be provided for function calls when using the native pipe, but they can be omitted when using the magrittr pipe

```
x %>% sum
[1] 15
x %>% sum()
[1] 15
x |> sum()
[1] 15
x |> sum # produces an error
Error: The pipe operator requires a function call
         as RHS (<text>:1:6)
```
#### **10 Lists**

Earlier, we already briefly mentioned lists in the context of data frames. Lists are one the most common types of data in R and they resemble basic vectors in many aspects. Like vectors, a list is an ordered sequence of elements, but unlike vectors, lists can contain elements of different types simultaneously and may even contain other lists. For example, we could construct a list that contains a logical value, a numeric value and a character value using the list() function as follows.

y <- list(TRUE, 7.2, "this is a list")

Subsetting a list object works slightly differently compared to vectors. When single brackets are used, a sublist is selected, i.e., a list object that contains the elements at the supplied indices, for example:

```
y[1:2]
[[1]]
[1] TRUE
[[2]]
[1] 7.2
```

```
typeof(y[1:2])
```
[1] "list"

To extract an actual element of a list, double brackets should be used:

y[[2]]

[1] 7.2

The elements of a list may also be named, which enables subsetting via the dollar sign operator similar to data frames, or by giving the element name in double brackets instead of the index:

```
z <- list(bool = TRUE, num = 7.2,
                         description = "this is another list")
z$bool
[1] TRUE
```
z[["description"]]

[1] "this is another list"

One benefit of using the dollar sign is that it is not necessary to provide the full element name, unlike when using the double brackets. It is sufficient to provide a prefix of the element name so that the full name can be uniquely determined from the prefix. Because all the names of our elements in the previous list z start with a different letter, the first letter of the name suffices as the prefix. The same functionality also applies when using the dollar sign to select columns of data frames.

z\$b [1] TRUE z\$n [1] 7.2 z\$d

[1] "this is another list"

#### **11 Functions**

A function is a set of statements that when organized together perform a specific task. Each function in R has a name, a set of arguments, a body and a return object. The name of the function usually describes the purpose of the function. For example, the base R function mean() computes the arithmetic mean of the argument vector. The result is returned as a numeric value.

```
x <- 1:5
mean(x)
```

```
[1] 3
```
Functions are often much more complicated, which is why it is often helpful to view the documentation of a function before using it in practice. To view the documentation pages of a function, one can simply write the name of the function prefixed by a question mark.

#### ?mean

In RStudio, the documentation will open in the "Help" tab in the bottom right pane by default. Functions will only be executed when they are called, i.e., when arguments are supplied to them. Simply writing the function name without parentheses will instead print the body of the function to the console, meaning the code that the function consists of and which is executed if the function is called.

In the previous sections, we've already familiarized ourselves with some commonly used basic functions such as c(), sd(), and summary(). Base R has a wide range of function to accomplish common tasks needed in data analysis, which is further extended by tidyverse and other R packages. This means that one does not typically have to write their own functions when programming in R. We will explore several of the functions provided by the tidyverse in later chapters.

#### **12 Conditional Statements**

Sometimes we may only wish to execute a piece of code when a certain condition is met. Conditional statements in R can be defined via the if and else clauses. The if clause evaluates a condition, which is an R expression that evaluates to a single logical value, and if this condition evaluates to TRUE, the expression following the clause if executed. Further, if an else clause is also provided, the expression following else will be executed instead if the condition evaluates to FALSE. As R code, the syntax for these clauses is

```
if (cond) expr
if (cond) expr else alt_expr
```
where cond is the condition being evaluated, expr is the expression that will be evaluated if cond == TRUE, and alt\_expr will be evaluated if cond == FALSE.

Note that if will only evaluate a single condition. If cond is a vector, an error will be produced:

```
cond <- c(TRUE, FALSE)
if (cond) {
  print("This will not be printed")
}
```

```
Error in if (cond) {: the condition has length > 1
```
As expected, the error message tells us that the condition contained more than one element when it was evaluated. However, there are often scenarios where we may wish to conditionally select or define values based on a vector of conditions. For such instances, the function ifelse() can be used. This function has three arguments: test, yes, and no. When called, the function will pick those elements of the vector yes for which the logical vector test evaluates to TRUE, and those elements of the vector no for which test evaluates to FALSE:

```
cond <- c(TRUE, FALSE, FALSE, TRUE)
x <- 1:4
y <- -(1:4)
ifelse(cond, x, y)
[1] 1 -2 -3 4
```
#### **13 Looping Constructs**

In some cases, we may wish to execute the same piece of code multiple times under varying conditions. Instead of writing the same code multiple times for each condition, we can use a looping construct. There are two types of loops in R: the for loop and the while loop. The main difference between the two loops is that for always executes the code associated with the loop a fixed number of times whereas while will continue executing the code as long as a specific condition remains satisfied. The syntax of these loops is

```
for (var in seq) expr
while (cond) expr
```
In other words, for will execute the expression expr for every element var in some object seq that can be indexed. For-loops are very general, and can be used to loop over most ordered structures such as vectors and lists. Similarly, while will execute the expression expr as long as the condition cond evaluates to TRUE. Care must be taken when using while-loops to ensure that the condition will eventually evaluate to FALSE, otherwise the loop will simply run indefinitely and the program will be stuck. As an example, we will print the number 1 through 5 to the console using both for and while loops:

```
x <- 1:5
for (i in x) {
  print(x[i])
}
```

```
[1] 1
[1] 2
[1] 3
[1] 4
[1] 5
i <- 0
while (i < length(x)) {
  i <- i + 1
  print(x[i])
}
[1] 1
[1] 2
[1] 3
[1] 4
```

```
[1] 5
```
In contrast to explicit for and while loops, the so-called apply family of functions can often be a simpler alternative (see ?apply). As the name suggests, these functions apply an operation to each element of a list or a vector (and other more general data structures). For example, the above loop example could also be accomplished with the lapply() function as follows:

```
y <- lapply(x, print)
[1] 1
[1] 2
[1] 3
```

```
[1] 4
[1] 5
```
#### **14 Discussion and Other Resources for Learning R**

The main aim behind this chapter was to introduce R to new users. This chapter is, of course, an initial step and can hardly cover all the basics. Interested users are advised to consult other resources e.g., open access books, tutorials, cheat sheets, and package manuals for more information. An introductory book "An Introduction to R" packaged with each R installation can be accessed from the RStudio Help panel by first selecting "Show R Help" and then selecting "An Introduction to R" under "Manuals". This book covers a wide range of topics on base R programming in great detail. The book "R for Data Science" by Hadley Wicham and Garret Grolemund provides a comprehensive tutorial on using R for data science under the tidyverse paradigm. The book is free to use and readily available online at https:// r4ds.had.co.nz/. As we go further, several questions will emerge and the reader will learn by doing and by consulting the literature and help files. In doing so, the reader will build knowledge and experience that helps advance their skills.

#### **References**

**Open Access** This chapter is licensed under the terms of the Creative Commons Attribution 4.0 International License (http://creativecommons.org/licenses/by/4.0/), which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons license and indicate if changes were made.

The images or other third party material in this chapter are included in the chapter's Creative Commons license, unless indicated otherwise in a credit line to the material. If material is not included in the chapter's Creative Commons license and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder.

# **An R Approach to Data Cleaning and Wrangling for Education Research**

**Juho Kopra, Santtu Tikka, Merja Heinäniemi, Sonsoles López-Pernas, and Mohammed Saqr**

#### **1 Introduction**

When analyzing data, it is crucial that the data is in a suitable format for the tools you will be using. This makes data wrangling essential. Data preparation and cleaning, such as extracting information from raw data or removing erroneous measurements, must also be done before data is ready for analysis. Data wrangling often takes up the majority of the time spent on analysis, sometimes up to 80%. To reduce the amount of work required, it is beneficial to use tools that follow the same design paradigm to minimize the time spent on data wrangling. The tidyverse [1] programming paradigm is currently the most popular approach for this in R.

The tidyverse has several advantages that make it preferable over other alternatives such as simply using base R for your data wrangling needs. All packages in the tidyverse follow a consistent syntax, making it intuitive to learn and use new tidyverse packages. This consistency also makes the code more easier to read, and maintain, and reduces the risk of errors. The tidyverse also has a vast range of readily available packages that are actively maintained, reducing the need for customized code for each new data wrangling task. Further, these packages integrate seamlessly with one another, facilitating a complete data analysis pipeline.

J. Kopra (-) · S. López-Pernas · M. Saqr

School of Computing, University of Eastern Finland, Kuopio, Finland e-mail<EMAIL_ADDRESS>
S. Tikka Department of Mathematics and Statistics, University of Jyväskylä, Finland

M. Heinäniemi Institute of Biomedicine, University of Eastern Finland, Kuopio, Finland

To fully realize the benefits of the tidyverse programming paradigm, one must first understand the key concepts of tidy data and pivoting. Tidy data follows three rules:

Let's consider examples of tidy data. For instance, if you have data from a Moodle course where two attempts of an exam for each student are located in a single column. This example data violates the first rule, because there are two variables in a single column instead of separate columns for each variable. What is needed to make this data tidy is to pivot it to a longer format. The pivoted data would have two rows for each student, both of which are different observations (exam attempts 1 and 2). Thus the pivoted data would not conflict with second rule.

Data can also be too long, but in practice, this is much more rare. This can occur if two or more variables are stored on a single column across multiple rows. A key indicators of this is if the different rows of the same column have different measurements units (e.g. lb vs. kg). It may also occur that your raw data has multiple values in a single cell. In these cases, it is necessary to split the cells to extract the necessary information. In a simple case, where you have two values systematically in a one cell, the values can be easily separated into their own columns.

Overall, using tidyverse and understanding the key concepts of tidy data and pivoting can streamline the data analysis process and make code easier to work with and maintain. The rest of this chapter will guide readers through the process of data cleaning and wrangling with R in the field of learning analytics. We demonstrate how data can be grouped and summarized, how to select and transform variables of interest, and how data can be rearranged, reshaped and joined with other datasets. We will strongly rely on the tidyverse programming paradigm for a consistent and coherent approach to data manipulation, with a focus on tidy data.

#### **2 Reading Data into R**

Data files come in many formats, and getting your data ready for analysis can often be a daunting task. The tidyverse offers much better alternatives to base R functions for reading data, especially in terms of simplicity and speed when reading large files. Additionally, most of the file input functions in the tidyverse follow a similar syntax, meaning that the user does not have to master every function for reading every type of data individually.

Often, just before the data can be read into R, user must specify the location of data files by setting a working directory. Perhaps most useful way to do that is to create a project in RStudio and then create a folder called "data" within the project folder. Data files can be put into that folder and user can refer to those files just by telling R-functions relative path of data file (e.g. "data/Final%20Dataset.csv") while the project takes care of the rest of the path. A more traditional way of setting this, which also works without RStudio, is by using a command such as setwd("/home/Projects/LAproject/data/Final%20Dataset.csv"). Here, a function called setwd() is used to set up a folder into location mentioned in a character string given as its argument. A getwd() lists current working directory, which can also be seen in RStudio just above the Console output.

Some of the most common text data formats are comma-separated files or semicolon-separated files, both of which typically have the file extension .csv. These files can be read into R using the readr [2] package and the functions read\_csv() and read\_csv2(), respectively. For instance, we can read a comma-separated file R as follows

```
library("readr")
url <- "https://github.com/lamethods/data/raw/main/2_moodleEdify/"
lms <- read_csv(paste0(url, "Final%20Dataset.csv"))
```
Functions in readr provide useful information about how the file was read into R, which can be used to assess if the input was successful and what assumptions about the data were made during the process. In the printout above, the read\_csv() function tells us the number of rows and columns in the data and the column specification, i.e., what type of data is contained within each column. In this case, we have 17 columns with character type of data, and 4 columns of double type of data. Functions in readr try to guess the column specification automatically, but it can also be specified manually when using the function. For more information about this dataset, please refer to Chapter 2 in this book [3].

Data from Excel worksheets can be read into R using the import() function from the rio [4] package. We will use synthetic data generated based on a real blended course of learning analytics for the remainder of this chapter. These data consist of three Excel files which we will first read into R.

```
library("rio")
url <- "https://github.com/lamethods/data/raw/main/1_moodleLAcourse/"
events <- import(paste0(url, "Events.xlsx"), setclass = "tibble")
results <- import(paste0(url, "Results.xlsx"), setclass = "tibble")
demographics <- import(paste0(url, "Demographics.xlsx"), setclass = "tibble")
```
The data files contain information on students' Moodle events, background information such as their name, study location and employment status, and various grades they've obtained during the course. For more information about the dataset, please refer to Chapter 2 in this book [3]. These data are read in the tibble [5] format, a special type of data.frame commonly used by tidyverse packages. We also load the dplyr [6] package which we will use for various tasks throughout this chapter.

library("dplyr")

#### **3 Grouping and Summarizing Data**

Instead of individual-level data metrics, we may be interested in specific groups as specified by some combination of data values. For example, we could compute the number of students studying in each location by gender. To accomplish this, we need to start by creating a grouped dataset with the function group\_by(). To compute the number of students, we can use the summarise() function which we already used previously in Chapter 1 and the function count(), which simply returns the number of observations in each category of its argument.

```
demographics |>
 group_by(Gender) |>
 count(Location)
# A tibble: 4 x 3
# Groups: Gender [2]
 Gender Location n
 <chr> <chr> <int>
1 F On campus 55
2 F Remote 10
3 M On campus 51
4 M Remote 14
```
The column n lists the number of students in each group. When a tibble that contains grouped data is printed into the console, the grouping variable and the number of groups will be displayed in the console below the dimensions. Next, we will compute the total number of Moodle events of each student, which we will also use in the subsequent sections.

```
events_summary <- events |>
  group_by(user) |>
  tally() |>
  rename(Frequency.Total = n)
events_summary
```

Here, the function tally() simply counts the number of number rows in the data related to each student, reported in the column n which we rename to Frequency.Total with the rename() function. We could also count the number of events by event type for each student

```
events |>
 group_by(user, Action) |>
 count(Action)
# A tibble: 1,439 x 3
# Groups: user, Action [1,439]
 user Action n
 <chr> <chr> <int>
1 00a05cc62 Applications 2
2 00a05cc62 Assignment 121
3 00a05cc62 Course_view 103
4 00a05cc62 Feedback 7
5 00a05cc62 General 10
# i 1,434 more rows
```
#### **4 Selecting Variables**

In the tidyverse paradigm, selecting columns, i.e., variables from data is done using the select() function. The select() function is very versatile, allowing the user to carry out selections ranging from simple selection of a single variable to highly complicated selections based on multiple criteria. The most basic selection selects only a single variable in the data based on its name. For example, we can select the employment statuses of students as follows

```
demographics |>
  select(Employment)
# A tibble: 130 x 1
  Employment
  <chr>
1 None
2 None
3 None
4 None
5 Part-time
# i 125 more rows
```
Note that using select() with a single variable is not the same as using the \$ symbol to select a variable, as the result is still a tibble, select() simply produces a subset of the data, where only the selected columns are present. Select is more similar to subset() in base R, which can accomplish similar tasks as select() and filter() in the tidyverse. However, we do not recommend using subset(), as it may not work correctly when the working environment has variables that have the same name as columns in the data, which can lead to undesired outcomes.

To extract the values of the selected column as a vector, we can use the function pull(). We use the head() function here to limit the console output to just the first few values of the vector (default is 6 values).

```
demographics |>
  pull(Employment) |>
  head()
```
[1] "None" "None" "None" "None" "Part-time" "Part-time"

The select() function syntax supports several operations that are similar to base R. We can select ranges of consecutive variables using :, complements using !, and combine selections using c(). The following selections illustrate some of these features:

```
demographics |>
    select(user:Origin)
# A tibble: 130 x 4
 user Name Surname Origin
 <chr> <chr> <chr> <chr>
1 6eba3ff82 Amanda Mora Costa Rica
2 05b604102 Lian Abdullah Yemen
3 111422ee7 Bekim Krasniqi Kosovo
4 b4658c3a9 Yusuf Kaya Turkey
5 e6ec47f29 Zoran Babić Serbia
# i 125 more rows
  demographics |>
    select(!Gender)
# A tibble: 130 x 7
 user Name Surname Origin Birthdate Location Employment
 <chr> <chr> <chr> <chr> <chr> <chr> <chr>
```

```
1 6eba3ff82 Amanda Mora Costa Rica 28.2.1998 On campus None
2 05b604102 Lian Abdullah Yemen 19.11.1996 On campus None
3 111422ee7 Bekim Krasniqi Kosovo 30.1.1999 On campus None
4 b4658c3a9 Yusuf Kaya Turkey 16.6.1998 On campus None
5 e6ec47f29 Zoran Babić Serbia 29.10.1998 On campus Part-time
# i 125 more rows
  demographics |>
    select(c(user, Surname))
# A tibble: 130 x 2
  user Surname
  <chr> <chr>
1 6eba3ff82 Mora
2 05b604102 Abdullah
3 111422ee7 Krasniqi
4 b4658c3a9 Kaya
5 e6ec47f29 Babić
# i 125 more rows
```
In the first selection, we select all variables starting from user on the left to Origin on the right. In the second, we select all variables except Gender. In the third, we select both user and Surname variables.

Sometimes, our selection might not be based directly on the variable names themselves as in the examples above but instead on vectors that contain the names of columns we wish to select. In such cases, we can use the function all\_of(). We can consider the intersections or unions of such selections using & and |, respectively.

```
cols_a <- c("user", "Name", "Surname")
cols_b <- c("Surname", "Origin")
demographics |>
  select(all_of(cols_a))
```

```
# A tibble: 130 x 3
 user Name Surname
 <chr> <chr> <chr>
1 6eba3ff82 Amanda Mora
2 05b604102 Lian Abdullah
3 111422ee7 Bekim Krasniqi
4 b4658c3a9 Yusuf Kaya
5 e6ec47f29 Zoran Babić
# i 125 more rows
```

```
demographics |>
    select(all_of(cols_a) & all_of(cols_b))
# A tibble: 130 x 1
 Surname
 <chr>
1 Mora
2 Abdullah
3 Krasniqi
4 Kaya
5 Babić
# i 125 more rows
  demographics |>
    select(all_of(cols_a) | all_of(cols_b))
# A tibble: 130 x 4
 user Name Surname Origin
 <chr> <chr> <chr> <chr>
1 6eba3ff82 Amanda Mora Costa Rica
2 05b604102 Lian Abdullah Yemen
3 111422ee7 Bekim Krasniqi Kosovo
4 b4658c3a9 Yusuf Kaya Turkey
5 e6ec47f29 Zoran Babić Serbia
# i 125 more rows
```
Often the names of data variables follow a similar pattern, and these patterns can be used to construct selections. Selections based on a prefix or a suffix in the variable name can be carried out with the functions starts\_with() and ends\_with(), respectively. The function contains() is used to look for a specific substring in the names of the variables, and more complicated search patterns can be defined with the function matches() that uses regular expressions (see ?tidyselect::matches for further information).

```
results |>
 select(starts_with("Grade"))
# A tibble: 130 x 13
 Grade.SNA_1 Grade.SNA_2 Grade.Review Grade.Group_self Grade.Group_All
    <dbl> <dbl> <dbl> <dbl> <dbl>
1 0 0 6.67 5 4
2 8 10 6.67 1 3
```

```
3 10 10 10 10 9.11
4 5 50 1 4
5 10 10 10 10 9.18
# i 125 more rows
# i 8 more variables: Grade.Excercises <dbl>, Grade.Project <dbl>,
# Grade.Literature <dbl>, Grade.Data <dbl>, Grade.Introduction <dbl>,
# Grade.Theory <dbl>, Grade.Ethics <dbl>, Grade.Critique <dbl>
results |>
 select(contains("Data"))
# A tibble: 130 x 1
 Grade.Data
     <dbl>
1 4
2 3
3 5
4 3
5 5
# i 125 more rows
```
So far, our selections have been based on variable names, but other conditions for selection are also feasible. The general-purpose helper function where() is used to select those variables for which a function provided to it returns TRUE. For example, we could select only those columns that contain character type data or double type data.

```
results |>
  select(where(is.character))
# A tibble: 130 x 1
  user
  <chr>
1 6eba3ff82
2 05b604102
3 111422ee7
4 b4658c3a9
5 e6ec47f29
# i 125 more rows
results |>
  select(where(is.double))
```

```
# A tibble: 130 x 14
 Grade.SNA_1 Grade.SNA_2 Grade.Review Grade.Group_self Grade.Group_All
    <dbl> <dbl> <dbl> <dbl> <dbl>
1 0 0 6.67 5 4
2 8 10 6.67 1 3
3 10 10 10 10 9.11
4 5 50 1 4
5 10 10 10 10 9.18
# i 125 more rows
# i 9 more variables: Grade.Excercises <dbl>, Grade.Project <dbl>,
# Grade.Literature <dbl>, Grade.Data <dbl>, Grade.Introduction <dbl>,
# Grade.Theory <dbl>, Grade.Ethics <dbl>, Grade.Critique <dbl>, Final_grade
<dbl>
```
#### **5 Filtering Observations**

In contrast to selection which relates to obtaining a subset of the columns of the data, filtering refers to obtaining a subset of the rows. In the tidyverse, data filtering is carried out with the dplyr package function filter(), which should not be confused with the base R filter() function in the stats package. As we have attached the dplyr package, the base R filter() function is masked, meaning that when we write code that uses filter(), the dplyr version of the function will automatically be called.

Filtering is often a much simpler operation than selecting variables, as the filtering conditions are based solely on the values of the data variables. Using filter() is analogous to the base R subset operator [, but the filtering condition is given as an argument to the filter() function instead. It is good to remind that in R a single equal sign (=) is merely for arguments of function calls, while double equal sign (==) is needed for comparison of two values. And example of filter:

```
demographics |>
    filter(Origin == "Bosnia") |>
    select(Name, Surname)
# A tibble: 2 x 2
  Name Surname
  <chr> <chr>
1 Hamza Hodžić
2 Davud Delić
```
The code above first filters our student demographics data to only those students whose country of origin is Bosnia. Then, we select their first and last names.

Multiple filtering conditions can be refined and combined using base R logical operators, such as & and |.

```
demographics |>
  filter(Gender == "F" & Location == "Remote")
# A tibble: 10 x 8
 user Name Surname Origin Gender Birthdate Location Employment
 <chr> <chr> <chr> <chr> <chr> <chr> <chr> <chr>
1 d93f7f0d3 Zahra Gul Afghanistan F 22.11.1999 Remote None
2 93d1f2e82 Louise Bernard France F 5.9.1998 Remote Part-time
3 417892918 Miora Rakotomalala Madagascar F 9.12.1995 Remote None
4 f98e6e2b8 Linda Mensah Ghana F 7.2.1991 Remote None
5 590846fe3 Lucija Horvat Croatia F 22.7.1998 Remote None
# i 5 more rows
```
Here, we filtered our data to female students who are studying remotely. The same result could also be obtained by using the filter() function two times

```
demographics |>
  filter(Gender == "F") |>
  filter(Location == "Remote")
```

This type of approach may improve the readability of your code especially when there are several independent filtering conditions to be applied simultaneously.

Filters can naturally be based on numeric values as well. For example, we could select those students whose final grade is higher than 8.

```
results |>
  filter(Final_grade > 8)
```

```
# i 9 more variables: Grade.Excercises <dbl>, Grade.Project <dbl>,
# Grade.Literature <dbl>, Grade.Data <dbl>, Grade.Introduction <dbl>,
# Grade.Theory <dbl>, Grade.Ethics <dbl>, Grade.Critique <dbl>, Final_grade <dbl>
```
Similarly, we could select students based on their total number of Moodle events.

```
events_summary |>
  filter(Frequency.Total > 100 & Frequency.Total < 500)
```

```
# A tibble: 44 x 2
 user Frequency.Total
 <chr> <int>
1 00a05cc62 417
2 046c35846 199
3 05b604102 199
4 0604ff3d3 436
5 0af619e4b 268
# i 39 more rows
```
#### **6 Transforming Variables**

In the best-case scenario, our data is already in the desired format after it has been read into R, but this is rarely the case with real datasets. We may need to compute new variables that were not present in the original data, convert measurements to different units, or transform text data into a numeric form. In the tidyverse, data transformations are carried out by the mutate() function of the dplyr package. This function can be used to transform multiple variables at the same time or to construct entirely new variables. The syntax of the function is the same in both cases: first, the name of the variable should be provided followed by an R expression that defines the variable. The transformed data is not automatically assigned to any variable, enabling transformations to be used as temporary variables within a chain of piped operations.

As a simple example, we could convert the students' locations into a factor variable.

```
demographics |>
    mutate(Location = factor(Location))
# A tibble: 130 x 8
 user Name Surname Origin Gender Birthdate Location Employment
 <chr> <chr> <chr> <chr> <chr> <chr> <fct> <chr>
1 6eba3ff82 Amanda Mora Costa Rica F 28.2.1998 On campus None
```

```
2 05b604102 Lian Abdullah Yemen F 19.11.1996 On campus None
3 111422ee7 Bekim Krasniqi Kosovo M 30.1.1999 On campus None
4 b4658c3a9 Yusuf Kaya Turkey M 16.6.1998 On campus None
5 e6ec47f29 Zoran Babić Serbia M 29.10.1998 On campus Part-time
# i 125 more rows
```
As we see from the tibble printout, the Location variable is a factor in the transformed data as indicated by the <fct> heading under the variable name. Note that the original demographics data was not changed, as we did not assign the result of the computation.

The gender and employment status of the students could also be used as factors, which we could do in a single mutate() call

```
demographics |>
  mutate(
    Gender = factor(Gender),
    Location = factor(Location),
    Employment = factor(Employment)
  )
```

```
# A tibble: 130 x 8
```

However, writing out individual identical transformations manually is cumbersome when the number of variables is large. For such cases, the across() function can be leveraged, which applies a function across multiple columns. This function uses the same selection syntax that we already familiarized ourselves with earlier to define the columns that will be transformed. To accomplish the same three transformations into a factor format, we could write

```
demographics |>
  mutate(across(c(Gender, Location, Employment), factor))
```

The first argument to the across() function is the selection that defines the variables to be transformed. The second argument defines the transformation, in this case, a function, to be used.

Working with dates can often be challenging. When we read the student demographic data into R, the variable Birthdate was assumed to be a character type variable. If we would like to use this variable to e.g., compute the ages of the students, we need to first convert it into a proper format using the as.Date function. Since the dates in the data are not in any standard format, we must provide the format manually. Afterwards, we can use the lubridate [7] package to easily compute the ages of the students, which we will save into a new variable called Age. We will also construct another variable called FullName which formats the first and last names of the students as "Last, First".

```
library("lubridate")
demographics |>
  mutate(
    Birthdate = as.Date(Birthdate, format = "%d.%m.%Y"),
    Age = year(as.period(interval(start = Birthdate, end = date("2023-03-12")))),
    FullName = paste0(Surname, ", ", Name)
  ) |>
  select(Age, FullName)
```

```
# A tibble: 130 x 2
   Age FullName
 <dbl> <chr>
1 25 Mora, Amanda
2 26 Abdullah, Lian
3 24 Krasniqi, Bekim
4 24 Kaya, Yusuf
5 24 Babić, Zoran
# i 125 more rows
```
The computation of the ages involves several steps. First, we construct a time interval object with the interval() function from the birthdate to the date for which we wish to compute the ages. Next, the as.period() function converts this interval into a time duration, from which we lastly get the number of years with the year() function.

Suppose that we would like to construct a new variable AchievingGroup that categorizes the students into top 50% achievers and bottom 50% achievers based on their final grade on the course. We leverage two functions from the dplyr package to construct this new variable: case\_when() and ntile(). The function case\_when() is used to transform variables based on multiple sequential conditions. The function ntile() has two arguments, a vector x and an integer n, and it splits x into n equal-sized groups based on the ranks of the values in x.

```
results <- results |>
  mutate(
    AchievingGroup = factor(
      case_when(
        ntile(Final_grade, 2) == 1 ~ "Low achiever",
        ntile(Final_grade, 2) == 2 ~ "High achiever"
      )
    )
  )
```
The syntax of case\_when() is very simple: we describe the condition for each case followed by ~ after which we define the value that the case should correspond to. We assign the result of the computation to the results data, as we will be using the AchievingGroup variable in later chapters.

We would also like to categorize the students based on their activity level, i.e., the number of total Moodle events. Our goal is to create three groups of equal size consisting of low activity, moderate activity and high activity students. The approach we applied to categorizing the achievement level of the students is also applicable for this purpose. We name our new variable as ActivityGroup, and we assign the result of the computation, as we will also be using this variable in later chapters.

```
events_summary <- events_summary |>
  mutate(
    ActivityGroup = factor(
      case_when(
        ntile(Frequency.Total, 3) == 1 ~ "Low activity",
        ntile(Frequency.Total, 3) == 2 ~ "Moderate activity",
        ntile(Frequency.Total, 3) == 3 ~ "High activity"
      )
    )
  )
```
#### **7 Rearranging Data**

Sometimes we may want to reorder the rows or columns of our data, for example in alphabetical order based on the names of students on a course. The arrange() function from the dplyr package orders the rows of the by the values of columns selected by the user. The values are sorted in ascending order by default, but the order can be inverted by using the desc() function if desired. The variable order in the selection defines how ties should be broken when duplicate values are encountered in the previous variables of the selection. For instance, the following code would arrange the rows of our demographics data by first comparing the surnames of the students, and then the given names for those students with the same surname. Missing values are placed last in the reordered data.

```
demographics |>
  arrange(Surname, Name)
# A tibble: 130 x 8
 user Name Surname Origin Gender Birthdate Location Employment
 <chr> <chr> <chr> <chr> <chr> <chr> <chr> <chr>
1 ba76ebfab Bismah Abbasi Pakistan F 2.4.1996 Remote Full-time
2 05b604102 Lian Abdullah Yemen F 19.11.1996 On campus None
3 d2c3ce9a4 Amir Ait Morocco M 19.6.1997 On campus None
4 68a377c82 Saliha Akmatova Kyrgyzstan F 19.5.1999 Remote None
5 7e2726f3c Kazi Akter Bangladesh M 22.12.1992 On campus None
# i 125 more rows
```
A descending order based on both names can be obtained by applying the desc() function.

```
demographics |>
  arrange(desc(Surname), desc(Name))
# A tibble: 130 x 8
 user Name Surname Origin Gender Birthdate Location Employment
 <chr> <chr> <chr> <chr> <chr> <chr> <chr> <chr>
1 a48165ad5 Liam Zambrano Ecuador M 4.4.1998 On campus None
2 1115dae61 Poema Wong Tahiti F 22.1.1999 Remote Part-time
3 0ef305578 Jack White Australia M 22.4.1995 Remote None
4 f753ce9bf Dechen Wangmo Bhutan F 29.4.1999 On campus None
5 f87eaa00c Prasert Wang Thailand M 9.4.1997 On campus None
# i 125 more rows
```
Column positions can be changed with the relocate() function of the dplyr package. Like arrange(), we first select the column or columns we wish to move into a different position in the data. Afterwards, we specify the position where the columns should be moved to in relation to positions of the other columns. In our demographics data, the user ID column user is the first column. The following code moves this column after the Employment column so that the user column becomes the last column in the data.

```
demographics |>
  relocate(user, .after = Employment)
# A tibble: 130 x 8
 Name Surname Origin Gender Birthdate Location Employment user
 <chr> <chr> <chr> <chr> <chr> <chr> <chr> <chr>
1 Amanda Mora Costa Rica F 28.2.1998 On campus None 6eba3ff82
2 Lian Abdullah Yemen F 19.11.1996 On campus None 05b604102
3 Bekim Krasniqi Kosovo M 30.1.1999 On campus None 111422ee7
4 Yusuf Kaya Turkey M 16.6.1998 On campus None b4658c3a9
5 Zoran Babić Serbia M 29.10.1998 On campus Part-time e6ec47f29
# i 125 more rows
```
The mutually exclusive arguments .before and .after of relocate() specify the new column position in relation to columns that were not selected. These arguments also support the select() function syntax for more general selections.

#### **8 Reshaping Data**

Broadly speaking, tabular data typically take one of two formats: wide or long. In the wide format, there is one row per subject, where the subjects are identified by an identifier variable, such as the user variable in our Moodle data, and multiple columns for each measurement. In the long format, there are multiple rows per subject, and the columns describe the type of measurement and its value. For example, the events data is in a long format containing multiple Moodle events per student, but the results and demographics data are in a wide format with one row per student.

In the previous section, we constructed a summary of the users' Moodle events in total and of different types. The latter data is also in a long format with multiple rows per subject, but we would instead like to have a column for each event type with one row per user, which means that we need to convert this data into a wide format. Conversion between the two tabular formats is often called *pivoting*, and the corresponding functions pivot\_wider() and pivot\_longer() from the tidyr [8] package are also named according to this convention. We will create a wide format data of the counts of different event types using the pivot\_wider() function as follows

```
library("tidyr")
events_types <- events |>
 group_by(user, Action) |>
 count(Action) |>
 pivot_wider(
   names_from = "Action",
   names_prefix = "Frequency.",
   values_from = "n",
   values_fill = 0
 )
events_types
# A tibble: 130 x 13
# Groups: user [130]
 user Frequency.Applications Frequency.Assignment Frequency.Course_view
 <chr> <int> <int> <int>
1 00a05cc62 2 121 103
2 042b07ba1 0 62 294
3 046c35846 0 41 53
4 05b604102 0 44 49
5 0604ff3d3 0 9 119
# i 125 more rows
# i 9 more variables: Frequency.Feedback <int>, Frequency.General <int>,
# Frequency.Group_work <int>, Frequency.Instructions <int>,
# Frequency.La_types <int>, Frequency.Practicals <int>, Frequency.Social <int>,
# Frequency.Ethics <int>, Frequency.Theory <int>
```
Here, we first specify the column name that the names of the wide format data should be taken from in the long format data with names\_from. In addition, we specify a prefix for the new column names using names\_prefix that helps to distinguish what these new columns will contain, but in general, the prefix is optional. Next, we specify the column that contains the values for the new columns with values\_from. Because not every student necessarily has events of every type, we also need to specify what the value should be in cases where there are no events of a particular type by using values\_fill. As we are considering the frequencies of the events, it is sensible to select 0 to be this value. We save the results to events\_types as we will use the event type data in later sections and chapters.

#### **9 Joining Data**

Now that we have computed the total number of events for each student and converted the event type data into a wide format, we still need to merge these new data with the demographics and results data. Data merges are also called *joins*, and the dplyr package provides several functions for different kinds of joins. Here, we will use the left\_join() function that will preserve all observations of the first argument.

```
left_join(demographics, events_summary, by = "user")
# A tibble: 130 x 10
 user Name Surname Origin Gender Birthdate Location Employment Frequency.Total
 <chr> <chr> <chr> <chr> <chr> <chr> <chr> <chr> <int>
1 6eba3ff~ Aman~ Mora Costa~ F 28.2.1998 On camp~ None 312
2 05b6041~ Lian Abdull~ Yemen F 19.11.19~ On camp~ None 199
3 111422e~ Bekim Krasni~ Kosovo M 30.1.1999 On camp~ None 532
4 b4658c3~ Yusuf Kaya Turkey M 16.6.1998 On camp~ None 246
5 e6ec47f~ Zoran Babić Serbia M 29.10.19~ On camp~ Part-time 356
# i 125 more rows
# i 1 more variable: ActivityGroup <fct>
```
In essence, the above left join adds all columns from events\_summary to demographics whenever there is a matching value in the by column. To continue, we can use additional left joins to add the remaining variables from the results data, and the Moodle event counts of different types from events\_types to have all the student data in a single object.

```
all_combined <- demographics |>
  left_join(events_types, by = "user") |>
  left_join(events_summary, by = "user") |>
  left_join(results, by = "user")
all_combined
```

We will use this combined dataset in the following chapters as well.

#### **10 Missing Data**

Sometimes it occurs that learning analytics data has cells for which the values are missing for some reason. The Moodle event data which we have utilized in this chapter does not naturally contain missing data. Thus, to have an example, we need to create a data which does. Second, handling of missing data is a vast topic of which we can only discuss some of the key points very briefly from a practical perspective. For a more comprehensive overview, we recommend reading [9] and [10] for a hands on approach. A short overview of missingness can be found in [11].

The code below will create missing values randomly to each column of events\_types data (user column is an exception). To do that, we use the mice [12] package which also has methods for the handling of missing data. Unfortunately, mice is not part of the tidyverse. For more information about mice, a good source is miceVignettes at https://www.gerkovink.com/ miceVignettes/. Now, let's create some missing data.

```
library("mice")
set.seed(44)
events_types <- events_types |>
  rename(
    "Ethics" = "Frequency.Ethics",
    "Social" = "Frequency.Social",
    "Practicals" = "Frequency.Practicals"
  )
ampute_list <- events_types |>
  ungroup(user) |>
  select(Ethics:Practicals)|>
  as.data.frame() |>
  ampute(prop = 0.3)
events_types_mis <- ampute_list$amp |>
  as_tibble()
events_types_mis[2, "Practicals"] <- NA
```
Above, we also rename the variables that contain the frequencies of Moodle events related to ethics, social and practicals into Ethics, Social and Practicals, respectively. Let's now see some of the values of events\_types\_mis

events\_types\_mis

```
# A tibble: 130 x 3
 Ethics Social Practicals
  <int> <int> <int>
1 NA 12 89
```

We can see that now the data contains NA values in some of the cells. These are the cells in which a missing value occurs, meaning that a value for those measurements has not been recorded. A missing data pattern, that is how missing of one variable affects missingness of other variables, can be show as:

```
md.pattern(events_types_mis, rotate.names = TRUE)
```
Above, each red square indicates a missing value while blue squares stand for observed ones. We can see that there are 95 complete rows, 10 for which Practicals are missing, 17 have missingness on Social and 9 are missing on Ethics. Also, one row has two red squares indicating a missing value on both Social and Practicals.

Let's now discuss options of handling missing data briefly. There are four classes of statistical methods for analyzing data with missing values: complete case (CC) methods, weighting methods, imputation methods, and model-based methods. The simplest of these is complete case analysis, which leaves missing values out of the analysis and only uses observations with all variables recorded. This can be done with the tidyr [8] package function drop\_na():

```
events_types_mis |>
 drop_na()
# A tibble: 95 x 3
 Ethics Social Practicals
  <int> <int> <int>
1 0 0 47
```

We can see that after using this method, our data has only 95 rows as those were the rows without any columns having missing values. This made our data much smaller! If there are a lot of missing values, the data may become too small to use for practical purposes.

A more novel group of methods are imputation methods. One of the options is using single imputation (SI) where the mean of each variable will determine the imputed value. The single mean imputation can be done as follows:

```
imp <- mice(events_types_mis, method = "mean",
m = 1, maxit = 1 , print = FALSE)
complete(imp) |>
  head()
```

We can see from above that the imputed values are not integers anymore. However, if we aim to estimate means or regression coefficients (see Chapter 5 [13] for details) that is not a problem. One of the problems with mean imputation is that the variance and standard error estimates will become downward biased. A mean of Ethics for mean imputation is:

```
fit <- with(imp, lm(Ethics ~ 1))
summary(fit)
# A tibble:1x6
 term estimate std.error statistic p.value nobs
 <chr> <dbl> <dbl> <dbl> <dbl> <int>
1 (Intercept) 7.55 0.811 9.32 4.20e-16 130
```
Next, let's briefly have a look at how we can utilize multiple imputation (MI) which is an improvement over single imputation. The multiple imputation approach generates more than one imputation thus creating many complete data sets for us. For each of these datasets, we can perform any analysis that we are interested in. After the analysis, one must pool the results from the impured datasets to get the final result. Here, we utilize a method called predictive mean matching (method = "pmm" in the code below), which uses the neighbour values of data as imputations.

```
imp2 <- mice(events_types_mis, method = "pmm",
m= 10, maxit = 100, print = FALSE)
fit2 <- with(imp2, lm(Ethics ~ Practicals))
pool_fit <- pool(fit2)
# Multiple imputation
summary(pool_fit)
```
term estimate std.error statistic df p.value 1 (Intercept) 1.62080372 2.18285149 0.7425167 101.7304 0.459485402 2 Practicals 0.08049328 0.02616765 3.0760606 106.0802 0.002668431 # Complete cases summary(lm(Ethics ~ Practicals, events\_types\_mis))["coefficients"] \$coefficients Estimate Std. Error t value Pr(>|t|) (Intercept) 2.15459162 2.12235051 1.015191 0.31226283 Practicals 0.06220884 0.02605001 2.388054 0.01865793 # Without missingness summary(lm(Ethics ~ Practicals, events\_types))["coefficients"] \$coefficients

Estimate Std. Error t value Pr(>|t|) (Intercept) 1.05409529 2.17821479 0.4839262 0.6292651258 Practicals 0.08891892 0.02590313 3.4327482 0.0008053447

From the results above, we can see that in this particular case the multiple imputation performs well in comparison to CC approach. The regression coefficient for full data without any missing values is 0*.*089, and it is 0*.*080 for multiple imputation, while complete case analysis gives 0*.*062. As all of them have very similar standard errors, this yields that MI and full data give statistically significant p-values for significance level 0.01, while CC does not.

#### **11 Correcting Erroneous Data**

Let's imagine that our data has an error on the surname variable Surname and that all the names ending with "sen" should end with "ssen". What we can do is that we can use regular expressions to detect the erroneous rows and we can also use them to replace the values. Let's first figure out which last names contain a name ending with "sen". We can use a function str\_detect() to return TRUE/FALSE for each row from stringr [14] package within a filter() function call. We define pattern = "sen\$" where \$ indicates the end of the string.

```
library("stringr")
demographics |>
  filter(str_detect(string = Surname, pattern = "sen$")) |>
  pull(Surname)
```

```
[1] "Nielsen" "Johansen" "Joensen" "Jansen" "Olsen"
```
After pulling the filtered surnames, there seems to be five surnames ending with "sen". Next, let's try to replace "sen" with "ssen". On the next row we filter just as previously to limit output.

```
demographics |>
  mutate(Surname = str_replace(
    string = Surname, pattern = "sen$", replacement = "ssen")
  ) |>
  filter(str_detect(string = Surname, pattern = "sen$")) |>
  pull(Surname)
```
[1] "Nielssen" "Johanssen" "Joenssen" "Janssen" "Olssen"

Thus, the following code updates the data so that all the surnames ending with "sen" now end with "ssen" instead.

```
demographics <- demographics |>
  mutate(Surname = str_replace(
    string = Surname, pattern = "sen$", replacement = "ssen")
  )
```
#### **12 Conclusion and Further Reading**

Data wrangling is one of the most important steps in any data analysis pipeline. This chapter introduced the tidyverse, tidy data, and several commonly used R packages for data manipulation and their use in basic scenarios in the context of learning analytics. However, the tidyverse is vast and can hardly be fully covered in a single chapter. We refer the reader to additional resources such as those found on the tidyverse website at https://www.tidyverse.org/learn/ and the book "R for Data Science" by Hadley Wicham and Garret Grolemund. The book is free to use and readily available online at https://r4ds.had.co.nz/.

#### **References**

**Open Access** This chapter is licensed under the terms of the Creative Commons Attribution 4.0 International License (http://creativecommons.org/licenses/by/4.0/), which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons license and indicate if changes were made.

The images or other third party material in this chapter are included in the chapter's Creative Commons license, unless indicated otherwise in a credit line to the material. If material is not included in the chapter's Creative Commons license and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder.

# **Introductory Statistics with R for Educational Researchers**

**Santtu Tikka, Juho Kopra, Merja Heinäniemi, Sonsoles López-Pernas, and Mohammed Saqr**

## **1 Introduction**

Learning analytics involves the practical application of statistical methods to quantitative data, which can represent various aspects of the learning process such as student engagement, progress, and outcomes. Thus, knowledge about basic statistical methods is essential. Let's start with how statistics connects with the research process and how it is also linked with philosophy of science. According to Niiniluoto [1], the research process can be described with the following eight steps:

Steps 1 and 2 require knowledge and skills related to the applied field, but also general scientific aptitude. Knowledge about statistics is central in steps 3, 4, 5, and

S. Tikka (-)

J. Kopra · S. López-Pernas · M. Saqr School of Computing, University of Eastern Finland, Joensuu, Finland

M. Heinäniemi Institute of Biomedicine, University of Eastern Finland, Kuopio, Finland

Department of Mathematics and Statistics, University of Jyväskylä, Jyväskylä, Finland e-mail<EMAIL_ADDRESS>
6. Finally, steps 7 and 8 mostly require skills in writing and communication. Overall, it can be argued that a solid understanding of statistics and statistical methods is crucial for anyone conducting research with quantitative data.

This chapter of the book concentrates on steps 4, 5, and 6 of the research process. We start with descriptive statistics, which are statistics that describe the overall features of the data. In contrast, inferential statistics are used to draw conclusions and make inferences about the population under study. Afterwards, we explain the basics of statistical hypothesis testing, which is the most common—although not the only—way to analyze data. The most common statistical tests, such as Student's t-test, Chi-squared test, Analysis of variance, Levene's test, and Shapiro-Wilk test are covered in this chapter. We also explain how to interpret the results of each test. We also present the linear regression model, which is not a statistical test but one of the most powerful statistical tools. Basic understanding of linear regression is essential for anyone interested in more advanced regression techniques, including logistic regression which is covered in the final section of this chapter. For a more in-depth view on the statistical tests covered in this chapter, we refer the reader to works such as [2, 3].

#### **2 Descriptive Statistics**

Descriptive statistics are used to provide a meaningful quantitative overview of data, and to summarize potentially vast amounts of information into more easily comprehensible and manageable quantities. In general, descriptive statistics are used as a first step in a data analysis workflow. In this section we will focus on numeric descriptors while visualizations are the topic of Chapter 6 [4]. For this chapter, we will use the combined Moodle data with students' demographics, results, and summarized Moodle event activity. For more information about the dataset, please refer to Chapter 2 in this book [5]. We begin by installing all R packages that we will use in this chapter.

```
install.packages(
  c("car", "rio", "see", "dplyr", "tidyr",
    "broom", "report", "correlation", "performance")
)
```
We use the rio [6] package to read the data files into R via the import() function. The Events dataset contains log data on the student's Moodle activity such as Moodle event types and names. The Demographics dataset contains background information on the students such as their gender and location of study (categorical variables). Finally, the Results data contains grades on various aspects of the Moodle course including the final course grade (numeric variables). We also create a new variable called AchievingGroup that categorizes the students

into bottom and top 50% of achievers in terms of the final grade. We will leverage the dplyr [7] and tidyr [8] packages to wrangle the data into a single combined dataset. We begin by reading the data files and by constructing the AchievingGroup variable.

```
library("rio")
library("dplyr")
library("tidyr")
url <- "https://github.com/lamethods/data/raw/main/1_moodleLAcourse/"
events <- import(paste0(url, "Events.xlsx"), setclass = "tibble")
demographics <- import(paste0(url, "Demographics.xlsx"), setclass = "tibble")
results <- import(paste0(url, "Results.xlsx"), setclass = "tibble") |>
  mutate(
    AchievingGroup = factor(
      case_when(
        ntile(Final_grade, 2) == 1 ~ "Low achiever",
        ntile(Final_grade, 2) == 2 ~ "High achiever"
      )
    )
  )
```
Next, we summarize the student's engagement based on their Moodle activity into three groups: Low activity, Moderate activity and High Activity.

```
events_summary <- events |>
  group_by(user) |>
  tally() |>
  rename(Frequency.Total = n) |>
  mutate(
    ActivityGroup = factor(
      case_when(
        ntile(Frequency.Total, 3) == 1 ~ "Low activity",
        ntile(Frequency.Total, 3) == 2 ~ "Moderate activity",
        ntile(Frequency.Total, 3) == 3 ~ "High activity"
      )
    )
  )
```
We also count the different types of Moodle events.

```
events_types <- events |>
  group_by(user, Action) |>
  count(Action) |>
  pivot_wider(
    names_from = "Action",
```

```
names_prefix = "Frequency.",
  values_from = "n",
  values_fill = 0
)
```
Finally, we combine the data.

```
all_combined <- demographics |>
  left_join(events_types, by = "user") |>
  left_join(events_summary, by = "user") |>
  left_join(results, by = "user")
```
The various steps of the combined data construction are discussed in greater detail in Chapter 4 [9].

#### *2.1 Measures of Central Tendency*

A typical way to summarize a univariate data sample is to describe its "middle point" using an appropriate statistic depending on the measurement scale of the data. The most common statistics to describe such a value are the *mean*, the *median*, and the *mode*.

For data on the interval or ratio scales (and sometimes also on the ordinal scale), the most common option is to use the arithmetic mean, which is available via the base R function mean(). This function takes a vector of values as its input.

```
all_combined |>
  summarise(
   mean_grade = mean(Final_grade),
   mean_total = mean(Frequency.Total)
  )
# A tibble: 1 x 2
 mean_grade mean_total
      <dbl> <dbl>
1 7.25 736.
```
The means are reported as a tibble with a single column for both variables.

For data on the ordinal scale (or interval or ratio scales), the median can be used which is defined as the value that separates the lower half from the bottom half of the data sample, i.e., the 50% quantile. The median can be computed in R using the built-in median() function. Similarly to mean(), this function also takes a vector of values as its input.

```
all_combined |>
  summarise(
   median_grade = median(Final_grade),
   median_total = median(Frequency.Total)
  )
# A tibble: 1 x 2
 median_grade median_total
        <dbl> <dbl>
1 7.95 606
```
Just like before, the medians are reported as a tibble with each value in its own column.

For data on the nominal or ordinal scale, the mode is a suitable choice as it describes the category with the highest number of observations. Unfortunately, there is no readily available function in R to compute the mode, and the reader should take care not to mistakenly use the mode() function, which is used to determine the internal storage mode of a variable (similar to the typeof() function). However, we can easily write our own function to compute the statistical mode as follows:

```
stat_mode <- function(x) {
  u <- unique(x)
  u[which.max(tabulate(match(x, u)))]
}
```
Functions in R are written using the following syntax. First, we define the name of the function, just like we would define the name of a variable when assigning data into it, in this case the name is stat\_mode. Next, we assign the function definition, which starts with the keyword function. Next, we describe the function arguments within the parentheses, in this case we call our argument x, which we assume contains the data vector we wish to compute the mode for. Next, we define the body of the function within the braces. The body determines what the function does and what its output should be. Within the body, we first determine the unique values in the data vector x, and assign the result to a variable u. Next, we need to count the number of occurrences of each unique value. To start, we first match each observed value in x to the unique values in u to get the corresponding indices, which we will then count using tabulate. We obtain the index of the value with the highest number of occurrences with the function which.max(), and finally the corresponding unique value by selecting it from u using the subset operator, i.e., the brackets. Our function will now work on all types of data.

```
all_combined |>
  summarise(
    mode_gender = stat_mode(Gender),
```

```
mode_location = stat_mode(Location)
  )
# A tibble: 1 x 2
 mode_gender mode_location
 <chr> <chr>
1 F On campus
```
The output is now similar to the mean and median functions that we used earlier, showing the modes of Gender and Location as a two-column tibble. For nominal variables, it is common to also compute the frequencies of each category. This can easily be done with the base R function table()

```
table(all_combined$Gender)
F M
65 65
table(all_combined$Location)
On campus Remote
     106 24
```
The function outputs the names of the categories and the frequency of each category as an integer vector.

#### *2.2 Measures of Dispersion*

For data on the interval and ratio scales (and sometimes also on the ordinal scale), it is also meaningful to describe how clustered or scattered the values in the sample are, i.e., how far apart the values are from one another. Commonly used measures of statistical dispersion include the *variance*, *standard deviation*, and the *interquartile range*. Typically, such measures have the value 0 when all values in the sample are identical, and the value increases as the dispersion in the data grows.

All three measures can be readily computed with built-in R functions var(), sd(), and IQR() respectively. Like mean() and median(), these functions accept a vector or numeric values as input.

```
all_combined |>
  summarise(
    var_grade = var(Final_grade),
```

```
sd_grade = sd(Final_grade),
   iqr_grade = IQR(Final_grade)
 )
# A tibble: 1 x 3
 var_grade sd_grade iqr_grade
     <dbl> <dbl> <dbl>
1 4.81 2.19 3.34
```
The variance, standard deviation and interquartile range of the final grade are returned as a tibble with three columns.

#### *2.3 Covariance and Correlation*

Covariance and correlation measure the linear dependence between two variables. Correlation is a unitless measure between −1 and 1, whereas covariance is not, and its scale depends on the scale of the variables. The sign of both measures indicates the tendency of the relationship. Positive sign means that as the value of one variable increases, the value of the other variable tends to increase as well. Conversely, a negative sign indicates that the value of the second variable tends to decrease as the value of the first variable increases.

Both covariance and correlation and be computed directly in R using the functions cov() and cor(), respectively.

```
all_combined |>
 summarise(
   cov_grade_total = cov(Final_grade, Frequency.Total),
   cor_grade_total = cor(Final_grade, Frequency.Total),
 )
# A tibble: 1 x 2
 cov_grade_total cor_grade_total
          <dbl> <dbl>
1 504. 0.513
```
We obtain the covariance and correlation between the final grade and total number of Moodle events. We will familiarize ourselves with correlations in greater depth in Sect. 4.

#### *2.4 Other Common Statistics*

The extreme values of a data sample can be found using the function range(), which computes the minimum and maximum values of the sample as a vector. These

values can also be computed individually with the corresponding functions min() and max(). Because summarise() only allows a single value as output per row, we use the reframe() function instead when computing the range.

```
all_combined |>
  reframe(
    range_grade = range(Final_grade)
  )
# A tibble: 2 x 1
 range_grade
       <dbl>
1 0
2 10
all_combined |>
  summarise(
    min = min(Final_grade),
    max = max(Final_grade)
  )
# A tibble: 1 x 2
   min max
 <dbl> <dbl>
1 0 10
```
With reframe(), we obtain a tibble with two rows, the first containing the minimum value and the second the maximum value of the final grade. If we instead use summarise() like before, we can only obtain one value per computed variable. The summary() function can also be used to quickly compute several of the most common descriptive statistics for all variables of a dataset.

```
results |>
  select(Grade.SNA_1:Grade.Group_self) |>
  summary()
 Grade.SNA_1 Grade.SNA_2 Grade.Review Grade.Group_self
Min. : 0.000 Min. : 0.000 Min. : 0.000 Min. : 0.000
1st Qu.: 8.000 1st Qu.: 9.000 1st Qu.: 6.670 1st Qu.: 9.000
Median : 9.000 Median :10.000 Median : 8.000 Median :10.000
Mean : 8.346 Mean : 9.262 Mean : 7.724 Mean : 8.085
3rd Qu.:10.000 3rd Qu.:10.000 3rd Qu.: 9.670 3rd Qu.:10.000
Max. :10.000 Max. :10.000 Max. :10.000 Max. :10.000
```
The output shows the minimum and maximum values, the quartiles, and the mean of each variable that we selected.

#### **3 Statistical Hypothesis Testing**

Statistical hypothesis testing aims to evaluate hypotheses about a population of interest using probabilistic inference. The starting point of any statistical test is a socalled *null hypothesis* (denoted by *H*0), which typically corresponds to a scenario, where evidence supporting a specific hypothesis is a result of pure chance. For example, when evaluating whether a new drug is an efficient form of treatment via a randomized controlled trial, the null hypothesis could be that the drug has no effect on the response. A null hypothesis is always associated with an *alternative hypothesis* (denoted by *H*1), which is typically the inverse of the null hypothesis and corresponds to the hypothesis of interest, e.g., the drug has an effect on the response.

Statistical tests operate by assuming that the null hypothesis is true, and highly unlikely events under this assumption are typically regarded as giving cause for rejecting the null hypothesis. A statistical test is associated with a *test statistic*, which is a measure of how much the observations deviate from the null hypothesis scenario. The distribution of the test statistic under the null hypothesis and the sample test statistic can be used to compute the probability of obtaining a test statistic as extreme or more extreme than the one observed, assuming that the null hypothesis is true. This probability is known as the *p-value*, which is often mischaracterized even in scientific literature. For instance, the p-value is not the probability that the null hypothesis is true or that the alternative hypothesis is false. The p-value also does not quantify the size of the observed effect, or its real-world importance.

Typically, a *confidence level* is decided before applying a statistical test (usually denoted by *α*), and the null hypothesis is rejected if the observed p-value is smaller than this confidence level. If the p-value is greater than the confidence level, the null hypothesis is not rejected. Traditionally, the confidence level is 0.05, but this convention varies by field, and should be understood as being arbitrary, i.e., there is nothing special about the value 0.05. If the p-value falls below the confidence level, the result is regarded as *statistically significant*.

Hypothesis testing is a powerful tool for drawing conclusions from data, but it is important to use it appropriately and to understand its limitations. Every statistical test is associated with a set of assumptions which are often related to the distribution of the data sample. If these assumptions are violated, then the results of the test may be unreliable. In the following sections, some of the most common statistical tests are introduced. We will take advantage of the report [10] package and the corresponding function report() to showcase the results of the various statistical tests. For more information on the concepts and principles related to statistical hypothesis testing, see e.g., [2, 3].

#### *3.1 Student's t-test*

Student's t-test [11] is one of the most well-known statistical tests. It compares the mean values of variables either between two populations or between a single population and a reference level and is thus applicable to continuous variables. The test assumes homogeneity of variance and that the data originates from a normal distribution. For nonhomogeneous data, the test can still be performed by using an approximation [12]. In R, all variants of the t-test test can be applied using the function t.test().

Our goal is to compare the students' Moodle activity with respect to their final grade. For this purpose, we use the binary variable called AchievingGroup which categorizes the students into top and bottom 50% achievers in terms of the final grade.

#### **3.1.1 One-Sample t-test**

The one-sample t-test compares the mean of a data sample against a reference value, typically defined by the null hypothesis. Let us begin by testing the hypothesis that the average number of Moodle events (Frequency.Total) is 600 (*H*<sup>0</sup> : *μ* = 600). The function t.test() can be used in various ways, but in this example we provide the function with a formula object Frequency.Total ~ 1 as the first argument. The formula syntax is a standard method for defining statistical models and other dependency structures in R. The formula defines that the left-hand side of the ~ symbol is a response variable which is explained by the terms on the right-hand side. Because we're not conducting the test with respect to any other variable, the right-hand side of the formula is simply 1, which means that it is a constant in the R formula syntax. This does not mean for example, that our null hypothesis would be that the number of Moodle events is 1. The expected value that the test is applied against (i.e., the value we assume *μ* to have under the null hypothesis) is defined via the argument mu, which by default has the value 0 for a one-sample t-test. Argument data defines in which environment the formula should be evaluated. By providing the all\_combined data, we do not have to explicitly extract the FrequencyTotal variable from the data in the formula by writing all\_combined\$Frequency.Total or by using pull(). This is especially useful when the formula contains several variables from the same data.

```
ttest_one <- t.test(Frequency.Total ~ 1, data = all_combined, mu = 600)
ttest_one
    One Sample t-test
data: Frequency.Total
t = 3.4511, df = 129, p-value = 0.0007553
```

```
alternative hypothesis: true mean is not equal to 600
95 percent confidence interval:
 657.8530 813.3163
sample estimates:
mean of x
 735.5846
```
As a result, we obtain the value of the test statistic (t = 3.4511), the degrees of freedom (df = 129), and the p-value of the test (p-value = 0.0007553). Because the p-value is very small (much smaller than the standard 0.05 confidence level), we reject the null hypothesis, which means that the average number of Moodle events is significantly different from 600. The output of the test result object also describes the alternative hypothesis *H*<sup>1</sup> under alternative hypothesis, and the confidence interval of the test statistic.

We produce a summarized report of the test results with the report() function.

```
Warning: Unable to retrieve data from htest object.
  Returning an approximate effect size using t_to_d().
Effect sizes were labelled following Cohen's (1988) recommendations.
The One Sample t-test testing the difference between Frequency.Total (mean =
735.58) and mu = 600 suggests that the effect is positive, statistically
significant, and small (difference = 135.58, 95% CI [657.85, 813.32], t(129) =
3.45, p < .001; Cohen's d = 0.30, 95% CI [0.13, 0.48])
```
This produces a description of the results of the test that is easier to read and interpret than the direct output of the test result object. We note that a warning is also produced which we can safely ignore in this case. The warning occurs because the test result object ttest\_one does not retain the original data all\_combined which we used to carry out the test. If non-approximate effect sizes are desired, the test should be carried out by supplying the variables being compared directly without using the formula interface of the t.test() function. For more information on the effect size, see e.g., [13, 14].

#### **3.1.2 Two-Sample t-test**

report(ttest\_one)

In contrast to the one-sample t-test, the two-sample t-test compares the means of two data samples against one another. For example, suppose that we're interested in the hypothesis that the average number of Moodle events is the same for the top and bottom 50% achievers (*H*<sup>0</sup> : *μ*<sup>1</sup> = *μ*2). We can once again leverage the formula syntax, but instead of the constant 1 on the right-hand side of the formula, we will now replace it with the variable Achievement which defines the achievement level.

```
ttest_two <- t.test(Frequency.Total ~ AchievingGroup, data = all_combined)
ttest_two
   Welch Two Sample t-test
data: Frequency.Total by AchievingGroup
t = 4.4749, df = 95.988, p-value = 2.102e-05
alternative hypothesis:
true difference in means between group 1 and group 2 is not equal to 0
95 percent confidence interval:
182.6427 473.8496
sample estimates:
mean in group High achiever mean in group Low achiever
                  899.7077 571.4615
```
The contents of the result object are mostly the same as in the case of the onesample t-test. The result is again statistically significant (using the 0.05 confidence level) meaning that according to the test, the average number of Moodle events is higher for the top 50% achievers. The report() function can be used to produce a similar summary as in the case of the one-sample t-test.

```
report(ttest_two)
Warning: Unable to retrieve data from htest object.
  Returning an approximate effect size using t_to_d().
Effect sizes were labelled following Cohen's (1988) recommendations.
The Welch Two Sample t-test testing the difference of Frequency.Total by
AchievingGroup (mean in group High achiever = 899.71, mean in group Low achiever =
571.46) suggests that the effect is positive, statistically significant, and large
(difference = 328.25, 95% CI [182.64, 473.85], t(95.99) = 4.47, p < .001; Cohen's d
= 0.91, 95% CI [0.49, 1.33])
```
This produces the same warning as before in the one-sample case, but we can safely ignore it again.

#### **3.1.3 Paired Two-Sample t-test**

Instead of directly comparing the means of two groups, it may sometimes be of interest to compare differences between pairs of measurements. Such a scenario typically arises in an experiment, where subjects are paired, or two sets of measurements are taken from the same subjects. While our Moodle event data does not contain such measurement pairs, we could still imagine that our data was organized such that each student in the bottom 50% achievers was paired with a

student in the top 50% achievers and that there is a one-to-one correspondence between the two achievement groups. This dependency between the two groups is the key difference between the paired test and the two-sample test.

A more suitable approach for paired data is to test the differences between the pairs, e.g., the differences between the number of Moodle events in our scenario. We supply the t.test() function with the argument paired = TRUE so that it will take the measurement pairs into account. In this test, the null hypothesis is that the mean difference between the student pairs is zero (*H*<sup>0</sup> : *μd* = 0).

```
ttest_paired <- t.test(
 Frequency.Total ~ AchievingGroup, data = all_combined, paired = TRUE
)
ttest_paired
    Paired t-test
data: Frequency.Total by AchievingGroup
t = 4.3733, df = 64, p-value = 4.599e-05
alternative hypothesis: true mean difference is not equal to 0
95 percent confidence interval:
 178.3014 478.1910
sample estimates:
mean difference
       328.2462
```
The result is once again statistically significant, and we reject the null hypothesis. Because the mean difference between the pairs is positive, this means average number of Moodle events is higher for the top 50% achievers of the pairs.

Paired two-sample t-test is also supported by report().

```
report(ttest_paired)
```
Warning: Unable to retrieve data from htest object. Returning an approximate effect size using t\_to\_d().

Effect sizes were labelled following Cohen's (1988) recommendations.

The Paired t-test testing the difference of Frequency.Total by AchievingGroup (mean difference = 328.25) suggests that the effect is positive, statistically significant, and medium (difference = 328.25, 95% CI [178.30, 478.19], t(64) = 4.37, p < .001; Cohen's d = 0.55, 95% CI [0.28, 0.81])

We can once again safely ignore the produced warning message.

#### *3.2 Chi-Squared Test*

The chi-squared test is used for the analysis of contingency tables; it tests whether two categorical variables are independent or not [15]. A typical use case for this test is to investigate differences between groups such as student attendance by location or gender. The basic idea of the test is to compare the observed contingency table to a table under the null hypothesis where the variables are independent. The chisquared test is based on the cell-specific differences between these two tables. As a general rule, the test assumes that the expected value is at least 5 in at least 80% of the cells, and that no expected values are below 1. If these assumptions are violated, the results of the test may not be reliable. In such cases, Fisher's exact test [16] can be used instead via the function fisher.test(), but it may be computationally slow for large contingency tables. Both the chi-squared test and Fisher's exact test assume that the data is a random sample from the population.

We will use the combined Moodle data to investigate whether the achievement level and the activity level of the students are independent. First, we must create the contingency table from the individual-level data. We use the table() function for this purpose.

```
tab <- table(all_combined$ActivityGroup, all_combined$AchievingGroup)
tab
```

The table shows the observed frequencies in each cell, i.e., for each combination of activity and achievement. Next, we apply the chi-squared test using the chisq.test() function.

```
Xsq_test <- chisq.test(tab)
Xsq_test
```
Pearson's Chi-squared test

data: tab X-squared = 9.2135, df = 2, p-value = 0.009984

Printing the test result object shows the test statistic (X-squared), the associated degrees of freedom (df) and the p-value (p-value). The p-value is very small, meaning that we reject the null hypothesis. In other words, the achievement and activity levels of the students are not independent. This is not a surprising result, as

more active students are more likely to engage with the course content and perform better in terms of the learning outcomes. We can confirm that the assumptions of the test related to the expected values of the cells were not violated by using the test result object, which contains the expected values of the cells in the element expected.

```
all(Xsq_test$expected >= 1)
[1] TRUE
mean(Xsq_test$expected >= 5) >= 0.80
[1] TRUE
```
All expected values were greater than one, and over 80% of the expected values were greater than 5. This means that the assumptions are satisfied for our data and thus the results are reliable. Here, we used the function all() which takes a logical vector as input and returns TRUE if all elements of the vector were TRUE. Otherwise, the function returns FALSE. Unfortunately, the report() function does not support the chi-squared test.

#### *3.3 Analysis of Variance*

Analysis of variance (ANOVA) [17] can be viewed as the generalization of Student's t-test, where instead of one or two groups, the means of a variable are compared across multiple groups simultaneously. The name of the method comes from its test statistic, which is based on a decomposition of the total variance of the variable into variance within the groups and between the groups. ANOVA makes several assumptions: the observations are independent, the residuals of the underlying linear model follow a normal distribution, and that the variance of the variable is the same across groups (homoscedasticity). If these assumptions are violated, the results of the test may not be reliable. One alternative in such instances is to use the non-parametric Kruskal-Wallis test [18] instead, which is available in R via the function kruskal.test(). This test uses the ranks of the observations, and the null hypothesis is that the medians are the same for each group.

We use our combined Moodle data to demonstrate ANOVA. Instead of comparing the total number of Moodle events between top and bottom 50% of achievers, this time we will compare the final grade of the students across three activity groups: low activity, moderate activity, and high activity, described by the variable ActivityGroup. Thus the null and alternative hypotheses are in this case:

• *H*0: The expected values of the final grade are the same across the three activity groups (*μ*<sup>1</sup> = *μ*<sup>2</sup> = *μ*3),

• *H*1: At least one activity group has a different expected final grade (*μi* /= *μj* for at least one pair *i* /= *j* ).

To carry out the analysis, we apply the aov function, which uses the same formula syntax to define the response variable and the groups as the t.test() function does. Next, we apply the summary() function to the aov() function return object fit, as the default output of aov() is not very informative.

```
fit <- aov(Final_grade ~ ActivityGroup, data = all_combined)
summary(fit)
              Df Sum Sq Mean Sq F value Pr(>F)
ActivityGroup 2 175.7 87.87 25.11 6.47e-10 ***
Residuals 127 444.4 3.50
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
```
The summary contains the following columns: Df describes the degrees of freedom of the *F*-distribution associated with the test, Sum Sq reports the sum of squares related to the groups and the residuals, Mean Sq reports the corresponding mean sum of squares, F value is the value of the test statistic, and finally Pr(>F) is the p-value of the test. For this example, the p-value is very small, which means that the null hypothesis is rejected, and there are statistically significant differences in the final grade between the groups according to the test. In the following sections we will learn how to test for the assumptions related to normality and homoscedasticity. The report() function can be used for the output of aov() as well.

```
report(fit)
```

```
The ANOVA (formula: Final_grade ~ ActivityGroup) suggests that:
  - The main effect of ActivityGroup is statistically significant and large (F(2,
127) = 25.11, p < .001; Eta2 = 0.28, 95% CI [0.17, 1.00])
Effect sizes were labelled following Field's (2013) recommendations.
```
This output also reports the degrees of freedom, the test statistic value and the p-value but in a more easily readable format.

Note that ANOVA simply measures if there are differences between the groups but does not provide information on how these differences emerge. For example, there could be a single group that is different from all the rest, or two subgroups where the means are similar within each group, but different between the subgroups. Visualizations can be a helpful tool for gaining more insight into the differences, and post-hoc pairwise tests can be carried out to compare the pairs of groups.

#### *3.4 Levene's Test*

Levene's test is used to investigate whether the variance of a variable is the same across two or more groups [19]. Compared to alternatives such as Bartlett's test [20], Levene's test is less sensitive to non-normal observations. The test is not available in base R, but it can be found in the car package as the function leveneTest(). The function uses the same formula syntax as t.test() and aov(). We will investigate the homogeneity of the variance of the final grade between the activity groups.

```
library("car")
leveneTest(Final_grade ~ ActivityGroup, data = all_combined)
Levene's Test for Homogeneity of Variance (center = median)
       Df F value Pr(>F)
group 2 5.7204 0.004181 **
      127
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
```
The output of leveneTest() is analogous to the output of the ANOVA summary, and it contains the degrees of freedom (Df), the value of the test statistic (F value) and the p-value of the test (Pr(>F)). The p-value is very small, so we reject the null hypothesis meaning that the variance of the final grade is not the same across the groups according to the test. This means that the assumption of homoscedasticity is violated for the analysis of variance of the final grade, and thus the results may not be reliable. The report() function is not supported for leveneTest().

#### *3.5 Shapiro-Wilk Test*

Shapiro-Wilk test tests the null hypothesis that a data sample originated from a normal distribution [21]. The test is available in base R as the function shapiro.test(). Unfortunately, this function does not support the formula syntax unlike the other test functions we have used thus far. The function only accepts a single numeric vector as its argument. Therefore, to test the normality of multiple groups simultaneously, the data must first be split into the groups to be tested. We apply the test to the final grade in each achievement group. With the help of the broom package [22], we wrap the test results into a tidy format.

```
library("broom")
all_combined |>
  # Performs the computations in each activity group
  group_by(ActivityGroup) |>
```

```
# Apply a function in each group
 group_modify(~{
   # Apply the Shapiro test in each group and create tidy output
   shapiro.test(.$Final_grade) |>
     tidy()
 }) |>
 # Selection of variables to keep in the output
 select(ActivityGroup, statistic, p.value)
# A tibble: 3 x 3
# Groups: ActivityGroup [3]
 ActivityGroup statistic p.value
 <fct> <dbl> <dbl>
1 High activity 0.918 0.00448
2 Low activity 0.909 0.00215
3 Moderate activity 0.862 0.000103
```
This is also a great example of the tidyverse paradigm. First, we group the data by ActivityGroup using group\_by(), and then apply a function in each group using group\_modify(). We apply the shapiro.test() function to the Final\_grade variable, and then we convert the test results into a tidy tibble using the tidy() function from the broom package. We also use the special dot notation . to select the final grade variable from the data in each group. Finally, we select the grouping variable (ActivityGroup), the test statistic (statistic) and the p-value (p.value) of each test using select() and print the results. The resulting object is a tibble with three columns: ActivityGroup, statistic and p.value, the last two of which give the test statistic and p-value of the test for the activity group of the first column.

We can see that according to the test, the Final\_grade variable is not normally distributed in any of the activity groups, as the p-values are very small. As a consequence, the results of the analysis of variance carried out earlier may not be reliable.

#### **4 Correlation**

In Sect. 2.3, we briefly described covariance and correlation, and showcased the base R functions to compute them. However, there are several powerful and user-friendly packages for the analysis and reporting of correlations, such as the correlation [23] package which we will demonstrate in this section.

```
library("correlation")
```
For example, the package can easily compute all pairwise correlations between the numeric variables of the data with the function correlation(). The argument select can be used to compute the correlations only for a subset of the variables.

```
corrs <- correlation(
  all_combined,
  select = c("Frequency.Total", "Grade.Theory", "Final_grade")
)
corrs
# Correlation Matrix (pearson-method)
Parameter1 | Parameter2 | r | 95% CI | t(128) | p
Frequency.Total | Grade.Theory | 0.31 | [0.15, 0.46] | 3.75 | < .001***
Frequency.Total | Final_grade | 0.51 | [0.37, 0.63] | 6.76 | < .001***
Grade.Theory | Final_grade | 0.45 | [0.30, 0.58] | 5.69 | < .001***
p-value adjustment method: Holm (1979)
Observations: 130
```
The columns Parameter1 and Parameter2 describe the variables that the correlation was computed for, r is the value of the sample correlation, and the remaining columns report the 95% confidence interval, the value of the test statistic, and the p-value of the test (a t-test for correlations) along with the statistical significance. By default, Pearson's correlation coefficient is calculated, but the package also supports many alternative correlation measures. The correlation coefficient to be computed can be selected with the argument method that has the value "pearson" by default. Selecting for example method = "spearman" would compute the Spearman correlation coefficient instead. We can also obtain a correlation matrix by using summary()

```
summary(corrs)
# Correlation Matrix (pearson-method)
Parameter | Final_grade | Grade.Theory
Frequency.Total | 0.51*** | 0.31***
Grade.Theory | 0.45*** |
```
p-value adjustment method: Holm (1979)

By default, redundant correlations are omitted, but they can be obtained by setting redundant = TRUE in the call to summary(). A plot of the correlation matrix can be produced with the help of the package see [24].

```
library("see")
corrs |>
  # Also include redundant correlations
  summary(redundant = TRUE) |>
  plot()
```
The plot shows the strength of the correlations where darker colors imply stronger correlations. Visualizations will be covered at greater length in Chapter 6 [4].

#### **5 Linear Regression**

Linear regression is a statistical tool where one continuous variable is explained by the values of other variables. The variable of interest is said to be a dependent, while the other variables are called predictors. Predictors may also be called explanatory variables, independent variables or covariates depending on the context, applied field, and perspective.

Consider a very simple case, where we only have one predictor, which happens to be a continuous variable. In this case, fitting a linear regression model is merely the same as fitting a straight line to a scatterplot. It is assumed that deviations from this line are simply a result of random variation.

Now, let's go through the formal definition of a linear regression model. Let *Y* be a dependent variable with measurements *y*<sup>1</sup> *...,yn*, and let *X*1*, X*2*,...,Xk* be predictor variables with measurements *x*1*i,...,xki* for all *i* = 1*,...,n* where *i* refers to an individual measurement. Then, the regression equation is

$$\mathbf{x}\_{l}\mathbf{y}\_{l} = \beta\_{0} + \beta\_{1}\mathbf{x}\_{ll} + \beta\_{2}\mathbf{x}\_{2l} + \dots + \beta\_{k}\mathbf{x}\_{kl} + \varepsilon\_{l}, \quad \varepsilon\_{l} \sim N(0, \sigma^{2}), \quad i = 1, \dots, n$$

where we have the regression coefficients *β*0*, β*1*,...βk* and the error variance *σ*2. The first parameter *β*<sup>0</sup> is called the intercept that models the conditional expectation of *Y* when all the predictors have the value 0. From the regression equation, several assumptions become apparent. First, as the name of the model suggests, a linear relationship is assumed between the response and the predictors. Next, the variance *σ*<sup>2</sup> of the errors *εi* is constant, and does not depend on the values of the predictors (homoscedasticity). The errors are also assumed independent. The predictor variables are assumed fixed and their values perfectly measured without error.

Let's fit a linear regression model that predicts the final grade with the number of Moodle events of different types. To simplify the exposition, we will only use three types of Moodle events as predictors. We use the lm() function which has the same formula interface that we are already familiar with. First, we must define the dependent variable on the left-hand side of the formula, followed by the predictors on the right-hand side separated by a + sign. We must also supply the data argument, which tells the function where the actual values of the variables can be accessed.

```
fit <- lm(
  Final_grade ~ Frequency.Applications + Frequency.Assignment +
    Frequency.La_types,
  data = all_combined
)
summary(fit)
Call:
lm(formula = Final_grade ~ Frequency.Applications +
   Frequency.Assignment + Frequency.La_types, data = all_combined)
Residuals:
   Min 1Q Median 3Q Max
Coefficients:
                      Estimate Std. Error t value Pr(>|t|)
(Intercept) 5.800211 0.405963 14.288 < 2e-16 ***
Frequency.Applications 0.076516 0.022294 3.432 0.000811 ***
Frequency.Assignment -0.005049 0.005734 -0.881 0.380225
Frequency.La_types 0.088252 0.027314 3.231 0.001574 **
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1''1
Residual standard error: 1.914 on 126 degrees of freedom
```

```
Multiple R-squared: 0.2559, Adjusted R-squared: 0.2382
F-statistic: 14.45 on 3 and 126 DF, p-value: 3.798e-08
```
The summary() function provides a compact overview of the model fit for lm objects. First, a summary of the residuals (i.e., the differences between the observed and predicted values) is provided under Residuals. Next, a summary of the regression coefficients *β* is provided under Coefficients, including their estimates (Estimate), standard errors (Std. Error), test statistics for t-tests that test whether the coefficients are significantly different from zero (t value), p-values of the tests (Pr(>|t|)) and statistical significance (indicated by the asterisks). For instance, we see that the number of group work events is statistically significant. The notation used for the significance levels of the tests is described following Signif. codes. Estimate of the square root of the error variance *σ*<sup>2</sup> is reported following Residual standard error. The two R-squared values are estimates of the proportion of variance of the data that is explained by the model. Finally, the F-statistic reports the results of ANOVA when applied with the same model formula that was used for the linear regression model.

The report() function provides a more comprehensive summary of the model fit and the regression coefficients:

#### report(fit)

We fitted a linear model (estimated using OLS) to predict Final\_grade with Frequency.Applications, Frequency.Assignment and Frequency.La\_types (formula: Final\_grade ~ Frequency.Applications + Frequency.Assignment + Frequency.La\_types). The model explains a statistically significant and moderate proportion of variance (R2 = 0.26, F(3, 126) = 14.45, p < .001, adj. R2 = 0.24). The model's intercept, corresponding to Frequency.Applications = 0, Frequency.Assignment = 0 and Frequency.La\_types = 0, is at 5.80 (95% CI [5.00, 6.60], t(126) = 14.29, p < .001). Within this model:

Standardized parameters were obtained by fitting the model on a standardized version of the dataset. 95% Confidence Intervals (CIs) and p-values were computed using a Wald t-distribution approximation.

Again, the report output has condensed the information about the model fit into a format that can be read in a straightforward manner.

The assumption of normality of the residuals can be assessed with a quantilequantile plot, or q-q plot for short. The residuals of the model fit can be accessed with the function resid(). The function qqnorm() draws the quantiles of the residuals against the quantiles of the normal distribution. The function qqline() adds a straight line through the plot that passes through the second and third quantiles, by default. Ideally, the residuals should fall on this line, and large deviations indicate that the normality assumption may not hold.

```
# Draw the quantiles of the residuals and the theoretical quantiles
qqnorm(resid(fit))
# Add a line through the theoretical quantiles
qqline(resid(fit))
```
**Normal Q−Q Plot**

```
Theoretical Quantiles
```
The vast majority of residuals fall nicely onto the line for our model. Besides the q-q plot, we can obtain more model diagnostics with the help of the performance [25] package. This package provides a wide array of tools to assess how well models fit to the data. The general-purpose function check\_model( provides a visual overview of the model fit using several metrics.

```
library("performance")
check_model(fit, theme = see::theme_lucid(base_size = 10))
```
The functions performs various tests related to the assumptions of the linear regression model. For example, the bottom right panel contains the same q-q plot that we previously constructed using the qqnorm() and qqline() functions. We refer the reader to the documentation of the performance package for more information on the remaining tests.

Low (< 5)

#### **6 Logistic Regression**

Logistic regression is a similar tool to linear regression, but with a binary outcome instead of a continuous one. Instead of modeling the outcome variable directly, a linear model is constructed for the logarithmic odds of the probability of "success" for the binary outcome, e.g., obtaining a passing grade. There is also no explicit error term *ε* in the model, as the uncertainty in the outcome is already captured by the success probability. Formally, the model is

$$\text{logit}\left(P(\mathbf{y}\_{l}=\mathbf{l})\right) = \beta\_{0} + \beta\_{1}\mathbf{x}\_{ll} + \beta\_{2}\mathbf{x}\_{2l} + \dots + \beta\_{k}\mathbf{x}\_{kl}, \quad l = 1, \dots, n, \quad l$$

where the logit-function is defined as logit*(x)* = log*(x/(*1 − *x))*. Here, the logitfunction serves as the so-called *link function* that connects the expected value of the response to the predictors.

We fit a logistic regression model where the outcome variable is the level of achievement (AchievingGroup) and the predictors are the Moodle event counts of each type. The logistic regression model is a *generalized linear model*: a class of models that extend the linear regression model and that can be fitted in R with the function glm(). The syntax of glm() is analogous to lm(), but we must also specify the distribution of the outcome and the link function via the family argument. We use the function binomial() and supply the argument link = "logit" to define that the model should be a logistic regression model (in this case the link argument is optional, as "logit" is the default value). Because AchievingGroup is a factor, we must first convert it into a binary response that attains values 1 and 0 (or TRUE and FALSE). We can do this within the formula via the I() function, so that we do not have to modify our data. When used in a formula, this function will first compute its argument expression when evaluated, so that the expression is not mistaken for a variable name in the data (that does not exist). We select the high achievers as the "success" category for the outcome.

```
fit_logistic <- glm(
  # Use the I() function to construct a binary response in the formula
  I(AchievingGroup == "High achiever") ~ Frequency.Applications +
    Frequency.Assignment + Frequency.La_types,
  data = all_combined,
  # Our response is binary, so we use the binomial family with logit link
  family = binomial(link = "logit")
)
summary(fit_logistic)
Call:
glm(formula = I(AchievingGroup == "High achiever") ~ Frequency.Applications +
   Frequency.Assignment + Frequency.La_types, family = binomial(link = "logit"),
   data = all_combined)
```
Coefficients:

```
Estimate Std. Error z value Pr(>|z|)
(Intercept) -0.66272 0.62914 -1.053 0.29217
Frequency.Applications 0.30443 0.07778 3.914 9.07e-05 ***
Frequency.Assignment -0.04477 0.01402 -3.193 0.00141 **
Frequency.La_types 0.12245 0.04710 2.600 0.00933 **
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1''1
(Dispersion parameter for binomial family taken to be 1)
   Null deviance: 180.22 on 129 degrees of freedom
Residual deviance: 120.21 on 126 degrees of freedom
AIC: 128.21
Number of Fisher Scoring iterations: 6
```
The summary of a glm() function output is very similar to the output of a lm() summary. First, the Call is reported, which simply restates how the model was fitted. Next, Coefficients reports the estimates of the regression coefficients *β*, and their standard errors and statistical significance. Lastly, two deviance measures and their degrees of freedom are reported. The null deviance is twice the difference between the log-likelihood of the saturated model and the null model, and residual deviance is twice the difference between the saturated model and the model that was fitted. In simpler terms, the saturated model is a perfect model in a sense that there is a parameter for each observation. Conversely, the null model only has the intercept term. The deviance serves as a generalization of the residual sum of squares of the linear regression model, and it can be used to assess the quality of the model fit [26].

The report() function is applicable to models fitted with glm().

#### report(fit\_logistic)

We fitted a logistic model (estimated using ML) to predict AchievingGroup with Frequency.Applications, Frequency.Assignment and Frequency.La\_types (formula: I(AchievingGroup == "High achiever") ~ Frequency.Applications + Frequency.Assignment + Frequency.La\_types). The model's explanatory power is substantial (Tjur's R2 = 0.38). The model's intercept, corresponding to Frequency.Applications = 0, Frequency.Assignment = 0 and Frequency.La\_types = 0, is at -0.66 (95% CI [-1.94, 0.56], p = 0.292). Within this model:

Standardized parameters were obtained by fitting the model on a standardized version of the dataset. 95% Confidence Intervals (CIs) and p-values were computed using a Wald z-distribution approximation.

The output describes succinctly the model that was fitted, and the effects of the predictors on the response. The performance package is also applicable to models fitted with the glm() function.

#### check\_model(fit\_logistic)

**Table 1** Summary of the statistical tests, their null hypotheses, and the number of groups they compare simultaneously

We note that a different set of diagnostic checks is carried out for the logistic regression model compared to the linear regression model. For example, there is no assumption of homoscedasticity of variance as there is no explicit error term *ε* in the model. Again, we refer the reader to the documentation of the performance package for more details on these checks.

#### **7 Conclusion**

Basic statistics are an essential component of learning analytics. Learning analytics involves the collection, analysis, and interpretation of data related to the learning process, and statistical methods are used to identify patterns and trends in this data and to draw conclusions. Basic descriptive statistics such as measures of central tendency, variability and correlation are crucial for analyzing, interpreting, and visualizing data. Understanding these concepts is important for anyone involved in conducting research with quantitative data in the field of learning analytics. Moreover, mastery of basic statistics can facilitate the comprehension of more advanced statistical methods that are commonly used in learning analytics, such as logistic regression and cluster analysis. Table 1 contains a summary of the statistical tests that were introduced in this chapter.

We emphasize that when using any statistical test or a statistical model, it is important to keep the various assumptions related to the chosen method in mind, and to assess them beforehand whenever possible. If the assumptions are violated, the results of the method may not be reliable, and thus suitable alternatives should be considered.

#### **8 Further Reading**

This chapter scratched the surface of the full features of packages such as correlation, report and performance that can streamline the statistical

analysis and reporting process. We refer the reader to the documentation of these packages to gain a more thorough understanding of their features. These packages are part of a package collection called *easystats* [27] (https://github.com/ easystats/easystats). There are several other packages in this collection that were not discussed in this chapter that can be useful for R users working with learning analytics. The book "Learning Statistics with R" by Danielle Navarro is freely available online and provides a comprehensive introduction to statistics using R (https://learningstatisticswithr.com/). For a general introductory text to statistical methods and inference, see e.g., [2].

#### **References**

**Open Access** This chapter is licensed under the terms of the Creative Commons Attribution 4.0 International License (http://creativecommons.org/licenses/by/4.0/), which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons license and indicate if changes were made.

The images or other third party material in this chapter are included in the chapter's Creative Commons license, unless indicated otherwise in a credit line to the material. If material is not included in the chapter's Creative Commons license and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder.

# **Visualizing and Reporting Educational Data with R**

**Sonsoles López-Pernas, Kamila Misiejuk, Santtu Tikka, Juho Kopra, Merja Heinäniemi, and Mohammed Saqr**

### **1 Introduction**

Data visualization can be defined as "the representation and presentation of data that exploits our visual perception abilities in order to amplify cognition" [1]. It has the power to transform complex information into stories that inform and inspire action. Data visualization is an effective tool for learning analytics, as it helps to present learners' data in a way that is easily understandable and intuitive for students, teachers, researchers, and other stakeholders. Through the use of graphs, charts, and other visual aids, it is possible to quickly identify patterns, trends, and relationships within data that may not be immediately apparent through purely numerical data analysis methods.

Visualization in learning analytics has two distinct applications. On the one hand, the use of visual dashboards has become the main vehicle for putting learning analytics into practice. Presenting data in visually appealing and intuitive ways can help promote data literacy among students and other stakeholders, encouraging

K. Misiejuk

S. Tikka

Department of Mathematics and Statistics, University of Jyväskylä, Jyväskylä, Finland

#### M. Heinäniemi

Institute of Biomedicine, University of Eastern Finland, Kuopio, Finland

M. Saqr School of Computing, University of Eastern Finland, Joensuu, Finland

S. López-Pernas (-) · J. Kopra School of Computing, University of Eastern Finland, Joensuu, Finland e-mail<EMAIL_ADDRESS>
Centre for the Science of Learning & Technology (SLATE), University of Bergen, Bergen, Norway

greater engagement with data and fostering a culture of continuous improvement. On the other hand, learning analytics scientific production heavily relies on data visualization to present research findings in a clear and accessible manner, making it easier for readers from different scholarly backgrounds to understand and act upon research insights. Regardless of the context, the power of visualization in learning analytics lies in its ability to take complex data and turn it into meaningful insights that support better decision-making and drive improvement.

In this chapter, the reader will be guided through the process of generating meaningful and aesthetically pleasing visualizations of different types of datasets using well-known R packages. Relevant plots and plot types will be demonstrated with an explanation of their usage and usage cases. Furthermore, learning-related examples will be discussed in detail. For instance, readers will learn how to visualize learners' logs extracted from learning management systems (LMSs) to show how trace data can be used to track students' learning activities. Other examples of common research scenarios in which learners' data are visualized will be illustrated throughout the chapter. In addition to creating compelling plots, readers will also be able to generate professional-looking tables with summary statistics to report descriptive statistics.

#### **2 Visualization in Learning Analytics**

Developing visualizations is a challenging task of balancing the cognitive load of users while not compromising on conveying specific insights from data [2]. Visualizations for practice in learning analytics are mostly developed for two main stakeholders: learners and instructors. Depending on the target group, a visualization or a dashboard (i.e., a collection of visualizations depicting multiple indicators) have different goals.

Learner-facing visualizations are meant to make learners aware of their own learning and to provide them with actionable feedback on their learning. Visualizations display learners' performance on a specific metric and compare it with a reference frame: other peers, desirable learning achievement, or their own progress over time [3]. Sense-making questions triggering reflection can be added to a visualization [4, 5], or some elements of the visualizations can be highlighted and described in words using layered storytelling [6, 7]. Another option is to gamify a dashboard, for example, by using badges [8]. To provide feedback to learners, visualizations can be augmented with links to recommended resources [9], information about specific topics to review to close the achievement gap [6], or explanations of the meaning of visualizations and their implications for the learner [10]. Current learner-facing dashboards mostly show resource use and assessment data [11], compare learners to their peers [12], display descriptive analytics rather than predictive or prescriptive analytics [10], and use self-regulated learning theory as their framework [12, 13]. Some reviews found a positive effect on student outcomes [10], while others reported mixed results [11, 14]. Showing visualizations to learners can change their behavior. For example, social network analysis visualizations have resulted in fewer cross-group commenting [15], while a visualization comparing individual submission patterns with the top 25% of students in a class led to earlier homework submissions [16].

In comparison, the goal of instructor-facing visualizations is to support teachers and their decision-making process by tracking student progress. Two main types can be distinguished. Mirroring or descriptive visualizations provide insights about the learners on an aggregated or an individual level using either descriptive or comparative data. Advising or prescriptive visualizations show not only information about the learners but also alert the instructor to undertake a pedagogical action [17, 18]. Current instructor-facing visualizations mostly display course-wide information about the learners or track group work [14]. These visualizations can support teachers in facilitating student collaboration [19], planning and collecting student feedback on learning activities [20], or obtaining insights into student interactions within an online environment, such as simulations, virtual labs or online games [21, 22]. However, interpreting dashboard information is a challenging task for instructors. Although some teachers use dashboards as complementary sources of information, others act based only on the dashboard information without further investigation [23].

A common point of criticism of learning analytics dashboards is that most of them are not grounded in learning theories [13, 14]. Data-driven evaluations of dashboards focused on dashboard acceptance, usefulness, or usability are more prevalent than pedagogically-focused evaluations [24]. Some approaches were developed to mitigate these issues. The model of user-centered learning analytics systems (MULAS) presents a set of recommendations on four interconnected dimensions: theory, design, evaluation, and feedback, and can be used to guide dashboard development [14]. Another approach is an iterative five-stage Learning Awareness Tools—User eXperience (LATUX) workflow, including problem identification, low-fidelity prototyping, high-fidelity prototyping, pilot studies, and classroom use, that can be used to develop visual analytics [25]. Finally, open learner model research could be used as a source of insights while developing learning analytics visualizations, such as dashboards [9].

#### **3 Generating Plots with ggplot2**

In the previous section, we have seen how central visualization is to learning analytics. In the remainder of the chapter, we will learn how to create different types of visualizations that are relevant to different types of data related to teaching and learning. We will mostly rely on ggplot2, a popular data visualization package in R that was developed by Hadley Wickham [26]. It is based on the grammar of graphics [27], which is a systematic way of thinking about and constructing visualizations. The ggplot2 library provides a flexible and intuitive framework for creating a wide range of graphics, from basic scatter plots to complex visualizations with multiple layers. It is known for its ability to produce visually appealing and informative graphics with relatively few lines of code. It enables users to define aesthetics, such as color and size, and add layers, such as points and lines, to create customized and interactive plots. In addition, ggplot2 allows for easy customization of plot features, such as titles, axis labels, and legends.

Overall, ggplot2 is a powerful and versatile tool for data visualization in R, and is widely used by data scientists, statisticians, and researchers in a variety of fields. In this chapter, we will cover the fundamental concepts and techniques of ggplot2, including how to create basic plots, and customize their appearance. We will start by introducing the building blocks of a ggplot2 plot, including aesthetics, layers, and scales. Then, we will create a plot from scratch step by step, showing how to customize its appearance, including how to change theme, colors, and scales. We will then explore the different types of plots that can be created with ggplot2, such as scatter plots, bar charts, and histograms.

Throughout this section, we will use datasets of students' learning data to demonstrate how to create effective visualizations for learning analytics with ggplot2. Please, refer to Chapter 2 of this book [28] to learn more about the datasets used. By the end of this section, you will have a solid foundation in ggplot2 and be able to create basic, yet compelling visualizations to explore your data.

#### *3.1 The ggplot2 Grammar*

The ggplot2 library is based on Wilkinson's grammar of graphics [27]. The main idea is that every plot can be broken down into a set of components, each of which can be customized and combined in a flexible way. These components are:

Through the combination and customization of these components, we can create a wide variety of complex and informative visualizations in ggplot2. The idea behind the graphics grammar is to provide a consistent framework for constructing plots, allowing users to focus on the data and the message they want to convey, rather than on the technical details of the visualization. In the following section, we will create a plot from scratch step by step to become familiar with the most relevant components.

#### *3.2 Creating Your First Plot*

We will now create our first plot using ggplot2. Our example deals with a widely studied matter in learning analytics, which is the relationship between online activity and achievement. We will use a bar chart to represent the number of students that have low, moderate and high activity levels in each achievement group (high achievers vs. low achievers). In order to become familiar with the syntax of ggplot2, we will recreate the plot step by step, explaining each of the elements in the plot. Below is the final result we aim at accomplishing (Fig. 1):

**Activity level by achievement group**

**Fig. 1** First plot with ggplot2

#### **3.2.1 Installing ggplot2**

Our first step is installing the ggplot2 library. This is usually the first step in any R script that makes use of external libraries.

```
install.packages("ggplot2")
```
To import ggplot2 we just need to use the library command and specify the ggplot2 library:

library(ggplot2)

#### **3.2.2 Downloading the Data**

Next, we need to import the data that we are going to plot. For this chapter, we are using synthetic data from a blended course on learning analytics. For more details about this dataset, refer to Chap. 2 in this book. The data is in Excel format. We can use the library rio since it makes it easy to read data in several formats. We first install the library:

```
install.packages("rio")
```
And import it so we can use its functions:

library(rio)

Now we can download the data using the import function from rio and assign it to a variable named df (short for dataframe).

```