This document contains the practical exercises for the course Basic Principles of Statistics. The exercises are designed to help you apply the statistical concepts and R commands introduced during the course to realistic datasets.

The aim is not simply to reproduce the code shown in the course, but to develop your ability to choose an appropriate analysis, perform it in R, and interpret the results.

For each exercise, you will be given a research question and a dataset. You will be asked to work through the analysis step by step. Depending on the exercise, this may involve:

  • Exploring and visualising the data
  • Calculating descriptive statistics
  • Choosing and performing an appropriate statistical test
  • Checking assumptions
  • Interpreting statistical output
  • Fitting a regression model
  • Making predictions from a statistical model

You can use the R Markdown guideline alongside these exercises. The guideline provides examples of the R commands and explains how to interpret the output. During the exercises, try to decide which commands and methods are appropriate for the question before looking at the examples.

There is often more than one way to perform an analysis in R. In these exercises, the focus is on learning a small number of reliable approaches that you can understand and apply to your own data.

Practical 1: Descriptive statistics

Exercise 1

In the demodata dataset you can find the participant list from 3 different StatUa workshops. The following variables are recorded: - The names of the participants (ID)
- The course they followed (workshop)
- Gender (gender)
- Exam result (exam) - The answer to an evaluation question at the end of the course (question 1: “Did you like the course ? “, with answers: 1=not at all, 2=not quite, 3=neutral, 4=yes, 5=very much).

Load the dataset into R and check the type of each variable. To read excel files into R, you can use the read_excel() function from the package readxl. Not that you first have to install and load the package into the R environment, before you can use the function:

install.packages('readxl')
library(readxl)
#> Warning: package 'readxl' was built under R version 4.5.3
demodata <- read_excel("demodata.xls")

Is each categorical variable recognized by R as the correct type? If not, change the type of the variable using the factor() function.

The variable gender was fully written out in the original data file as “female” and “male”. For categorical (factor) variables, it is good practice and necessary for certain statistical tests to encode the different factor levels as numbers. For instance, female = 0 and male = 1. Make a new variable sex that contains the gender variable encoded by numbers.

The new variable sex now has two possible values 0 and 1. In the output of our analyses (plots, tables), this is not very descriptive. To have the names of the factor levels in our output, we attach a label to the factor levels, using the function levels().

The variable exam is recognized by R as a ‘character’. Can you explain why this is the case?

Exercise 2

The dataset morphometrics contains several morphometric measurements on students from the University of Antwerp. Load the dataset into R and check what the data look like:

morphometrics <- read.table("morphometric.txt",header=T)
head(morphometrics)
#>   smoke length weight sex    oz_date     b_date haircolor eyecolor heeltoe bal
#> 1     0   1.77     69   0 2002-02-24 1981-08-07     bruin    bruin     240  95
#> 2     0   1.78     66   0 2002-02-25 1981-03-30       ros    bruin     240  85
#> 3     1   1.59     56   1 2002-02-22 1979-08-10     blond    blauw     235  92
#> 4     0   1.65     62   0 2002-02-22 1981-05-22     bruin    blauw     231  91
#> 5     0   1.64     60   0 2002-02-24 1980-07-23     blond    blauw     234  86
#> 6     0   1.63     52   0 2002-02-25 1981-09-14     bruin    blauw     230  82

In the current dataset, the variables length and weight are expressed in meters and kilo, respectively. However, we want to make two new variables where length and weight are expressed in centimeters and gram, respectively. This can be done using the $ operator. For example:

morphometrics$length_cm <- morphometrics$length*100

Create a new variable yourself that contains the variable weight recorded in gram.

Using the functions mean() and var() we can calculate the mean and the variance of continuous variables:

mean(morphometrics$length)
#> [1] NA
var(morphometrics$length)
#> [1] NA

When you do this, however, you see that the results is ‘NA’. This is due to the fact that there are missing observations in the dataset. You can tell R to ignore these missing observations, by adding the option na.rm = TRUE to the mean and variance statement. This option stands for ‘remove the missing values from the variable’.

mean(morphometrics$length, na.rm=TRUE)
#> [1] 1.701604
var(morphometrics$length, na.rm=TRUE)
#> [1] 0.006329308

Calculate the mean and variance for both the original length and weight variables, as well as for the newly created length_cm and weight_gram variables. How do the mean and the variances from the original and the transformed variables relate?

Exercise 3

Load the dataset bldpres2 into R and explore the data. This dataset represents the results of a random sample (n=328) collected in 2 villages in Pakistan. The aim of the study was to elucidate the influence of environmental risk factors on blood pressure. The variables that are recorded are:

  • Age (in years)
  • Sys (systolic blood pressure in mmHg)
  • Dias (diastolic blood pressure in mmHg)
  • Weight (in kg)
  • Pulse (heart rate in beats/min)
  • Sex (0 = male, 1 = female )
  • SES (socioeconomic status: 1 = low ; 2 = middle ; 3 = high)
bldpres2 <- read.table("bldpres2.txt", header=T)

Set all variables to the correct type.

Make the following graphs to explore the variables in the dataset:

  • Barplot for variable sex and for SES.
  • Stacked barplot showing the info for sex and SES simultaneously.
  • Boxplot for systolic blood pressure.
  • Boxplots for systolic blood pressure for each category of SES

Make a frequency table for the variable SES.

By using the function summary() you can get a quick overview of some summary statistics fro the variables in a dataset:

summary(bldpres2)

Look at the variable age:

  • What is the smallest resp. largest observed value for age?
  • What is the age range ?
  • What is the median?
  • What are the values for p25 resp p75 for age ?
  • What is the interquartile range?
  • Look at the distribution of the variable age using a histogram. Is the distribution symmetric or skewed?
  • Look at the distribution of the variable age using a boxplot. Are there any outliers?

Make a scatterplot of systolic blood pressure (sys) versus body weight (weight). Do the variables seem positively, negatively or non-correlated?

Calculate the correlation between systolic blood pressure and body weight. Does this correspond to what you see in the scatterplot?

Exercise 4

The grrtx2 dataset contains data from 52 patients who underwent a first kidney transplant before age 15. At the end of the study their final height was recorded. The aim of the study was looking for factors influencing height:

  • rangnr: patient ID
  • finalht: final body length
  • sex: 0 = man ; 1 = woman
  • tarht: target height
  • diagn: primary kidney disease (1=urinairy tract abnormality , 2=glomerulopathy)
  • htsdshd: standard deviation of body length at the start of hemodialysis (accounting for age and sex)
  • durhd: duration (in years) of hemodialysis
  • agetx: age at first kidney transplant
  • htsdtx: standard deviation of body length at first kidney transplant (accounting for age and sex)
  • txnum: total number of kidney transplants

Create a boxplot for the variable finalht, grouped for each level of sex. Something is clearly wrong. One of the individuals is obviously very very small. The outlier is individual nr. 50

Look at the outlying variable in the data. If you know that the data were manually entered, it is quite obvious what typing error was made here. What would be the most likely (true) value?

Make a (simple) scatterplot for finalht (Y-axis) versus tarht (X-axis). Can you explain what is going on here?

Practical 2: Parameter estimation

Exercise 1

The aim of parameter estimation is to estimate the properties of an unknown population based upon a limited sample that was randomly chosen from that population. In this practical, we will simulate random sampling using a large (artificial) dataset of 10,000 persons, which represents the total population. We will then take limited random samples from that population, and illustrate how we can draw conclusions about the large population based upon a limited sample, accounting for the uncertainty caused by the random sampling.

Let’s assume that we have a population of 10,000 people that all take an exam. We know that the scores of the exam are approximately normally distributed with a mean equal to 11 and standard deviation equal to 6. We can generate a random population with this distribution by using the rnorm() function:

set.seed(123)
population <- rnorm(n=10000, mean=11, sd=6)

Unlike in typical research, here we do know the properties of the entire population. We can calculate the mean and the variance of our population:

mean(population)
#> [1] 10.98577
var(population)
#> [1] 35.9019

These are the true, underlying population parameters, denoted with the Greek letters mu (for mean) and sigma (for standard deviation). If we assume that the score is normally distributed in the population, this distribution is completely described using these two parameters.

Usually, these parameters are unknown as it is practically impossible to analyze the entire population. Researchers will try to estimate the parameters of the entire population based upon a random sample. And most researchers can only draw one single sample, but since we’re working on simulated data, we can draw as many samples as we want.

You can draw a random sample from exactly 50 individuals from the total population by using the sample() function.

sample1 <- sample(population, size=50)

Calculate the mean and the variance of the sample.

The mean and variance we’ve calculated here, are not the true population mean and variance, but estimates based upon a random sample. A random sample is subject to statistical fluctuation: if you take several random samples, your estimates won’t always be exactly the same. To illustrate this, generate 25 other random samples and calculate the mean and the variance. Not that you can automate this by adapting the code below:


sample_size <- 50
number_of_samples <- 5

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

Repeat these 25 sampling experiments but now taking a random sample of 12 individuals instead of 50.

Compare the results of the sampling experiment with sample size 12 or 50, and compare these to the true population mean. Check the single experiment where you have the largest error (i.e. the worst estimate) compared to the population mean?

Generate a histogram of the sample distribution of means for sample size 12, and one for sample size 50. Calculate the mean and variance for the both sample sizes (i.e. the mean of means and variance of means).

  • What is the mean of means for sample size 12 and sample size 50?
  • What is the variance around the mean of means for sample size 12 resp sample size 50?
  • Check whether the mean and variance of your sampling distribution of means are indeed the population mean mu and the population variance divided by the sample size.
  • Unlike us, who can draw several samples from the entire population, researchers typically can only draw one sample and have to infer the population mean from that one sample. Would you trust a researcher who drew one large sample, or the one who drew a small sample? Why?

Exercise 2

The Central Limit Theorem states that when we sample many times from a population (with any kind of distribution), and each time calculate the mean, these means form a normal distribution, with mean mu (true population mean) and sigma^2/n as variance. We will demonstrate this with the datasets sampling2 and sampling2b.

sampling2 <- read.table("sampling2.txt")
sampling2b <- read.table("sampling2b.txt")

The dataset sampling2 contains a variable finalexam. If we make a histogram of this variable, the results does not look normal at all. Yet you can estimate the mean using the normal approximation, when n is >25.

hist(sampling2$finalexam)

The dataset sampling2b contains the results of 200 samples taken from the population in the dataset sampling2.
- samplenumber: number of the sample (1 to 200)
- meanof30: mean of the sample, each sample containing 30 observations from the entire population
- meanof100: mean of samples with 100 observations
- meanof500: mean of samples with 500 observations
- sdof500: standard deviation of samples with 500 observations

We will demonstrate that with increasing sample size, the distribution of the means will be increasingly normal. Also, as the sample size increases, the estimated means will be increasingly close to the true population mean. The chance that your one sample makes a wrong estimate of the population mean, gets increasingly smaller with increasing sample size.

  • Compute mean and standard deviation of finalexam.
  • Using the dataset sampling2b, compute the mean and standard deviation of meanof30, meanof100 and meanof500. How do they compare to the values of the mean and variance of finalexam?
  • Make a histogram of meanof30, meanof100 and meanof500. What can you conclude?

For the samples of size 500 we have the mean and standard deviation of each sample so we can compute a 95% confidence interval for the population mean using these values. We can use these then to check if the true population mean calculated is included in the 95% confidence interval:

sampling2b$LCL = data.frame(sampling2b$meanof500-1.96*sampling2b$sdof500/sqrt(500))
sampling2b$UCL = data.frame(sampling2b$meanof500+1.96*sampling2b$sdof500/sqrt(500))

colnames(sampling2b$LCL) <- c("LCL")
colnames(sampling2b$UCL) <- c("UCL")

sampling2b$check = ifelse(((sampling2b$LCL$LCL <= 10.19) & (10.19 <= sampling2b$UCL$UCL)), 1, 0)

Make a frequency table of the variable Check. What are your findings? Does this make sense?

Practical 3: Hypothesis testing

Exercise 1

The dataset sampling3 contains measurements of height and IQ of Dutch citizens. We will take random samples out of this dataset and use a one-sample T-test to test the following null hypotheses:

  • The population mean of the length of the Dutch people is 185 cm
  • The population mean of the IQ of the Dutch people is 105.51

Using the entire population, test both hypotheses.

Take a random sample of 12 individuals. Using this sample, test both hypotheses. Repeat this process 15 times and record the p-values each time.

Count the number of times you find a significant difference in the random sample of 12 individuals. Compare this conclusion to the properties of the complete population. How many times did you make a correct conclusion based upon the sample of 12 individuals ?

On the complete population, you have analyzed whether the H0 is true for length and IQ, respectively. Type 1 and type 2 error describe how well the conclusions from your statistical test in the random sample, match the properties of the complete (underlying) population.

  • Type 1 error : Rejecting H0 when it’s true (or : Finding a significant difference when there actually is no difference)
  • Type 2 error : Accepting H0 when it’s not true (or: Failing to find a significant difference when in reality there is a difference)
  • Power : 1 – probability of a type 2 error (or: Your chance to detect a significant difference that is present in your population).

Now turn to the results of your samplings. Fill out the percentage of Type 1 or Type 2 error.

Exercise 2

The company AGEMO FARMA charges 1.23 EUR fixed costs per order, but the manager suspects after a while that this amount is not correct. To make an informed decision about increasing the fixed cost or not the manager will check the fixed cost of a number of orders and then perform a statistical test.

  • Which test would you perform? Formulate null and alternative hypothesis.
  • When would you decide to increase the costs?
  • How many orders does the manager have to check if she only wants 5% risk of increasing the cost if this is not necessary and not more than 10% risk to keep the 1.23 EUR if the mean cost is 1.36 EUR?

You may assume that the fixed cost per order is normally distributed with variance 0.04.

Practical 4: Association between continuous and categorical variable with 2 levels

The scopinaro dataset contains data from 344 persons with overweight and obesity (BMI>25), that were treated using one of two possible surgical interventions (variable surgery): either a gastric band(=0) or a gastric bypass(=1). Both interventions reduce the active volume of the stomach, but the latter one is much more invasive. BMI had been measured before (variable BMIpre) and after (BMIpost) the intervention.

The T-test tests whether a numeric variable differs between two groups. In this dataset, you can answer the following research questions using one of the t-tests we’ve seen.

Important! For each test we’ve seen, there is a parametric and a non-parametric variant. Before you start testing, always check whether the conditions to work parametrically are fulfilled.

Question 1

Is the BMI before the operation (BMIpre) different between males and females?

Question 2

Did the intervention lead to a significant difference in BMI in males? How large is the difference in BMI after operation (use variables BMIpre-BMIpost)? Give a 95% confidence interval. Perform the same analysis for females. Do you notice something unusual?

Question 3

Define the variable BMIdiff= BMIpre-BMIpost giving the difference in BMI before and after the operation. Using this variable offers an alternative for Question 2. Use the variable BMIdiff to test whether the difference in BMI before and after the operation differs from zero. Perform this test for the two genders separately, and compare the results with Question 2.

Question 4

Use the variable BMIdiff to test whether the drop in BMI differs between the two surgical techniques (variable surgery).

Question 5

Select only the persons with diabetes (variable diabetes=1). Test if the BMIpre in diabetic patients differs significantly between males and females.

Practical 5: Association between continuous and categorical variable with >2 levels

The database voicedata contains vocal characteristics of children up to age 12. It consists of healthy children (controls=0) and toddlers(=5), as well as children with diverse conditions (deaf=1, cerebral palsy=2, mental retardation=3, learning disorder=4). The research question is to what extent the vocal characteristics differ between the different groups of children (variable group). Vocal characteristics include :

  • The maximal phoneme time (variable mpt), being the time one can say “aaaaaaaaaaaahhhh”.
  • Jitter and Shimmer are measures of how hoarse the voice sounds
  • Harmonics to noise ratio (hn) is a measure of vocal purity

Answer the following questions using one of the hypothesis tests we’ve seen so far. Don’t forget to check whether you’re allowed to work with parametric tests!!

Question 1

Create a subset of the data that only includes children that are deaf, have mental retardation or learning disorders:

voicedata <- read.table("voicedata.txt", header=T)
data <- subset(voicedata, (voicedata$group ==1 | voicedata$group ==3 | voicedata$group ==4))

Analyze whether the maximal phoneme time mpt is different between the 3 groups of children.

Question 2

Analyze whether the jitter is different for different between the 3 groups of children.

In case you find an overall difference, find out what groups are different from each other.

Question 3

Include all children (i.e. use the original dataset).

Test whether the maximal phoneme time mpt differs between girls and boys. You may assume normality is met for both girls and boys.

Question 4

Analyze the difference in fundamental frequency (f0) according to group. You may assume that the normality and variance assumptions are met. Try different correction methods for multiple testing and check whether your conclusions change based on the chosen method.

Practical 6: Association between categorical variables

Going back to the scopinaro dataset, answer the following questions.

Question 1

Is there a relationship between comorbidity (variable comorb) and the occurrence of complications (variable complic)?

Question 2

Does the occurrence of complications differ between the two operation techniques (variable surgery)?

Question 3

Are there more complications (variable complic) in men compared to women?

Question 4

Researchers looked at 20 children of 14 years old and recorded if they have a mobile phone and if they have a Facebook account. The data can be summarized using the following table:

a <- c(10, 8)
b <- c(0, 2)
tab <- data.frame(a,b)
rownames(tab) <- c("facebook_yes","facebook_no")
colnames(tab) <- c("mobile_yes","mobile_no")
print(tab)
#>              mobile_yes mobile_no
#> facebook_yes         10         0
#> facebook_no           8         2

Test whether there is a relation between having a mobile phone and having a Facebook account.

Practical 7: Simple linear regression

A criminologist studying the relationship between level of education and crime rate in medium-sized US counties collected the following data for a random sample of 84 counties:

  • countyid: id of the county
  • education: percentage of individuals in the county having at least a high-school diploma
  • crimerate: crimes reported per 100 000 residents last year
  • income: 1=if the median income of that specific county is below the median over all county’s; 2=if the median income of that specific county is above the median over all county’s

The dataset is calles crimedata. In this exercise, we focus on the relation between the variables education level and crime rate.

Graphical analysis

Make a simple scatter plot of crime rate against education level. Add a regression line to the plot.

Fitting the model

Fit a regression model with crim rate as the dependent variable and education level as the independent variable. Look at the summary of the model.

  • Look at the p-value that tells you whether the regression model is significant or not. What is your conclusion?
  • Give the equation of the regression line.
  • What is the mean decrease in reported crimes per 100 000 residents per percentage increase of education level? Give the estimate and the corresponding 95% confidence interval.
  • What is the estimate of the residual standard deviation?
  • What is the percentage of variability in crime rate explained by education level?

Analysis of residuals

The model fit can be analyzed by looking at the residuals. Make a QQ-plot of the residuals and a residual plot. Discuss the assumptions of a linear regression model.

Formally test whether Normality of the residuals is fulfilled.

Prediction

Predict the crimerate in a county where the level of education is equal to 50%. Give a 99% confidence interval for this prediction.

Plot the 99% prediction interval of the regression line on a scatterplot.