The purpose of this lab is to have you investigate more complex statistical models, ones that include multiple variables!

NOTE: there is code in the hypothesis testing guide document to help create code to run these types of tests, but use this link for more detail and for other methods of graphing: https://www.datanovia.com/en/lessons/ancova-in-r/

PART 1: ANCOVA: flower characteristics

The “iris.csv” dataset is commonly used for examples in R. The dataset utilizes measurements of flower parts across different species of irises. Please download and bring in the “iris.csv” file.

For this dataset, I would like you to do the following: 1. Code: run the ANCOVA to determine how ‘variety’ and ‘sepal.width’ influence ‘sepal.length’ NOTE: don’t worry about testing assumptions for this (you would need to do this normally)

library(ggplot2)
## Warning: package 'ggplot2' was built under R version 4.4.3
iris <- read.csv("~/Desktop/BIN510-assignments/iris.csv")
model <- aov(sepal.length ~ variety + sepal.width, data = iris)
summary(model)
##              Df Sum Sq Mean Sq F value   Pr(>F)    
## variety       2  63.21  31.606   164.8  < 2e-16 ***
## sepal.width   1  10.95  10.953    57.1 4.19e-12 ***
## Residuals   146  28.00   0.192                     
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
  1. Question: how do you interpret the results of ANCOVA? The ANCOVA results show that both variety and sepal width significantly influence sepal length. Variety had a significant effect on sepal length (p < 0.001), as well as sepal width (p < 0.001). Both length and width differ amoung the iris varieties.

PART 2: Two Way ANOVA: tilling methods and plant yield

Download and bring in the “biomassTill.csv” file. This file contains data on how different soil tilling methods impact the growth (biomass) of the plants. But the researchers also manipulated the amount of water received by each plant, which would obviously influence any measurement of growth.

For this dataset, I would like you to do the following:

  1. Code: change the “DVS” variable to a factor
  2. Code: Run a two-way ANOVA that looks at how “Tillage” and “DVS” (watering regime) influence “Biomass” NOTE: don’t worry about testing assumptions for this (you would need to do this normally)
biomass <- read.csv("~/Desktop/BIN510-assignments/biomassTill.csv")
biomass$DVS <- as.factor(biomass$DVS)
model <- aov(Biomass ~ Tillage * DVS, data = biomass)
summary(model)
##             Df Sum Sq Mean Sq F value   Pr(>F)    
## Tillage      2  36447   18224   2.686   0.0796 .  
## DVS          4 711964  177991  26.238 4.71e-11 ***
## Tillage:DVS  8  34574    4322   0.637   0.7422    
## Residuals   43 291703    6784                     
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
  1. Question: how do you interpret the results of the two way ANOVA? The two-way ANOVA shows that DVS has a significant effect on biomass (p < 0.001). Also Tillage does not have a significant effect on biomass (p = 0.0796). There is no clear interaction between Tillage and DVS (p = 0.7422). This means that watering regime significantly affects biomass.
  2. Code: running a Tukey test for a two way ANOVA can be confusing if there are many iterations, the easiest way to look at comparisons is through a visual. Create a visual that shows the two way model. Note: you can use code from my example on the guide OR use code from the website linked at the top.
interaction.plot(biomass$DVS, biomass$Tillage, biomass$Biomass,
                 xlab = "Watering Regime (DVS)",
                 ylab = "Mean Biomass",
                 trace.label = "Tillage")

PART 3: Logistic Regression: Non-restorative sleep and age

Download and bring in the “sleep.csv” file. This file contains data on participants in a sleep study and documenting the occurance of non-restorative sleep. Here, data in the “NRS” column are either a 0 or a 1, 0 meaning NRS not observed, 1 meaning NRS observed.

For this dataset, I would like you to do the following:

  1. Code: run a logistic regression that looks at how “Age_2013” correlates with NRS.
sleep <- read.csv("~/Desktop/BIN510-assignments/sleep.csv")
model <- glm(NRS ~ Age_2013, data = sleep, family = binomial)
summary(model)
## 
## Call:
## glm(formula = NRS ~ Age_2013, family = binomial, data = sleep)
## 
## Coefficients:
##              Estimate Std. Error z value Pr(>|z|)    
## (Intercept)  0.868077   0.068995   12.58   <2e-16 ***
## Age_2013    -0.031965   0.001079  -29.63   <2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for binomial family taken to be 1)
## 
##     Null deviance: 98557  on 90188  degrees of freedom
## Residual deviance: 97701  on 90187  degrees of freedom
## AIC: 97705
## 
## Number of Fisher Scoring iterations: 4
  1. Question: how would you interpret the results of the logistic regression function? The logistic regression shows that age is significantly realted to NRS. The coefficient for Age2013 is negative (-0.031965) and statistically significant (p < 0.001). This can mean that as age increases the probability of experiencing NRS decreases.
  2. Code: create a logistic regression curve for the model
ggplot(sleep, aes(x = Age_2013, y = NRS)) +
  geom_jitter(height = 0.05, width = 0) +
  stat_smooth(method = "glm", method.args = list(family = "binomial"),
              se = TRUE) +
  labs(x = "Age", y = "Probability of Non-Restorative Sleep") +
  theme_classic()
## `geom_smooth()` using formula = 'y ~ x'

  1. Question: using the curve, how would you explain patterns of NRS across age? The curve shows that the probability of NRS decreases as age increases. Younger participants have a higher probability of experiencing NRS, while the older have a lower probability.