Statistical Inference Final Assignment

================================================================================ ### Part 1 - Questions

Part1A: - Exponential Distribution can be simulating using rexp(n,lambda)
- Mean and SD of exp dist is 1/lambda
- Set lambda = 0.2
- Investigate distribution of 40 exopnentials
- Need 1000 simulations

Part1B: - Show the sample mean and compare it to the theoretical mean of the distribution
- Show how variable the sample is (via variance) and compare it to the theoretical variance of the distribution
- Show that distribution is approximately normal

================================================================================

Sample vs. Theoretical Comparisons of the Exponential Distribution

Introduction

In this project, we will be comparing the theoretical mean and standard deviation of the exponential distribution to 40 simulated samples using R programming. The theoretical mean and standard deviation of an exponential distribution are 1/lambda, where lambda is the rate. Thus, the variance of an exponential distribution would be (1/lambda)^2. In our examples, we will be using a rate of lambda=0.2 which would make a theoretical mean and SD of 5 and variance of 25.

Sample vs. Theoretical Mean Comparison

We will start with the easier of the two, comparing the means of 1000 simulations of 40 random exponentials.

## Set seed for reproducibility and define our variables
set.seed(13245)
lambda <- 0.2
sim <- 1000
n <- 40

## Run a simulation and put the results into a matrix 
## with 'sim' cols and 'n' rows
simulation<-matrix(rexp(sim*n,lambda),sim,n)

## Let's find the variable we want to compare 'mean' of all rows
sim_mean<-apply(simulation,1,mean)
mean(sim_mean)
## [1] 5.03555
## Plot the results for visual comparison
hist(sim_mean, breaks = 25, col = 'green',
     main='Distribution of the Mean of 40 Random Exponentials', 
     xlab='Mean',
     ylab='Frequency')
abline(v=1/lambda,lty = 1, lwd=5, col='red')
legend('topright',lty=1,lwd=5,col='red',legend= 'Theoretical Mean')

Summary
Here, we can see that the frequency histogram of the means follows a gaussian distribution as we would expect with it’s mean near 5 (5.04), which is very close to the theoretical mean of 5.

Sample vs. Theoretical Variance Comparison

We will now replicate the previous code with a slight alteration in oreder to compare the variance of the simulations.

## Set seed for reproducibility and define our variables
set.seed(13245)
lambda <- 0.2
sim <- 1000
n <- 40

## Run a simulation and put the results into a matrix 
## with 'sim' cols and 'n' rows
simulation<-matrix(rexp(sim*n,lambda),sim,n)

## Let's find the variable we want to compare 'variance' of all rows
sim_var<-apply(simulation,1,var)
mean(sim_var)
## [1] 25.05544
## Plot the results for visual comparison
hist(sim_var, breaks = 25, col = 'magenta',
     main='Distribution of the Variance of 40 Random Exponentials', 
     xlab='Variance',
     ylab='Frequency')
abline(v=(1/lambda)^2,lty = 1, lwd=5, col='blue')
legend('topright',lty=1,lwd=5,col='blue',legend= 'Theoretical Variance')

Summary Here, we can see that the frequency histogram of the sample variance (mean = 25.06) shows a distribution centered around the theoretical variance (mean = 25.00) of an exponential distribution.

Show that the Sample Means are Normally Distributed (Compare to the CLT)

## Set seed for reproducibility and define our variables
library(rcompanion)
## Warning: package 'rcompanion' was built under R version 4.3.3
set.seed(13245)
lambda <- 0.2
sim1 <- 10
sim3 <- 100
sim2 <- 1000
n <- 40
## Run a simulation and put the results into a matrix 
## with 'sim' cols and 'n' rows
simulation1<-matrix(rexp(sim1*n,lambda),sim1,n)
simulation2<-matrix(rexp(sim2*n,lambda),sim2,n)
simulation3<-matrix(rexp(sim3*n,lambda),sim3,n)
## Let's find the variable we want to compare 'mean' of all rows
sim_mean1<-apply(simulation1,1,mean)
sim_mean2<-apply(simulation2,1,mean)
sim_mean3<-apply(simulation3,1,mean)

## Plot the results for visual comparison with an appropriate normal Curve

## Set our layout for our charts
par(mfrow=c(3,1))

## Plot a normal curve for 10 simulations
plotNormalHistogram(sim_mean1, prob=FALSE, breaks = 25, 
     main='Ten Simulations', 
     xlab='Mean',
     ylab='Frequency',
     length=10)
abline(v=1/lambda,lty = 1, lwd=5, col='red')
legend('topright',lty=1,lwd=5,col='red',legend= 'Theoretical Mean')

## Plot a normal curve for 100 simulations
plotNormalHistogram(sim_mean3, prob=FALSE, breaks = 25, 
     main='Ten Thousand Simulations', 
     xlab='Mean',
     ylab='Frequency',
     length=100)
abline(v=1/lambda,lty = 1, lwd=5, col='red')

##  Plot a normal curve for 1000 simulations
plotNormalHistogram(sim_mean2, prob=FALSE, breaks = 25, col = 'green',
     main='One Thousand Simulations (Our Study)', 
     xlab='Mean',
     ylab='Frequency',
     length=1000)
abline(v=1/lambda,lty = 1, lwd=5, col='red')

Summary According to the central limit theory, the more simulations we perform (i.e the bigger the sample), the closer we will get to a normal distribution around the mean of the population (exponential distribution). Here, I have demonstrated that as our simulations increase from 10 to 100 to 1000 our mean distribution appears more Gaussian around the mean by the typical ‘bell shaped curve’ and 1000 iterations appears accurate enough for analysis.

================================================================================ ### Part 2 - Questions
Analyze ToothGrowth Data in R datasets package
- Load and perform basic explotratory analysis
- Provide a basic summary of the data
- Use confidence intervals and/or hypothesis tests to compare tooth growth by supp and dose (only techniques from class)
- State your conclusion and the assumptions needed for your conclusions ================================================================================

Exploratory Analysis of Tooth Growth

Introduction

We will be performing an exploratory analysis of tooth growth data from the R data sets package. This data set includes data based on the length of teeth in Guinea Pigs that have been given three vitamin C doses (0.5, 1.0 and 2.0mg). The guinea pigs received either Vitamin C or Orange Juice according to the data.

library(datasets)
data("ToothGrowth")
head(ToothGrowth)
dim(ToothGrowth)
## [1] 60  3
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 ...

There are 60 total observations with 3 variables recorded. ‘len’ refers to the length of the teeth of each guinea pig, ‘supp’ refers to the supplement type (VC-vitamin C, OJ-Orange Juice), and ‘dose’ refers to the amount given in milligrams (mg).

Let’s start by looking at all of the data as a whole separated by dose and by supplement type.

library(ggplot2)
g<-ggplot(ToothGrowth,aes(dose,len,color=supp))
g<-g+theme_bw()
g<-g+geom_point(alpha=1/2)
g<-g+facet_grid(.~supp)
g<-g+geom_smooth(method="lm",alpha=1/2,se=FALSE)
g<-g+labs(title='Teeth Length Based on Dose and Route of Administration')
g<-g+labs(x='Dose of Vitamin C',y='Length of Teeth',legend='supplement')
g<-g+guides(color=guide_legend(title="Supplement"))
g
## `geom_smooth()` using formula = 'y ~ x'

Summary Right away, we can see some interesting information. It appears as though vitamin C administration leads to increased length of Teeth. It also appears as though Lower doses of Vitamin C work better in Orange Juice rather than Ascorbic Acid. However, there appears to be more variation in this group. Higher doses of Vitamin C can be delivered either way for a similar result, however there is more variation in the Ascorbic acid group. Let’s create a hypothesis to test.

Hypothesis

It’s easy to see that we have improved growth with incrased dosages of vitamin C, however, we don’t have a control group here (no exposure to Vitamin C). It is for this reason that I am going to compare the route of administration in my hypothesis testing.

Let’s assume that H0 (our null hypothesis) is that consuming vitamin C via Oranje Juice or Ascorbic Acid results in no difference in tooth growth. Therefore, our Ha (Alternate Hypothesis) would be that there is more or less growth when Vitamin C is delivered in a different form. H0!=Ha

Assumptions

This will be a two sided test because we are looking at H0!=Ha (vs H0<Ha). I assume that the guinea pigs are not paired (i.e. they didn’t get all three doses at different times). The variables are iid (independent and identically distributed). Because tooth growth is nondiscrete we assume a normal distribution. We should also assume (based on our graph above) that the different methods have different variances.

Let’s perform a basic t.test to determine if there is significance between the two forms of administration.

OrangeJuice<-ToothGrowth$len[ToothGrowth$supp == "OJ"]
AscorbicAcid<-ToothGrowth$len[ToothGrowth$supp == "VC"]
t.test(OrangeJuice,AscorbicAcid,alternative='two.sided',paired=FALSE,var.equal=FALSE,conf.level=0.975)
## 
##  Welch Two Sample t-test
## 
## data:  OrangeJuice and AscorbicAcid
## t = 1.9153, df = 55.309, p-value = 0.06063
## alternative hypothesis: true difference in means is not equal to 0
## 97.5 percent confidence interval:
##  -0.7510928  8.1510928
## sample estimates:
## mean of x mean of y 
##  20.66333  16.96333

Summary We can see that taking the data aggregated as a whole, our P-Value is not significant and thus we fail to reject the null hypothesis. We can also see that our 95% confidence interval on a two.sided test (-.751 to 8.15), contains 0.0 which means there is a chance of no difference under our H0. I think a lot of this is due to the higher doses of Vitamin C being so similar, so let’s look at the data for each dose (with lower degrees of freedom).

OJlow<-ToothGrowth$len[ToothGrowth$supp == "OJ" & ToothGrowth$dose == 0.5]
AAlow<-ToothGrowth$len[ToothGrowth$supp == "VC" & ToothGrowth$dose == 0.5] 
t.test(OJlow,AAlow,alternative='two.sided',paired=FALSE,var.equal=FALSE,conf.level=0.975)
## 
##  Welch Two Sample t-test
## 
## data:  OJlow and AAlow
## t = 3.1697, df = 14.969, p-value = 0.006359
## alternative hypothesis: true difference in means is not equal to 0
## 97.5 percent confidence interval:
##  1.125088 9.374912
## sample estimates:
## mean of x mean of y 
##     13.23      7.98

Summary At lower doses (0.5mg), we reject our null hypothesis and determine that there is a difference in tooth growth based on route of delivery with a P-value of 0.006 and a confidence interval of 1.12 - 9.37

OJmed<-ToothGrowth$len[ToothGrowth$supp == "OJ" & ToothGrowth$dose == 1.0]
AAmed<-ToothGrowth$len[ToothGrowth$supp == "VC" & ToothGrowth$dose == 1.0] 
t.test(OJmed,AAmed,alternative='two.sided',paired=FALSE,var.equal=FALSE,conf.level=0.975)
## 
##  Welch Two Sample t-test
## 
## data:  OJmed and AAmed
## t = 4.0328, df = 15.358, p-value = 0.001038
## alternative hypothesis: true difference in means is not equal to 0
## 97.5 percent confidence interval:
##  2.2781 9.5819
## sample estimates:
## mean of x mean of y 
##     22.70     16.77

Summary At a dose of 1.0mg, we also reject the null hypothesis and determine that there is a difference in tooth growth based on route of delivery with a P-value of 0.001 and a confidence interval of 2.28 - 9.58.

OJhigh<-ToothGrowth$len[ToothGrowth$supp == "OJ" & ToothGrowth$dose == 2.0]
AAhigh<-ToothGrowth$len[ToothGrowth$supp == "VC" & ToothGrowth$dose == 2.0] 
t.test(OJhigh,AAhigh,alternative='two.sided',paired=FALSE,var.equal=FALSE,conf.level=0.975)
## 
##  Welch Two Sample t-test
## 
## data:  OJhigh and AAhigh
## t = -0.046136, df = 14.04, p-value = 0.9639
## alternative hypothesis: true difference in means is not equal to 0
## 97.5 percent confidence interval:
##  -4.430131  4.270131
## sample estimates:
## mean of x mean of y 
##     26.06     26.14

Summary As expected, we fail to reject the null hypothesis with a p-value of 0.96 at higher doses of vitamin C where it appears as though the route of delivery is inconsequential.

Multiple Testing

We only performed a total of 4 comparisons, however, the more hypothesis that we test, the higher the possibility that we will have a false positive. We can control for this by using the FWER or FDR technique. I will use the latter.

pValues<-c(0.061,0.0064,.0010,.9693)
sum(p.adjust(pValues,method='BH')<0.025)
## [1] 2

Summary So even after adjusting for multiple comparisons we can see that we have two values of significance which are our lower dose treatments.

Conclusions

In conclusion, we have demonstrated with significance that delivery of low dose vitamin C through Oranje Juice is significantly more effective than Ascorbic Acid. However, at higher doses of vitamin C, the method of delivery becomes irrelevant.