BIN510 R Lab 7

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.

irises <- read.csv("iris.csv", header = TRUE)

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)

ancovaIris <- lm(sepal.length ~ variety + sepal.width, data = irises)
summary(ancovaIris)

## 
## Call:
## lm(formula = sepal.length ~ variety + sepal.width, data = irises)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -1.30711 -0.25713 -0.05325  0.19542  1.41253 
## 
## Coefficients:
##                   Estimate Std. Error t value Pr(>|t|)    
## (Intercept)         2.2514     0.3698   6.089 9.57e-09 ***
## varietyVersicolor   1.4587     0.1121  13.012  < 2e-16 ***
## varietyVirginica    1.9468     0.1000  19.465  < 2e-16 ***
## sepal.width         0.8036     0.1063   7.557 4.19e-12 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.438 on 146 degrees of freedom
## Multiple R-squared:  0.7259, Adjusted R-squared:  0.7203 
## F-statistic: 128.9 on 3 and 146 DF,  p-value: < 2.2e-16

anova(ancovaIris)

## Analysis of Variance Table
## 
## Response: sepal.length
##              Df Sum Sq Mean Sq F value    Pr(>F)    
## variety       2 63.212 31.6061 164.781 < 2.2e-16 ***
## sepal.width   1 10.953 10.9525  57.102 4.187e-12 ***
## Residuals   146 28.004  0.1918                      
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

2. Question: how do you interpret the results of ANCOVA?

• The ANCOVA reports that both variety and sepal.width were significantly associated with sepal length (p < 0.001). Therefore, both of these variables significantly influence sepal.length.

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.

biomass <- read.csv("biomassTill.csv", header = TRUE)

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

1. Code: change the “DVS” variable to a factor

biomass$DVS <- as.factor(
  biomass$DVS
)

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)

twoway_biomass <- aov(Biomass~Tillage*DVS, data = biomass)
summary(twoway_biomass)

##             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

3. Question: how do you interpret the results of the two way ANOVA?

• We fail to reject the null hypothesis for tillage as there is no evidence of tillage unfluencing biomass in this model. However, there is strong evidence that biomass differs among DVS watering regimes (p < 0.001).

4. 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.

library(ggplot2)

ggplot(data = biomass,
       aes(x = Tillage,
           y = Biomass)) +
  geom_boxplot(
    aes(fill = DVS)
  ) 

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.

sleep <- read.csv("sleep.csv", header = TRUE)

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.

log_model <- glm(
  NRS ~ Age_2013, 
  family = binomial(link='logit'),
  data = sleep)
summary(log_model)

## 
## Call:
## glm(formula = NRS ~ Age_2013, family = binomial(link = "logit"), 
##     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

anova(log_model, test = "Chisq")

## Analysis of Deviance Table
## 
## Model: binomial, link: logit
## 
## Response: NRS
## 
## Terms added sequentially (first to last)
## 
## 
##          Df Deviance Resid. Df Resid. Dev  Pr(>Chi)    
## NULL                     90188      98557              
## Age_2013  1    855.7     90187      97701 < 2.2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

2. Question: how would you interpret the results of the logistic regression function?

• The logistic regression illustrates a significant relationship between age and NGS (p < 0.001).

3. Code: create a logistic regression curve for the model

ggplot(
  data = sleep,
  aes(
    x = Age_2013,
    y = NRS
  )
) +
  geom_point() +
  stat_smooth(
    method = "glm",
    method.args = list(
      family = "binomial"
    ),
    se = FALSE
  )

## `geom_smooth()` using formula = 'y ~ x'

4. Question: using the curve, how would you explain patterns of NRS across age? - The figure shows that the predicted probability of NRS decreases as age increases. Therefore, younger participants are more likely to exhibit NRS than older patients.