Now in the second portion of the project, we’re going to analyze the ToothGrowth data in the R datasets package. 1. Load the ToothGrowth data and perform some basic exploratory data analyses. 2. Provide a basic summary of the data. 3. Use confidence intervals and/or hypothesis tests to compare tooth growth by supp and dose. (Only use the techniques from class, even if there’s other approaches worth considering) 4. State your conclusions and the assumptions needed for your conclusions.

library(tidyverse)
library(grid)
library(gridExtra)
library(flextable)

1. Loading the ToothGrowth data and performing some basic exploratory data analyses

Description of the data set

We can see through the ?ToothGrouthcommand that the response is the length of odontoblasts (cells responsible for tooth growth) in 60 guinea pigs. Each animal received one of three dose levels of vitamin C (0.5, 1, and 2 mg/day) by one of two delivery methods, orange juice (OJ) or ascorbic acid (a form of vitamin C and coded as VC).

It is a data frame with 60 observations on 3 variables.

[,1] len Tooth length [,2] supp Supplement type (VC or OJ). [,3] dose Dose in milligrams/day

Format of the variables: len is numeric, supp is a factor and dose is numeric.

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

Description of the variable supp: in this variable we have 2 levels: OJ and VC, both of 30 observations each.

summary(ToothGrowth$supp)
OJ VC 
30 30

Description of the variable dose: in this variable we have 3 levels: 0.5, 1.0 and 2.0, all of 20 observations each.

ToothGrowth %>% group_by(dose) %>% count()
# A tibble: 3 x 2
# Groups:   dose [3]
   dose     n
  <dbl> <int>
1   0.5    20
2   1      20
3   2      20

Description of the variable len: this is a response numeric variable, with mean = 18.8 and median = 19.2.

summary(ToothGrowth$len)
   Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
   4.20   13.07   19.25   18.81   25.27   33.90

So, the very small difference between median and mean indicates the lack of skewness in the sample distribution, as we can see in the following histogram:

ggplot(ToothGrowth, aes(x=len))+
    geom_histogram(bins = 6, fill= "steelblue1", colour="black")+
    geom_vline(xintercept = mean(ToothGrowth$len),lwd=1,colour="red")+
    geom_vline(xintercept = median(ToothGrowth$len),lwd=1,colour="green")+
    annotate("text",c(17.5,20.5),c(15,16),
             label=c("Mean","Median"),
             colour= c("red","green4"))+
    ylab("Frequency")+
    ggtitle("Frequency of len (Tooth length)")+
    theme(plot.title = element_text(hjust = 0.5))

2. Provide a basic summary of the data

So we have a response variable, len, and two explanatory variables, app and dose, which try to explain the response.

So we will build a summary, grouping len in relation to doseand supp.

We obtain three different tables that summarize the averages of the variable len (Tooth lengths) in relation to supp (Supplement type) and dose.

ToothGrowth %>% group_by(supp) %>% summarise(n=n(),
                                             mean=round(mean(len),2),
                                             sd=round(sd(len),2)) %>% 
    flextable() %>%
    align(align = "center",part = "all")%>%
    set_caption(caption = "Average tooth length (len) and standard deviation by supplement type (supp)")%>%
    width(width = 1.5)
ToothGrowth %>% group_by(dose) %>% summarise(n=n(),
                                             mean=round(mean(len),2),
                                             sd=round(sd(len),2)) %>% 
    flextable() %>%
    align(align = "center",part = "all")%>%
    set_caption(caption = "Average tooth length (len) and standard deviation by dose")%>%
    width(width = 1.5)
ToothGrowth %>% group_by(supp,dose) %>% summarise(n=n(),
                                             mean=round(mean(len),2),
                                             sd=round(sd(len),2)) %>% 
    flextable() %>%
    align(align = "center",part = "all")%>%
    set_caption(caption = "Average tooth length (len) and standard deviation by supplement type (supp) and dose")%>%
    width(width = 1.5)
a<-ToothGrowth %>% group_by(supp) %>% summarise(meanlen=mean(len)) %>% 
    ggplot(aes(x=supp,y=meanlen, fill= supp))+
    geom_bar(stat = "identity")+
    geom_text(aes(label=round(meanlen, digits=1)), vjust=4,colour="white",size=5)+
    ylab("Means of len")+
    theme(legend.position = "none")

c<-ToothGrowth %>% group_by(dose) %>% summarise(meanlen=mean(len)) %>% 
    ggplot(aes(x=dose,y=meanlen, fill= dose))+
    geom_bar(stat = "identity")+
    geom_text(aes(label=round(meanlen, digits=1)), vjust=4,colour="white",size=4)+
    ylab("Means of len")+
    theme(legend.position = "none")

b<-ToothGrowth %>% group_by(supp,dose) %>% summarise(meanlen=mean(len)) %>% 
    ggplot(aes(x=dose,y=meanlen, fill=supp))+
    geom_bar(stat = "identity")+
    facet_wrap(~supp)+
    geom_text(aes(label=round(meanlen, digits=1)),
              vjust=2,colour="white",size=3)+
    ylab("Means of len")+
    theme(legend.position = "none")

grid.arrange(a,c,b, nrow = 2, top= textGrob("Summary of means",
                                            gp= gpar(fontsize=16)))

The summaries of the data indicate a possible correlations of both variables supp and dose with the response variable len. Which will be our alternative hypothesis in the hypothesis test of the next chapter.

3.1 Confidence intervals and/or hypothesis tests to compare tooth growth by supp

Confidence interval

As you can see in the summaries above, there is an important difference in the variances between the groups. So we think it is better to work by the assumption of unequal variances: a. We apply the confidence interval general formula: \[\bar{X_{OJ}}-\bar{X_{CV}}\pm t_{df}\times SE => \bar{X_{OJ}}-\bar{Y_{CV}}\pm t_{df}\times \sqrt{\frac{S^2_{OJ}}{n_{OJ}}+\frac{S^2_{CV}}{n_{CV}}}\] b. To determine \(df\), we will use the Welch method, as it is explained in our textbook (p. 83). \[df=\frac{(S^2_x/n_x+S^2_y/n_y)^2}{(\frac{S^2_x}{n_x})^2/(n_x-1)+(\frac{S^2_y}{n_y})^2/(n_y-1)}\] We prepare the data generating the correspondent R objects;

sdOJ <- sd(ToothGrowth$len[ToothGrowth$supp == "OJ"])
sdVC <- sd(ToothGrowth$len[ToothGrowth$supp == "VC"])
se <- sqrt((sdOJ^2/30)+(sdVC^2/30))
meanOJ <- mean(ToothGrowth$len[ToothGrowth$supp == "OJ"])
meanVC <- mean(ToothGrowth$len[ToothGrowth$supp == "VC"])
difOJ_VC <- meanOJ-meanVC

and we create the specific function freedom (nice name!) in order to easily calculate de \(df\):

freedom <- function(s_x=1, s_y=1, n_x, n_y){
    p <- ((s_x^2/n_x) + (s_y^2/n_y))^2 / (((s_x^2/n_x)^2 /(n_x-1)) + ((s_y^2/n_y)^2 /(n_y-1)))
    print(p)
}
dfsupp <- freedom(sdOJ,sdVC,30,30)
[1] 55.30943

Now we can compute the CI:

difOJ_VC+c(-1,1)*qt(.975,df = dfsupp)*se
[1] -0.1710156  7.5710156

We can see in that the confidence interval of \(\bar{X_{OJ}}-\bar{X_{CV}}\) also includes the point \(0\).

Hypothesis test

We can execute an hypothesis test: \(H_0: \mu_1=\mu_2 => \mu_1-\mu_2 = 0; H_A: \mu_1 \neq \mu_2 => \mu_1-\mu_2 \neq 0\)

The test statistic can be calculated es follows: \[ TS=\frac{\bar{X_{OJ}}-\bar{X_{CV}}-(\mu_1-\mu_2)}{SE};(\mu_1-\mu_2)=0=>\frac{\bar{X_{OJ}}-\bar{X_{CV}}}{SE} \] We will reject if: \(\mid TS\mid \geq t_{1-\alpha/2}\) OR \(p.value\times2\leq \alpha\)

The p-value obtained will be:

pt(difOJ_VC/se,58,lower.tail = FALSE)*2
[1] 0.06039337

The p-value is bigger than \(\alpha\). Which means that we fail to reject the null hypothesis \(H_0\).

T.test

Anyhow, we can perform the same inference analysis you have seen above using the function t.test(), much faster and easier to use.

t.test(ToothGrowth$len[ToothGrowth$supp== "OJ"],
       ToothGrowth$len[ToothGrowth$supp== "VC"],
       paired = FALSE,
       alternative = "two.sided",
       var.equal = FALSE,
       conf.level = .95)

    Welch Two Sample t-test

data:  ToothGrowth$len[ToothGrowth$supp == "OJ"] and ToothGrowth$len[ToothGrowth$supp == "VC"]
t = 1.9153, df = 55.309, p-value = 0.06063
alternative hypothesis: true difference in means is not equal to 0
95 percent confidence interval:
 -0.1710156  7.5710156
sample estimates:
mean of x mean of y 
 20.66333  16.96333

Giving to the function the parameters paired = FALSE, alternative = "two.sided", var.equal = FALSE(see above) and conf.level = .95, we get the same outcomes obtained above.

3.2 Confidence intervals and/or hypothesis tests to compare tooth growth by dose

We do not dispose of a control group, so we can’t compare the impact of different doses comparing with a \(dose=0\). So we will compare: \(dose=0.5\) with \(dose=1.0\) and \(dose=1.0\) with \(dose=2.0\).

Comparing \(dose=0.5\) and \(dose=1.0\)

Our hypothesis test is:

\(H_0: \mu_{0.5}=\mu_{1.0} => \mu_{0.5}-\mu_{1.0} = 0\)
\(H_A: \mu_{0.5} \neq \mu_{1.0} => \mu_{0.5}-\mu_{1.0} \neq 0\)

We have already showed the different options to get IC and p-values. So, in this case we will obtain the p-values using the t.test() function, which we have seen es being the quickest option. We will use the same option, var.equal=FALSE, which computes the \(df\) using the Welch option. And the confidence level is 95%.

t1<-t.test(ToothGrowth$len[ToothGrowth$dose== 0.5],
       ToothGrowth$len[ToothGrowth$dose== 1],
       paired = FALSE,
       alternative = "two.sided",
       var.equal = FALSE,
       conf.level = .95)$p.value

In this case the p-value is equal to 0.0000001. Very tiny. Which allows us to reject the null hypothesis.

The power of this test

We set delta equal to the difference of means between dose=0.5 and dose=1.0 groups; we set sd equal to pooled s of both groups; and the size of each group equal to 20.

std0.5<- sd(ToothGrowth$len[ToothGrowth$dose== 0.5])
std1.0<- sd(ToothGrowth$len[ToothGrowth$dose== 1])
spool0.5_1.0 <- sqrt((std0.5^2*20-1+std1.0^2*20-1)/(20+20-2))
delta0.5_1.0 <- mean(ToothGrowth$len[ToothGrowth$dose==0.5])-mean(ToothGrowth$len[ToothGrowth$dose== 1])
p1 <- power.t.test(n=20,
             delta = delta0.5_1.0,
             sd=spool0.5_1.0,
             sig.level = 0.05,
             type = "two.sample",
             alternative = "two.sided")$power

Power of this test is 0.9999863, very high. Which indicates that we could obtain the rejection of the null hypothesis even with a sample much smaller or with a \(\alpha\) smaller than the current one.

Comparing \(dose=1.0\) and \(dose=2.0\)

Our hypothesis test is:
\(H_0: \mu_{1.0}=\mu_{2.0} => \mu_{1.0}-\mu_{2.0} = 0\)
\(H_A: \mu_{1.0} \neq \mu_{2.0} => \mu_{1.0}-\mu_{2.0} \neq 0\)

t2<-t.test(ToothGrowth$len[ToothGrowth$dose== 1],
       ToothGrowth$len[ToothGrowth$dose== 2],
       paired = FALSE,
       alternative = "two.sided",
       var.equal = FALSE,
       conf.level = .95)$p.value

Using the same methodology as above we obtain a p-value of 0.0000191, very tiny. So we reject the null hypothesis.

The power of this test

We set delta equal to the difference of means between dose=0.5 and dose=1.0 groups; we set sd equal to pooled s of both groups; and the size of each group equal to 20.

std1.0<- sd(ToothGrowth$len[ToothGrowth$dose== 1])
std2.0<- sd(ToothGrowth$len[ToothGrowth$dose== 2])
spool1.0_2.0 <- sqrt((std1.0^2*20-1+std2.0^2*20-1)/(20+20-2))
delta1.0_2.0 <- mean(ToothGrowth$len[ToothGrowth$dose==1])-mean(ToothGrowth$len[ToothGrowth$dose== 2])
p2 <- power.t.test(n=20,
             delta = delta1.0_2.0,
             sd=spool1.0_2.0,
             sig.level = 0.05,
             type = "two.sample",
             alternative = "two.sided")$power

Power of this test is 0.9965287, very high. Which indicates that, as in the previous case, we could obtain the rejection of the null hypothesis even with a sample much smaller or with a \(\alpha\) smaller than the current one.

4. Conclusions and the assumptions

  1. The difference observed in the sample between the averages of Tooth length (len) by Supplement type (supp) is not enough to reject the null hypothesis and to be considered statistically significant. So we can’t say that the supplement type has any correlation with tooth length.
  2. The difference observed in the sample between the averages of Tooth length (len) by Dose (dose) is enough to reject the null hypothesis in the case of doses of 0.5 and 1.0. So we can infer that to increase the dose from 0.5 to 1.0 has a correlation with tooth length. On the other hand we have to note tha this correlation has an high power that could allow us to get the same outcome even with smaller samples.
  3. The difference observed in the sample between the averages of Tooth length (len) by Dose (dose) is enough to reject the null hypothesis in the case of doses of 1.0 and 2.0. So we can infer that to increase the dose from 1.0 to 2.0 has a correlation with tooth length. As in the previous case, we have to note tha this correlation has an high power that could allow us to get the same outcome even with smaller samples.
  4. As a final conclusion, we have we can state that it exist a very robust correlation between dose of the supplement and tooth length; while the correlation between supplement type and tooth length is unclear.