This solution manual provides example solutions to the practical exercises in the course. It is intended as a guide to the analysis process, rather than as a set of definitive answers.

In R, there are often several ways to approach the same problem. Different commands, functions, or workflows may lead to the same conclusion, and there may be valid alternatives to the approach shown here.

The solutions should therefore be used to understand the reasoning behind an analysis, how to implement it in R, and how to interpret the results. They are not intended to represent the only correct way of solving each exercise.

When working on the exercises yourself, focus on the statistical question, the choice of method, and the interpretation of the results rather than trying to reproduce the solution exactly.

1 Practical 1: Descriptive statistics

1.1 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.

class(demodata$ID)
#> [1] "character"
class(demodata$Workshop)
#> [1] "numeric"
class(demodata$Gender)
#> [1] "character"
class(demodata$Question1)
#> [1] "numeric"
class(demodata$Exam)
#> [1] "character"

demodata$Workshop <- factor(demodata$Workshop)
demodata$Gender <- factor(demodata$Gender)

The variable gender was fully written out in the original data file as “male” and “female”. 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.

demodata$Sex <- demodata$Gender
levels(demodata$Sex) <- c(0,1)

#Another way of doing this, is by using the ifelse() function:

demodata$Sex <- ifelse(demodata$Gender=="male", 1, 0)

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().

demodata$Sex <- as.factor(demodata$Sex)
levels(demodata$Sex) <- c("female", "male")

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

demodata$Exam
#>  [1] "9"    "16.5" "14"   "12"   "?"    "10"   "14.5" "15"   "9"    "8"   
#> [11] "16"   "15"

Looking at the variable Exam we can see that one of the values is equal to ‘?’. This denotes a missing value. However, missing values in R are normally coded as ‘NA’ (Not Available). For R to recognize this question mark as a missing observation, rather than a value for the variable Exam, we can change this value to NA:

demodata$Exam[demodata$Exam == "?"] <- NA

To have R recognize the variable Exam as a continuous (numerical) variable, we can specify:

demodata$Exam <- as.numeric(demodata$Exam)

1.2 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.

morphometrics$weight_gram <- morphometrics$weight*1000

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?

mean(morphometrics$length, na.rm  = TRUE)
#> [1] 1.701604
mean(morphometrics$length_cm, na.rm  = TRUE)
#> [1] 170.1604
mean(morphometrics$weight, na.rm  = TRUE)
#> [1] 60.05566
mean(morphometrics$weight_gram, na.rm  = TRUE)
#> [1] 60055.66

var(morphometrics$length, na.rm  = TRUE)
#> [1] 0.006329308
var(morphometrics$length_cm, na.rm  = TRUE)
#> [1] 63.29308
var(morphometrics$weight, na.rm  = TRUE)
#> [1] 86.23563
var(morphometrics$weight_gram, na.rm  = TRUE)
#> [1] 86235634

The mean value of the new variable length_cm is equal to 100 times the mean of the original variable length. The mean value of the new variable weight_gram is equal to 1000 times the mean of the original variable weight.

The variance of the new variable length_cm is equal to 10000 (100^2) times the variance of the original variable length. The variance of the new variable weight_gram is equal to 1000000 (1000^2) times the variance of the original variable weight.

1.3 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.

str(bldpres2)
#> 'data.frame':    328 obs. of  7 variables:
#>  $ age   : int  28 24 65 32 23 24 39 22 30 31 ...
#>  $ sys   : num  130 118 115 125 128 ...
#>  $ dias  : num  75 72.5 70 87.5 87.5 72.5 80 62.5 87.5 87.5 ...
#>  $ weight: int  49 54 50 61 67 40 44 47 62 78 ...
#>  $ pulse : int  72 80 80 88 76 80 88 68 88 96 ...
#>  $ sex   : chr  "women" "women" "women" "women" ...
#>  $ SES   : chr  "high" "high" "high" "high" ...

bldpres2$sex <- factor(bldpres2$sex)
levels(bldpres2$sex) = c("men", "women")

bldpres2$SES <- factor(bldpres2$SES)
levels(bldpres2$SES) = c("low", "middle", "high")

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

# Barplot for sex
barplot(table(bldpres2$sex), xlab="sex", ylab="Frequency", main="Barplot for sex")

# Barplot for SES
barplot(table(bldpres2$SES), xlab="SES", ylab="Frequency", main="Barplot for SES")

# Stacked barplot showing the info for `sex` and `SES` simultaneously
barplot(table(bldpres2$sex,bldpres2$SES), ylab="Frequency", xlab="SES", legend=c("female","male"),args.legend = list(x = "topleft"))

# Boxplot for systolic blood pressure.
boxplot(bldpres2$sys, ylab="Systolic Blood Pressure", main="Boxplot for systolic blood pressure")

# Boxplots for systolic blood pressure for each category of SES
boxplot(bldpres2$sys~bldpres2$SES, ylab="Systolic Blood Pressure", xlab="SES", main="Boxplots for systolic blood pressure")

Make a frequency table for the variable SES.

freq <- table(bldpres2$SES)
percent <- prop.table(table(bldpres2$SES)) 
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
#> low           63  0.1926606             0.1926606
#> middle       138  0.4220183             0.6146789
#> high         126  0.3853211             1.0000000

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?

# What is the smallest resp. largest observed value for age? 
min(bldpres2$age)
max(bldpres2$age)

# What is the age range ?   
max(bldpres2$age)-min(bldpres2$age)

# What is the median?
median(bldpres2$age)

# What are the values for p25 resp p75 for age ? 
quantile(bldpres2$age, probs = c(0.25, 0.75))

# What is the interquartile range?   
IQR(bldpres2$age)

Look at the distribution of the variable age using a histogram. Is the distribution symmetric or skewed?

hist(bldpres2$age, main="Histogram of age")

The distribution of age is very skewed, with a right tail.

Look at the distribution of the variable age using a boxplot. Are there any outliers?

boxplot(bldpres2$age, main="Boxplot for variable age")

There are no outliers.

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

plot(bldpres2$sys,bldpres2$weight, main="Systolic blood pressure vs Weight", xlab="Systolic blood pressure", ylab="Weight", pch=19, col="blue")

There is no clear (linear) pattern between Weight and systolic blood pressure.

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

cor(bldpres2$sys,bldpres2$weight)
#> [1] 0.1994665

The correlation is relatively small (close to zero), indicating that there is no strong (linear) relationship between the variables Weight and Systolic blood pressure.

1.4 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.

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

grrtx2$sex <- factor(grrtx2$sex)

boxplot(grrtx2$finalht~grrtx2$sex, main="Boxplot of final body length")

which(grrtx2$finalht<100)
#> [1] 50

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?

grrtx2$finalht[50]
#> [1] 73

The value that was entered is equal to 73. Most likely this is a typo and the actual value should have been 173. To fix it, we can manually change the value of this observation:

grrtx2$finalht[50] <- 173

boxplot(grrtx2$finalht~grrtx2$sex, main="Adjusted boxplot of final body length")

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


plot(grrtx2$finalht~grrtx2$tarht, main="Final height vs target height")

range(grrtx2$tarht)
#> [1] 160.5 999.0

There are a couple of values for variable tarht that are very large. If we look into the data, these values are equal to 999. This is sometimes used to denote missing values in a dataset. However, R recognizes this as a value for the variable tarht, rather than missing observations (NA). We can change these observations to be recognized by R as missing observations:


grrtx2$tarht[grrtx2$tarht==999.0] <- NA

plot(grrtx2$finalht~grrtx2$tarht, main="Adjusted plot for final height vs target height")

2 Practical 2: Standardization & Parameter estimation

2.1 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.

set.seed(25)
sample1 <- sample(population, size=50)

Calculate the mean and the variance of the sample.

mean(sample1)
#> [1] 10.30701
var(sample1)
#> [1] 42.47738

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.

set.seed(23)

sample_size <- 50
number_of_samples <- 25

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

sample_variances <- replicate(
  number_of_samples,
  var(sample(population, sample_size)))

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

set.seed(23)

sample_size <- 12
number_of_samples <- 25

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

sample_variances12 <- replicate(
  number_of_samples,
  var(sample(population, sample_size)))

Compare the results of the sampling experiment with sample size 12 or 50, and compare these to the true population mean.

mean(population)
#> [1] 10.98577

sample_means
#>  [1] 11.730389 10.666103 10.586724 10.781876 12.733569 11.146130 12.381214
#>  [8] 11.619996 12.487757  9.915755 10.737076 12.334567 11.144813 11.201867
#> [15] 10.960572 11.284609 11.166399 10.594232 10.839012  9.979151 11.240029
#> [22] 11.976789 12.077147  9.037001 10.303638
sample_means12
#>  [1] 11.473224 13.081854  7.905699 13.546069 11.699813  9.878658 12.716780
#>  [8] 10.155676  7.682037 12.140980 12.322018  8.934542  9.108420  8.622693
#> [15] 12.007219 10.584252 14.937343 12.008024 11.952820 13.154770 12.470660
#> [22] 10.582680 11.041953 11.818626 10.621686

var(population)
#> [1] 35.9019

sample_variances
#>  [1] 35.74572 46.99762 59.67489 32.75696 45.71006 49.13731 37.77806 37.12525
#>  [9] 25.76864 50.87203 41.22843 32.91100 34.34365 38.23959 35.31664 32.22290
#> [17] 29.04526 42.39794 38.06533 43.89108 31.80376 38.12590 29.36895 38.95019
#> [25] 48.34383
sample_variances12
#>  [1] 42.10960 17.62389 17.76187 26.62455 18.87393 19.67393 24.37621 27.27519
#>  [9] 28.80198 24.71522 30.08312 52.07771 28.61768 60.83244 15.44987 49.47518
#> [17] 55.21169 85.88281 19.06387 32.66654 22.75701 33.80871 28.18831 30.07393
#> [25] 37.57994

Check the single experiment where you have the largest error (i.e. the worst estimate) compared to the population mean?

#Absolute value of the difference between the sample means and the population mean for sample size n=50
difference_50 <- abs(sample_means-mean(population))
range(difference_50)
#> [1] 0.0251976 1.9487692

#Absolute value of the difference between the sample means and the population mean for sample size n=12
difference_12 <- abs(sample_means12-mean(population))
range(difference_12)
#> [1] 0.0561831 3.9515729

The largest ‘error’ (i.e. difference between the true population mean and the sample means) is made when taking a sample of only 12 observations.

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).

# To put two figures next to each other we can specify:
par(mfrow=c(1,2))

hist(sample_means, main="Histogram of the sample means with sample size n=50")
hist(sample_means12, main="Histogram of the sample means with sample size n=12")

mean(sample_means)
#> [1] 11.15706
mean(sample_means12)
#> [1] 11.21794

var(sample_means)
#> [1] 0.7799813
var(sample_means12)
#> [1] 3.283267

What is the mean of means for sample size 12 and sample size 50?

mean(sample_means)
#> [1] 11.15706
mean(sample_means12)
#> [1] 11.21794

What is the variance around the mean of means for sample size 12 resp sample size 50?

var(sample_means)
#> [1] 0.7799813
var(sample_means12)
#> [1] 3.283267

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.

mean(sample_means)
#> [1] 11.15706
mean(sample_means12)
#> [1] 11.21794
mean(population)
#> [1] 10.98577

var(sample_means)
#> [1] 0.7799813
var(population)/50
#> [1] 0.7180381

var(sample_means12)
#> [1] 3.283267
var(population)/12
#> [1] 2.991825

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?

A larger sample leads to a smaller variance and thus less uncertainty.

2.2 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, main="Histogram of final scores")

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.

mean(sampling2$finalexam)
#> [1] 10.18657
sd(sampling2$finalexam)
#> [1] 5.995861

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?

mean(sampling2b$meanof30)
#> [1] 10.20303
sd(sampling2b$meanof30)
#> [1] 1.087197

mean(sampling2b$meanof100)
#> [1] 10.15277
sd(sampling2b$meanof100)
#> [1] 0.5899379

mean(sampling2b$meanof500)
#> [1] 10.22395
sd(sampling2b$meanof500)
#> [1] 0.2665706

Make a histogram of meanof30, meanof100 and meanof500. What can you conclude?

# To put three figures next to each other we can specify:
par(mfrow=c(1,3))

hist(sampling2b$meanof30, main="Histogram of Meanof30")
hist(sampling2b$meanof100, main="Histogram of Meanof100")
hist(sampling2b$meanof500, main="Histogram of Meanof500")

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?

freq <- table(sampling2b$check)
percent <- prop.table(table(sampling2b$check)) 
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
#> 0        10       0.05                  0.05
#> 1       190       0.95                  1.00

We find that 5% of the confidence intervals does not contain the true underlying population mean. This is in accordance to our type I error: alpha of 5%.

3 Practical 3: Hypothesis testing

3.1 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.

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

str(sampling3)
#> 'data.frame':    5000 obs. of  3 variables:
#>  $ IQ         : num  123.9 96.2 86.1 120.1 125 ...
#>  $ length     : num  189 182 177 181 179 ...
#>  $ nationality: int  1 1 1 1 1 1 1 1 1 1 ...

# Test for normality of the variables:
shapiro.test(sampling3$length)
#> 
#>  Shapiro-Wilk normality test
#> 
#> data:  sampling3$length
#> W = 0.99974, p-value = 0.831
shapiro.test(sampling3$IQ)
#> 
#>  Shapiro-Wilk normality test
#> 
#> data:  sampling3$IQ
#> W = 0.9995, p-value = 0.2181

par(mfrow=c(2,2))
hist(sampling3$length)
qqnorm(sampling3$length)
qqline(sampling3$length)

hist(sampling3$IQ)
qqnorm(sampling3$IQ)
qqline(sampling3$IQ)

t.test(sampling3$length, mu=185)
#> 
#>  One Sample t-test
#> 
#> data:  sampling3$length
#> t = -65.58, df = 4999, p-value < 2.2e-16
#> alternative hypothesis: true mean is not equal to 185
#> 95 percent confidence interval:
#>  177.0635 177.5243
#> sample estimates:
#> mean of x 
#>  177.2939
t.test(sampling3$IQ, mu=105.51)
#> 
#>  One Sample t-test
#> 
#> data:  sampling3$IQ
#> t = -0.020872, df = 4999, p-value = 0.9833
#> alternative hypothesis: true mean is not equal to 105.51
#> 95 percent confidence interval:
#>  105.0927 105.9185
#> sample estimates:
#> mean of x 
#>  105.5056

The population mean of the length of the Dutch people is 185 cm

The p-value for the Shapiro-Wilk test is larger than 0.05, so we do not reject the null hypothesis of normality. Therefore, the variable length follows a normal distribution. The p-value for the t-test is equal to p-value < 2.2e-16, indicating that we reject the null hypothesis. This means that the the population mean of the length of the Dutch people is significantly different from 185 cm.

The population mean of the IQ of the Dutch people is 105.51

The p-value for the Shapiro-Wilk test is larger than 0.05, so we do not reject the null hypothesis of normality. Therefore, the variable IQ follows a normal distribution. The p-value for the t-test is equal to p-value = 0.9833, indicating that we do not reject the null hypothesis. This means that the the population mean of the IQ of the Dutch people is not significantly different from 105.51

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.


# Initiate a variable to store the p-values in
pvalues_test_length <- NULL

# For loop that goes through the process 15 times:
# 1. set the seed, so your generated results for the sampling are reproducible
# 2. take a random sample of size twelve from the length variable
# 3. Perform the t-test and store the p-value in the variable

for(i in 1:15){
  set.seed(i)
  sample <- sample(sampling3$length, size=12)
  pvalues_test_length[i] <- t.test(sample, mu=185)$p.value
}

# Initiate a variable to store the p-values in
pvalues_test_IQ <- NULL

# For loop that goes through the process 15 times:
# 1. set the seed, so your generated results for the sampling are reproducible
# 2. take a random sample of size twelve from the IQ variable
# 3. Perform the t-test and store the p-value in the variable

for(i in 1:15){
  set.seed(i)
  sample <- sample(sampling3$IQ, size=12)
  pvalues_test_IQ[i] <- t.test(sample, mu=105.51)$p.value
}

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 ?

# Test which p-values are smaller than 0.05
pvalues_test_length < 0.05
#>  [1]  TRUE  TRUE  TRUE  TRUE  TRUE  TRUE FALSE  TRUE  TRUE  TRUE FALSE  TRUE
#> [13]  TRUE  TRUE  TRUE

# Count how many p-values are smaller than 0.05
sum(pvalues_test_length < 0.05)
#> [1] 13

# Test which p-values are smaller than 0.05
pvalues_test_IQ < 0.05
#>  [1] FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE
#> [13] FALSE FALSE FALSE

# Count how many p-values are smaller than 0.05
sum(pvalues_test_IQ < 0.05)
#> [1] 0

For the variable length we find that there is a significant difference in 13 out of 15 random samples. Based on the entire population, we found that the length was significantly different from 185 cm. Therefore we make a correct conclusion in 13 out of the 15 random samples we’ve drawn.

For the variable IQ we find that there is a significant difference in 0 out of 15 random samples. Based on the entire population, we found that the IQ was not significantly different from 105.51. Therefore we make a correct conclusion based on 15 out of the 15 random samples we’ve drawn.

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.

3.2 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.

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.

power.t.test(
  power = 0.90,
  delta = 1.36-1.23,
  sd = sqrt(0.04),
  sig.level = 0.05,
  type = "one.sample",
  alternative = "one.sided"
)
#> 
#>      One-sample t test power calculation 
#> 
#>               n = 21.69627
#>           delta = 0.13
#>              sd = 0.2
#>       sig.level = 0.05
#>           power = 0.9
#>     alternative = one.sided

4 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.

scopinaro <- read.table("scopinaro.txt", header=T)
str(scopinaro)
#> 'data.frame':    344 obs. of  24 variables:
#>  $ surgery : int  0 0 0 0 0 0 0 0 0 0 ...
#>  $ sex     : int  0 0 0 0 0 0 0 0 0 0 ...
#>  $ age     : int  26 26 53 18 37 43 36 43 41 24 ...
#>  $ weight  : num  100 71 97 109 100 94 105 107 103 93 ...
#>  $ length  : num  1.57 1.48 1.6 1.68 1.63 1.56 1.64 1.57 1.6 1.69 ...
#>  $ BMIpre  : num  40.6 32.4 37.9 38.6 37.6 ...
#>  $ ideal   : num  54.2 48.2 56.3 62.1 58.5 ...
#>  $ excess  : num  45.8 22.8 40.7 46.9 41.5 ...
#>  $ current : num  76.5 65 NA 88 73 85 66 69 96 78 ...
#>  $ weightlo: num  23.5 6 NA 21 27 9 39 38 7 15 ...
#>  $ excess1 : num  51.3 26.3 NA 44.8 65 ...
#>  $ excess2 : num  23.5 8.45 NA 19.27 27 ...
#>  $ BMIpost : num  31 29.7 NA 31.2 27.5 ...
#>  $ chir    : int  0 0 0 0 0 0 0 0 0 0 ...
#>  $ complic : int  0 0 0 0 0 0 1 1 1 0 ...
#>  $ comorb  : int  0 0 0 1 0 0 1 0 0 0 ...
#>  $ diabetes: int  0 0 0 1 0 0 0 0 0 0 ...
#>  $ hypercho: int  0 0 0 0 0 0 0 0 0 0 ...
#>  $ aht     : int  0 0 0 0 0 0 0 0 0 0 ...
#>  $ copd    : int  0 0 0 0 0 0 0 0 0 0 ...
#>  $ sleepapn: int  0 0 0 0 0 0 0 0 0 0 ...
#>  $ dispneu : int  0 0 0 0 0 0 0 0 0 0 ...
#>  $ locomot : int  0 0 0 0 0 0 1 0 0 0 ...
#>  $ depressi: int  0 0 0 0 0 0 0 0 0 0 ...
head(scopinaro)
#>   surgery sex age weight length   BMIpre   ideal  excess current weightlo
#> 1       0   0  26    100   1.57 40.56960 54.2278 45.7722    76.5     23.5
#> 2       0   0  26     71   1.48 32.41417 48.1888 22.8112    65.0      6.0
#> 3       0   0  53     97   1.60 37.89062 56.3200 40.6800      NA       NA
#> 4       0   0  18    109   1.68 38.61961 62.0928 46.9072    88.0     21.0
#> 5       0   0  37    100   1.63 37.63785 58.4518 41.5482    73.0     27.0
#> 6       0   0  43     94   1.56 38.62590 53.5392 40.4608    85.0      9.0
#>    excess1   excess2  BMIpost chir complic comorb diabetes hypercho aht copd
#> 1 51.34121 23.500000 31.03574    0       0      0        0        0   0    0
#> 2 26.30287  8.450704 29.67495    0       0      0        0        0   0    0
#> 3       NA        NA       NA    0       0      0        0        0   0    0
#> 4 44.76925 19.266055 31.17914    0       0      1        1        0   0    0
#> 5 64.98476 27.000000 27.47563    0       0      0        0        0   0    0
#> 6 22.24375  9.574468 34.92768    0       0      0        0        0   0    0
#>   sleepapn dispneu locomot depressi
#> 1        0       0       0        0
#> 2        0       0       0        0
#> 3        0       0       0        0
#> 4        0       0       0        0
#> 5        0       0       0        0
#> 6        0       0       0        0

4.1 Question 1

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


# Check normality of BMIpre for females
shapiro.test(scopinaro$BMIpre[scopinaro$sex==0])
#> 
#>  Shapiro-Wilk normality test
#> 
#> data:  scopinaro$BMIpre[scopinaro$sex == 0]
#> W = 0.92408, p-value = 1.761e-10

par(mfrow=c(1,2))
hist(scopinaro$BMIpre[scopinaro$sex==0])
qqnorm(scopinaro$BMIpre[scopinaro$sex==0])
qqline(scopinaro$BMIpre[scopinaro$sex==0])

# Normality of BMIpre for females is not fulfilled

# Check normality of BMIpre for males
shapiro.test(scopinaro$BMIpre[scopinaro$sex==1])
#> 
#>  Shapiro-Wilk normality test
#> 
#> data:  scopinaro$BMIpre[scopinaro$sex == 1]
#> W = 0.96435, p-value = 0.06491

par(mfrow=c(1,2))
hist(scopinaro$BMIpre[scopinaro$sex==1])
qqnorm(scopinaro$BMIpre[scopinaro$sex==1])
qqline(scopinaro$BMIpre[scopinaro$sex==1])

# Normality of BMIpre for males is fulfilled

# We have to use a nonparametric test

wilcox.test(BMIpre~factor(sex), data=scopinaro)
#> 
#>  Wilcoxon rank sum test with continuity correction
#> 
#> data:  BMIpre by factor(sex)
#> W = 6747, p-value = 0.01184
#> alternative hypothesis: true location shift is not equal to 0

The p-value for the Wilcoxon Rank test is equal to 0.01184, indicating that we reject the null hypothesis of equal distributions. Therefore we can conclude that there is a significant difference in the BMI before the operation (BMIpre) between males and females.

4.2 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?

# Create a dataset that only contains the male observations
scopinaro_males <- scopinaro[scopinaro$sex==1,]

# Create a variable that calculates the difference in BMI post and pre intervention
scopinaro_males$DIFF <- scopinaro_males$BMIpost-scopinaro_males$BMIpre

# Check normality of the difference 
shapiro.test(scopinaro_males$DIFF)
#> 
#>  Shapiro-Wilk normality test
#> 
#> data:  scopinaro_males$DIFF
#> W = 0.90299, p-value = 0.001771

par(mfrow=c(1,2))
hist(scopinaro_males$DIFF)
qqnorm(scopinaro_males$DIFF)
qqline(scopinaro_males$DIFF)


# Normality of DIFF for males is not fulfilled
# We have to use a nonparametric test

wilcox.test(scopinaro_males$BMIpost,scopinaro_males$BMIpre, paired=T, conf.int = T)
#> 
#>  Wilcoxon signed rank exact test
#> 
#> data:  scopinaro_males$BMIpost and scopinaro_males$BMIpre
#> V = 0, p-value = 4.547e-13
#> alternative hypothesis: true location shift is not equal to 0
#> 95 percent confidence interval:
#>  -9.120525 -6.346310
#> sample estimates:
#> (pseudo)median 
#>      -7.526576

The p-value for the Wilcoxon test is equal to p-value = 4.547e-13, indicating that we reject the null hypothesis. There is a significant difference in BMI post and pre intervention for males. The median difference is equal to -7.5, indicating that the median drop in BMI after intervention is equal to 7.5 units. The 95% confidence interval is equal to [-9.120525 ; -6.346310].

# Create a dataset that only contains the male observations
scopinaro_females <- scopinaro[scopinaro$sex==0,]

# Create a variable that calculates the difference in BMI post and pre intervention
scopinaro_females$DIFF <- scopinaro_females$BMIpost-scopinaro_females$BMIpre

# Check normality of the difference 
shapiro.test(scopinaro_females$DIFF)
#> 
#>  Shapiro-Wilk normality test
#> 
#> data:  scopinaro_females$DIFF
#> W = 0.169, p-value < 2.2e-16

par(mfrow=c(1,2))
hist(scopinaro_females$DIFF)
qqnorm(scopinaro_females$DIFF)
qqline(scopinaro_females$DIFF)


# Normality of DIFF for females is not fulfilled
# We have to use a nonparametric test

wilcox.test(scopinaro_females$BMIpost,scopinaro_females$BMIpre, paired=T, conf.int = T)
#> 
#>  Wilcoxon signed rank test with continuity correction
#> 
#> data:  scopinaro_females$BMIpost and scopinaro_females$BMIpre
#> V = 236, p-value < 2.2e-16
#> alternative hypothesis: true location shift is not equal to 0
#> 95 percent confidence interval:
#>  -7.672179 -6.521152
#> sample estimates:
#> (pseudo)median 
#>      -7.087582

The p-value for the Wilcoxon test is equal to p-value < 2.2e-16, indicating that we reject the null hypothesis. There is a significant difference in BMI post and pre intervention for females. The median difference is equal to -7.1, indicating that the median drop in BMI after intervention is equal to 7.1 units. The 95% confidence interval is equal to [-7.672179 ; -6.521152].

4.3 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.

# Create a variable that calculates the difference in BMI post and pre intervention
scopinaro_males$BMIdiff <- scopinaro_males$BMIpost-scopinaro_males$BMIpre
scopinaro_females$BMIdiff <- scopinaro_females$BMIpost-scopinaro_females$BMIpre

# Check normality of the difference 
shapiro.test(scopinaro_males$BMIdiff) # Normality not fulfilled
#> 
#>  Shapiro-Wilk normality test
#> 
#> data:  scopinaro_males$BMIdiff
#> W = 0.90299, p-value = 0.001771
shapiro.test(scopinaro_females$BMIdiff) # Normality not fulfilled
#> 
#>  Shapiro-Wilk normality test
#> 
#> data:  scopinaro_females$BMIdiff
#> W = 0.169, p-value < 2.2e-16

# We have to use a nonparametric test

wilcox.test(scopinaro_males$BMIdiff)
#> 
#>  Wilcoxon signed rank exact test
#> 
#> data:  scopinaro_males$BMIdiff
#> V = 0, p-value = 4.547e-13
#> alternative hypothesis: true location is not equal to 0
wilcox.test(scopinaro_females$BMIdiff)
#> 
#>  Wilcoxon signed rank test with continuity correction
#> 
#> data:  scopinaro_females$BMIdiff
#> V = 236, p-value < 2.2e-16
#> alternative hypothesis: true location is not equal to 0

4.4 Question 4

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

# Create a variable that calculates the difference in BMI post and pre intervention
scopinaro$BMIdiff <- scopinaro$BMIpost-scopinaro$BMIpre

# Check normality of the difference for the two surgical techniques
shapiro.test(scopinaro$BMIdiff[scopinaro$surgery==0]) # Normality not fulfilled
#> 
#>  Shapiro-Wilk normality test
#> 
#> data:  scopinaro$BMIdiff[scopinaro$surgery == 0]
#> W = 0.21321, p-value = 1.873e-14
shapiro.test(scopinaro$BMIdiff[scopinaro$surgery==1]) # Normality not fulfilled
#> 
#>  Shapiro-Wilk normality test
#> 
#> data:  scopinaro$BMIdiff[scopinaro$surgery == 1]
#> W = 0.9485, p-value = 2.759e-07

# We have to use a nonparametric test

wilcox.test(BMIdiff~surgery,data=scopinaro)
#> 
#>  Wilcoxon rank sum test with continuity correction
#> 
#> data:  BMIdiff by surgery
#> W = 6033, p-value = 0.2099
#> alternative hypothesis: true location shift is not equal to 0

The p-value for the Wilcoxon test is equal to p-value = 0.2099, indicating that there is no significant difference in the change in BMI between the two types of surgery.

4.5 Question 5

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

# Create a subset of the dataset that only contains observations of diabetic patients
diabetic <- scopinaro[scopinaro$diabetes==1,]

# Check normality of the variable BMIpre for males and females
shapiro.test(diabetic$BMIpre[diabetic$sex==0]) # Normality fulfilled
#> 
#>  Shapiro-Wilk normality test
#> 
#> data:  diabetic$BMIpre[diabetic$sex == 0]
#> W = 0.90039, p-value = 0.05836
shapiro.test(diabetic$BMIpre[diabetic$sex==1])# Normality fulfilled
#> 
#>  Shapiro-Wilk normality test
#> 
#> data:  diabetic$BMIpre[diabetic$sex == 1]
#> W = 0.86359, p-value = 0.2775

# We can do a parametric t-test

t.test(diabetic$BMIpre~factor(diabetic$sex))
#> 
#>  Welch Two Sample t-test
#> 
#> data:  diabetic$BMIpre by factor(diabetic$sex)
#> t = -1.3935, df = 2.4573, p-value = 0.2761
#> alternative hypothesis: true difference in means between group 0 and group 1 is not equal to 0
#> 95 percent confidence interval:
#>  -28.94620  12.85365
#> sample estimates:
#> mean in group 0 mean in group 1 
#>        41.98710        50.03337

The p-value for the t-test is equal to p-value = 0.2761, indicating that we do not reject the null hypothesis. There is no proof in the data that there is a significant difference in BMIpre for male and female diabetic patients.

5 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!!

5.1 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.

data$gender <- factor(data$gender)
data$group <- factor(data$group)

group1 <- subset(data, data$group == 1)
group3 <- subset(data, data$group == 3)
group4 <- subset(data, data$group == 4)

# normality test
shapiro.test(group1$mpt) # Normality not fulfilled
#> 
#>  Shapiro-Wilk normality test
#> 
#> data:  group1$mpt
#> W = 0.9542, p-value = 0.002061
shapiro.test(group3$mpt) # Normality not fulfilled
#> 
#>  Shapiro-Wilk normality test
#> 
#> data:  group3$mpt
#> W = 0.94065, p-value = 0.0002112
shapiro.test(group4$mpt) # Normality not fulfilled
#> 
#>  Shapiro-Wilk normality test
#> 
#> data:  group4$mpt
#> W = 0.96942, p-value = 0.02659

# We need to perform a non-parametric test

kruskal.test(mpt ~ group, data = data)
#> 
#>  Kruskal-Wallis rank sum test
#> 
#> data:  mpt by group
#> Kruskal-Wallis chi-squared = 48.531, df = 2, p-value = 2.894e-11

The p-value = 2.894e-11 indicating that we reject the null hypothesis. The maximal phoneme time mpt is different between the 3 groups of children.

5.2 Question 2

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

# normality test
shapiro.test(group1$jitter)  # Normality not fulfilled
#> 
#>  Shapiro-Wilk normality test
#> 
#> data:  group1$jitter
#> W = 0.69523, p-value = 6.856e-13
shapiro.test(group3$jitter)  # Normality not fulfilled
#> 
#>  Shapiro-Wilk normality test
#> 
#> data:  group3$jitter
#> W = 0.66325, p-value = 8.127e-14
shapiro.test(group4$jitter)  # Normality not fulfilled
#> 
#>  Shapiro-Wilk normality test
#> 
#> data:  group4$jitter
#> W = 0.7394, p-value = 1.273e-11

# We need to perform a non-parametric test
kruskal.test(jitter ~ group, data = data)
#> 
#>  Kruskal-Wallis rank sum test
#> 
#> data:  jitter by group
#> Kruskal-Wallis chi-squared = 6.842, df = 2, p-value = 0.03268

The p-value = 0.03268 indicating that we reject the null hypothesis. The jitter is different between the 3 groups of children.

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

library(FSA)
#> Warning: package 'FSA' was built under R version 4.5.3
#> ## FSA v0.10.1. See citation('FSA') if used in publication.
#> ## Run fishR() for related website and fishR('IFAR') for related book.

dunnTest(jitter ~ group, data = data, method="bh")
#>   Kruskal-Wallis rank sum test
#> 
#> data: x and g
#> Kruskal-Wallis chi-squared = 6.842, df = 2, p-value = 0.03
#> 
#>                       Dunn's Pairwise Comparison of x by g                      
#>                               (Benjamini-Hochberg)                              
#> 
#> Col Mean-│
#> Row Mean │          1          3
#> ─────────┼──────────────────────
#>        3 │  -2.097045
#>          │     0.0540 
#>          │
#>        4 │  -2.402305  -0.339977
#>          │     0.0489*    0.7339 
#> 
#> FDR = 0.05
#> Reject Ho if adjusted p ≤ FDR with stopping rule, where (unadjusted) p = Pr(|Z| ≥ |z|)
#> Dunn (1964) Kruskal-Wallis multiple comparison
#> 
#>   p-values adjusted with the Benjamini-Hochberg method.
#>   Comparison          Z    P.unadj      P.adj
#> 1      1 - 3 -2.0970453 0.03598956 0.05398435
#> 2      1 - 4 -2.4023059 0.01629208 0.04887624
#> 3      3 - 4 -0.3399766 0.73387417 0.73387417

There is only a significant difference between group 1 and group 4.

5.3 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.

library(car)
#> Loading required package: carData
#> Registered S3 methods overwritten by 'car':
#>   method       from
#>   hist.boot    FSA 
#>   confint.boot FSA
#> 
#> Attaching package: 'car'
#> The following object is masked from 'package:FSA':
#> 
#>     bootCase

# Test for equality of variances
leveneTest(voicedata$mpt~factor(voicedata$gender))
#> Levene's Test for Homogeneity of Variance (center = median)
#>        Df F value  Pr(>F)  
#> group   1  3.1476 0.07634 .
#>       987                  
#> ---
#> Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

# Parametric ANOVA

model <- aov(mpt ~ factor(gender), data = voicedata)

summary(model)
#>                 Df Sum Sq Mean Sq F value Pr(>F)  
#> factor(gender)   1    102  102.46   3.903 0.0485 *
#> Residuals      987  25912   26.25                 
#> ---
#> Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
#> 1 observation deleted due to missingness

Based on the Levene’s test, we can conclude that the variances of mpt are equal in the boys and girls group.

The ANOVA indicates that there is a significant difference in the mean mpt value for girls and boys.

Not that there are only two groups here. We could have also performed a t-test, which is equivalent:

t.test(mpt ~ factor(gender), data = voicedata, var.equal=T)
#> 
#>  Two Sample t-test
#> 
#> data:  mpt by factor(gender)
#> t = -1.9756, df = 987, p-value = 0.04848
#> alternative hypothesis: true difference in means between group 0 and group 1 is not equal to 0
#> 95 percent confidence interval:
#>  -1.283383805 -0.004312568
#> sample estimates:
#> mean in group 0 mean in group 1 
#>        11.12379        11.76763

5.4 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.

model <- aov(voicedata$f0~factor(voicedata$group))
summary(model)
#>                          Df Sum Sq Mean Sq F value  Pr(>F)    
#> factor(voicedata$group)   5  12530  2506.0   4.281 0.00074 ***
#> Residuals               973 569540   585.3                    
#> ---
#> Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
#> 11 observations deleted due to missingness

pairwise.t.test(voicedata$f0,voicedata$group, p.adjust.method = "bonferroni")
#> 
#>  Pairwise comparisons using t tests with pooled SD 
#> 
#> data:  voicedata$f0 and voicedata$group 
#> 
#>   0       1       2       3       4      
#> 1 0.00036 -       -       -       -      
#> 2 1.00000 1.00000 -       -       -      
#> 3 1.00000 0.75191 1.00000 -       -      
#> 4 1.00000 0.10642 1.00000 1.00000 -      
#> 5 1.00000 0.00727 1.00000 1.00000 1.00000
#> 
#> P value adjustment method: bonferroni
pairwise.t.test(voicedata$f0,voicedata$group, p.adjust.method = "BH")
#> 
#>  Pairwise comparisons using t tests with pooled SD 
#> 
#> data:  voicedata$f0 and voicedata$group 
#> 
#>   0       1       2       3       4      
#> 1 0.00036 -       -       -       -      
#> 2 0.29143 0.48022 -       -       -      
#> 3 0.25099 0.18798 0.77339 -       -      
#> 4 0.57747 0.03547 0.54161 0.57747 -      
#> 5 0.87276 0.00363 0.29143 0.29143 0.57747
#> 
#> P value adjustment method: BH
pairwise.t.test(voicedata$f0,voicedata$group, p.adjust.method = "holm")
#> 
#>  Pairwise comparisons using t tests with pooled SD 
#> 
#> data:  voicedata$f0 and voicedata$group 
#> 
#>   0       1       2       3       4      
#> 1 0.00036 -       -       -       -      
#> 2 1.00000 1.00000 -       -       -      
#> 3 0.92029 0.60153 1.00000 -       -      
#> 4 1.00000 0.09223 1.00000 1.00000 -      
#> 5 1.00000 0.00678 1.00000 1.00000 1.00000
#> 
#> P value adjustment method: holm

TukeyHSD(model)
#>   Tukey multiple comparisons of means
#>     95% family-wise confidence level
#> 
#> Fit: aov(formula = voicedata$f0 ~ factor(voicedata$group))
#> 
#> $`factor(voicedata$group)`
#>            diff         lwr       upr     p adj
#> 1-0 -11.3232254 -18.9415630 -3.704888 0.0003456
#> 2-0  -6.2529574 -18.3026583  5.796743 0.6760912
#> 3-0  -4.5608374 -12.0816775  2.960003 0.5110869
#> 4-0  -1.8754255  -9.5961324  5.845281 0.9826721
#> 5-0   0.4009364  -6.7455691  7.547442 0.9999854
#> 2-1   5.0702680  -8.5517146 18.692251 0.8959665
#> 3-1   6.7623880  -3.0827201 16.607496 0.3655546
#> 4-1   9.4478000  -0.5508217 19.446422 0.0764953
#> 5-1  11.7241618   2.1619623 21.286361 0.0064257
#> 3-2   1.6921200 -11.8755761 15.259816 0.9992497
#> 4-2   4.3775319  -9.3019658 18.057030 0.9430476
#> 5-2   6.6538938  -6.7099336 20.017721 0.7138174
#> 4-3   2.6854119  -7.2391234 12.609947 0.9720739
#> 5-3   4.9617738  -4.5229311 14.446479 0.6683955
#> 5-4   2.2763619  -7.3675952 11.920319 0.9847682

The ANOVA indicates that there is a significant difference between the groups (p-value=0.00074)

Based on the Bonferroni correction:

There is a significant difference between group 0 and 1, and between group 1 and 5.

Based on the Benjamini-Hochberg correction:

There is a significant difference between group 0 and 1, between group 1 and 4, and between group 1 and 5.

Based on the Holm correction:

There is a significant difference between group 0 and 1, and between group 1 and 5.

Based on Tukey correction:

There is a significant difference between group 0 and 1, and between group 1 and 5.

6 Practical 6: Association between categorical variables

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

6.1 Question 1

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

scopinaro$comorb <- as.factor(scopinaro$comorb)
scopinaro$complic <- as.factor(scopinaro$complic)

test <- chisq.test(table(scopinaro$comorb, scopinaro$complic))
test$expected # All expected counts are larger than 5
#>    
#>             0         1
#>   0  55.09145  60.90855
#>   1 105.90855 117.09145

test
#> 
#>  Pearson's Chi-squared test with Yates' continuity correction
#> 
#> data:  table(scopinaro$comorb, scopinaro$complic)
#> X-squared = 2.8844, df = 1, p-value = 0.08944

The p-value = 0.08944 indicates that we do not reject the null hypothesis. The variables comorbidity (variable comorb) and the occurrence of complications (variable complic) are independent (i.e. there is no relationship between the variables).

6.2 Question 2

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

scopinaro$surgery <- as.factor(scopinaro$surgery)

test2 <- chisq.test(table(scopinaro$surgery, scopinaro$complic))
test2$expected # All expected counts are larger than 5
#>    
#>             0         1
#>   0  30.87021  34.12979
#>   1 130.12979 143.87021

test2
#> 
#>  Pearson's Chi-squared test with Yates' continuity correction
#> 
#> data:  table(scopinaro$surgery, scopinaro$complic)
#> X-squared = 3.355, df = 1, p-value = 0.067

The p-value = 0.067 indicates that we do not reject the null hypothesis. The occurrence of complications does not differ between the two operation techniques.

6.3 Question 3

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

scopinaro$sex <- as.factor(scopinaro$sex)

test3 <- chisq.test(table(scopinaro$sex, scopinaro$complic))
test3$expected # All expected counts are larger than 5
#>    
#>             0         1
#>   0 129.65487 143.34513
#>   1  31.34513  34.65487

test3
#> 
#>  Pearson's Chi-squared test with Yates' continuity correction
#> 
#> data:  table(scopinaro$sex, scopinaro$complic)
#> X-squared = 0.10063, df = 1, p-value = 0.7511

The p-value = 0.7511 indicates that we do not reject the null hypothesis. There is no relationship between the occurence of complications and sex.

6.4 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.

test4 <- chisq.test(tab)
#> Warning in chisq.test(tab): Chi-squared approximation may be incorrect
test4$expected # Not all expected counts are larger than 5
#>              mobile_yes mobile_no
#> facebook_yes          9         1
#> facebook_no           9         1

fisher.test(tab)
#> 
#>  Fisher's Exact Test for Count Data
#> 
#> data:  tab
#> p-value = 0.4737
#> alternative hypothesis: true odds ratio is not equal to 1
#> 95 percent confidence interval:
#>  0.1911327       Inf
#> sample estimates:
#> odds ratio 
#>        Inf

The p-value = 0.4737 indicates that we do not reject the null hypothesis. There is no relationship between having a mobile phone and having a Facebook account.

7 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.

crimedata <- read.table("crimedata.txt", header=T)
head(crimedata)
#>   countyid education crimerate income
#> 1        1        74      8487      1
#> 2        3        81      8362      1
#> 3        6        66      9100      1
#> 4        8        81      5873      1
#> 5        9        74      7993      1
#> 6       14        84      4595      1
str(crimedata)
#> 'data.frame':    84 obs. of  4 variables:
#>  $ countyid : int  1 3 6 8 9 14 16 18 19 20 ...
#>  $ education: int  74 81 66 81 74 84 79 73 77 65 ...
#>  $ crimerate: int  8487 8362 9100 5873 7993 4595 4427 10768 8335 12311 ...
#>  $ income   : int  1 1 1 1 1 1 1 1 1 1 ...

crimedata$income <- factor(crimedata$income)

7.1 Graphical analysis

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


model <- lm(crimerate ~ education, data = crimedata)

plot(crimerate ~ education, data = crimedata)
abline(model, col="red")

7.2 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.


model <- lm(crimerate ~ education, data = crimedata)

summary(model)
#> 
#> Call:
#> lm(formula = crimerate ~ education, data = crimedata)
#> 
#> Residuals:
#>     Min      1Q  Median      3Q     Max 
#> -5278.3 -1757.5  -210.5  1575.3  6803.3 
#> 
#> Coefficients:
#>             Estimate Std. Error t value Pr(>|t|)    
#> (Intercept) 20517.60    3277.64   6.260 1.67e-08 ***
#> education    -170.58      41.57  -4.103 9.57e-05 ***
#> ---
#> Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
#> 
#> Residual standard error: 2356 on 82 degrees of freedom
#> Multiple R-squared:  0.1703, Adjusted R-squared:  0.1602 
#> F-statistic: 16.83 on 1 and 82 DF,  p-value: 9.571e-05

confint(model)
#>                  2.5 %      97.5 %
#> (Intercept) 13997.3245 27037.87538
#> education    -253.2798   -87.87061

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?

The p-value for the slope is equal to 9.57e-05, indicating that we reject the null hypothesis. This means that there is a significant effect of the education level on the crimerate. The equation of the regression line is: crimerate=20517.50-170.58 x education. On average the crimerate goes down by 170.58 units when education level increases by one unit. The 95% confidence interval for the effect of education is equal to [-253.2798 ; -87.87061]. The percentage of variability in the crimerate that is explained by the education percentage is equal to 17%.

7.3 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.


# Residual plot to check homoscedasticity and linearity

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

abline(h = 0)

# Checking normality

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

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

Linearity: There is no clear pattern in the residual plot. The points are randomly scattered around 0. Linearity is ok.
Homoscedasticity: There is a slight decrease in the variability for larger values of the response, i.e. the points are less widely scattered around 0, however the trend is not very strong. Therefore, constant variance can be assumed ok.
Normality: The p-value of the Shapiro-Wilk test indicates that the residuals do not deviate from normality. The QQ-plot also shows that normality is ok.

7.4 Prediction

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

new_data <- data.frame(education = 50)

predict(
  model,
  newdata = new_data,
  interval = "prediction",
  level = 0.99
)
#>        fit      lwr     upr
#> 1 11988.84 4995.977 18981.7

The predicted crimerate for a county where the education level is equal to 50%, is equal to 11988.84. The 99% prediction interval is equal to [4995.977 ; 18981.7].

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

new_data <- data.frame(
  education = seq(
    min(crimedata$education, na.rm = TRUE),
    max(crimedata$education, na.rm = TRUE),
    length.out = 100
  )
)

confidence <- predict(
  model,
  newdata = new_data,
  interval = "prediction",
  level = 0.99
)

plot(
  crimedata$education,
  crimedata$crimerate,
  main = "Regression with 99% confidence interval",
  xlab = "Education",
  ylab = "Crime Rate",
  pch = 19
)

# Draw the shaded confidence interval
polygon(
  c(new_data$education, rev(new_data$education)),
  c(confidence[, "lwr"], rev(confidence[, "upr"])),
  col = "lightblue",
  border = NA
)

# Add the observed data and regression line on top
points(crimedata$education, crimedata$crimerate, pch = 19)
lines(new_data$education, confidence[, "fit"], lwd = 2)

Basic Principles of Statistics - Practical Session Solutions

dr. Annelies Agten

2026-10-04

1 Practical 1: Descriptive statistics

1.1 Exercise 1

1.2 Exercise 2

1.3 Exercise 3

1.4 Exercise 4

2 Practical 2: Standardization & Parameter estimation

2.1 Exercise 1

2.2 Exercise 2

3 Practical 3: Hypothesis testing

3.1 Exercise 1

3.2 Exercise 2

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

4.1 Question 1

4.2 Question 2

4.3 Question 3

4.4 Question 4

4.5 Question 5

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

5.1 Question 1

5.2 Question 2

5.3 Question 3

5.4 Question 4

6 Practical 6: Association between categorical variables

6.1 Question 1

6.2 Question 2

6.3 Question 3

6.4 Question 4

7 Practical 7: Simple linear regression

7.1 Graphical analysis

7.2 Fitting the model

7.3 Analysis of residuals

7.4 Prediction