Testing Assumptions, Non-parametric, Regression and Correlation
Gordon Ober
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.
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.
- Question = What type of test would be appropriate to run to compare the means of ‘WidthChange’ by ‘Tide.Height’? (explain why)
- Code = Create code to test the data for normality (for this, create a histogram AND run a Shapiro-Wilks test)
- Question = what do the tests for normality tell you about the data? Is this assumption met?
- Code = Create code to test the data for equal variance (use a Levene Test)
- Question = what does the test for equal variance tell you about the data? Is this assumption met?
- Code = log transform the WidthChange column
- Code = run a shaprio-wilks and levene test on the transformed data
- Question = what did the transformation do to your assumptions? Have they all been met?
- Question = given the outcome, what type of non-parametric test should be run?
- Code = code for the test you decided in Q9
- Question = interpret the outcome of the code from Q10
- Code = create an appropriate visual to show WidthChange by Tide.Height
barnacle <- read.csv("barnacle_data_2.csv", header = TRUE)
str(barnacle)## 'data.frame': 191 obs. of 4 variables:
## $ Tide.Height : chr "H" "H" "H" "H" ...
## $ OperculumAug: num 5.26 5.62 5.48 5.07 5.23 ...
## $ BasalAug : num 7.1 9.26 8.29 11.58 9.56 ...
## $ WidthChange : num 68.019 18.658 0.977 6.366 12.214 ...table(barnacle$Tide.Height)##
## H L
## 94 97Part 1 Answers
- A two sample t-test would be appropriate to compare WidthChange (a numerical reponse value) with Tide.Height (a categorical explanatory value), two indpenedent groups. A two sample t-test compares the means of a quantitative variable between two independent groups.
- Test Normality
hist(barnacle$WidthChange)
shapiro.test(barnacle$WidthChange)##
## Shapiro-Wilk normality test
##
## data: barnacle$WidthChange
## W = 0.4299, p-value < 2.2e-16- Both the histogram as well as the shapiro test indicate that the data is not normally distributed. This does not meet the assumption of normally distributed data.
- Levene Test
library(car)## Loading required package: carDataleveneTest(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- The p-value < 0.05, therefore we reject the null. This indicates that the data does not have equal variance and the assumption is therefore not met.
- Log transform
barnacle$logWidthChange <- log(barnacle$WidthChange)
str(barnacle)## 'data.frame': 191 obs. of 5 variables:
## $ Tide.Height : chr "H" "H" "H" "H" ...
## $ OperculumAug : num 5.26 5.62 5.48 5.07 5.23 ...
## $ BasalAug : num 7.1 9.26 8.29 11.58 9.56 ...
## $ WidthChange : num 68.019 18.658 0.977 6.366 12.214 ...
## $ logWidthChange: num 4.2198 2.9263 -0.0231 1.8509 2.5026 ...- Retest
shapiro.test(barnacle$logWidthChange)##
## Shapiro-Wilk normality test
##
## data: barnacle$logWidthChange
## W = 0.95978, p-value = 2.894e-05leveneTest(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- The transformation took the natural log of the WidthChange numerical category. While it did bring the data closer to normality, the p-value is still less than 0.05 so the assumption of normal distribution is still not met. The levene test, however, showed improvement with a p-value = 0.6358 > 0.05. Therefore we fail to reject the null and the assumption of equal variance is met.
- Based on the results from the log transformation, we should run the Wilcox rank sum test because the normality assumption is not met after log transformation. This test provides a non-parametric alternative to the two sample t-test.
- Wilcox rank sum test
wilcox.test(WidthChange ~ Tide.Height, data = barnacle)##
## Wilcoxon rank sum test with continuity correction
##
## data: WidthChange by Tide.Height
## W = 2200, p-value = 6.623e-10
## alternative hypothesis: true location shift is not equal to 0- The p-value < 0.05, therefore we reject the null hypothesis. Based on this test, the distribution of WidthChange differed between high and low Tide.Height groups.
- Visual
boxplot(WidthChange ~ Tide.Height,
data = barnacle,
main = "Barnacle Width Change by Tide Height",
xlab = "Tide Height",
ylab = "Width Change",
col = c("blue4", "lightblue"))
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
- Code: create a scatter plot of OperculumAug vs. BasalAug
- Question: what do you suspect the correlation is? Does it look positive or negative? Does it look strong or weak?
- Code: Run a correlation test on this relationship
- Question: interpret the p-value and tell me what this means for this relationship.
Part 2 Answers
- Scatter Plot
plot(x = barnacle$BasalAug,
y = barnacle$OperculumAug,
pch = 16,
col = "blue",
xlab = "Basal Width",
ylab = "Operculum Length",
main = "Operculum Length vs. Basal Width")
2. I am expecting a strong positive correlation because there seems to
be a somewhat clear trend in the data from the scatter plot. 3.
Correlation Test
cor(x = barnacle$BasalAug, y = barnacle$OperculumAug) ## [1] 0.7152863cor.test(x = barnacle$BasalAug, y = barnacle$OperculumAug, alternative = "two.sided")##
## 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- The correlation result was 0.7153, this indicates a positive correlation between the two variables, so I reject the null hypothesis that states there is no linear correlatoin. This provides evidence of a positive relationship between basal width and operculum lenght.
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
- Code: create a scatter plot of ‘radiation’ vs. ‘temperature’ (x axis should be radiation)
enviro <- read.csv("environmental.csv", header = TRUE)
str(enviro)## 'data.frame': 111 obs. of 5 variables:
## $ rownames : int 1 2 3 4 5 6 7 8 9 10 ...
## $ ozone : int 41 36 12 18 23 19 8 16 11 14 ...
## $ radiation : int 190 118 149 313 299 99 19 256 290 274 ...
## $ temperature: int 67 72 74 62 65 59 61 69 66 68 ...
## $ wind : num 7.4 8 12.6 11.5 8.6 13.8 20.1 9.7 9.2 10.9 ...plot(x = enviro$radiation,
y = enviro$temperature,
pch = 16,
col = "black",
xlab = "Radiation",
ylab = "Temperature",
main = "Relationship between Radiation and Temperature")
- Code: create and run a linear regression model for ‘radiation’ vs. ‘temperature’
linreg <- lm(temperature ~ radiation,
data = enviro)
summary(linreg)##
## Call:
## lm(formula = temperature ~ radiation, data = enviro)
##
## 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- Question: how do you interpret the p-value? How do you interpret the R^2 value?
- The p-value is 0.00173, which is less than 0.05 so, therefore, I reject the null hypothesis that the radiation slope is zero. The R square value is 0.0865, or 8.6%. This indicates that 8.6% of the variation in temperature is explalained by its relationship with radiation, a very small portion of total variation in temperature.
- Code: add a least squares regression line to your plot
plot(enviro$radiation,
enviro$temperature,
main = "Temperature vs. Radiation",
xlab = "Radiation",
ylab = "Temperature",
pch = 19)
abline(linreg,
col = "red",
lwd = 2)
- Code: create code that predicts what the temperature in NYC will be with a radiation value of 180
new_radiation <- data.frame(radiation = 180)
predict(linreg,
newdata = new_radiation)## 1
## 77.64515