The purpose of this R lab is to have you practice testing assumptions, interpretting output, practicing ways to deal with violations of assumptions and to dip into regression and correlation

NOTE: please refer to the hypothesis testing guide document for code examples!

Part 1: Testing assumptions

Create a code chunk(s) below to bring in “barnacle_data.csv” and complete the following actions. There are questions that will follow your coding. This data comes from some of my research looking at how things like microhabitat and shore height influence the growth and abundance of barnacles.

barnacle <- read.csv("~/Desktop/BIN510-files/barnacle_data_2.csv")

Actions/Questions: You are interested in seeing if there is an effect of shore height (‘Tide.Height’) on the growth of individual barnacles over a 6 month period.

  1. Question = What type of test would be appropriate to run to compare the means of ‘WidthChange’ by ‘Tide.Height’? (explain why) Answer: A one-way ANOVA would be best to compare the WidthChange (numerical) and Tide.Height(categorical) because one is a numerical variable and the other is a categorical value.

  2. Code = Create code to test the data for normality (for this, create a histogram AND run a Shapiro-Wilks test)

hist(barnacle$WidthChange,
     main = "Histogram of Barnacle WidthChange",
     xlab = "WidthChange",
     col= c("red","orange","yellow","limegreen","skyblue"))

shapiro.test(barnacle$WidthChange)
## 
##  Shapiro-Wilk normality test
## 
## data:  barnacle$WidthChange
## W = 0.4299, p-value < 2.2e-16
  1. Question = what do the tests for normality tell you about the data? Is this assumption met? Answer: Since the p-vaule was 2.2e-16 this means that the data is signficantly different from normal.

  2. Code = Create code to test the data for equal variance (use a Levene Test)

#install.packages("car")
library(car)
## Warning: package 'car' was built under R version 4.4.3
## Loading required package: carData
## Warning: package 'carData' was built under R version 4.4.3
leveneTest(WidthChange ~ Tide.Height, data = barnacle)
## 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.1557 0.0243 *
##       189                 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
  1. Question = what does the test for equal variance tell you about the data? Is this assumption met? Answer: The Levene test had 0.0243 for a p-value, this means that the variances differ signficantly in the Tide.Height groups.

  2. Code = log transform the WidthChange column

summary(barnacle$WidthChange)
##     Min.  1st Qu.   Median     Mean  3rd Qu.     Max. 
##   0.0688   6.3909  13.4567  22.4644  26.3543 447.3620
barnacle$logWidthChange <- log(barnacle$WidthChange)
hist(barnacle$logWidthChange,
     main = "Histogram of Log Width Change",
     xlab = "Log Width Change",
      col= c("red","orange","yellow","limegreen","skyblue","royalblue","violet","purple"))

  1. Code = run a shaprio-wilks and levene test on the transformed data
shapiro.test(barnacle$logWidthChange)
## 
##  Shapiro-Wilk normality test
## 
## data:  barnacle$logWidthChange
## W = 0.95978, p-value = 2.894e-05
leveneTest(logWidthChange ~ Tide.Height, data = barnacle)
## 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   0.225 0.6358
##       189
  1. Question = what did the transformation do to your assumptions? Have they all been met? Answer: The log transformation improved the equal variance assumption. The Levene test changed from 0.0243 to 0.6358.The Shapiro-Wilk test was still significant (p = 2.894e-05), meaning the data were still not normally distributed.
  2. Question = given the outcome, what type of non-parametric test should be run? Answer:A Kruskal-Wallis test should be run.
  3. Code = code for the test you decided in Q9
kruskal.test(WidthChange ~ Tide.Height, data = barnacle)
## 
##  Kruskal-Wallis rank sum test
## 
## data:  WidthChange by Tide.Height
## Kruskal-Wallis chi-squared = 38.145, df = 1, p-value = 6.568e-10
  1. Question = interpret the outcome of the code from Q10 Answer: Since the p-value is much smaller than 0.05, this indicates there is significant difference in WidthChange and Tide.Height.

  2. Code = create an appropriate visual to show WidthChange by Tide.Height

boxplot(WidthChange ~ Tide.Height,
        data = barnacle,
        main = "Barnacle Width Change by Tide Height",
        xlab = "Tide Height",
        ylab = "Width Change",
        col= c("pink","gold"))


Part 2: Correlation

The barnacle data set had measured key size metric for this species over time. In Part 1, you explored how tide height would impact the growht of one of these metrics, basal width. However, we collected two types of size data on the barnacles, their basal width and their operculum length. Much like we see body and brain size correlate, do we see basal width and operculum length do the same thing?

Actions/Questions: Complete the following using one or several code chunks with the ‘barnacle_data_2.csv’ file

  1. Code: create a scatter plot of OperculumAug vs. BasalAug
plot(barnacle$BasalAug,
     barnacle$OperculumAug,
     main = "Operculum Length versus Basal Width",
     xlab = "Basal Width",
     ylab = "Operculum Length")

  1. Question: what do you suspect the correlation is? Does it look positive or negative? Does it look strong or weak? Answer: It seems like there is a positive correlation between basal width and operculum length. The relationship seems strong since it seems to be continuing to rise in a pattern.
  2. Code: Run a correlation test on this relationship
cor.test(barnacle$BasalAug,
         barnacle$OperculumAug)
## 
##  Pearson's product-moment correlation
## 
## data:  barnacle$BasalAug and barnacle$OperculumAug
## t = 14.071, df = 189, p-value < 2.2e-16
## alternative hypothesis: true correlation is not equal to 0
## 95 percent confidence interval:
##  0.6381109 0.7782316
## sample estimates:
##       cor 
## 0.7152863
  1. Question: interpret the p-value and tell me what this means for this relationship. Answer: The correlation showed a positive and strong relationship between basal width and operculum length, this means that as one measurment gets bigger so does the other. The P-value was less than 0.05.

Part 3: Regression

Load the data ‘environmental.csv’, this dataset has tracked environmental parameters in NYC throughout the course of a summer. You are interested developing a predictive model that tracks how radiation influences temperature. This would allow you to better hone in on weather forecasting.

Actions/Questions: Complete the following using one or several code chunks

  1. Code: create a scatter plot of ‘radiation’ vs. ‘temperature’ (x axis should be radiation)
environmental <- read.csv("~/Desktop/BIN510-files/environmental.csv")
plot(environmental$radiation,
     environmental$temperature,
     main = "Radiation versus Temperature",
     xlab = "Radiation",
     ylab = "Temperature")

  1. Code: create and run a linear regression model for ‘radiation’ vs. ‘temperature’
model <- lm(temperature ~ radiation, data = environmental)

summary(model)
## 
## Call:
## lm(formula = temperature ~ radiation, data = environmental)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -19.735  -6.292   1.080   6.231  18.648 
## 
## Coefficients:
##              Estimate Std. Error t value Pr(>|t|)    
## (Intercept) 72.110720   1.970502  36.595  < 2e-16 ***
## radiation    0.030747   0.009571   3.212  0.00173 ** 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 9.15 on 109 degrees of freedom
## Multiple R-squared:  0.08649,    Adjusted R-squared:  0.07811 
## F-statistic: 10.32 on 1 and 109 DF,  p-value: 0.001731
  1. Question: how do you interpret the p-value? How do you interpret the R^2 value? Answer: The p-value for radiation was 0.00173, which is less than 0.05. This shows that radiation is statistically significant of temperature. The R2 value was 0.08649, which means that approximately 8.65% of the variation in temperature can be explained by radiation.
  2. Code: add a least squares regression line to your plot
plot(environmental$radiation,
     environmental$temperature,
     main = "Radiation versus Temperature",
     xlab = "Radiation",
     ylab = "Temperature")

abline(model)

  1. Code: create code that predicts what the temperature in NYC will be with a radiation value of 180
predict(model, newdata = data.frame(radiation = 180))
##        1 
## 77.64515