## -----------------------------------------------------------------------------
#| warning: false
library(ggplot2)
theme_set(theme_classic())
diarrhea <- read.delim("https://raw.githubusercontent.com/IowaBiostat/data-sets/main/diarrhea/diarrhea.txt")


## -----------------------------------------------------------------------------
ggplot(diarrhea, aes(x=Stool, y = Group)) +
  labs(
    title = "Effect of Bismuth Salicylate", 
    x = "Infant Stool Output (ml/kg)"
  )+
  geom_boxplot(fill = c("grey", "snow2"))


## -----------------------------------------------------------------------------
ggplot(diarrhea, aes(x=log(Stool), y = Group)) +
  labs(
    title = "Effect of Bismuth Salicylate", 
    x = "Infant Stool Output (ml/kg)\nLog Scale"
  )+
  geom_boxplot(fill = c("grey", "snow2"))


## -----------------------------------------------------------------------------
ggplot(diarrhea, aes(x = Stool)) + 
  geom_histogram(fill = "snow2", 
                 color = "black", bins = 20) + 
  facet_wrap(~Group)


## -----------------------------------------------------------------------------
ggplot(diarrhea, aes(x = log(Stool))) + 
  geom_histogram(fill = "snow2", 
                 color = "black", bins = 20) + 
  facet_wrap(~Group)


## -----------------------------------------------------------------------------
t.test(log(Stool) ~ Group, data = diarrhea)


## -----------------------------------------------------------------------------
wilcox.test(Stool ~ Group, data = diarrhea)


## -----------------------------------------------------------------------------
oatbran <- read.delim('https://raw.githubusercontent.com/IowaBiostat/data-sets/main/oatbran/oatbran.txt')
wilcox.test(oatbran$CornFlakes, oatbran$OatBran, paired=TRUE)


## -----------------------------------------------------------------------------
t.test(oatbran$CornFlakes, oatbran$OatBran, paired = TRUE)$p.value
t.test(oatbran$CornFlakes, oatbran$OatBran, paired = TRUE)$conf.int


## -----------------------------------------------------------------------------
#| layout-ncol: 2
#| warning: false
#| code-fold: true

# create an outlier and attach it to oatbran dataset
outlier <- list(CornFlakes = c(2.5), OatBran = c(6))
oatbran_outlier <- rbind(oatbran, outlier)

# generate plots with and without outlier
ggplot(oatbran, aes(x=CornFlakes, y  = OatBran)) +
  geom_point(size = 3)

ggplot(oatbran_outlier, aes(x = CornFlakes, y = OatBran)) + 
  geom_point(size = 3) +
  geom_point(aes(x=outlier$CornFlakes, y = outlier$OatBran), 
             color = "red", size = 3) +
  geom_text(
    aes(x=outlier$CornFlakes, y = outlier$OatBran - .25), 
    label = "Outlier", 
    color = "red", 
    size = 7
  )


## -----------------------------------------------------------------------------
pears <- cor.test(oatbran_outlier$CornFlakes, oatbran_outlier$OatBran, method = "pearson")
pears$estimate
pears$p.value


## -----------------------------------------------------------------------------
#| warning: false
spear <- cor.test(oatbran_outlier$CornFlakes, oatbran_outlier$OatBran, method = "spearman")
spear$estimate
spear$p.value


## -----------------------------------------------------------------------------
flicker <- read.delim('https://raw.githubusercontent.com/IowaBiostat/data-sets/main/flicker/flicker.txt')

model <- aov(Flicker ~ Color , data = flicker)
summary(model)


## -----------------------------------------------------------------------------
null_model <- aov(Flicker ~ 1, data = flicker)

# Now calculate the residual sums of squares for each model and plug them into the formula written above
RSS_0 <- sum(null_model$residuals^2)
RSS_1 <- sum(model$residuals^2)
(RSS_0 - RSS_1) / RSS_0


## -----------------------------------------------------------------------------
summary(lm(Flicker ~ Color , data = flicker))$r.squared


## -----------------------------------------------------------------------------
#| code-fold: true
ggplot(flicker, aes(Color, Flicker, fill = Color)) +
  scale_fill_manual( # this maps my chosen colors to the variable-level names
    values = c("Blue" = "steelblue3", 
               "Brown" = "brown4", 
               "Green" = "darkgreen")
    )+
  geom_boxplot() +
  guides(fill = "none") # removes my legend, since the labels are on x axis


## -----------------------------------------------------------------------------
TukeyHSD(model)


## -----------------------------------------------------------------------------
pairwise.t.test(flicker$Flicker, flicker$Color, p.adj = "bonf")


## -----------------------------------------------------------------------------
pairwise.t.test(flicker$Flicker, flicker$Color, p.adj = "none")