Case studies of inferential analysis

1. Exponential distribution and the CLT

\[\begin{aligned} &X \text{~}Exponential(\lambda) \\ &f(x)=\lambda e^{-\lambda x} \ \ \text{ for }\ \ x>0 \\ &F(x) = 1- e^{-\lambda x} \\ &E(X)=\frac{1}{\lambda} \\ &VAR(X)=\frac{1}{\lambda^2} \end{aligned}\]

# Load of necessary libraries
library(lattice)
require(gridExtra)
library(ggplot2)
library(cowplot)
library(patchwork)



# initialization of values
lambda = 0.2
expo_values = rexp(40000, lambda)
lambdas = apply(matrix(expo_values, 1000), 1, mean)


# Calculation of statistics and parameters
sample_mean = round(mean(lambdas), 2)
theoretical_mean = round(1/lambda, 2)
sample_variance = round(var(expo_values), 2)
theoretical_variance = round(1/lambda^2, 2)


## Plotting

# Histogram of means of 1000 random exponential of size 40
class_width = round(3.5*sd(lambdas)*length(lambdas)^(-1/3), 1)
hist_means = ggplot(data.frame(lambdas), aes(x=lambdas)) 
hist_means = hist_means + geom_histogram(aes(y=after_stat(density)), binwidth = class_width, fill="#94c6ff") 
hist_means = hist_means + geom_density(kernel="gaussian", aes(y=after_stat(density))) 
hist_means = hist_means + geom_vline(xintercept = theoretical_mean, colour="red")
hist_means = hist_means + labs(title="Sample means", x="", y="Density")
hist_means = hist_means + annotate("text", x=theoretical_mean, y=0.6, parse=TRUE, label=paste("~mu ==", theoretical_mean))
hist_means = hist_means + theme_cowplot()

# Histogram of 40_000 random exponential 
class_width = round(3.5*sqrt(sample_variance)*length(expo_values)^(-1/3), 1)
hist_raw = ggplot(data.frame(expo_values), aes(x=expo_values))
hist_raw = hist_raw + geom_histogram(aes(y=after_stat(density)), binwidth=class_width, fill="#94c6ff")
hist_raw = hist_raw + geom_density(aes(y=after_stat(density)))
hist_raw = hist_raw + geom_vline(xintercept = sample_mean, colour="red")
hist_raw = hist_raw + labs(title="Raw exponential", x="", y="Density")
hist_raw = hist_raw + annotate("text", x=sample_mean, y=0.19, parse=TRUE, label=paste("bar(x)==", sample_mean))
hist_raw = hist_raw + theme_cowplot()

# print both plots
layout = c(area(1,1,1,2), area(2,3,2,4))
hist_raw + hist_means + plot_layout(design=layout)

\[\textbf{Mean comparison} \\ \bar{x}=4.99 \ \ \mu=5 \\ \textbf{Variance comparison} \\ S^2=24.88 \ \ \sigma^2=25\]

Given that \(n \rightarrow \infty\) from the CLT and LLN is safe to assume that both sample statistics (mean and variance) are normally distributed and centered at what they try to estimate (the population parameters, in this case: \(1/\lambda\) and \(1/\lambda^2\)), this is further confirmed when observing the shape of the distribution of sample means of 40 random exponential draws, as well as the closeness of the observed statistics and theoretical parameters.

2. Guinea pigs’ tooth growth analysis

Distribution of tooth length by group

# Getting first group (VC) from the ToothGrowth dataframe
vc_supp = ToothGrowth[ToothGrowth$supp=="VC",c("len", "dose")]

# Calculating sample statistics
mean_vc = round(mean(vc_supp$len), 2)
sd_vc = round(sd(vc_supp$len), 2)

# Building histogram for tooth length
class_width = 3.5*sd(vc_supp$len)*length(vc_supp$len)^(-1/3)
hist_vc = ggplot(data=vc_supp, aes(x=len))
hist_vc = hist_vc + geom_histogram(aes(y=after_stat(density)), binwidth=class_width, fill="#94c6ff")
hist_vc = hist_vc + geom_density(aes(y=after_stat(density)))
hist_vc = hist_vc + geom_vline(xintercept = mean_vc, colour="red")
hist_vc = hist_vc + labs(title = "Supplement VC", y="Density", x="Tooth length")
hist_vc = hist_vc + annotate("text", x=37, y=.04, parse=TRUE, label=paste("bar(x)==", mean_vc)) + annotate("text", x=36.5, y=.03, parse=TRUE, label=paste("S==", sd_vc))
hist_vc = hist_vc + theme_cowplot()

# Getting second group (OJ) from the ToothGrowth dataframe
oj_supp = ToothGrowth[ToothGrowth$supp=="OJ",c("len", "dose")]

# Calculating sample statistics
mean_oj = round(mean(oj_supp$len), 2)
sd_oj = round(sd(oj_supp$len), 2)

# Building histogram for tooth length
class_width = 3.5*sd(oj_supp$len)*length(oj_supp$len)^(-1/3)
hist_oj = ggplot(data=oj_supp, aes(x=len))
hist_oj = hist_oj + geom_histogram(aes(y=after_stat(density)), binwidth=class_width, fill="#94c6ff")
hist_oj = hist_oj + geom_density(aes(y=after_stat(density)))
hist_oj = hist_oj + geom_vline(xintercept = mean_oj, colour="red")
hist_oj = hist_oj + labs(title = "Supplement OJ", y="Density", x="Tooth length")
hist_oj = hist_oj + annotate("text", x=10, y=.06, parse=TRUE, label=paste("bar(x)==", mean_oj)) + annotate("text", x=9.7, y=.045, parse=TRUE, label=paste("S==", sd_oj))
hist_oj = hist_oj + theme_cowplot()

# print both plots
hist_oj + hist_vc + plot_layout(design=layout)

Boxplot by supplement and dose

box_plot = ggplot(data = ToothGrowth, aes(x = factor(dose)))
box_plot = box_plot + geom_boxplot(aes(y = len, fill = supp))
box_plot = box_plot + labs(x = "Dose [mg/day]", y = "Tooth length")
box_plot = box_plot + guides(fill=guide_legend(title="Supplement"))
box_plot = box_plot + stat_summary(fun=mean, shape=18,size=2, col="red", geom = "point", show_guide=FALSE, aes(y=len, fill = supp)) #show mean as a red dot
box_plot = box_plot + theme_cowplot() + scale_fill_brewer(palette="Set2")
box_plot

Hypothesis testing for differences in means

It is safe to assume that for any kind of drug or substance the greater the dose the higher the effect, for simplicity of the current analysis a far more interesting question to ask is whether or not the supplement type affects in any way the absorption rate of vitamin C by the guinea pigs and therefore tooth growth. Formally: \[H_0: \mu_{oj}-\mu_{vc}=0 \\ H_a:\mu_{oj}-\mu_{vc}>0 \]

t.test(oj_supp$len ,vc_supp$len, var.equal = TRUE, alternative = "greater")

## 
##  Two Sample t-test
## 
## data:  oj_supp$len and vc_supp$len
## t = 1.9153, df = 58, p-value = 0.0302
## alternative hypothesis: true difference in means is greater than 0
## 95 percent confidence interval:
##  0.4708204       Inf
## sample estimates:
## mean of x mean of y 
##  20.66333  16.96333

Due to \(P\text{-}value<\alpha\) the null hypothesis is rejected in favor of the alternative, concluding that there’s enough statistical evidence to claim (with 95% confidence) that when guinea pigs are given vitamin C via orange juice the tooth length increases significantly more than when given ascorbic acid.

The previous T-test was performed under the following assumptions:

\(\bullet\) Independence: Tooth growth in one group is not affected by the other and viceversa.

\(\bullet\) Normality: Because \(n\geq30\) the use of t-test is still available.