For Female Mice

# Box plot of BMP_2 for females grouped by disease stage.
female_data <- data[data$Sex == "Female",]
ggplot(female_data, aes(x=Age_group, y=BMP_2, fill=Category)) + 
    geom_boxplot() +
    facet_wrap(~Age_group, scale="free") 

When controlling for Female mice we see that difference between End Stage levels of BMP_1 may be statistically significant.

# Box plot of BMP_1 for females grouped by disease stage.
female_data <- data[data$Sex == "Female",]
ggplot(female_data, aes(x=Age_group, y=BMP_1, fill=Category)) + 
    geom_boxplot() +
    facet_wrap(~Age_group, scale="free")

When controlling for female mice, we see that there is a strong case that there is a statistically significant difference in the mean BMP_2 level for the four disease stages given.

For Male Mice

# This is for BMP_1
male_data <- data[data$Sex == "Male",]
ggplot(male_data, aes(x=Age_group, y=BMP_1, fill=Category)) + 
    geom_boxplot() +
    facet_wrap(~Age_group, scale="free")

When controlling for male mice we do not see promising data as both control groups have higher mean BMP_1 levels than their test counterparts.

# This is for BMP_2
male_data <- data[data$Sex == "Male",]
ggplot(male_data, aes(x=Age_group, y=BMP_2, fill=Category)) + 
    geom_boxplot() +
    facet_wrap(~Age_group, scale="free")

When controlling for male mice, we see some promising results with the mean BMP_2 level being higher in the test groups, when compared to their respective control groups for the disease stages Early, End, and Intermediate Stages.

While these box-whisker plots have mixed results, we will do a Three-Way ANOVA to determine whether the means of each sub-population are equivalent.

# Three-Way ANOVA 
# Source: https://stats.oarc.ucla.edu/stata/faq/how-can-i-understand-a-three-way-interaction-in-anova/
# These are great guys ^^^
res.aov2 <- aov(BMP_2 ~ Age_group * Sex * Category, data = data)
#print(summary(res.aov2)[1])
pander(summary(res.aov2),caption="Three-way ANOVA of factors Test/Control, Disease Stage, Sex on BMP_2 levels")

The results of the F-test following the ANOVA, we see that the Age Group is a significant factor \(p=0.000426< 0.001\), while neither Sex nor Category are significant factors. However there is a strong interaction effect between Sex and Category, and a small interaction effect between Sex, Age_group and Category. Lets do more ANOVA’s to understand these effects. We will first group by Age_group and run ANOVAs to understand the interaction effect of Sex*Category.

# iterate through the age groups and conduct a 2-way ANOVA for each group
# where the factors are Test/Control and Male/Female, and the output variable is BMP_2 level
for (group in unique(data$Age_group)){
  temp_data <- data[data$Age_group == group,]
  tryCatch({
  temp_anova <- aov(BMP_2 ~ Category*Sex, data=temp_data)
  pander(summary(temp_anova), caption=group)
  }, error=function(e){})
  }

At this stage we run a Tukey Honest Significant Difference test to compare the four groups and determine which groups are different.

# iterate through the age groups and conduct a 2-way ANOVA for each group
# where the factors are Test/Control and Male/Female, and the output variable is BMP_1 level
for (group in unique(data$Age_group)){
  temp_data <- data[data$Age_group == group,]
  tryCatch({
  temp_anova <- aov(BMP_2 ~ Category*Sex, data=temp_data)
  temp_tukey <- TukeyHSD(temp_anova, conf.level=.95)
  pander(temp_tukey,caption=group)
  }, error=function(e){})
  }

## Warning in pander.default(temp_tukey, caption = group): No pander.method for
## "TukeyHSD", reverting to default.No pander.method for "multicomp", reverting to
## default.

• Category:

  	diff	lwr	upr	p adj
  Test-Control	-41.93	-210.1	126.3	0.6045

• Sex:

  	diff	lwr	upr	p adj
  Male-Female	-83.08	-253.9	87.72	0.3178

• Category:Sex:

  	diff	lwr	upr	p adj
  Test:Female-Control:Female	251.6	-105.5	608.8	0.2232
  Control:Male-Control:Female	174.5	-164.3	513.3	0.4751
  Test:Male-Control:Female	-57.36	-373.9	259.2	0.9535
  Control:Male-Test:Female	-77.16	-416	261.6	0.9134
  Test:Male-Test:Female	-309	-625.6	7.551	0.05696
  Test:Male-Control:Male	-231.8	-527.6	63.88	0.1539

## Warning in pander.default(temp_tukey, caption = group): No pander.method for
## "TukeyHSD", reverting to default.No pander.method for "multicomp", reverting to
## default.

• Category:

  	diff	lwr	upr	p adj
  Test-Control	36.33	-50.94	123.6	0.3792

• Sex:

  	diff	lwr	upr	p adj
  Male-Female	38.78	-48.48	126	0.349

• Category:Sex:

  	diff	lwr	upr	p adj
  Test:Female-Control:Female	38.9	-124.1	201.9	0.8879
  Control:Male-Control:Female	38.77	-124.3	201.8	0.8888
  Test:Male-Control:Female	78.12	-97.96	254.2	0.5613
  Control:Male-Test:Female	-0.13	-163.2	162.9	1
  Test:Male-Test:Female	39.22	-136.9	215.3	0.9061
  Test:Male-Control:Male	39.35	-136.7	215.4	0.9053

## Warning in pander.default(temp_tukey, caption = group): No pander.method for
## "TukeyHSD", reverting to default.No pander.method for "multicomp", reverting to
## default.

• Category:

  	diff	lwr	upr	p adj
  Test-Control	35.72	-65.18	136.6	0.4757

• Sex:

  	diff	lwr	upr	p adj
  Male-Female	-26.98	-128.2	74.25	0.5906

• Category:Sex:

  	diff	lwr	upr	p adj
  Test:Female-Control:Female	63.25	-119.2	245.7	0.7831
  Control:Male-Control:Female	4.577	-188.3	197.5	0.9999
  Test:Male-Control:Female	6.153	-186.8	199.1	0.9998
  Control:Male-Test:Female	-58.67	-247	129.7	0.8323
  Test:Male-Test:Female	-57.1	-245.4	131.2	0.8432
  Test:Male-Control:Male	1.576	-196.9	200.1	1

## Warning in pander.default(temp_tukey, caption = group): No pander.method for
## "TukeyHSD", reverting to default.No pander.method for "multicomp", reverting to
## default.

• Category:

  	diff	lwr	upr	p adj
  Test-Control	131.5	-7.123	270.1	0.06178

• Sex:

  	diff	lwr	upr	p adj
  Male-Female	27.47	-111.2	166.1	0.6837

• Category:Sex:

  	diff	lwr	upr	p adj
  Test:Female-Control:Female	184.7	-90.26	459.6	0.2678
  Control:Male-Control:Female	80.44	-194.5	355.4	0.8448
  Test:Male-Control:Female	167.4	-98.43	433.3	0.3196
  Control:Male-Test:Female	-104.2	-366.3	157.9	0.686
  Test:Male-Test:Female	-17.25	-269.8	235.3	0.9974
  Test:Male-Control:Male	86.97	-165.6	339.6	0.7711

While this may seem like a lot of information Lets make a table for the adjusted p-value for the groups we are actually interested in comparing.

res.aov2 <- aov(BMP_1 ~ Age_group * Sex * Category, data = data)
pander(summary(res.aov2),caption="Three-way ANOVA where the output variable is BMP_1 level, and the factors are Disease stage, Sex, Test/Control")

# iterate through the age groups and conduct a 2-way ANOVA for each group
# where the factors are Test/Control and Male/Female, and the output variable is BMP_1 level
for (group in unique(data$Age_group)){
  temp_data <- data[data$Age_group == group,]
  tryCatch({
  temp_anova <- aov(BMP_1 ~ Category*Sex, data=temp_data)
  pander(summary(temp_anova), caption=group)
  }, error=function(e){})
  }

# source: https://www.r-bloggers.com/2021/08/how-to-perform-tukey-hsd-test-in-r/
for (group in unique(data$Age_group)){
  temp_data <- data[data$Age_group == group,]
  tryCatch({
  temp_anova <- aov(BMP_1 ~ Category*Sex, data=temp_data)
  temp_tukey <- TukeyHSD(temp_anova, conf.level=.95)
  pander(temp_tukey,caption=group)
  }, error=function(e){})
  }

## Warning in pander.default(temp_tukey, caption = group): No pander.method for
## "TukeyHSD", reverting to default.No pander.method for "multicomp", reverting to
## default.

## Warning in pander.default(temp_tukey, caption = group): No pander.method for
## "TukeyHSD", reverting to default.No pander.method for "multicomp", reverting to
## default.

Multiple Regression

A linear model can developed to associate the disease stage with the level of BMP_2 and BMP_1. This is resolved in the ANOVA, but not explicitly investigated.

linearmodel <- lm(BMP_2 ~ Age_group, data = data)
plot(linearmodel)

pander(summary(linearmodel))

T-tests if you need them

T-tests are included but not recommended because of significant increase in the Family Wise error rate. Meaning you are more likely to have a false positive in your results.

# t tests by sex
pander(t.test (BMP_2 ~ Category , var.equal=TRUE, alternative="less", data = subset(data, Sex == "Male")))

## 
## --------------------------------------------------------------------------------
##  Test statistic   df   P value   Alternative hypothesis   mean in group Control 
## ---------------- ---- --------- ------------------------ -----------------------
##     0.09206       54   0.5365             less                    249.1         
## --------------------------------------------------------------------------------
## 
## Table:  Two Sample t-test: `BMP_2` by `Category` (continued below)
## 
##  
## --------------------
##  mean in group Test 
## --------------------
##        245.6        
## --------------------

t.test (BMP_2 ~ Category , var.equal=TRUE, alternative="less", data = subset(data, Sex == "Female"))

## 
##  Two Sample t-test
## 
## data:  BMP_2 by Category
## t = -2.242, df = 44, p-value = 0.01503
## alternative hypothesis: true difference in means between group Control and group Test is less than 0
## 95 percent confidence interval:
##       -Inf -30.31352
## sample estimates:
## mean in group Control    mean in group Test 
##              199.3309              320.3100

t.test (BMP_1 ~ Category , var.equal=TRUE, alternative="less", data = subset(data, Sex == "Male"))

## 
##  Two Sample t-test
## 
## data:  BMP_1 by Category
## t = 0.80213, df = 16, p-value = 0.7829
## alternative hypothesis: true difference in means between group Control and group Test is less than 0
## 95 percent confidence interval:
##      -Inf 26.20367
## sample estimates:
## mean in group Control    mean in group Test 
##             100.43500              92.18595

t.test (BMP_1 ~ Category , var.equal=TRUE, alternative="less", data = subset(data, Sex == "Female"))

## 
##  Two Sample t-test
## 
## data:  BMP_1 by Category
## t = -1.2369, df = 16, p-value = 0.117
## alternative hypothesis: true difference in means between group Control and group Test is less than 0
## 95 percent confidence interval:
##      -Inf 4.777702
## sample estimates:
## mean in group Control    mean in group Test 
##              59.19129              70.80095

ttest_table <- function(data,bmp1or2){
  dstage <-degfree <-tstat <- pval <- c()
  for (stage in unique(data$Age_group)){
    tryCatch({
      if (bmp1or2){
        ttest <- t.test(BMP_1 ~ Category , var.equal=TRUE, alternative="less", data = subset(data, Age_group == stage))
      }else{
        ttest <- t.test(BMP_2 ~ Category , var.equal=TRUE, alternative="less", data = subset(data, Age_group == stage))
      }
      degfree <- c(degfree,ttest$parameter)
      tstat <- c(tstat,ttest$statistic)
      pval <- c(pval,ttest$p.value)
      dstage <- c(dstage,stage)
      }, error=function(e){})
  }
  return(data.frame(dstage,degfree,tstat,pval))
  }
kable(ttest_table(data,bmp1or2=TRUE),caption="bmp1 t-test results") 

kable(ttest_table(data,bmp1or2=FALSE),caption="bmp2 t-test results") 

female_data <- data[data$Sex == "Female",]
kable(ttest_table(female_data,bmp1or2=TRUE),caption="bmp1 t-test results") 

kable(ttest_table(female_data,bmp1or2=FALSE),caption="bmp2 t-test results") 

male_data <- data[data$Sex == "Male",]
kable(ttest_table(male_data,bmp1or2=TRUE),caption="bmp1 t-test results") 

kable(ttest_table(male_data,bmp1or2=FALSE),caption="bmp2 t-test results") 

Appendix

This section contains the QQ plots and tests of normality for each variable.

# The kinds of QQ plots for each ANOVA are

# creating Q-Q plot for each group.
QQplotter_BMP_1 <- function(dataset, s, ag){
  test <-  dataset[dataset$Category == "Test",]
  control <- dataset[dataset$Category == "Control", ]
  qqnorm(test$BMP_1, pch=1, frame=FALSE, main=sprintf("Normal Q-Q Plot for Test group, Biomarker: BMP_1, Sex: %s, Age_group: %s", s, ag))
  qqline(test$BMP_1, col = "steelblue", lwd = 2)
  qqnorm(control$BMP_1,pch=1, frame=FALSE, main=sprintf("Normal Q-Q Plot for Control group, Biomarker: BMP_1, Sex: %s, Age_group: %s", s, ag))
  qqline(control$BMP_1, col = "steelblue", lwd = 2)
}

# Total 
QQplotter_BMP_1(data, s="All", ag="All")

# By Sex 
for (sex in unique(data$Sex)){
  temp_data <-  data[data$Sex == sex,]
  QQplotter_BMP_1(temp_data,s=sex, ag="All")
}

# By Group 
for (group in unique(data$Age_group)){
  tryCatch({
    temp_data <-  data[data$Sex == group,]
  QQplotter_BMP_1(temp_data,s="All", ag=group)
  }, error=function(e){})
}
# By Group Given Sex is Male
male_data <- data[data$Sex == "Male",]
for (group in unique(male_data$Age_group)){
  tryCatch({
    temp_data <-  male_data[male_data$Sex == group,]
  QQplotter_BMP_1(temp_data,s="All", ag=group)
  }, error=function(e){})
}
# By Group Given Sex is Female 
male_data <- data[data$Sex == "Female",]
for (group in unique(male_data$Age_group)){
  tryCatch({
    temp_data <-  male_data[male_data$Sex == group,]
  QQplotter_BMP_1(temp_data,s="All", ag=group)
  }, error=function(e){})
}