Statistical Inference Course Project

- A Simulation Exercise and A Basic Inferential Data Analysis

Overview

This project consists of 2 parts:

A simulation exercise
A basic inference data analysis

For the simulation exercise, we will investigate the exponential distribution in R and compare it with the Central Limit Theorem. The Sample Mean versus Theoretical Mean and Sample Variance versus Theoretical Variance will be discussed. In addition, we will show that the distribution is approximately normal.

For the basic inference data analysis, we will analyse the ToothGrowth data in the R datasets package. We will use hypothesis tests,confidence intervals and p-values to compare tooth growth by supplement types and dosage.

PART 1: A SIMULATION EXERCISE

For this simulation exercise, we will simulate the exponential distribution in R with rexp(n, lambda) where lambda is the rate parameter. We will set lambda = 0.2 for all of the simulations, and investigate the distribution of averages of 40 exponentials. We will do a thousand simulations.

Simulations

First, let’s simulate some random exponentials for the simulation exercise.

## Setting the seed of the simulations to be 123456 so that the simulations can be reproduced.
set.seed(123456)
n = 40        ## number of exponentials = 40
lambda = 0.2  ## set lambda = 0.2
nosim = 1000  ## number of simulations = 1000

## Generating 1000 simulations, based on the averages of 40 exponentials
for (i in 1:1000) {
    temp <- matrix(sample(rexp(n, lambda), n*nosim, replace = TRUE), ncol = n)
    temp <- apply(temp, 1, mean)  ## finding the mean of each 40 exponentials
}

Sample Mean versus Theoretical Mean

Let’s calculate the sample mean and compare it with the theoretical mean of the distribution.

The mean of exponential distribution is 1/lambda. As lambda = 0.2, the theoretical mean of the distribution is 1/0.2 = 5.

sample_mean <- mean(temp)
sample_mean

[1] 4.949389

theoretical_mean <- 1/lambda
theoretical_mean

[1] 5

From the calculations, we can see that the sample mean is very close to the theoretical mean. Let’s plot a histogram to look at the distribution. The red vertical line represents the sample mean while the blue vertical line represents the theoretical mean.

hist(temp, main = "1000 simulations based on the averages of 40 exponentials", xlab = "Average of 40 exponentials")
abline(v = sample_mean, col = "red")
abline(v = theoretical_mean, col = "blue")

From the histogram, we can see that the sample mean of 4.949 is very close to the theoretical mean of 5.

Sample Variance versus Theoretical Variance

Next, let’s calculate the sample variance and compare it with the theoretical variance of the distribution.

The standard deviation of the exponential distribution is equal to (1/lambda)/sqrt(n). Square the standard deviation to calculate the variance.

sample_variance <- sd(temp)^2
sample_variance

[1] 0.6161016

theoretical_variance <- ((1/lambda)/sqrt(n))^2
theoretical_variance

[1] 0.625

From the calculations, we can see that the sample variance of 0.616 is very close to the theoretical variance of 0.625.

Distribution

Finally, we will investigate if the distribution is approximately normal. Based on Central Limit Theorem, the distribution of averages of iid variables, properly normalised, becomes that of a standard normal as the sample size increases. Then $\frac{\bar X_n - \mu}{\sigma / \sqrt{n}}$ has a distribution like that of a standard normal for large n.

temp <- (temp - 5)/(sqrt(0.625) / sqrt(n))  ## μ = 5, σ = sqrt(0.625)
dty <- density(temp)
plot(dty$x, dty$y, xlab = "", ylab = "density", type = "n")
title(paste("Sample mean of", n, "observations"))
lines(dty$x, dty$y, lwd = 2)

Based on the plot, the distribution of the sample mean of 40 observations appears to follow a normal distribution. However, if we increase the number of observations (n), the distribution would align more closely to a normal distribution.

PART 2: A BASIC INFERENTIAL DATA ANALYSIS

For the basic inference data analysis, we will analyse the ToothGrowth data in the R datasets package.

Load datasets

First, let’s load the ToothGrowth data in the R datasets package and perform some basic exploratory data analyses by using the “head”, “tail” and “str” functions. These functions perform the following exploratory data analyses:

“head”: to look at the top 6 rows of the dataset
“tail”: to look at the bottom 6 rows of the dataset
“str”: to look at the structure of the dataset

## Load the dataset "ToothGrowth"
data("ToothGrowth")

library(xtable)

## Look at the top 6 rows of the dataset
head <- xtable(head(ToothGrowth))
print(head, type = "html")

	len	supp	dose
1	4.20	VC	0.50
2	11.50	VC	0.50
3	7.30	VC	0.50
4	5.80	VC	0.50
5	6.40	VC	0.50
6	10.00	VC	0.50

## Look at the bottom 6 rows of the dataset
tail <- xtable(tail(ToothGrowth))
print(tail, type = "html")

	len	supp	dose
55	24.80	OJ	2.00
56	30.90	OJ	2.00
57	26.40	OJ	2.00
58	27.30	OJ	2.00
59	29.40	OJ	2.00
60	23.00	OJ	2.00

## Look at the structure of the dataset
str(ToothGrowth)

‘data.frame’: 60 obs. of 3 variables: $ len : num 4.2 11.5 7.3 5.8 6.4 10 11.2 11.2 5.2 7 … $ supp: Factor w/ 2 levels “OJ”,“VC”: 2 2 2 2 2 2 2 2 2 2 … $ dose: num 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 …

The dataset has 60 rows and 3 columns.

Basic summary of the data

Let’s look at some basic summary of the data.

summary <- xtable(summary(ToothGrowth))
print(summary, type = "html")

	len	supp	dose
X	Min. : 4.20	OJ:30	Min. :0.500
X.1	1st Qu.:13.07	VC:30	1st Qu.:0.500
X.2	Median :19.25		Median :1.000
X.3	Mean :18.81		Mean :1.167
X.4	3rd Qu.:25.27		3rd Qu.:2.000
X.5	Max. :33.90		Max. :2.000

unique(ToothGrowth$dose)

[1] 0.5 1.0 2.0

As there are 3 unique values for dose, I will convert it to a factor type for easier analysis.

ToothGrowth$dose <- as.factor(ToothGrowth$dose)

str(ToothGrowth)

‘data.frame’: 60 obs. of 3 variables: $ len : num 4.2 11.5 7.3 5.8 6.4 10 11.2 11.2 5.2 7 … $ supp: Factor w/ 2 levels “OJ”,“VC”: 2 2 2 2 2 2 2 2 2 2 … $ dose: Factor w/ 3 levels “0.5”,“1”,“2”: 1 1 1 1 1 1 1 1 1 1 …

Let’s plot a boxplot to look at the length of tooth growth based on supp and dose.

library(ggplot2)

ggplot(ToothGrowth, aes(x = dose, y = len, fill = supp)) +
    geom_boxplot() +
    ggtitle("Length of tooth growth based on supp and dose") +
    theme_classic()

From the boxplot, we can observe the following:

As dosage increases, the length of tooth growth increases as well.
There is a larger differences in length between supplement OJ and VC for doses 0.5mg/day and 1mg/day, as compared to dose 2mg/day.
Supplement OJ is more effective for doses 0.5mg/day and 1mg/day, as compared to supplement VC, but there is not much difference between supplement OJ and VC for dose 2mg/day.

Based on these observations, let’s perform some statistical tests to see if the differences in the length of tooth growth is significant by supplement types and dosage.

Using hypothesis tests,confidence intervals and p-values to compare tooth growth by supplement types and dosage

Our null hypothesis is that there is no difference in the length of tooth growth between supplements OJ and VC
Our alternate hypothesis is that there is a difference in the length of tooth growth between supplements OJ and VC

We will test our hypothesis on all 3 dosages.

First, let’s filter the dataset based on dosage so that it’s easier to perform t-test.

dose_0.5 <- subset(ToothGrowth, ToothGrowth$dose == 0.5)
dose_1 <- subset(ToothGrowth, ToothGrowth$dose == 1)
dose_2 <- subset(ToothGrowth, ToothGrowth$dose == 2)

T-Test for 0.5mg/day

t.test(len ~ supp, dose_0.5, paired = FALSE)

Welch Two Sample t-test

data: len by supp t = 3.1697, df = 14.969, p-value = 0.006359 alternative hypothesis: true difference in means between group OJ and group VC is not equal to 0 95 percent confidence interval: 1.719057 8.780943 sample estimates: mean in group OJ mean in group VC 13.23 7.98

T-Test for 1mg/day

t.test(len ~ supp, dose_1, paired = FALSE)

Welch Two Sample t-test

data: len by supp t = 4.0328, df = 15.358, p-value = 0.001038 alternative hypothesis: true difference in means between group OJ and group VC is not equal to 0 95 percent confidence interval: 2.802148 9.057852 sample estimates: mean in group OJ mean in group VC 22.70 16.77

T-Test for 2mg/day

t.test(len ~ supp, dose_2, paired = FALSE)

Welch Two Sample t-test

data: len by supp t = -0.046136, df = 14.04, p-value = 0.9639 alternative hypothesis: true difference in means between group OJ and group VC is not equal to 0 95 percent confidence interval: -3.79807 3.63807 sample estimates: mean in group OJ mean in group VC 26.06 26.14

Based on the three t-tests,

0.5mg/day:
- As the p-value is 0.006359 (< 0.05), we will reject the null hypothesis that there is no difference in the length of tooth growth between supplements OJ and VC.
1mg/day:
- As the p-value is 0.001038 (< 0.05), we will reject the null hypothesis that there is no difference in the length of tooth growth between supplements OJ and VC.
2mg/day:
- As the p-value is 0.9639(> 0.05), we fail to reject the null hypothesis that there is no difference in the length of tooth growth between supplements OJ and VC.

Conclusions and Assumptions

Based on the result of the t-tests, we conclude that supplement OJ is significantly better than supplement VC in enhancing tooth growth at dosage 0.5mg/day and 1mg/day. However, there is no significant difference between both supplements at dosage 2mg/day.

However, do note that these conclusions are based on the following assumption that the sample is representative of the population, and the variables are IID random variables.

For the t-test performed, we assumed that the data are independent (unpaired) and that there is a difference in the variance of the data.

Setting Global Options