R Notebook

library(readr)
library(ggplot2)
library(dplyr)
library(lattice)
library(car)
library(knitr)
library(car)
library(cowplot)
library(datasets)

barley <- read_delim("http://sciences.ucf.edu/biology/d4lab/wp-content/uploads/sites/125/2017/01/barley.txt", "\t", escape_double = FALSE, trim_ws = TRUE)
zinc <- read_delim("http://sciences.ucf.edu/biology/d4lab/wp-content/uploads/sites/125/2017/01/zinc.txt", "\t", escape_double = FALSE, trim_ws = TRUE)
possum <- read_delim("http://sciences.ucf.edu/biology/d4lab/wp-content/uploads/sites/125/2017/01/possum.txt", "\t", escape_double = FALSE, trim_ws = TRUE)

attach(zinc)

The following objects are masked from zinc (pos = 5):

    DIVERSITY, STREAM, ZINC

The following objects are masked from zinc (pos = 9):

    DIVERSITY, STREAM, ZINC

The following objects are masked from zinc (pos = 13):

    DIVERSITY, STREAM, ZINC

The following objects are masked from zinc (pos = 20):

    DIVERSITY, STREAM, ZINC

str(zinc)

Classes ‘tbl_df’, ‘tbl’ and 'data.frame':   34 obs. of  3 variables:
 $ STREAM   : chr  "Arkan" "Arkan" "Arkan" "Arkan" ...
 $ ZINC     : chr  "LOW" "LOW" "LOW" "LOW" ...
 $ DIVERSITY: num  1.4 2.18 1.83 1.88 2.1 2.02 1.94 1.7 2.05 1.98 ...
 - attr(*, "spec")=List of 2
  ..$ cols   :List of 3
  .. ..$ STREAM   : list()
  .. .. ..- attr(*, "class")= chr  "collector_character" "collector"
  .. ..$ ZINC     : list()
  .. .. ..- attr(*, "class")= chr  "collector_character" "collector"
  .. ..$ DIVERSITY: list()
  .. .. ..- attr(*, "class")= chr  "collector_double" "collector"
  ..$ default: list()
  .. ..- attr(*, "class")= chr  "collector_guess" "collector"
  ..- attr(*, "class")= chr "col_spec"

fZinc <- as.factor(ZINC)
leveneTest(DIVERSITY ~ ZINC)

group coerced to factor.

Levene's Test for Homogeneity of Variance (center = median)
      Df F value Pr(>F)
group  3  0.0195 0.9962
      30               

Zin2 <- aov(DIVERSITY ~ ZINC)
Zin3 <- lm(DIVERSITY ~ ZINC)
summary(Zin2)

            Df Sum Sq Mean Sq F value Pr(>F)  
ZINC         3  2.567  0.8555   3.939 0.0176 *
Residuals   30  6.516  0.2172                 
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

summary(Zin3)

Call:
lm(formula = DIVERSITY ~ ZINC)

Residuals:
     Min       1Q   Median       3Q      Max 
-1.03750 -0.22896  0.07986  0.33222  0.79750 

Coefficients:
            Estimate Std. Error t value Pr(>|t|)    
(Intercept)  1.79750    0.16478  10.909 5.81e-12 ***
ZINCHIGH    -0.51972    0.22647  -2.295   0.0289 *  
ZINCLOW      0.23500    0.23303   1.008   0.3213    
ZINCMED     -0.07972    0.22647  -0.352   0.7273    
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 0.4661 on 30 degrees of freedom
Multiple R-squared:  0.2826,    Adjusted R-squared:  0.2108 
F-statistic: 3.939 on 3 and 30 DF,  p-value: 0.01756

interaction.plot(fZinc, trace.factor=ZINC, response=DIVERSITY, 
fun=mean, type="b")

posthoc <- TukeyHSD(x=Zin2, ordered = FALSE, conf.level=0.95)
plot(posthoc)

posthoc

  Tukey multiple comparisons of means
    95% family-wise confidence level

Fit: aov(formula = DIVERSITY ~ ZINC)

$ZINC
                 diff        lwr        upr     p adj
HIGH-BACK -0.51972222 -1.1355064 0.09606192 0.1218677
LOW-BACK   0.23500000 -0.3986367 0.86863665 0.7457444
MED-BACK  -0.07972222 -0.6955064 0.53606192 0.9847376
LOW-HIGH   0.75472222  0.1389381 1.37050636 0.0116543
MED-HIGH   0.44000000 -0.1573984 1.03739837 0.2095597
MED-LOW   -0.31472222 -0.9305064 0.30106192 0.5153456

1. The levene’s test is there and shows that the data set does have a normal distribution (we can’t reject the null hypothesis) so it’s safe to use the dataset without doing transformations to the data.

1. In the one-way ANOVA statistically we see that Zinc has a significant affect on the data set, meaning that within the dataset there is something that is causing the means to be statistically significant. The LM shows us that between all of the populations within the dataset, the mean that is significantly different is only High levels of Zinc. This means that something is occuring in the High groups that’s causing the means to shift in a significant way in that data group. Biologically High levels of zinc is causing a statistically significant decrease in diversity.

attach(barley)
barley$fblock <- factor(barley$Block)
barley$fvariety <- factor(barley$variety)
leveneTest(yield ~ barley$fvariety)

Levene's Test for Homogeneity of Variance (center = median)
      Df F value Pr(>F)
group  3  0.2868  0.834
      12               

leveneTest(yield ~ barley$fblock )

Levene's Test for Homogeneity of Variance (center = median)
      Df F value  Pr(>F)  
group  3  3.3031 0.05756 .
      12                  
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

shapiro.test(yield)

    Shapiro-Wilk normality test

data:  yield
W = 0.84163, p-value = 0.01028

barley$YIELD<- log(barley$yield)
shapiro.test(barley$YIELD)

    Shapiro-Wilk normality test

data:  barley$YIELD
W = 0.90372, p-value = 0.09219

modelA <- aov(barley$YIELD ~ barley$fvariety+ barley$fblock)
modelL <- lm(barley$YIELD ~ barley$fvariety+ barley$fblock)
modelT <- aov(barley$YIELD ~ barley$fblock)
modelR <- aov(barley$YIELD ~ barley$fvariety)
summary(modelA)

                Df Sum Sq Mean Sq F value   Pr(>F)    
barley$fvariety  3 0.0538  0.0179   1.714    0.233    
barley$fblock    3 1.0817  0.3606  34.451 2.91e-05 ***
Residuals        9 0.0942  0.0105                     
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

summary(modelL)

Call:
lm(formula = barley$YIELD ~ barley$fvariety + barley$fblock)

Residuals:
      Min        1Q    Median        3Q       Max 
-0.156441 -0.040672  0.006526  0.030688  0.163595 

Coefficients:
                    Estimate Std. Error t value Pr(>|t|)    
(Intercept)          4.51009    0.06767  66.650 1.95e-13 ***
barley$fvarietyvar2 -0.03993    0.07234  -0.552 0.594418    
barley$fvarietyvar3 -0.14896    0.07234  -2.059 0.069580 .  
barley$fvarietyvar4 -0.01724    0.07234  -0.238 0.816992    
barley$fblock2      -0.38596    0.07234  -5.335 0.000471 ***
barley$fblock3      -0.55392    0.07234  -7.657 3.14e-05 ***
barley$fblock4      -0.69470    0.07234  -9.603 5.01e-06 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 0.1023 on 9 degrees of freedom
Multiple R-squared:  0.9234,    Adjusted R-squared:  0.8723 
F-statistic: 18.08 on 6 and 9 DF,  p-value: 0.0001496

summary(modelT)

              Df Sum Sq Mean Sq F value   Pr(>F)    
barley$fblock  3  1.082  0.3606   29.23 8.45e-06 ***
Residuals     12  0.148  0.0123                     
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

summary(modelR)

                Df Sum Sq Mean Sq F value Pr(>F)
barley$fvariety  3 0.0538 0.01794   0.183  0.906
Residuals       12 1.1759 0.09800               

posthoc1 <- TukeyHSD(x=modelT, ordered = FALSE, conf.level=0.95)
posthoc2 <- TukeyHSD(x=modelR, ordered = FALSE, conf.level=0.95)
posthoc1

  Tukey multiple comparisons of means
    95% family-wise confidence level

Fit: aov(formula = barley$YIELD ~ barley$fblock)

$`barley$fblock`
          diff        lwr         upr     p adj
2-1 -0.3859587 -0.6191245 -0.15279288 0.0017495
3-1 -0.5539242 -0.7870900 -0.32075839 0.0000683
4-1 -0.6946958 -0.9278616 -0.46153002 0.0000069
3-2 -0.1679655 -0.4011313  0.06520030 0.1960145
4-2 -0.3087371 -0.5419029 -0.07557133 0.0093122
4-3 -0.1407716 -0.3739374  0.09239418 0.3231282

posthoc2

  Tukey multiple comparisons of means
    95% family-wise confidence level

Fit: aov(formula = barley$YIELD ~ barley$fvariety)

$`barley$fvariety`
                 diff        lwr       upr     p adj
var2-var1 -0.03992891 -0.6971075 0.6172497 0.9977970
var3-var1 -0.14895801 -0.8061366 0.5082206 0.9053096
var4-var1 -0.01723823 -0.6744168 0.6399403 0.9998205
var3-var2 -0.10902910 -0.7662077 0.5481495 0.9592583
var4-var2  0.02269068 -0.6344879 0.6798692 0.9995915
var4-var3  0.13171979 -0.5254588 0.7888983 0.9316379

2. What we saw here was that there was no difference between the cultivares of barley but rather the significance that we had observed in the ANOVA was due to the block designs that the barley were placed in. We see that because when you do the LM and then when you run a tukey post-hoc we see that there’s no variation between the yields but that there’s significant variation between the block designs. So we can essentially conclude that there were no discernable difference in yield between the different barley cultivars but any significant differencein yield was due to different environmental factors in the different blocks.

attach(possum)
possum <-na.omit(possum)
leveneTest(totalL ~ sex)

Levene's Test for Homogeneity of Variance (center = median)
       Df F value Pr(>F)
group   1  0.5079 0.4777
      100               

shapiro.test(totalL)

    Shapiro-Wilk normality test

data:  totalL
W = 0.98521, p-value = 0.3154

shapiro.test(skullW)

    Shapiro-Wilk normality test

data:  skullW
W = 0.93754, p-value = 0.0001169

PopulationWhat <- aov(totalL ~ pop)
modelA <- aov(totalL ~ sex*age)
totallength <- aov(totalL ~ pop + Error(sex+age))
summary(PopulationWhat)

             Df Sum Sq Mean Sq F value Pr(>F)
pop           1     26   25.97   1.481  0.226
Residuals   100   1753   17.53               

summary(modelA)

            Df Sum Sq Mean Sq F value Pr(>F)  
sex          1   34.3   34.33   2.090 0.1514  
age          1  112.9  112.87   6.872 0.0101 *
sex:age      1   22.7   22.69   1.382 0.2427  
Residuals   98 1609.5   16.42                 
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

summary(totallength)

Error: sex
    Df Sum Sq Mean Sq
pop  1  34.33   34.33

Error: age
    Df Sum Sq Mean Sq
pop  1  112.9   112.9

Error: Within
          Df Sum Sq Mean Sq F value Pr(>F)
pop        1    9.8   9.771    0.59  0.444
Residuals 98 1622.4  16.555               

PopulationWha <- aov(skullW ~ pop)
model <- aov(skullW ~ sex*age)
skulllength <- aov(skullW ~ pop + Error(sex+age))
summary(PopulationWha)

             Df Sum Sq Mean Sq F value Pr(>F)
pop           1    1.5   1.475   0.153  0.696
Residuals   100  961.2   9.612               

summary(model)

            Df Sum Sq Mean Sq F value  Pr(>F)   
sex          1   10.3   10.31   1.222 0.27174   
age          1   82.3   82.30   9.752 0.00236 **
sex:age      1   43.0   43.01   5.096 0.02620 * 
Residuals   98  827.0    8.44                   
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

summary(skulllength)

Error: sex
    Df Sum Sq Mean Sq
pop  1  10.31   10.31

Error: age
    Df Sum Sq Mean Sq
pop  1   82.3    82.3

Error: Within
          Df Sum Sq Mean Sq F value Pr(>F)
pop        1    1.5   1.456   0.164  0.686
Residuals 98  868.6   8.863               

3. So with the total length what we see is that between the popoulations there’s no significant difference between the two populations when you’re controlling for sex and age, but there is a significant difference among ages, with possums getting longer as they get older independent of their sex or the population. Regarding the skull width, I ran the tests on the popoulations but it doesn’t particularly matter because no amount of transforming the data could make the data pass the shapiro or levenes’ tests. So the tests are unreliable and since we don’t know any non parametric tests, I’m going to forgoe attempting to infer the data and assume that it’s not dependable. (Should I attempt to discern the data, there’s an interaction between sex and age on skull width).