title: “BIO 415/514 Module 4 Report” author: “Alaina Breakiron” output: html_document —
R Markdown
This is an R Markdown document. Markdown is a simple formatting syntax for authoring HTML, PDF, and MS Word documents. For more details on using R Markdown see http://rmarkdown.rstudio.com.
When you click the Knit button a document will be generated that includes both content as well as the output of any embedded R code chunks within the document. R code chunks look like this:
# This is a comment line, which we will often use as a question prompt.
# Do not change the name of the chunk because it is used when grading.
note <- "This is a line of R code that sets a variable to a value"
note # This line will print the value to the knitted document.## [1] "This is a line of R code that sets a variable to a value"We will also use multiple choice questions within this assignment.
You will use an X to indicate you choice, as shown
below.
Which of the following is a shade of red?
Gradescope Unit Tests
Gradescope is the website that we will be using to grade your RMarkdown exercises, and you will access it via Canvas. To submit, you will upload your RMarkdown file to Gradescope, and Gradescope will run your file and its output through many different unit tests to validate your code. Each unit test will check something specific about your code and if something fails, you will received a report about the failure. You may submit as many times as you want up until the deadline for the assignment.
There is a learning curve attached to Gradescope, but once you get the hang of using unit testing to fix your code, it will improve your productivity. Due to limitations of the unit tests, some unit tests will be difficult to pass on first submission, but in these cases the failure messages will help you fix your code to match what the unit test expects.
Unit testing is new to BIO 415/514 this semester. There may be errors with the assignments or the unit tests, please reach out to us with any issues or feedback. We have tried our best to make the unit test output understandable by novice R uses, if you have problems understanding your unit test reports, please reach out to us for help.
Gradescope Tips
- Treat unit testing as a form of real-time feedback that allows you to improve your code.
- Pay attention to the failure reports, as they often contain information about what you need to fix.
- Do not upload data files, as they will be available on available in the working directory on Gradscope’s server.
- Do not call
setwd()in your RMarkdown, as your local directory structure will not match the server’s directory structure. - Do not try to install packages or other interactive operations as Gradescope’s interface does not support that.
Introduction
In this assignment, you will test hypotheses using three data sets, described below. For each data set, you will evaluate whether the data meet the assumptions for parametric tests, transform the data if necessary, and carry out an appropriate test. The three different data sets all require a different approach: one is best analyzed with a Mann-Whitney test, one with a Welch’s t-test, and one with an ordinary two-sample t-test after log transformation. Your task is to figure out which approach is most appropriate for which data set, explain why, and perform the analysis.
For each data set, do the following:
- Make a box plot showing the data from the two groups that are being compared.
- Test the assumption of normality, both graphically (using normal
probability plots) and through hypothesis tests (Shapiro-Wilk).
- Test the assumption of equal variance by carrying out Levene’s test.
- If the assumptions are adequately met for a Welch’s t-test, then perform it in R.
- If the assumptions are not adequately met for any kind of t-test, try log transforming the data. Then re-evaluate to see if it now conforms to the assumptions well enough for a t-test. If so, then perform the test in R.
- If log transformation does not yield data that meet t-test assumptions, determine whether the assumptions of the Mann-Whitney test are met. If so, perform this test in R.
- Make a clear and concise statement of the result of your hypothesis test.
Question 1: Carotenoids and immune function.
Carotenoids are pigments responsible for much of the red we see in nature, including the bright red beak color of the male zebra finch. Carotenoids also act as antioxidants, suggesting that carotenoids in birds may affect immune system function. To test this, a group of zebra finches was randomly divided into two groups. Ten finches received supplemental carotenoids in their diet, and ten individuals did not. All 20 birds were then measured using an assay for cell-mediated immunocompetence (PHA). The data are given in file “finches.csv”. Use these results to test whether there is a difference in PHA, and thus in immune system function, between the two groups.
# Read the data into R. You can also use this chunk as a scratch area to explore
# your data.
finch.data <- read.csv("finches.csv")
finch.data## diet PHA
## 1 carot 2.3
## 2 carot 2.1
## 3 carot 2.7
## 4 carot 2.6
## 5 carot 2.2
## 6 carot 2.9
## 7 carot 3.0
## 8 carot 2.5
## 9 carot 2.4
## 10 carot 2.8
## 11 no.carot 1.2
## 12 no.carot 1.4
## 13 no.carot 1.3
## 14 no.carot 1.1
## 15 no.carot 1.5
## 16 no.carot 1.6
## 17 no.carot 1.0
## 18 no.carot 1.3
## 19 no.carot 1.2
## 20 no.carot 1.4# Generate a box-and-whisker plot of birds with and without supplemental
# carotenoids.
boxplot(PHA ~ diet, data = finch.data,
xlab = "Diet Group",
ylab = "PHA (immune response)",
main = "PHA by Carotenoid Supplementation")
carot_group <- finch.data$PHA[finch.data$diet == "carot"]
z_carot <- scale(carot_group)[, 1]
carot_df <- data.frame(x = z_carot, y = carot_group)
library(ggplot2)
plot_carot <- ggplot(carot_df, aes(x = x, y = y)) +
geom_point() +
labs(x = "Standardized PHA (Carot)", y = "Raw PHA")
question_1c <- list(value = plot_carot)
question_1c## $value
no_carot_group <- finch.data$PHA[finch.data$diet == "no.carot"]
z_no_carot <- scale(no_carot_group)[, 1]
nocarot_df <- data.frame(x = z_no_carot, y = no_carot_group)
plot_nocarot <- ggplot(nocarot_df, aes(x = x, y = y)) +
geom_point() +
labs(x = "Standardized PHA (No Carot)", y = "Raw PHA")
question_1d <- list(value = plot_nocarot)
question_1d## $value
# Perform a Shapiro-Wilk test for birds with supplemental carotenoids.
carot_group <- finch.data$PHA[finch.data$diet == "carot"]
shapiro.test(carot_group)##
## Shapiro-Wilk normality test
##
## data: carot_group
## W = 0.97016, p-value = 0.8924# Perform a Shapiro-Wilk test for birds without supplemental carotenoids.
no_carot_group <- finch.data$PHA[finch.data$diet == "no.carot"]
shapiro.test(no_carot_group)##
## Shapiro-Wilk normality test
##
## data: no_carot_group
## W = 0.98372, p-value = 0.9819# Test for equal variances between birds with and without supplemental
# carotenoids.
library(car)
leveneTest(PHA ~ diet, data = finch.data)## Warning in leveneTest.default(y = y, group = group, ...): group coerced to
## factor.## Levene's Test for Homogeneity of Variance (center = median)
## Df F value Pr(>F)
## group 1 3.5822 0.0746 .
## 18
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1# If necessary, use this chunk to evaluate to if log transformed data conforms
# well enough to the assumptions for a t-test.
# Create log-transformed PHA column
finch.data$logPHA <- log(finch.data$PHA)
# Subset the logPHA values by diet group
logPHA_carot <- subset(finch.data, diet == "carot")$logPHA
logPHA_nocarot <- subset(finch.data, diet == "no.carot")$logPHA
# Plot QQ plots
par(mfrow = c(1, 2)) # side-by-side plots
qqnorm(logPHA_carot, main = "QQ Plot: logPHA (Carotenoids)")
qqline(logPHA_carot)
qqnorm(logPHA_nocarot, main = "QQ Plot: logPHA (No Carotenoids)")
qqline(logPHA_nocarot)
par(mfrow = c(1, 1)) # reset layoutfinch.data$logPHA <- log(finch.data$PHA)
t_res <- t.test(logPHA ~ diet, data = finch.data)
question_1i <- list(value = list(
method = t_res$method,
statistic = unname(t_res$statistic),
p.value = t_res$p.value
))
question_1i## $value
## $value$method
## [1] "Welch Two Sample t-test"
##
## $value$statistic
## [1] 11.47919
##
## $value$p.value
## [1] 1.44492e-09Summary
Which sentence provides the most accurate conclusion from the results of your hypothesis test? (Fill in the blanks.)
At a 0.05 level of significance, I conclude that immunocompetence (PHA level) is {QUESTION 1J} in birds fed supplemental carotenoids than in birds that were not ({QUESTION 1K}).
Question 1j
Question 1k
Question 2: Tobacco smoke exposure
Exposure to tobacco smoke can be measured by the urinary cotinine-to-creatinine ratio (CCR). (Cotinine is formed in the body by the breakdown of nicotine.) Scientists measured this in infants from smoking households. The houses were divided into two groups, based on their previous behavior: ones with strict controls to prevent exposure of the infant to smoke, and ones with looser controls. Their data are in the file “smoking.csv”. You must test whether CCR levels differ between the two types of household.
# Read the data into R. You can also use this chunk as a scratch area to explore
# your data.
smoking.data <- read.csv("smoking.csv")
smoking.data## behavior CCR
## 1 loose 30.2
## 2 loose 45.2
## 3 loose 28.9
## 4 loose 60.3
## 5 loose 40.1
## 6 strict 10.3
## 7 strict 12.2
## 8 strict 14.5
## 9 strict 9.9
## 10 strict 13.7# Generate a box-and-whisker plot of infant CCR in both types of households.
boxplot(CCR ~ behavior, data = smoking.data,
xlab = "Household Smoking Behavior",
ylab = "Cotinine-Creatinine Ratio (CCR)",
main = "Infant CCR by Household Smoking Behavior")
loose_ccr <- smoking.data$CCR[smoking.data$behavior == "loose"]
z_loose <- scale(loose_ccr)[,1]
loose_df <- data.frame(x = z_loose, y = loose_ccr)
plot_loose <- ggplot(loose_df, aes(x = x, y = y)) +
geom_point() +
labs(x = "Standardized CCR (Loose)", y = "Raw CCR")
question_2c <- list(value = plot_loose)
question_2c## $value
strict_ccr <- smoking.data$CCR[smoking.data$behavior == "strict"]
z_strict <- scale(strict_ccr)[,1]
strict_df <- data.frame(x = z_strict, y = strict_ccr)
plot_strict <- ggplot(strict_df, aes(x = x, y = y)) +
geom_point() +
labs(x = "Standardized CCR (Strict)", y = "Raw CCR")
question_2d <- list(value = plot_strict)
question_2d## $value
loose_ccr <- smoking.data$CCR[smoking.data$behavior == "loose"]
shapiro_loose <- shapiro.test(loose_ccr)
question_2e <- list(value = shapiro.test(loose_ccr))
question_2e## $value
##
## Shapiro-Wilk normality test
##
## data: loose_ccr
## W = 0.91711, p-value = 0.5115strict_ccr <- smoking.data$CCR[smoking.data$behavior == "strict"]
shapiro_strict <- shapiro.test(strict_ccr)
question_2f <- list(value = shapiro.test(strict_ccr))
question_2f## $value
##
## Shapiro-Wilk normality test
##
## data: strict_ccr
## W = 0.91654, p-value = 0.5079# Test for equal variances in infant CCRR in both types of households.
library(car)
leveneTest(CCR ~ behavior, data = smoking.data)## Warning in leveneTest.default(y = y, group = group, ...): group coerced to
## factor.## Levene's Test for Homogeneity of Variance (center = median)
## Df F value Pr(>F)
## group 1 5.1167 0.05355 .
## 8
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1# If necessary, use this chunk to evaluate to if log transformed data conforms
# well enough to the assumptions for a t-test.
# Create a new column for log-transformed CCR
smoking.data$logCCR <- log(smoking.data$CCR)
# Subset log-transformed CCR by group
log_loose <- smoking.data$logCCR[smoking.data$behavior == "loose"]
log_strict <- smoking.data$logCCR[smoking.data$behavior == "strict"]
# Check normality with QQ plots
qqnorm(log_loose, main = "QQ Plot: logCCR in Loose Households")
qqline(log_loose)
qqnorm(log_strict, main = "QQ Plot: logCCR in Strict Households")
qqline(log_strict)
# Optional: Shapiro-Wilk test on log-transformed data
shapiro.test(log_loose)##
## Shapiro-Wilk normality test
##
## data: log_loose
## W = 0.93611, p-value = 0.6386shapiro.test(log_strict)##
## Shapiro-Wilk normality test
##
## data: log_strict
## W = 0.91246, p-value = 0.4825# Optional: Levene's test on log-transformed data
leveneTest(logCCR ~ behavior, data = smoking.data)## Warning in leveneTest.default(y = y, group = group, ...): group coerced to
## factor.## Levene's Test for Homogeneity of Variance (center = median)
## Df F value Pr(>F)
## group 1 1.3126 0.285
## 8# Perform the appropriate test between the distributions.
t.test(logCCR ~ behavior, data = smoking.data, var.equal = TRUE)##
## Two Sample t-test
##
## data: logCCR by behavior
## t = 7.6808, df = 8, p-value = 5.847e-05
## alternative hypothesis: true difference in means between group loose and group strict is not equal to 0
## 95 percent confidence interval:
## 0.8335442 1.5487877
## sample estimates:
## mean in group loose mean in group strict
## 3.674698 2.483532Summary
Which sentence provides the most accurate conclusion from the results of your hypothesis test? (Fill in the blanks.)
At a 0.05 level of significance, I conclude that mean CCR level is {QUESTION 2J} in households with loose smoking controls than in households with stict controls ({QUESTION 2K}).
Question 2j
Question 2k
Question 3: Wolbachia and dengue fever
Dengue fever is caused by an RNA virus that is transmitted by mosquitoes. Infection by the bacterium Wolbachia confers immunity to RNA viruses, at least in fruit flies. Researchers are trying to determine if a similar infection could be spread among mosquitoes, and thus prevent them from acquiring and spreading the dengue fever virus. They created a strain of mosquitoes (called WB1) infected with Wolbachia, which is normally not found in mosquitoes. Then they infected a random sample of these WB1 mosquitoes with dengue fever virus, and did the same to a random sample of wild type mosquitoes that lacked Wolbachia. Fourteen days later, they allowed the mosquitoes to feed on artificial food solution, and then they measured dengue virus titers in the food. The results are given in file “dengue.csv”. Use these data to determine whether WB1 differs from the wild strain in virus titer.
# Read the data into R. You can also use this chunk as a scratch area to explore
# your data.
dengue.data <- read.csv("dengue.csv")
dengue.data## strain titer
## 1 WB1 100
## 2 WB1 120
## 3 WB1 95
## 4 WB1 105
## 5 WB1 90
## 6 WB1 130
## 7 WB1 125
## 8 WB1 140
## 9 Wild 110
## 10 Wild 115
## 11 Wild 108
## 12 Wild 112
## 13 Wild 117
## 14 Wild 113
## 15 Wild 111
## 16 Wild 116# Generate a box-and-whisker plot of virus titer in both strains.
boxplot(titer ~ strain, data = dengue.data,
main = "Dengue Virus Titers by Mosquito Strain",
xlab = "Mosquito Strain",
ylab = "Virus Titer",
col = c("skyblue", "orange"),
border = "black")
wb1_titer <- dengue.data$titer[dengue.data$strain == "WB1"]
z_wb1 <- scale(wb1_titer)[,1]
wb1_df <- data.frame(x = z_wb1, y = wb1_titer)
plot_wb1 <- ggplot(wb1_df, aes(x = x, y = y)) +
geom_point() +
labs(x = "Standardized Titer (WB1)", y = "Raw Titer")
question_3c <- list(value = plot_wb1)
question_3c## $value
wild_titer <- dengue.data$titer[dengue.data$strain == "Wild"]
z_wild <- scale(wild_titer)[,1]
wild_df <- data.frame(x = z_wild, y = wild_titer)
plot_wild <- ggplot(wild_df, aes(x = x, y = y)) +
geom_point() +
labs(x = "Standardized Titer (Wild)", y = "Raw Titer")
question_3d <- list(value = plot_wild)
question_3d## $value
# Perform a Shapiro-Wilk test for virus titer in WB1 mosquitoes.
# Subset virus titer values for WB1 mosquitoes (case-sensitive match)
wb1.titer <- na.omit(dengue.data$titer[dengue.data$strain == "WB1"])
# Perform Shapiro-Wilk test
shapiro.test(wb1.titer)##
## Shapiro-Wilk normality test
##
## data: wb1.titer
## W = 0.94336, p-value = 0.6444# Perform a Shapiro-Wilk test for virus titer in wild mosquitoes.
# Subset virus titer values for Wild mosquitoes (case-sensitive match)
wild.titer <- na.omit(dengue.data$titer[dengue.data$strain == "Wild"])
# Perform Shapiro-Wilk test
shapiro.test(wild.titer)##
## Shapiro-Wilk normality test
##
## data: wild.titer
## W = 0.97207, p-value = 0.9137library(car)
levene_result <- leveneTest(titer ~ strain, data = dengue.data)## Warning in leveneTest.default(y = y, group = group, ...): group coerced to
## factor.question_3g <- list(value = levene_result)
question_3g## $value
## Levene's Test for Homogeneity of Variance (center = median)
## Df F value Pr(>F)
## group 1 26.439 0.0001497 ***
## 14
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1# If necessary, use this chunk to evaluate to if log transformed data conforms
# well enough to the assumptions for a t-test.
# Log-transform virus titer data (adding a small constant if needed to avoid log(0))
dengue.data$log_titer <- log10(dengue.data$titer + 1e-6) # Add tiny value to avoid log(0)
# Subset log-transformed titers
wb1.log <- na.omit(dengue.data$log_titer[dengue.data$strain == "WB1"])
wild.log <- na.omit(dengue.data$log_titer[dengue.data$strain == "Wild"])
# Q-Q plots to check normality
par(mfrow = c(1, 2)) # Side-by-side plots
qqnorm(wb1.log, main = "Q-Q Plot (Log Titer, WB1)")
qqline(wb1.log, col = "red")
qqnorm(wild.log, main = "Q-Q Plot (Log Titer, Wild)")
qqline(wild.log, col = "blue")
# Shapiro-Wilk tests for normality
shapiro.wb1.log <- shapiro.test(wb1.log)
shapiro.wild.log <- shapiro.test(wild.log)
# F-test for equal variances
var.test(wb1.log, wild.log)##
## F test to compare two variances
##
## data: wb1.log and wild.log
## F = 33.958, num df = 7, denom df = 7, p-value = 0.0001394
## alternative hypothesis: true ratio of variances is not equal to 1
## 95 percent confidence interval:
## 6.798511 169.616845
## sample estimates:
## ratio of variances
## 33.95794dengue.data$log_titer <- log10(dengue.data$titer + 1e-6)
wb1_log <- dengue.data$log_titer[dengue.data$strain == "WB1"]
wild_log <- dengue.data$log_titer[dengue.data$strain == "Wild"]
t_result <- t.test(wb1_log, wild_log, var.equal = TRUE)
question_3i <- list(value = t.test(log_titer ~ strain, data = dengue.data, var.equal = TRUE))
question_3i## $value
##
## Two Sample t-test
##
## data: log_titer by strain
## t = -0.13212, df = 14, p-value = 0.8968
## alternative hypothesis: true difference in means between group WB1 and group Wild is not equal to 0
## 95 percent confidence interval:
## -0.05699946 0.05038458
## sample estimates:
## mean in group WB1 mean in group Wild
## 2.048665 2.051972Summary
Which sentence provides the most accurate conclusion from the results of your hypothesis test? (Fill in the blanks.)
At a 0.05 level of significance, I conclude that the wild strain has {QUESTION 3J} dengue titer than the WB1 strain ({QUESTION 3K}).