DiaData <- read.delim("https://raw.githubusercontent.com/IowaBiostat/data-sets/main/diarrhea/diarrhea.txt")

boxplot(DiaData$Stool ~ DiaData$Group, 
        col = c("burlywood4", "darkgoldenrod3"), 
        ylab = "Drug Group",
        xlab = "Stool Output",
        main = "Effect of bismuth salicylate", 
        horizontal = T)

#Find n in each group
table(DiaData$Group)

#Find mean of each group
by(DiaData$Stool, DiaData$Group, mean)

#Find standard deviation of each group 
by(DiaData$Stool, DiaData$Group, sd)

x1 <- 260.2976 # Mean of control group
x2 <- 181.8706 # Mean of treatment group

sd1 <- 253.6135 # SD of control group
sd2 <- 197.3636 # SD of treatment group

n1 <- 84 # n for control group
n2 <- 85 # n for treatment group

sd_pooled <- sqrt(((n1 - 1)*sd1^2 + (n2 - 1)*sd2^2) / (n1 + n2 - 2))
SE <- sd_pooled * sqrt((1 / n1) + (1 / n2))
df <- n1 + n2 - 2

t <- (x1 - x2) / (SE)
t
2*pt(t, df = df, lower.tail = F)

(x1 - x2) + qt(c(0.025, 0.975), df = df)*SE

t.test(DiaData$Stool~DiaData$Group, var.equal = TRUE)

by(DiaData$Stool, DiaData$Group, sd)

t.test(DiaData$Stool~DiaData$Group, var.equal = FALSE)

heart <- read.delim("https://raw.githubusercontent.com/IowaBiostat/data-sets/main/infant-heart/infant-heart.txt")

by(heart$PDI, heart$Treatment, sd)
by(heart$MDI, heart$Treatment, sd)

t.test(heart$PDI ~ heart$Treatment, var.equal = T)

t.test(heart$MDI ~ heart$Treatment, var.equal = T)