diff --git a/karch_06.Rmd b/karch_06.Rmd new file mode 100644 index 0000000..82e0b83 --- /dev/null +++ b/karch_06.Rmd @@ -0,0 +1,138 @@ +--- +title: 'Homework #6' +author: "Jessica Karch" +date: "October 17, 2016" +output: html_document +--- + +```{r setup, include=FALSE} +knitr::opts_chunk$set(echo = TRUE) +``` + +##19.## +**a.** Draw a scatter plot of these data. Which variable should be the explanatory variable (*X*), and wich should be the response variable (*Y*)? +```{r = 19a} +# Load libraries +library(readr) +library(dplyr) +library(ggplot2) +# Load data +plants <- read_csv("data/chap17q19GrasslandNutrientsPlantSpecies.csv") +# Inspect data +plants +# Draw a scatter plot +plants_plot <- ggplot(data = plants, mapping = aes(x = nutrients, y = species)) + + geom_point() + + labs(title = "The Effect of Different Nutrients in Fertilizer on the Number of Plant Species") +# Look at plot +plants_plot +``` + +The number of nutrients added is the explanatory variable *X* and the number of plant species is the response variable *Y*. + +**b.** What is the rate of change in the number of plant species supported per nutrient type added? Provide a standard error for your estimate. +```{r = 19b} +# Fit a linear model +plants_lm <- lm(species ~ nutrients, data = plants) +# Look at the slope, standard error +summary(plants_lm) +``` + +The rate of change is -3.3 +/- 1.1 plant species supported per nutrient type added. + +**c.** Add the least-squares regression line to your scatter plot. What fraction of the variation in the number of plant species is "explained" by the number of nutrients added? +```{r = 19c} +# Make a regression line to the scatter plot +plants_plot + + stat_smooth(method = "lm") + + theme_bw() +``` + +2/3 of the variation is "explained" by the number of nutrients added (6 of the 9 points fall within the confidence interval of the linear model). + +**d.** Test the null hypothesis of no treatment effect on the number of plant species. +```{r = 19d} +# Do an f test +anova(plants_lm) +# Look at the R^2 values, p-values +summary(plants_lm) +``` + +Our F-value is 9.24, and we have a pretty small p-value (0.016), especially taking into consideration the fact we only have 9 data points which bypasses the problem of the "obese n". Therefore we can feel fairly confident rejecting the null hypothesis that there is no treatment effect on the plant species. + +##24.## +**a.** Calculate the mean, standard deviation, and sample size of the slope for penguins in each of the three groups. Display your results in a table. +```{r = 24a} +# Load data +penguins <- read_csv("data/chap17q24PenguinTreadmill.csv") +# Inspect data +penguins +# Group data by group +penguins_table <- penguins %>% + group_by(group) %>% + # calculate the mean and standard deviation + mutate(mean = mean(slope), sd = sd(slope)) %>% + # find the number of points in a group + group_by(group, mean, sd) %>% + summarise(n = n()) %>% + ungroup() +# look at table +penguins_table +``` + +**b.** Test whether the means of the slopes are equal between the three groups. +```{r = 24b} +# Test null hypothesis, that means of the slopes are equal +# Ha: the means of the slopes are not equal +# Set up dummy vectors for each group w/ mean, sd, samp size from part a, assuming a normal distribution +# use the means and sd calculated in part a +set.seed(42) +BF <- rnorm(9, mean = 0.3155556, sd = 0.05725188) +BM <- rnorm(9, mean = 0.3355556, sd = 0.02697736) +MF <- rnorm(6, mean = 0.3033333, sd = 0.06683313) +# Set up a data frame with response variable and group +dat = data.frame(factor = c("BF","BF","BF","BF","BF","BF","BF","BF","BF","BM", "BM","BM","BM","BM","BM", "BM", "BM", "BM", "MF", "MF", "MF", "MF", "MF", "MF"), + response = c(BF, BM, MF)) +# Run ANOVA analysis +summary(aov(response ~ factor, data = dat)) +# compare each mean to each other with a t test +t.test(BF, BM) +t.test(BF, MF) +t.test(BM, MF) +``` +We have failed to reject the null hypothesis that the means are equal. The F value from our ANOVA test is very small (2.417) and the p value is large (0.114). Of the three t-tests, each of them has a 95% confidence interval that includes 0. Therefore we have failed to rject the null. +##25.## +**a.** Use these results to calculate the residuals. +```{r = 25a} +# Load data +beetles <- read_csv("data/chap17q25BeetleWingsAndHorns.csv") +# Inspect data +beetles +# Fit a linear model +beetles_lm <- lm(wingMass ~ hornSize, data = beetles) +# Calculate residuals +res_beetles <- residuals(beetles_lm) +res_beetles +``` + +**b.** Use your results from part (a) to produce a residual plot. +```{r = 25b} +# Load library +library(modelr) +# Add residuals column to main tibble +res_beetles_tb <- beetles %>% add_residuals(beetles_lm) +ggplot(data = res_beetles_tb, mapping = aes(x = hornSize, y = resid)) + + geom_point() + + labs(title = "Residual plot") + + theme_bw() +``` + +**c.** Use the graph provided and your residual plot to evaluate the main assumptions of linear regression. + +We see that the residual plot graph has a funnel effect, so that at larger relative horn sizes (>0.0) there is a wider range of residual values. The plot of the data also shows that larger relative horn sizes there is a wider range of relative wing masses. The linear regression assumes that there is a direct causal relationship between relative horn size and relative wing mass. Due to the funnel effect in both the residual plot and the main graph and the very large residuals (on the same order of magnitude as the original data), I am not sure we can assume a causal relationship between the two variables, or at least we cannot assume a linear relationship. + +**d.** In light of your conclusions in part (c), what steps should be taken? + +First, run other statistical tests, such as f test and a t test. Then try log-transforming the data, as it looks like there might be an exponential relationship in the original graph (at a certain point, the relative wing mass values get smaller very quickly, whereas the head of the graph is fairly flat). + +The statistical tests of the log-transformed data can be compared to that of the linear regression to see if that relationship would better fit the data. Then evaluate the results of the log and the linear regressions in light of what is known from literature about these beetles and evaluate what could have been the source of so much variation in the population (samples from different areas, etc.). If there is no satisfactory explanation, it might be justified to collect more data, especially because the authors only evaluated a sample size of 19. \ No newline at end of file diff --git a/karch_06.html b/karch_06.html new file mode 100644 index 0000000..eaaef3c --- /dev/null +++ b/karch_06.html @@ -0,0 +1,420 @@ + + + + + + + + + + + + + + + +Homework #6 + + + + + + + + + + + + + + + + + + + + + + + + + + +
+ + + + + + + + + + + + + + + +
+

19.

+

a. Draw a scatter plot of these data. Which variable should be the explanatory variable (X), and wich should be the response variable (Y)?

+
# Load libraries
+library(readr)
+library(dplyr)
+
## 
+## Attaching package: 'dplyr'
+
## The following objects are masked from 'package:stats':
+## 
+##     filter, lag
+
## The following objects are masked from 'package:base':
+## 
+##     intersect, setdiff, setequal, union
+
library(ggplot2)
+# Load data
+plants <- read_csv("data/chap17q19GrasslandNutrientsPlantSpecies.csv")
+
## Parsed with column specification:
+## cols(
+##   nutrients = col_integer(),
+##   species = col_integer()
+## )
+
# Inspect data
+plants
+
## # A tibble: 10 × 2
+##    nutrients species
+##        <int>   <int>
+## 1          0      36
+## 2          0      36
+## 3          0      32
+## 4          1      34
+## 5          2      33
+## 6          3      30
+## 7          1      20
+## 8          3      23
+## 9          4      21
+## 10         4      16
+
# Draw a scatter plot
+plants_plot <- ggplot(data = plants, mapping = aes(x = nutrients, y = species)) +
+  geom_point() +
+  labs(title = "The Effect of Different Nutrients in Fertilizer on the Number of Plant Species")
+# Look at plot
+plants_plot
+

+

The number of nutrients added is the explanatory variable X and the number of plant species is the response variable Y.

+

b. What is the rate of change in the number of plant species supported per nutrient type added? Provide a standard error for your estimate.

+
# Fit a linear model
+plants_lm <- lm(species ~ nutrients, data = plants)
+# Look at the slope, standard error
+summary(plants_lm)
+
## 
+## Call:
+## lm(formula = species ~ nutrients, data = plants)
+## 
+## Residuals:
+##     Min      1Q  Median      3Q     Max 
+## -10.771  -1.856   1.068   2.894   5.907 
+## 
+## Coefficients:
+##             Estimate Std. Error t value Pr(>|t|)    
+## (Intercept)   34.110      2.599   13.12 1.08e-06 ***
+## nutrients     -3.339      1.098   -3.04   0.0161 *  
+## ---
+## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
+## 
+## Residual standard error: 5.336 on 8 degrees of freedom
+## Multiple R-squared:  0.536,  Adjusted R-squared:  0.478 
+## F-statistic: 9.241 on 1 and 8 DF,  p-value: 0.01607
+

The rate of change is -3.3 +/- 1.1 plant species supported per nutrient type added.

+

c. Add the least-squares regression line to your scatter plot. What fraction of the variation in the number of plant species is “explained” by the number of nutrients added?

+
# Make a regression line to the scatter plot
+plants_plot +
+  stat_smooth(method = "lm") +
+  theme_bw()
+

+

2/3 of the variation is “explained” by the number of nutrients added (6 of the 9 points fall within the confidence interval of the linear model).

+

d. Test the null hypothesis of no treatment effect on the number of plant species.

+
# Do an f test
+anova(plants_lm)
+
## Analysis of Variance Table
+## 
+## Response: species
+##           Df Sum Sq Mean Sq F value  Pr(>F)  
+## nutrients  1 263.11 263.112  9.2406 0.01607 *
+## Residuals  8 227.79  28.474                  
+## ---
+## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
+
# Look at the R^2 values, p-values
+summary(plants_lm)
+
## 
+## Call:
+## lm(formula = species ~ nutrients, data = plants)
+## 
+## Residuals:
+##     Min      1Q  Median      3Q     Max 
+## -10.771  -1.856   1.068   2.894   5.907 
+## 
+## Coefficients:
+##             Estimate Std. Error t value Pr(>|t|)    
+## (Intercept)   34.110      2.599   13.12 1.08e-06 ***
+## nutrients     -3.339      1.098   -3.04   0.0161 *  
+## ---
+## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
+## 
+## Residual standard error: 5.336 on 8 degrees of freedom
+## Multiple R-squared:  0.536,  Adjusted R-squared:  0.478 
+## F-statistic: 9.241 on 1 and 8 DF,  p-value: 0.01607
+

Our F-value is 9.24, and we have a pretty small p-value (0.016), especially taking into consideration the fact we only have 9 data points which bypasses the problem of the “obese n”. Therefore we can feel fairly confident rejecting the null hypothesis that there is no treatment effect on the plant species.

+
+
+

24.

+

a. Calculate the mean, standard deviation, and sample size of the slope for penguins in each of the three groups. Display your results in a table.

+
# Load data
+penguins <- read_csv("data/chap17q24PenguinTreadmill.csv")
+
## Parsed with column specification:
+## cols(
+##   group = col_character(),
+##   slope = col_double()
+## )
+
# Inspect data
+penguins
+
## # A tibble: 24 × 2
+##    group slope
+##    <chr> <dbl>
+## 1     BM  0.31
+## 2     BM  0.34
+## 3     BM  0.30
+## 4     BM  0.38
+## 5     BM  0.35
+## 6     BM  0.33
+## 7     BF  0.30
+## 8     BF  0.32
+## 9     BF  0.23
+## 10    BF  0.38
+## # ... with 14 more rows
+
# Group data by group
+penguins_table <- penguins %>%
+  group_by(group) %>%
+  # calculate the mean and standard deviation
+  mutate(mean = mean(slope), sd = sd(slope)) %>%
+  # find the number of points in a group
+  group_by(group, mean, sd) %>%
+  summarise(n = n()) %>%
+  ungroup()
+# look at table
+penguins_table
+
## # A tibble: 3 × 4
+##   group      mean         sd     n
+##   <chr>     <dbl>      <dbl> <int>
+## 1    BF 0.3155556 0.05725188     9
+## 2    BM 0.3355556 0.02697736     9
+## 3    MF 0.3033333 0.06683313     6
+

b. Test whether the means of the slopes are equal between the three groups.

+
# Test null hypothesis, that means of the slopes are equal
+# Ha: the means of the slopes are not equal
+# Set up dummy vectors for each group w/ mean, sd, samp size from part a, assuming a normal distribution
+# use the means and sd calculated in part a
+set.seed(42)
+BF <- rnorm(9, mean = 0.3155556, sd = 0.05725188)
+BM <- rnorm(9, mean = 0.3355556, sd = 0.02697736)
+MF <- rnorm(6, mean = 0.3033333, sd = 0.06683313)
+# Set up a data frame with response variable and group
+dat = data.frame(factor = c("BF","BF","BF","BF","BF","BF","BF","BF","BF","BM", "BM","BM","BM","BM","BM", "BM", "BM", "BM", "MF", "MF", "MF", "MF", "MF", "MF"),
+                 response = c(BF, BM, MF))
+# Run ANOVA analysis
+summary(aov(response ~ factor, data = dat))
+
##             Df  Sum Sq  Mean Sq F value Pr(>F)
+## factor       2 0.01920 0.009601   2.417  0.114
+## Residuals   21 0.08343 0.003973
+
# compare each mean to each other with a t test
+t.test(BF, BM)
+
## 
+##  Welch Two Sample t-test
+## 
+## data:  BF and BM
+## t = 0.81371, df = 15.181, p-value = 0.4284
+## alternative hypothesis: true difference in means is not equal to 0
+## 95 percent confidence interval:
+##  -0.02739273  0.06127988
+## sample estimates:
+## mean of x mean of y 
+## 0.3507698 0.3338263
+
t.test(BF, MF)
+
## 
+##  Welch Two Sample t-test
+## 
+## data:  BF and MF
+## t = 1.5971, df = 6.5563, p-value = 0.1572
+## alternative hypothesis: true difference in means is not equal to 0
+## 95 percent confidence interval:
+##  -0.03586072  0.17897715
+## sample estimates:
+## mean of x mean of y 
+## 0.3507698 0.2792116
+
t.test(BM, MF)
+
## 
+##  Welch Two Sample t-test
+## 
+## data:  BM and MF
+## t = 1.2507, df = 5.9679, p-value = 0.2578
+## alternative hypothesis: true difference in means is not equal to 0
+## 95 percent confidence interval:
+##  -0.05237394  0.16160322
+## sample estimates:
+## mean of x mean of y 
+## 0.3338263 0.2792116
+

We have failed to reject the null hypothesis that the means are equal. The F value from our ANOVA test is very small (2.417) and the p value is large (0.114). Of the three t-tests, each of them has a 95% confidence interval that includes 0. Therefore we have failed to rject the null. ##25.##
+a. Use these results to calculate the residuals.

+
# Load data
+beetles <- read_csv("data/chap17q25BeetleWingsAndHorns.csv")
+
## Parsed with column specification:
+## cols(
+##   hornSize = col_double(),
+##   wingMass = col_double()
+## )
+
# Inspect data
+beetles
+
## # A tibble: 19 × 2
+##    hornSize wingMass
+##       <dbl>    <dbl>
+## 1     0.074    -42.8
+## 2     0.079    -21.7
+## 3     0.019    -18.8
+## 4     0.017    -16.0
+## 5     0.085    -12.8
+## 6     0.081     11.6
+## 7     0.011      7.6
+## 8     0.023      1.6
+## 9     0.005      3.7
+## 10    0.007      1.1
+## 11    0.004     -0.8
+## 12   -0.002     -2.9
+## 13   -0.065     12.1
+## 14   -0.065     20.1
+## 15   -0.014     21.2
+## 16   -0.014     22.2
+## 17   -0.132     20.1
+## 18   -0.143     12.5
+## 19   -0.177      7.0
+
# Fit a linear model
+beetles_lm <- lm(wingMass ~ hornSize, data = beetles)
+# Calculate residuals
+res_beetles <- residuals(beetles_lm)
+res_beetles
+
##           1           2           3           4           5           6 
+## -32.8573453 -11.0942627 -16.1512536 -13.6164867  -1.3985636  22.4709703 
+##           7           8           9          10          11          12 
+##   9.1878142   4.7792124   4.4921151   2.1573482  -0.1405014  -3.0362005 
+##          13          14          15          16          17          18 
+##   3.6089591  11.6089591  19.4724013  20.4724013   2.7236525  -6.3351291 
+##          19 
+## -16.3440907
+

b. Use your results from part (a) to produce a residual plot.

+
# Load library
+library(modelr)
+# Add residuals column to main tibble
+res_beetles_tb <- beetles %>% add_residuals(beetles_lm)
+ggplot(data = res_beetles_tb, mapping = aes(x = hornSize, y = resid)) +
+  geom_point() +
+  labs(title = "Residual plot") +
+  theme_bw()
+

+

c. Use the graph provided and your residual plot to evaluate the main assumptions of linear regression.

+

We see that the residual plot graph has a funnel effect, so that at larger relative horn sizes (>0.0) there is a wider range of residual values. The plot of the data also shows that larger relative horn sizes there is a wider range of relative wing masses. The linear regression assumes that there is a direct causal relationship between relative horn size and relative wing mass. Due to the funnel effect in both the residual plot and the main graph and the very large residuals (on the same order of magnitude as the original data), I am not sure we can assume a causal relationship between the two variables, or at least we cannot assume a linear relationship.

+

d. In light of your conclusions in part (c), what steps should be taken?

+

First, run other statistical tests, such as f test and a t test. Then try log-transforming the data, as it looks like there might be an exponential relationship in the original graph (at a certain point, the relative wing mass values get smaller very quickly, whereas the head of the graph is fairly flat).

+

The statistical tests of the log-transformed data can be compared to that of the linear regression to see if that relationship would better fit the data. Then evaluate the results of the log and the linear regressions in light of what is known from literature about these beetles and evaluate what could have been the source of so much variation in the population (samples from different areas, etc.). If there is no satisfactory explanation, it might be justified to collect more data, especially because the authors only evaluated a sample size of 19.

+
+ + + + +
+ + + + + + + +