The Effect of Vitamin C on Tooth Growth in Guinea Pigs

The response is the length of odontoblasts (teeth) in each of 10 guinea pigs at each of three dose levels of Vitamin C (0.5, 1, and 2 mg) with each of two delivery methods (orange juice or ascorbic acid).We are interested to know the effect of dosage and delivery method on tooth growth

Data Loading,Preproccessing and some exploratory Tables

Load Tooth Data and create the required Factors for T-tests

data(ToothGrowth)
ToothGrowth$dose<-as.factor(ToothGrowth$dose)
colnames(ToothGrowth)<-c("ToothLength","Supplement","Dose")
set.seed(2017)
par(mfrow = c(1, 2))
kable(list(summary(ToothGrowth),head(ToothGrowth)),format="html")%>%
      kable_styling(bootstrap_options = "striped", full_width = T,position="center")
ToothLength Supplement Dose
Min. : 4.20 OJ:30 0.5:20
1st Qu.:13.07 VC:30 1 :20
Median :19.25 NA 2 :20
Mean :18.81 NA NA
3rd Qu.:25.27 NA NA
Max. :33.90 NA NA
ToothLength Supplement Dose
4.2 VC 0.5
11.5 VC 0.5
7.3 VC 0.5
5.8 VC 0.5
6.4 VC 0.5
10.0 VC 0.5 
ToothGrowth$DoseGroup<-ToothGrowth$Dose
ToothGrowth$DoseGroup<-mapvalues(ToothGrowth$DoseGroup, from = c("0.5", "1","2"), to = c("1", "2","3"))
ToothGrowth$DoseSuppl <- with(ToothGrowth, interaction(Dose,Supplement))
ToothGrowth$DoseSuppl<- factor(ToothGrowth$DoseSuppl, levels = c("0.5.VC", "0.5.OJ","1.VC","1.OJ","2.VC","2.OJ"))
ToothGrowth$DoseSupplGroup<-mapvalues(ToothGrowth$DoseSuppl, from = c("0.5.VC", "0.5.OJ","1.VC","1.OJ","2.VC","2.OJ"), to = c("1", "2","3","4","5","6"))

Examine the case if Tooth Length is influenced by different Supplements

g1=ggplot(ToothGrowth,aes(x=Supplement, y=ToothLength))+ geom_boxplot(aes(fill=Supplement))+scale_fill_brewer(palette="Set1")+labs(y="")
title1=textGrob("Figure1.ToothLength per Supplement", gp=gpar(fontface="bold",fontsize=12))
grid.arrange(g1,top=title1)

### Examine ToothLength Variances per Supplement

tapply(ToothGrowth$ToothLength,ToothGrowth$Supplement,var)
##       OJ       VC 
## 43.63344 68.32723

Differentiation in Variances may be assumed for the t-test

T-test Results Toothlength per Supplement

t.test(ToothGrowth$ToothLength~ToothGrowth$Supplement,paired = F,var.equal = F)
## 
##  Welch Two Sample t-test
## 
## data:  ToothGrowth$ToothLength by ToothGrowth$Supplement
## 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 in group OJ mean in group VC 
##         20.66333         16.96333

Conclussion:

We assume that different supplements have no effect on Tooth Length because p-value is above 0.05

Examine the case if Tooth Length is influenced by different Dose

Examine ToothLength Variances per Dose

tapply(ToothGrowth$ToothLength,ToothGrowth$Dose,var)
##      0.5        1        2 
## 20.24787 19.49608 14.24421

Equality of the Variances may be assumed for the t-test

Running t-test and store pvalues,p-Adusted(Bonferroni), p-Adjusted(BH) in a database

PVal<-as.vector("numeric")
PValID<-as.vector("Character")
a1<-1
for (i1 in seq(1, 2, 1)){
for (j1 in seq(i1+1, 3, 1)){
PVal[a1]<-t.test(ToothGrowth$ToothLength[ToothGrowth$DoseGroup==i1],                   ToothGrowth$ToothLength[ToothGrowth$DoseGroup==j1],paired = F,var.equal = T)$p.value
PValID[a1]<- paste(i1, j1, sep="-")
a1<-a1+1
}
}
PValDose<-as.data.frame(cbind(PValID,PVal))
names(PValDose)<-c("PValID","PVal")
PValDose$PVal<-as.numeric(paste(PValDose$PVal))
PValDose$PValBonferoni<-round(p.adjust(PValDose$PVal, method = "bonferroni"),3)
PValDose$PValBH<-round(p.adjust(PValDose$PVal, method = "BH"),3)
PValDose$PVal<-round(PValDose$PVal,3)

Ploting the results

g2=ggplot(ToothGrowth,aes(x=Dose, y=ToothLength))+ geom_boxplot(aes(fill=Dose))+scale_fill_brewer(palette="Set2")
title2=textGrob("Figure2.ToothLength per Dose", gp=gpar(fontface="bold",fontsize=12))
grid.arrange(g2,top=title2)

kable(PValDose,format="html")%>%
  kable_styling(bootstrap_options = "striped", full_width = F)
PValID PVal PValBonferoni PValBH
1-2 0 0 0
1-3 0 0 0
2-3 0 0 0 
par(mfrow=c(1,2),oma = c(0, 0, 1, 0))
plot(PValDose$PVal,PValDose$PValBonferoni,ylim=c(0, 1), xlim=c(0,1),xlab = "P Values", ylab="P.Adj Bonferroni",pch = 19)
plot(PValDose$PVal,PValDose$PValBH,ylim=c(0, 1), xlim=c(0,1),xlab = "P Values",ylab="P.Adj BH", pch = 19)
mtext("Ploting P-values against P-Adjusted",outer = TRUE,cex = 1)

Conclussion:

We see that different Doses influence the Tooth Length in all tests we ran. There is a difference in all pairs compared

Examine the case if Tooth Length is influenced by combined Dose-Supplement

Examine ToothLength Variances per Dose-Supplement

tapply(ToothGrowth$ToothLength,ToothGrowth$DoseSuppl,var)
##    0.5.VC    0.5.OJ      1.VC      1.OJ      2.VC      2.OJ 
##  7.544000 19.889000  6.326778 15.295556 23.018222  7.049333

We assume that the Variances are different for the t-test

Running t-test and store pvalues,p-Adusted(Bonferroni), p-Adjusted(BH) in a database

PVal2<-as.vector("numeric")
PValID2<-as.vector("Character")
a2<-1
for (i2 in seq(1, 5, 1)){
for (j2 in seq(i2+1, 6, 1)){
PVal2[a2]<-t.test(ToothGrowth$ToothLength[ToothGrowth$DoseSupplGroup==i2],                   ToothGrowth$ToothLength[ToothGrowth$DoseSupplGroup==j2],paired = F,var.equal = F)$p.value
PValID2[a2]<- paste(i2, j2, sep="-")
a2<-a2+1
}
}
PValDoseSuppl<-as.data.frame(cbind(PValID2,PVal2))
names(PValDoseSuppl)<-c("PValID2","PVal2")
PValDoseSuppl$PVal2<-as.numeric(paste(PValDoseSuppl$PVal2))
PValDoseSuppl$PVal2Bonferoni<-round(p.adjust(PValDoseSuppl$PVal2, method = "bonferroni"), 3)
PValDoseSuppl$PVal2BH<-round(p.adjust(PValDoseSuppl$PVal2, method = "BH"),3)
PValDoseSuppl$PVal2<-round(PValDoseSuppl$PVal2,3)

Ploting the results

g3=ggplot(ToothGrowth,aes( x=DoseSupplGroup,y=ToothLength,label=DoseSupplGroup))+ geom_boxplot(aes(fill=Supplement))+scale_fill_brewer(palette="Set3")+labs(x="Dose-Supplement Group")+facet_grid(~Dose)
title3=textGrob("Figure2. ToothLength per Dose-Supplement Combined", gp=gpar(fontface="bold",fontsize=14))
grid.arrange(g3,top=title3)

kable(PValDoseSuppl,format="html")%>%
  kable_styling(bootstrap_options = "striped", full_width = F)
PValID2 PVal2 PVal2Bonferoni PVal2BH
1-2 0.006 0.095 0.009
1-3 0.000 0.000 0.000
1-4 0.000 0.000 0.000
1-5 0.000 0.000 0.000
1-6 0.000 0.000 0.000
2-3 0.046 0.690 0.053
2-4 0.000 0.001 0.000
2-5 0.000 0.000 0.000
2-6 0.000 0.000 0.000
3-4 0.001 0.016 0.002
3-5 0.000 0.001 0.000
3-6 0.000 0.000 0.000
4-5 0.097 1.000 0.103
4-6 0.039 0.588 0.049
5-6 0.964 1.000 0.964 
par(mfrow=c(1,2),oma = c(0, 0, 1, 0))
plot(PValDoseSuppl$PVal2,PValDoseSuppl$PVal2Bonferoni,ylim=c(0, 1), xlim=c(0,1),xlab = "P Values", ylab="P.Adj Bonferroni",pch = 19)
plot(PValDoseSuppl$PVal2,PValDoseSuppl$PVal2BH,ylim=c(0, 1), xlim=c(0,1),xlab = "P Values",ylab="P.Adj BH", pch = 19)
mtext("Ploting P-values against P-Adjusted",outer = TRUE,cex = 1)

Conclussion:

After completed pairwise comparisons we are having quite similar results from each method

  1. In individual pair comparison all pairs seem to have significant differences except pairs 5-6, 4-5, and 2-3.

  2. Applying Bonferroni method, again all pairs seem to have significant differences except pairs 5-6, 4-6,4-5,2-3 and 1-2.

  3. Applying BH method, again all pairs seem to have significant differences except pairs 5-6, 4-6,4-5 and 2-3.