This document provides a practical introduction to performing statistical analyses in R. Throughout the course, we discuss a range of statistical methods and concepts. Here, we will translate these concepts into practice by looking at the R code used to perform the corresponding analyses.

The analyses will be demonstrated using an example dataset: apartments. For each analysis, we will walk through the relevant R code step by step, explain what the different commands do, and look at the resulting output. We will also discuss how to interpret the results and how the output relates to the statistical concepts covered in the course.

The goal is not to memorise the code, but to understand the general workflow and learn how R can be used as a tool for statistical analysis. The code examples can serve as a reference when you work with your own data. In many cases, the same code can be adapted to a different dataset by changing the relevant variable names or other arguments.

As you work through the examples, it is useful to keep the statistical reasoning in mind: What question are we trying to answer? Which analysis is appropriate? How do we perform it in R? And how do we interpret the results? These questions are at least as important as the R code itself.

The examples in this document are intended as a practical guide that you can return to when performing similar analyses on your own data.

1 Getting started in R

Throughout this course, the APARTMENT dataset will be used to show how to perform certain statistical analyses in R. The APARTMENT dataset contains 117 observations on apartments sold in Antwerp. The variables recorded in the dataset are:

PRICE: selling price in euro
PRICECLASS: 1 (less than 85000 euro), 2 (more than 85000 and less than 125000 euro), 3 (more than 125 000 euro)
AGE: age of apartment at time of sale
FEATS: number of features out of a list of 11 luxury articles (e.g. dishwasher, micro-wave,…)
AREA: residential area (1=yes, 0=no)
TAX: cadastrale income
SQM: living area in square meters

1.1 Reading in a txt file

Since the APARTMENT dataset is recorded in a text file (extension .txt), we can use the read.table() command in R to load the dataset into the R environment:

apartment <- read.table("apartment.txt", header=T)

Here, apartment is the name we give to our dataset. The argument header = TRUE tells R that the first row of the file contains the variable names.

Once the data have been read into R, the dataset is stored as an object called apartment. We can then use this object in our analyses.

1.2 Viewing the data

There are several simple ways to inspect a dataset.

The View() function opens the dataset in a spreadsheet-like window:

View(apartment)

This is useful when you want to look through the data and get a general idea of its structure.

You can also use head() to look at the first rows of the dataset:

head(apartment)
#>   price priceclass age feats area tax    sqm feats2
#> 1 90000          2  16     2    1 731 63.000      1
#> 2 86000          2  11     2    1 653 62.325      1
#> 3 92500          2  45     2    1 553 47.250      1
#> 4 89900          2   9     2    1 566 65.880      1
#> 5 85000          1  41     1    1 600 53.550      1
#> 6 87600          2   8     1    1 650 52.020      1

By default, head() shows the first six rows. You can specify the number of rows you want to see:

head(apartment,10)
#>    price priceclass age feats area tax    sqm feats2
#> 1  90000          2  16     2    1 731 63.000      1
#> 2  86000          2  11     2    1 653 62.325      1
#> 3  92500          2  45     2    1 553 47.250      1
#> 4  89900          2   9     2    1 566 65.880      1
#> 5  85000          1  41     1    1 600 53.550      1
#> 6  87600          2   8     1    1 650 52.020      1
#> 7  89000          2  12     2    1 591 78.570      1
#> 8  87000          2  23     1    1 599 57.600      1
#> 9  72000          1   7     1    1 436 47.250      1
#> 10 81000          1  33     2    1 673 61.425      1

Similarly, tail() shows the last rows of the dataset:

tail(apartment)
#>      price priceclass age feats area  tax     sqm feats2
#> 112 215000          3   6     5    1 1193 119.880      2
#> 113 215000          3   3     6    1 1635 131.445      2
#> 114 215000          3   4     6    1 1487 128.160      2
#> 115 159900          3  24     5    1 1265 109.800      2
#> 116 129900          3  25     5    1 1232 123.435      2
#> 117 184400          3  40     6    0  915 101.250      2

These commands are particularly useful for quickly checking whether the data have been read into R correctly.

1.3 Extracting a variable

Variables (columns) can be extracted from a dataset using the $ operator.

For example, the apartment dataset contain the variable ‘price’, so we can extract it using:

apartment$price
#>   [1]  90000  86000  92500  89900  85000  87600  89000  87000  72000  81000
#>  [11]  79900  75000  69000  67000  61900  93900  82000  70000  54000 107000
#>  [21]  66000  58000  69900 103000 158000 144900 115500 111000 113900  99500
#>  [31]  99500  97500  97500 102000 102000  92200  70000  72000  72500  67000
#>  [41] 123900 112500 110000 105000  95500  93400  87500  88900  80500  75000
#>  [51]  75900  73000  71000  97500  78000  77000  62000  72500  60000 116000
#>  [61] 110900 112900 105000 104500 102000  97500  95000  94000  92000  94500
#>  [71]  87400  87200  87000  86900  76600  73900 145000 125000 118000 155300
#>  [81] 130000 125000 120000 108000 104900  85500 210000 105000 208000 199900
#>  [91] 190000 180000 169500 125000 135000 129500  96000  74900  73100  83500
#> [101]  75500  72900  77300 100000 156000 137500 127000 123500 117000 133000
#> [111] 205000 215000 215000 215000 159900 129900 184400

This returns the values of the variable price.

The $ operator can also be used when performing calculations or creating graphs. For example, we can calculate the mean age with:

mean(apartment$price)
#> [1] 106273.5

In general, the syntax for extracting a variable is:

Syntax:
data$variable

Where:

data: Is the name of your dataset

variable: Is the name of the variable you want to access

These basic commands are enough to start exploring a dataset in R. Once the data are loaded and we know how to access the variables, we can move on to describing and visualising the data.

2 Descriptive statistics

The first step in any statistical analysis is to explore and describe the data. Before applying statistical tests or fitting statistical models, it is important to develop a good understanding of the data you are working with.

This initial exploration helps us to understand the structure and characteristics of the dataset. We can, for example, look at the distributions of variables, identify patterns or differences between groups, and detect unusual observations or potential data problems. It also helps us to decide which statistical analyses may be appropriate for the data.

To explore our data, we use a combination of visualisations and summary statistics. Graphs allow us to see patterns and distributions in the data, while numerical summaries provide useful information about the central tendency, variability, and other characteristics of our variables.

In this section, we will use R to create several commonly used visualisations and calculate descriptive statistics. The specific methods we use will depend on the type of variable we are interested in and the question we want to answer.

2.1 Checking and changing the type of a variable

Before calculating summary statistics or creating graphs, it is important to check what type of variable we are working with. The type of variable determines which analyses and visualisations are appropriate.

We can use the class() function to check what type of variable we are working with:

class(apartment$price)
#> [1] "numeric"

In this example, R returns: [1] "numeric". This tells us that price is a numerical variable.

For a categorical variable, we generally want R to recognise the variable as a factor. We can check this in the same way:

class(apartment$area)
#> [1] "integer"

If the variable is not recognised as a factor by R, we can convert it using factor(), and check the results again:

apartment$area <- factor(apartment$area)
class(apartment$area)
#> [1] "factor"

It is important to make sure that variables have the correct type before starting the analysis. For example, a variable representing groups or categories should generally be treated as a factor rather than as a numerical variable.

When we convert a variable to a factor, R uses the values in the dataset as the levels of the factor. We can inspect these levels using levels():

levels(apartment$area)
#> [1] "0" "1"

Sometimes the levels are not very informative or are difficult to understand. We can change them to more meaningful labels using the same function. For example, the variable area in the dataset contains the levels 0 and 1, representing non-residential and residential areas. We can replace these numbers with more understandable labels:

levels(apartment$area) <- c("Non-residential","Residential")
levels(apartment$area)
#> [1] "Non-residential" "Residential"

The order of the labels corresponds to the order of the existing factor levels. Therefore, it is important to check the original levels before changing them.

In order to avoid making mistakes and changing variables permanently in your dataset, it is best to create a new variable in your dataset that contains the levels that you want. For example:

apartment$priceclass2 <- factor (apartment$priceclass,
                            levels = c(1,2,3),
                            labels = c("low", "medium", "high"))
apartment$priceclass2
#>   [1] medium medium medium medium low    medium medium medium low    low   
#>  [11] low    low    low    low    low    medium low    low    low    medium
#>  [21] low    low    low    medium high   high   medium medium medium medium
#>  [31] medium medium medium medium medium medium low    low    low    low   
#>  [41] medium medium medium medium medium medium medium medium low    low   
#>  [51] low    low    low    medium low    low    low    low    low    medium
#>  [61] medium medium medium medium medium medium medium medium medium medium
#>  [71] medium medium medium medium low    low    high   medium medium high  
#>  [81] high   medium medium medium medium medium high   medium high   high  
#>  [91] high   high   high   medium high   high   medium low    low    low   
#> [101] low    low    low    medium high   high   high   medium medium high  
#> [111] high   high   high   high   high   high   high  
#> Levels: low medium high

This creates a new categorical variable called priceclass2 in the dataset apartment, which contains the levels low, medium, and high as defined by the categories 1,2, and 3 in the original variable priceclass.

2.2 Frequency tables

For categorical variables, a useful first step is to create a frequency table. The table() function shows how many observations belong to each category:

table(apartment$area)
#> 
#> Non-residential     Residential 
#>              39              78

To calculate the percentages, we can divide the frequencies by the total number of observations:

table(apartment$area) / length(apartment$area)
#> 
#> Non-residential     Residential 
#>       0.3333333       0.6666667

Alternatively, we can use prop.table():

prop.table(table(apartment$area))
#> 
#> Non-residential     Residential 
#>       0.3333333       0.6666667

The result gives the percentage of observations in each category.

We can also calculate the cumulative percentage using cumsum():

cumsum(prop.table(table(apartment$area)))
#> Non-residential     Residential 
#>       0.3333333       1.0000000

We can combine the frequency, percentage, and cumulative percentage into one table:

freq <- table(apartment$area)
percent <- prop.table(table(apartment$area)) 
cumulative_percent <- cumsum(percent)

freq_table <- data.frame(
  Frequency = as.vector(freq),
  Percentage = as.vector(percent),
  Cumulative_Percentage = cumulative_percent
)
freq_table
#>                 Frequency Percentage Cumulative_Percentage
#> Non-residential        39  0.3333333             0.3333333
#> Residential            78  0.6666667             1.0000000

The cumulative percentage is particularly useful when the categories have a meaningful order.

2.3 Barplot

A barplot is useful for visualising the frequencies of a categorical variable. In base R, we can create one directly from a frequency table:

barplot(table(apartment$area))

The height of each bar represents the number of observations in that category.

We can customise the appearance of the barplot by adding different arguments to the barplot() function:

barplot(
  table(apartment$area),
  main = "Barplot of the variable AREA",
  xlab = "AREA",
  ylab = "Frequency",
  col="darkred"
)

The different arguments control different aspects of the graph:

main adds a title to the graph

xlab specifies the label of the x-axis.

ylab specifies the label of the y-axis.

col specifies the colour of the bars.

For example, if we want to use a different colour:

barplot(
  table(apartment$area),
  main = "Barplot of the variable AREA",
  xlab = "AREA",
  ylab = "Frequency",
  col="lightblue"
)

You can use many different colour names in R, such as “lightblue”, “lightgreen”, “orange”, “pink”, “grey”, or “steelblue”. You find lists online of all the colors that are recognized by R.

You can also specify a different colour for each bar by providing a vector of colours:

barplot(
  table(apartment$area),
  main = "Barplot of the variable AREA",
  xlab = "AREA",
  ylab = "Frequency",
  col=c("lightblue","lightgreen")
)

The first colour is used for the first bar, the second colour for the second bar, and so on.

You do not need to specify all of these arguments. If you leave an argument out, R will use its default setting. As you become more familiar with R, you can explore additional options for controlling the appearance of your graphs. For now, the arguments introduced above are sufficient for creating clear and informative barplots.

2.4 Histogram

For a numerical variable, we are often interested in its distribution. A histogram shows how the observations are distributed across different values or intervals.

For example:

hist(apartment$price)

We can also add labels to the histogram, change the color, …:

hist(apartment$price,
     main="Distribution of Price",
     xlab="Price",
     ylab="Frequency",
     col="darkred")

A histogram can help us assess the shape of a distribution and identify features such as skewness, gaps, or unusually high or low values.

2.5 Boxplot

Another useful way to visualise a numerical variable is with a boxplot:

boxplot(apartment$price)

We can add labels and colours in the same way:

boxplot(apartment$price,
        main = "Distribution of Price",
        ylab = "Price",
        col="yellow")

A boxplot provides a compact summary of the distribution. It shows the median, the quartiles, and the overall spread of the data. Individual observations that are far from the rest of the data may also be displayed as potential outliers.

Boxplots can also be used to compare a numerical variable between groups. For example, if we want to compare price between residential and non-residential areas:

boxplot(apartment$price~apartment$area,
        main = "Distribution of Price by Area",
        ylab = "Price",
        xlab="Area",
        col="darkgreen")

Here, the ~ notation tells R that we want to examine price by area.

2.6 Summary statistics

Graphs provide a visual overview of our data, but we can also describe a numerical variable using summary statistics. These statistics give us information about the centre, spread, and range of the distribution.

2.6.1 Measures of central tendency

Measures of central tendency describe the centre of a distribution. The two most commonly used measures are the mean and the median.

The mean is calculated by adding all observations and dividing by the number of observations. In R, we use the mean() function:

mean(apartment$price)
#> [1] 106273.5

The mean can be sensitive to extreme values. If a variable contains very large or very small observations, the mean can be pulled towards these values.

The median is the middle value when the observations are ordered from smallest to largest. In R, we use the median() function:

median(apartment$price)
#> [1] 96000

The median is less affected by extreme values than the mean and can therefore be a useful measure of the centre when a distribution is skewed.

2.6.2 Measures of spread

Measures of spread tell us how much the observations differ from each other.

The variance measures how much the observations vary around the mean. It is based on the squared deviations from the mean. In R, we use var():

var(apartment$price)
#> [1] 1447322999

Because the deviations are squared, the variance is expressed in squared units of the original variable. This can make it difficult to interpret directly. The standard deviation is closely related to the variance, i.e. the square root of the variance, and is often easier to interpret because it is expressed in the same units as the original variable. The standard deviation gives an indication of the typical distance of observations from the mean. In R, we use sd():

sd(apartment$price)
#> [1] 38043.7

The minimum and maximum give us the smallest and largest observed values.

min(apartment$price)
#> [1] 54000
max(apartment$price)
#> [1] 215000

We can also obtain both values at the same time using range():

range(apartment$price)
#> [1]  54000 215000

The result gives the minimum followed by the maximum. The range describes the distance between the minimum and maximum values. It can be calculated by subtracting the minimum from the maximum:

max(apartment$price)-min(apartment$price)
#> [1] 161000

The range gives a simple indication of the total spread of the data, but it can be strongly affected by extreme observations.

Quartiles divide the ordered data into four parts.

The first quartile (Q1) is the value below which approximately 25% of the observations fall.
The second quartile (Q2) is the median, with approximately 50% of observations below it.
The third quartile (Q3) is the value below which approximately 75% of the observations fall.

In R, we can calculate the quartiles using quantile():

quantile(apartment$price)
#>     0%    25%    50%    75%   100% 
#>  54000  78000  96000 120000 215000

By default, R returns the minimum, Q1, median, Q3, and maximum. We can also request specific quantiles. For example:

quantile(apartment$price, probs = c(0.25, 0.50, 0.75))
#>    25%    50%    75% 
#>  78000  96000 120000

Here, probs specifies which quantiles we want to calculate. 0.25 corresponds to the 25th percentile (Q1), 0.50 to the median (Q2), and 0.75 to the 75th percentile (Q3).

The interquartile range (IQR) describes the spread of the middle 50% of the observations. It is calculated as IQR = Q3 − Q1. In R, we can calculate it directly using:

IQR(apartment$price)
#> [1] 42000

The IQR is less affected by extreme values than the range and is therefore a useful measure of spread, particularly for skewed distributions.

Instead of calculating each statistic separately, we can use the summary() function:

summary(apartment$price)
#>    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
#>   54000   78000   96000  106274  120000  215000

For a numerical variable, summary() provides the minimum, first quartile, median, mean, third quartile, and maximum. This is a convenient way to get a quick numerical overview of a variable.

When we have two numerical variables, we can investigate whether they are related using a correlation. The Pearson correlation coefficient, usually denoted by r, measures the strength and direction of a linear relationship between two numerical variables.

The correlation coefficient ranges from -1 to +1:

A value close to +1 indicates a strong positive linear relationship.
A value close to -1 indicates a strong negative linear relationship.
A value close to 0 indicates little or no linear relationship.

For example, we can calculate the correlation between two numerical variables using cor():

cor(apartment$price,apartment$tax)
#> [1] 0.8768248

We can also visualise the relationship between two numerical variables using a scatterplot:

plot(apartment$price,
     apartment$tax,
     col="darkred",
     pch=19,
     xlab="PRICE",
     ylab="TAX")

Here, pch = 19 specifies that the observations should be displayed as solid circles. There are many different plotting symbols available in base R. For example:

plot(apartment$price,apartment$tax,pch=1)
plot(apartment$price,apartment$tax,pch=23)
plot(apartment$price,apartment$tax,pch=17)

You can find an overview online of the different values for different shapes of points.

The scatterplot is an important complement to the correlation coefficient because it allows us to see the relationship between the variables. A correlation coefficient alone does not tell us everything about the relationship and can, for example, hide non-linear patterns or the influence of unusual observations.

3 Parameter estimation

In statistics, we are often interested in learning something about a population, but we usually do not have access to data from the entire population. Instead, we collect a sample and use this sample to learn about the population.

For example, suppose we are interested in the average value of a variable in a population. We could calculate the mean of a sample and use this sample mean as an estimate of the population mean.

An important question is: How much would the sample mean vary if we repeatedly took different samples from the same population? The Central Limit Theorem (CLT) helps us answer this question.

The Central Limit Theorem states, roughly, that when we repeatedly take random samples of a sufficiently large size from a population, the distribution of the sample means will be approximately normal, even if the original population itself is not normally distributed. There are therefore two different distributions to keep in mind:

The population distribution: the distribution of the individual observations.
The sampling distribution of the mean: the distribution we would obtain if we calculated the mean for many different samples from the population.

The CLT is about the second distribution. We can illustrate this with a simulation in R. The details of the simulation are not important to remember. The purpose is simply to visualise what happens when we repeatedly take samples and calculate their means.

3.1 Simulating samples

Let’s start with a population that is clearly not normally distributed. For example, we can generate observations from an exponential distribution:

set.seed(123)

population <- rexp(100000, rate = 1)

We can visualise this population:

hist(
  population,
  main = "Population distribution",
  xlab = "Value",
  col = "lightblue"
)

As you can see, the population distribution is strongly right-skewed rather than normally distributed.

Now, imagine that we repeatedly take samples from this population. For each sample, we calculate the mean. We can simulate this process using the following code:

sample_size <- 30
number_of_samples <- 10000

sample_means <- replicate(
  number_of_samples,
  mean(sample(population, sample_size))
)

Here, we have taken 10,000 samples, each containing 30 observations, and calculated the mean of each sample. The object sample_means therefore contains 10,000 sample means.

We can now look at the distribution of these means:

hist(
  sample_means,
  main = "Sampling distribution of the mean",
  xlab = "Sample mean",
  col = "lightblue"
)

Although the original population was strongly skewed, the distribution of the sample means is approximately normal.

This is the central idea behind the Central Limit Theorem.

What happens with smaller samples? The sample size is important. Let’s repeat the simulation using a much smaller sample size, for example 5 observations per sample:

sample_size <- 5

sample_means_small <- replicate(
  number_of_samples,
  mean(sample(population, sample_size)))

We can again look at the distribution of the sample means:

hist(
  sample_means_small,
  main = "Sampling distribution (n = 5)",
  xlab = "Sample mean",
  col = "lightblue"
)

With such a small sample size, the distribution of the sample means is less clearly normal and retains more of the skewness of the original population.

If we increase the sample size:

sample_size <- 100

sample_means_large <- replicate(
  number_of_samples,
  mean(sample(population, sample_size))
)

and plot the resulting means:

hist(
  sample_means_large,
  main = "Sampling distribution (n = 100)",
  xlab = "Sample mean",
  col = "lightblue"
)

The distribution of the sample means becomes even more clearly approximately normal.

4 Confidence interval for the mean

A sample mean gives us an estimate of the population mean. However, because we usually work with a sample rather than the entire population, our estimate is subject to sampling variability. If we were to take a different sample, we would generally obtain a somewhat different sample mean.

Rather than reporting only the sample mean, we can therefore also give an interval of plausible values for the population mean. This is called a confidence interval.

A commonly used confidence level is 95%. A 95% confidence interval gives us an interval constructed from the sample data that, under repeated sampling and the assumptions of the method, would contain the true population mean in approximately 95% of repeated samples.

Suppose we want to estimate the population mean of price based on the observations in our apartment dataset. R provides a convenient way to calculate a confidence interval for a mean using the t.test() function:

t.test(apartment$price)
#> 
#>  One Sample t-test
#> 
#> data:  apartment$price
#> t = 30.216, df = 116, p-value < 2.2e-16
#> alternative hypothesis: true mean is not equal to 0
#> 95 percent confidence interval:
#>   99307.36 113239.65
#> sample estimates:
#> mean of x 
#>  106273.5

The relevant part of the output for our purposes is:

95 percent confidence interval: 99307.36 113239.65

sample estimates: mean of x 106273.5

We estimate the population mean to be 106273.5, with a 95% confidence interval from 99307.36 to 113239.65. In other words, the interval gives us a range of plausible values for the population mean based on our sample.

It is important to be precise about what the 95% refers to. The confidence level describes the method used to construct the interval, not the probability that the particular interval we have calculated contains the population mean. In repeated sampling, approximately 95% of intervals constructed using this method would contain the true population mean.

A 95% confidence interval is common, but we can choose other confidence levels. For example, t.test() can calculate a 99% confidence interval:

t.test(apartment$price, conf.level = 0.99)
#> 
#>  One Sample t-test
#> 
#> data:  apartment$price
#> t = 30.216, df = 116, p-value < 2.2e-16
#> alternative hypothesis: true mean is not equal to 0
#> 99 percent confidence interval:
#>   97062.54 115484.47
#> sample estimates:
#> mean of x 
#>  106273.5

A higher confidence level produces a wider interval. For example, a 99% confidence interval will generally be wider than a 95% confidence interval because we want to be more confident that the interval captures the population mean.

Similarly, a 90% confidence interval will generally be narrower:

t.test(apartment$price, conf.level = 0.90)
#> 
#>  One Sample t-test
#> 
#> data:  apartment$price
#> t = 30.216, df = 116, p-value < 2.2e-16
#> alternative hypothesis: true mean is not equal to 0
#> 90 percent confidence interval:
#>  100441.7 112105.3
#> sample estimates:
#> mean of x 
#>  106273.5

The width of a confidence interval is influenced by the sample size. Larger samples generally produce more precise estimates and therefore narrower confidence intervals.

This connects directly to the Central Limit Theorem and sampling variability discussed above: as the sample size increases, the sample mean becomes less variable, which results in a smaller standard error and a more precise estimate of the population mean.

5 Hypothesis testing

In the previous section, we used a sample to estimate a population mean and quantified the uncertainty around this estimate using a confidence interval. Sometimes, however, our research question is not simply “What is the population mean?” but rather “Is the population mean different from a particular value?” This is where hypothesis testing comes in.

The one-sample t-test can be used to test whether the mean of a population differs from a specific value.

For example, suppose we want to investigate whether the mean price in our population differs from 100.000 euro.

We start by defining two hypotheses:

Null hypothesis ($H_0$): the population mean price is equal to 100.000.
Alternative hypothesis ($H_1$): the population mean price is different from 100.000.

We can perform this test in R using t.test():

t.test(apartment$price, mu=100000)
#> 
#>  One Sample t-test
#> 
#> data:  apartment$price
#> t = 1.7837, df = 116, p-value = 0.07709
#> alternative hypothesis: true mean is not equal to 1e+05
#> 95 percent confidence interval:
#>   99307.36 113239.65
#> sample estimates:
#> mean of x 
#>  106273.5

The argument mu = 40 tells R that we want to test the null hypothesis that the population mean is 100.000.

The output from t.test() contains several pieces of information. The p-value is particularly important for hypothesis testing. The p-value tells us how compatible our observed data are with the null hypothesis. More specifically, it is the probability of obtaining a result at least as extreme as the one observed, assuming that the null hypothesis is true.

A common significance level is alpha=0.05 We compare the p-value with this significance level:

If p < 0.05, we reject the null hypothesis.
If p >= 0.05, we do not reject the null hypothesis.

For the example above, the p-value is 0.07709. Since this is larger than 0.05, we do not reject the null hypothesis and conclude that the data do not provide evidence that the population mean price is different from 100000.

It is important to distinguish between “rejecting the null hypothesis” and “proving that the alternative hypothesis is true”. A hypothesis test provides evidence against the null hypothesis; it does not prove a hypothesis with certainty.

5.1 One-sided versus two-sided tests

The test above is a two-sided test because our alternative hypothesis is “not equal to”.

If we have a directional research question (meaning “larger than”, or “smaller than”), we can perform a one-sided test. For example, if we specifically want to test whether the population mean is greater than 100000:

t.test(apartment$price, mu=100000, alternative = "greater")
#> 
#>  One Sample t-test
#> 
#> data:  apartment$price
#> t = 1.7837, df = 116, p-value = 0.03854
#> alternative hypothesis: true mean is greater than 1e+05
#> 95 percent confidence interval:
#>  100441.7      Inf
#> sample estimates:
#> mean of x 
#>  106273.5

For a test of whether the mean is less than 100000:

t.test(apartment$price, mu=100000, alternative = "less")
#> 
#>  One Sample t-test
#> 
#> data:  apartment$price
#> t = 1.7837, df = 116, p-value = 0.9615
#> alternative hypothesis: true mean is less than 1e+05
#> 95 percent confidence interval:
#>      -Inf 112105.3
#> sample estimates:
#> mean of x 
#>  106273.5

The default in t.test() is: alternative = "two.sided"

Therefore, if you do not specify alternative, R performs a two-sided test.

A one-sided test should only be used when a directional hypothesis was specified before looking at the results. We should not choose a one-sided test simply because it produces a smaller p-value.

5.2 Statistical power and sample size

When planning a study, we face an important question: How large should our sample be?

If our sample is too small, we may have little ability to detect an effect that is genuinely present in the population. This ability to detect an effect is called statistical power. Statistical power is the probability of rejecting the null hypothesis when the alternative hypothesis is actually true.

Power tells us how likely our study is to detect an effect of a particular size if that effect really exists. A commonly used target for statistical power is 80%.

Power depends on several factors, including:

The sample size
The size of the effect we want to detect
The variability in the data
The chosen significance level

Generally, power increases when the sample size increases. A larger sample provides more information and therefore makes it easier to distinguish a real effect from random sampling variation.

R provides the power.t.test() function to perform power and sample size calculations for t-tests.

For example, suppose we want to know how much power we would have with a sample of 50 observations if we expect an effect size of 0.5.

power.t.test(
  n = 50,
  delta = 0.5,
  sd = 1,
  sig.level = 0.05,
  type = "one.sample",
  alternative = "two.sided"
)
#> 
#>      One-sample t test power calculation 
#> 
#>               n = 50
#>           delta = 0.5
#>              sd = 1
#>       sig.level = 0.05
#>           power = 0.9338976
#>     alternative = two.sided

Here:

n is the sample size
delta is the difference between the mean under the alternative hypothesis and the mean under the null hypothesis
sd is the assumed standard deviation
sig.level is the significance level alpha
type = “one.sample” specifies that we are dealing with a one-sample t-test;
alternative = “two.sided” specifies a two-sided test.

The output includes the estimated power for these assumptions.

We can also use power.t.test() to determine the sample size needed to achieve a desired level of power.

Instead of specifying n, we specify the desired power:

power.t.test(
  power = 0.80,
  delta = 0.5,
  sd = 1,
  sig.level = 0.05,
  type = "one.sample",
  alternative = "two.sided"
)
#> 
#>      One-sample t test power calculation 
#> 
#>               n = 33.3672
#>           delta = 0.5
#>              sd = 1
#>       sig.level = 0.05
#>           power = 0.8
#>     alternative = two.sided

R will calculate the required sample size to achieve approximately 80% power, assuming the other values are correct.

6 Statistical testing

In this section, we will look at how to perform several of the statistical tests discussed in the course using R. The focus here is mainly on the practical implementation: how to specify the test in R, which arguments are important, and how to find the relevant information in the output.

The choice of statistical test depends on the type of variables, the number of groups, and whether the observations are independent or paired. The theoretical background and assumptions of these tests are discussed in the course slides. Here, we focus on translating these analyses into R.

6.1 Comparing two means

When we want to compare the mean of a numerical variable between two groups, we can use a t-test. The appropriate version depends on whether the observations in the two groups are independent or paired.

6.1.1 Independent samples t-test

An independent samples t-test is used when the two groups consist of different, independent observations.

For example, suppose we want to compare the mean price between two groups defined by area.

t.test(price~area, data=apartment)
#> 
#>  Welch Two Sample t-test
#> 
#> data:  price by area
#> t = -1.9675, df = 92.859, p-value = 0.05211
#> alternative hypothesis: true difference in means between group Non-residential and group Residential is not equal to 0
#> 95 percent confidence interval:
#>  -27100.3514    125.9925
#> sample estimates:
#> mean in group Non-residential     mean in group Residential 
#>                      97282.05                     110769.23

The general structure is:

t.test(numerical_variable ~ grouping_variable, data = data)

Here:

price is the numerical variable we want to compare.

area is the categorical variable defining the two groups.

data = apartment tells R that these variables can be found in the dataset called apartment.

The ~ symbol can be read as “by”. Thus, price ~ area means that we want to compare price by area.

By default, t.test() performs Welch’s two-sample t-test, which does not assume that the variances in the two groups are equal. If we specifically want to perform the version that assumes equal variances, we can use:

t.test(price~area, data=apartment, var.equal = TRUE)
#> 
#>  Two Sample t-test
#> 
#> data:  price by area
#> t = -1.8258, df = 115, p-value = 0.07048
#> alternative hypothesis: true difference in means between group Non-residential and group Residential is not equal to 0
#> 95 percent confidence interval:
#>  -28119.51   1145.15
#> sample estimates:
#> mean in group Non-residential     mean in group Residential 
#>                      97282.05                     110769.23

The var.equal argument therefore specifies whether we want to assume equal variances: var.equal = TRUE or var.equal = FALSE. The default is FALSE.

The output contains the t-statistic, degrees of freedom, p-value, confidence interval, and the estimated mean in each group.

We can also specify whether we want a two-sided or one-sided test using the alternative argument:

# For one-sided test:
t.test(price~area, data=apartment, alternative="greater")
#> 
#>  Welch Two Sample t-test
#> 
#> data:  price by area
#> t = -1.9675, df = 92.859, p-value = 0.9739
#> alternative hypothesis: true difference in means between group Non-residential and group Residential is greater than 0
#> 95 percent confidence interval:
#>  -24876.47       Inf
#> sample estimates:
#> mean in group Non-residential     mean in group Residential 
#>                      97282.05                     110769.23
t.test(price~area, data=apartment, alternative="less")
#> 
#>  Welch Two Sample t-test
#> 
#> data:  price by area
#> t = -1.9675, df = 92.859, p-value = 0.02606
#> alternative hypothesis: true difference in means between group Non-residential and group Residential is less than 0
#> 95 percent confidence interval:
#>       -Inf -2097.892
#> sample estimates:
#> mean in group Non-residential     mean in group Residential 
#>                      97282.05                     110769.23

# For two-sided test:

t.test(price~area, data=apartment, alternative="two.sided")
#> 
#>  Welch Two Sample t-test
#> 
#> data:  price by area
#> t = -1.9675, df = 92.859, p-value = 0.05211
#> alternative hypothesis: true difference in means between group Non-residential and group Residential is not equal to 0
#> 95 percent confidence interval:
#>  -27100.3514    125.9925
#> sample estimates:
#> mean in group Non-residential     mean in group Residential 
#>                      97282.05                     110769.23

The default is “two.sided”.

6.1.2 Non-parametric independent t-test: Wilcoxon rank-sum test

If the assumptions of the independent samples t-test are not appropriate, we can use the Wilcoxon rank-sum test as a non-parametric alternative. The Wilcoxon test works with the ranks of the observations rather than directly comparing their means. It is therefore often used when the assumptions required for the t-test are not met. In R, this is performed using wilcox.test():

wilcox.test(price~area, data=apartment)
#> 
#>  Wilcoxon rank sum test with continuity correction
#> 
#> data:  price by area
#> W = 1240.5, p-value = 0.1054
#> alternative hypothesis: true location shift is not equal to 0

6.1.3 Paired samples t-test

Sometimes the two groups are not independent. Instead, each observation in one group is directly related to an observation in the other group. For example, suppose we measure the blood pressure of the same individuals before and after an intervention. Because each person’s measurements are paired, we need a paired t-test.

For this purpose, we will use a different dataset: bloodpressure.txt: we measure the blood pressure of a group of nonpregnant, premenopausal women of age 16-49 who do not use oral contraceptives (OC). After this the women start to take OC and one 1 year later we measure the blood pressure again. The dataset has 4 variables:

Index: indicator for the women
SBPstart: Systolic blood pressure before taking OC
SBPend: Systolic blood pressure after taking OC for a year
Difference: Difference in systolic blood pressure (end-start)

bp <- read.table("bloodpressure.txt", header=T)
head(bp)
#>   index SBPstart SBPend Difference
#> 1     1      115    128         13
#> 2     2      112    115          3
#> 3     3      107    106         -1
#> 4     4      119    128          9
#> 5     5      115    122          7
#> 6     6      138    145          7

We can perform a paired t-test using the t.test() function, specifying the paired = TRUE option.

t.test(bp$SBPstart,bp$SBPend, paired=T)
#> 
#>  Paired t-test
#> 
#> data:  bp$SBPstart and bp$SBPend
#> t = -3.3247, df = 9, p-value = 0.008874
#> alternative hypothesis: true mean difference is not equal to 0
#> 95 percent confidence interval:
#>  -8.066013 -1.533987
#> sample estimates:
#> mean difference 
#>            -4.8

The paired = TRUE option tells R that the two measurements are paired and should not be treated as independent observations.

6.1.4 Non-parametric dependent t-test: Wilcoxon signed-rank test

The non-parametric alternative to the paired t-test is the Wilcoxon signed-rank test. We use the same wilcox.test() function, but specify: paired = TRUE

wilcox.test(bp$SBPstart,bp$SBPend, paired=T)
#> Warning in wilcox.test.default(bp$SBPstart, bp$SBPend, paired = T): cannot
#> compute exact p-value with ties
#> 
#>  Wilcoxon signed rank test with continuity correction
#> 
#> data:  bp$SBPstart and bp$SBPend
#> V = 3.5, p-value = 0.01646
#> alternative hypothesis: true location shift is not equal to 0

6.2 Comparing more than two means: ANOVA

When we want to compare the means of more than two groups, we can use an analysis of variance (ANOVA).

In our apartment data, we have a numerical variable age and a categorical variable priceclass with three levels. If we want to compare the mean age of the buyers in the three different price classes, we can perform a one-way ANOVA using aov():


model <- aov(age ~ factor(priceclass), data = apartment)

summary(model)
#>                     Df Sum Sq Mean Sq F value  Pr(>F)   
#> factor(priceclass)   2   1373   686.7   5.367 0.00592 **
#> Residuals          114  14588   128.0                   
#> ---
#> Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

The general structure is:

model <- aov(numerical_variable ~ factor(grouping_variable), data = data)

summary(model)

The first line creates the ANOVA model and stores it in an object called model.

The second line displays the results of the ANOVA.

The output contains an F-statistic and a p-value. The p-value tells us whether there is evidence that the group means are not all equal.

6.2.1 Pairwise comparisons

When comparing more than two groups with an ANOVA, a significant result tells us that at least one group differs from another, but it does not tell us which groups are different.

We can perform pairwise comparisons between the groups using pairwise.t.test(). For example, if we want to compare age between the different levels of priceclass:


pairwise.t.test(apartment$priceclass, apartment$priceclass)
#> 
#>  Pairwise comparisons using t tests with pooled SD 
#> 
#> data:  apartment$priceclass and apartment$priceclass 
#> 
#>   1      2     
#> 2 <2e-16 -     
#> 3 <2e-16 <2e-16
#> 
#> P value adjustment method: holm

The first argument is the numerical variable we want to compare, and the second argument is the grouping variable. R will perform a separate t-test for each possible pair of groups and provide the p-value for each comparison.

When we perform several statistical tests, the probability of obtaining a significant result by chance increases. We therefore need to adjust the p-values for multiple comparisons. We can specify the method used for this adjustment with the p.adjust.method argument.


pairwise.t.test(apartment$priceclass, apartment$priceclass,
  p.adjust.method = "bonferroni")
#> 
#>  Pairwise comparisons using t tests with pooled SD 
#> 
#> data:  apartment$priceclass and apartment$priceclass 
#> 
#>   1      2     
#> 2 <2e-16 -     
#> 3 <2e-16 <2e-16
#> 
#> P value adjustment method: bonferroni

The Bonferroni correction is a simple and conservative way to adjust the p-values.

Another commonly used adjustment is Holm’s method:


pairwise.t.test(apartment$priceclass, apartment$priceclass,
  p.adjust.method = "holm")
#> 
#>  Pairwise comparisons using t tests with pooled SD 
#> 
#> data:  apartment$priceclass and apartment$priceclass 
#> 
#>   1      2     
#> 2 <2e-16 -     
#> 3 <2e-16 <2e-16
#> 
#> P value adjustment method: holm

Holm’s method is generally less conservative than Bonferroni while still controlling the overall error rate.

The output is a table showing the adjusted p-value for each pairwise comparison. We can then see which specific groups differ from each other.

Another commonly used post-hoc test after a one-way ANOVA is Tukey’s HSD (Honestly Significant Difference). Tukey’s HSD compares all pairs of group means while adjusting for multiple comparisons.

model <- aov(age ~ factor(priceclass), data = apartment) 
TukeyHSD(model)
#>   Tukey multiple comparisons of means
#>     95% family-wise confidence level
#> 
#> Fit: aov(formula = age ~ factor(priceclass), data = apartment)
#> 
#> $`factor(priceclass)`
#>          diff        lwr       upr     p adj
#> 2-1 -6.938596 -12.657394 -1.219799 0.0130439
#> 3-1 -8.291667 -15.370646 -1.212687 0.0172859
#> 3-2 -1.353070  -7.889662  5.183522 0.8754681

The output contains a comparison for each pair of groups. For each comparison, you will see:

diff: the difference between the two group means
lwr: lower limit of the confidence interval
upr: upper limit of the confidence interval
p adj: p-value adjusted for multiple comparisons

Important: Pairwise comparisons should generally be performed as follow-up analyses after an overall test such as ANOVA, rather than simply performing many separate t-tests without adjusting for multiple comparisons.

6.2.2 Non-parametric anova: Kruskal-Wallis test

The non-parametric alternative to a one-way ANOVA is the Kruskal-Wallis test.


kruskal.test(age~priceclass, data=apartment)
#> 
#>  Kruskal-Wallis rank sum test
#> 
#> data:  age by priceclass
#> Kruskal-Wallis chi-squared = 10.602, df = 2, p-value = 0.004987

The Kruskal-Wallis test is based on the ranks of the observations rather than directly comparing the means.

As with ANOVA, a significant Kruskal-Wallis test tells us that there is evidence of a difference between the groups, but does not by itself tell us which particular groups differ.

6.3 Testing assumptions

Many statistical tests make certain assumptions about the data. Before performing a statistical test, it is therefore important to check whether these assumptions are reasonably met. The assumptions we need to check depend on the statistical test. Two common assumptions we encounter are:

Normality: the data, or more specifically the relevant residuals/differences for a particular test, are approximately normally distributed.
Homogeneity of variances: the variability is approximately equal across the groups being compared.

In R, we can assess these assumptions using both graphical methods and statistical tests. Graphical methods are particularly useful because they allow us to see the actual distribution of the data, while statistical tests provide a formal test of a specific assumption.

6.3.1 Checking normality

One assumption of several parametric statistical tests is that the relevant data are approximately normally distributed. There are several ways to assess normality. Two common approaches are:

a QQ plot, which provides a visual assessment;
a Shapiro-Wilk test, which provides a formal statistical test.

A QQ plot compares the distribution of our data with what we would expect if the data came from a normal distribution. In base R, we can create a QQ plot using qqnorm(). A reference line can be added using qqline().

qqnorm(apartment$price)
qqline(apartment$price)

The interpretation is primarily visual. If the observations approximately follow the reference line, the data are reasonably consistent with a normal distribution. The points do not need to lie exactly on the line. In real data, some deviation from the line is expected. We are interested in whether there are substantial or systematic deviations from normality. A QQ plot is particularly useful for identifying features such as strong skewness or unusually extreme observations.

We can also perform a formal test of normality using the Shapiro-Wilk test:

shapiro.test(apartment$price)
#> 
#>  Shapiro-Wilk normality test
#> 
#> data:  apartment$price
#> W = 0.86152, p-value = 4.561e-09

The null hypothesis of the Shapiro-Wilk test is that the data are normally distributed. We can interpret the p-value in the usual way:

p < 0.05: there is evidence against normality.
p >= 0.05: there is not sufficient evidence to conclude that the data deviate from normality.

However, the Shapiro-Wilk test should not be used on its own. Its result is strongly influenced by the sample size. With a large sample, even small and practically unimportant deviations from normality can result in a significant p-value. With a small sample, substantial deviations may not be detected.

For this reason, it is useful to consider the QQ plot together with the Shapiro-Wilk test.

When comparing two groups with an independent samples t-test, we are interested in the distribution of the numerical variable within each group.

For example, if we compare price between areas, we can create separate QQ plots:

qqnorm(apartment$price[apartment$area=="Non-residential"], main="QQplot Price for Non-Residential area")
qqline(apartment$price[apartment$area=="Non-residential"])

qqnorm(apartment$price[apartment$area=="Residential"], main="QQplot Price for Residential area")
qqline(apartment$price[apartment$area=="Residential"])

We can also perform a Shapiro-Wilk test separately for each group:

shapiro.test(apartment$price[apartment$area=="Non-residential"])
#> 
#>  Shapiro-Wilk normality test
#> 
#> data:  apartment$price[apartment$area == "Non-residential"]
#> W = 0.85483, p-value = 0.0001416
shapiro.test(apartment$price[apartment$area=="Residential"])
#> 
#>  Shapiro-Wilk normality test
#> 
#> data:  apartment$price[apartment$area == "Residential"]
#> W = 0.86271, p-value = 5.789e-07

For a paired t-test, the relevant assumption concerns the differences between the paired observations, rather than the two variables separately. For example:

difference <- bp$SBPend-bp$SBPstart

qqnorm(difference)
qqline(difference)

shapiro.test(difference)
#> 
#>  Shapiro-Wilk normality test
#> 
#> data:  difference
#> W = 0.97703, p-value = 0.9473

6.3.2 Checking equality of variances

When comparing the means of two or more independent groups, we may also need to consider whether the groups have similar variances. This is referred to as homogeneity of variances or equal variances.

We can visually investigate the spread of the groups using a boxplot:

boxplot(
  price ~ area,
  data = apartment,
  main = "Price by Area",
  xlab = "Area",
  ylab = "Price"
)

The boxplots allow us to compare the spread of the observations between the groups.

We can also formally test the equality of variances using Levene’s test. Levene’s test is available in the car package. In order to be able to use the car package, it first needs to be installed and then loaded into the R environment:

install.packages('car')

We can then perform the test:

library(car)
#> Loading required package: carData
leveneTest(price ~ area, data = apartment)
#> Levene's Test for Homogeneity of Variance (center = median)
#>        Df F value Pr(>F)
#> group   1    2.38 0.1256
#>       115

The general syntax is:

leveneTest(numerical_variable ~ grouping_variable, data = data)

The null hypothesis of Levene’s test is that the variances are equal across the groups. As usual, we can use the p-value to assess the evidence against the null hypothesis:

p < 0.05: evidence that the variances differ between the groups.
p >= 0.05: not enough evidence to conclude that the variances differ.

Remember that a non-significant result does not prove that the variances are identical. It means that the test did not find sufficient evidence of a difference.

6.4 Association between categorical variables

The tests above are mainly used when we have a numerical outcome and want to compare groups. When both variables are categorical, we can instead investigate whether there is an association between the variables.

For example, suppose we want to investigate whether the priceclass is associated with area. We can first create a contingency table:

table(apartment$area, apartment$priceclass)
#>                  
#>                    1  2  3
#>   Non-residential 13 21  5
#>   Residential     23 36 19

This gives us the number of observations in each combination of the two variables.

6.4.1 Chi-square test

The chi-square test of independence can be used to test whether two categorical variables are associated.

chisq.test(table(apartment$area, apartment$priceclass))
#> 
#>  Pearson's Chi-squared test
#> 
#> data:  table(apartment$area, apartment$priceclass)
#> X-squared = 2.1283, df = 2, p-value = 0.345

The table() function creates the contingency table, which is then supplied to chisq.test().

The output contains the chi-square statistic, degrees of freedom, and p-value.

The chi-square test is based on comparing the observed frequencies with the frequencies we would expect if the two variables were independent. R allows us to inspect these expected frequencies:

test <- chisq.test(table(apartment$area, apartment$priceclass))
test$expected
#>                  
#>                    1  2  3
#>   Non-residential 12 19  8
#>   Residential     24 38 16

This can be useful for checking whether the assumptions of the chi-square test are reasonable.

6.4.2 Fisher’s exact test

When the sample size is small, or when some of the expected frequencies in the contingency table are very small, the chi-square approximation may not be appropriate. In this situation, we can use Fisher’s exact test.

fisher.test(table(apartment$area, apartment$priceclass))
#> 
#>  Fisher's Exact Test for Count Data
#> 
#> data:  table(apartment$area, apartment$priceclass)
#> p-value = 0.3689
#> alternative hypothesis: two.sided

Fisher’s exact test calculates the exact probability based on the observed table rather than relying on the chi-square approximation.

7 Simple Linear Regression

Linear regression is used to describe the relationship between two numerical variables and to make predictions. We use one variable, the predictor or explanatory variable, to predict another variable, the outcome or response variable.

For example, we might want to investigate whether price is related to tax:

price = predictor variable
tax = outcome variable

7.1 Scatterplot

The first step is to visualise the relationship using a scatterplot. A scatterplot allows us to see whether there appears to be a linear relationship between two numerical variables.

plot(
  apartment$price,
  apartment$tax,
  main = "Tax vs Price",
  xlab = "Price",
  ylab = "Tax",
  pch = 19
)

Here, each point represents one observation.When looking at a scatterplot, consider:

Is there a roughly linear relationship?
Is the relationship positive or negative?
Are there any unusual observations or outliers?

7.2 Performing linear regression

In R, linear regression can be performed using the lm() function.

model <- lm(tax ~ price, data = apartment)

The ~ symbol can be read as “is modelled by” or “as a function of”.

To see the results of the regression model, use:

summary(model)
#> 
#> Call:
#> lm(formula = tax ~ price, data = apartment)
#> 
#> Residuals:
#>     Min      1Q  Median      3Q     Max 
#> -420.52  -79.35    3.27   77.71  460.82 
#> 
#> Coefficients:
#>              Estimate Std. Error t value Pr(>|t|)    
#> (Intercept) 2.229e+01  4.109e+01   0.542    0.589    
#> price       7.122e-03  3.642e-04  19.556   <2e-16 ***
#> ---
#> Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
#> 
#> Residual standard error: 149.2 on 115 degrees of freedom
#> Multiple R-squared:  0.7688, Adjusted R-squared:  0.7668 
#> F-statistic: 382.5 on 1 and 115 DF,  p-value: < 2.2e-16

The output contains several important pieces of information.

The Coefficients section of the output contains the intercept and slope.

For example:

Coefficients:
             Estimate Std. Error t value Pr(>|t|)    
(Intercept) 2.229e+01  4.109e+01   0.542    0.589    
price       7.122e-03  3.642e-04  19.556   <2e-16 ***

The regression equation can be written as:

predicted tax = intercept + slope × price

In this example:

predicted tax = 2.229e+01 + 7.122e-03 × height

The intercept is the predicted value of the outcome when the predictor is equal to zero. In this example, the intercept is 2.229e+01. This means that the model predicts a tax of 2.229e+01 when the price is zero.

The intercept is necessary for defining the regression line, but its practical interpretation is not always meaningful. For example, a price of zero may be outside the range of realistic observations.

The slope describes how much the predicted outcome changes for a one-unit increase in the predictor.

Here, the slope is 7.122e-03. This means that for every one-unit increase in price, the predicted tax increases by 7.122e-03 units, on average. A positive slope indicates a positive relationship, while a negative slope indicates a negative relationship.

7.3 Hypothesis test for the slope

The regression output also provides a hypothesis test for each coefficient. For the slope, the null and alternative hypotheses are:

H0: slope = 0
H1: slope is not equal to 0

The Pr(>|t|) column gives the p-value for this test.

For example:

Coefficients:
             Estimate Std. Error t value Pr(>|t|)    
(Intercept) 2.229e+01  4.109e+01   0.542    0.589    
price       7.122e-03  3.642e-04  19.556   <2e-16 ***

A small p-value provides evidence that the slope differs from zero. In this case, there is evidence of a linear association between price and age, because the p-value for the slope is <2e-16.

7.4 R-squared

The regression output also contains the R-squared value.

For example:

Multiple R-squared:  0.7688

R-squared describes the proportion of variation in the outcome that is explained by the linear regression model. An R-squared of 0.77 means that approximately 77% of the variation in tax is explained by price in this linear model.

The R-squared is always a value between 0 and 1:

values closer to 0 indicate that the model explains relatively little of the variation
values closer to 1 indicate that the model explains more of the variation

A high R-squared does not necessarily mean that the model is appropriate, and a low R-squared does not necessarily mean that a model is useless. The context and the purpose of the analysis are important.

7.5 Plotting the regression line

The regression line can be added to the scatterplot using abline().

plot(
  apartment$price,
  apartment$tax,
  main = "Tax vs Price",
  xlab = "Price",
  ylab = "Tax",
  pch = 19
)

abline(model)

7.6 Checking the assumptions

Before interpreting a linear regression model, it is important to check whether its assumptions are reasonable.

Important assumptions include:

Linearity – the relationship between predictor and outcome is approximately linear.
Normality of residuals – the residuals are approximately normally distributed.
Homoscedasticity – the variability of the residuals is approximately constant across the fitted values.
Observations should be independent.

Regression assumptions are mainly assessed using the residuals.

A residual is the difference between the observed value and the value predicted by the regression model.

A residual plot can be created using:

plot(
  model$fitted.values,
  model$residuals,
  xlab = "Fitted values",
  ylab = "Residuals",
  main = "Residual plot",
  pch = 19
)

abline(h = 0)

The residuals should be scattered randomly around zero.

A reasonable residual plot should not show a clear pattern or systematic shape.

For example:

a curved pattern can indicate that the relationship is not linear
a funnel-shaped pattern can indicate that the variability changes across fitted values
individual points far from the others may indicate unusual observations

The residual plot therefore helps us assess both linearity and homoscedasticity.

The normality assumption concerns the residuals, rather than the original outcome variable.

A Q-Q plot can be used to assess normality:

qqnorm(model$residuals)
qqline(model$residuals)

If the points follow the reference line reasonably closely, the residuals are approximately normally distributed.

The Shapiro-Wilk test can also be used:

shapiro.test(model$residuals)
#> 
#>  Shapiro-Wilk normality test
#> 
#> data:  model$residuals
#> W = 0.97332, p-value = 0.01963

The null hypothesis of this test is that the residuals are normally distributed.

As with other normality tests, the result should not be interpreted in isolation. The Q-Q plot provides useful visual information, particularly because formal normality tests can be sensitive to sample size.

7.7 Making predictions

Once a regression model has been fitted, it can be used to predict the outcome for a new observation.

Suppose we want to predict the tax of an appartment with a price of 250.000 euro:

new_data <- data.frame(price = 250000)

predict(model, newdata = new_data)
#>        1 
#> 1802.695

This gives the predicted mean outcome for an observation with a price of 250.000 euro.

When predicting an individual new observation, there is uncertainty around the prediction. A prediction interval gives a range of plausible values for the outcome of a new individual observation.

For example, to obtain a 95% prediction interval:

predict(
  model,
  newdata = new_data,
  interval = "prediction",
  level = 0.95
)
#>        fit     lwr      upr
#> 1 1802.695 1488.29 2117.101

The output contains:

```
`fit` – predicted value  
```

`lwr` – lower limit of the prediction interval

`upr` – upper limit of the prediction interval

For example:

      fit     lwr      upr
1 1802.695 1488.29 2117.101

This means that the predicted tax is 1802.695 with a 95% prediction interval from 1488.29 to 2117.101.

Basic Principles of Statistics

dr. Annelies Agten

2026-09-30

1 Getting started in R

1.1 Reading in a txt file

1.2 Viewing the data

1.3 Extracting a variable

2 Descriptive statistics

2.1 Checking and changing the type of a variable

2.2 Frequency tables

2.3 Barplot

2.4 Histogram

2.5 Boxplot

2.6 Summary statistics

2.6.1 Measures of central tendency

2.6.2 Measures of spread

3 Parameter estimation

3.1 Simulating samples

4 Confidence interval for the mean

5 Hypothesis testing

5.1 One-sided versus two-sided tests

5.2 Statistical power and sample size

6 Statistical testing

6.1 Comparing two means

6.1.1 Independent samples t-test

6.1.2 Non-parametric independent t-test: Wilcoxon rank-sum test

6.1.3 Paired samples t-test

6.1.4 Non-parametric dependent t-test: Wilcoxon signed-rank test

6.2 Comparing more than two means: ANOVA

6.2.1 Pairwise comparisons

6.2.2 Non-parametric anova: Kruskal-Wallis test

6.3 Testing assumptions

6.3.1 Checking normality

6.3.2 Checking equality of variances

6.4 Association between categorical variables

6.4.1 Chi-square test

6.4.2 Fisher’s exact test

7 Simple Linear Regression

7.1 Scatterplot

7.2 Performing linear regression

7.3 Hypothesis test for the slope

7.4 R-squared

7.5 Plotting the regression line

7.6 Checking the assumptions

7.7 Making predictions