Load packages

library(ggplot2)
library(dplyr)
library(tidyr)
library(readxl)

Question 1: Squirrels

Load and review data

squirrels <- read_excel("C:/Users/charl/OneDrive - Nanyang Technological University/ES Mods/ES3307 Experimental Design/Assignment 1/Assignment 1 data.xls")

# Check data structure and summary
str(squirrels)
## tibble [50 × 2] (S3: tbl_df/tbl/data.frame)
##  $ MALE  : num [1:50] 0.41 0.38 1.13 0.61 0.73 0.56 0.62 0.75 0.5 0.75 ...
##  $ FEMALE: num [1:50] 0.4 0.39 0.59 0.62 0.53 0.34 0.41 0.59 0.46 0.63 ...
summary(squirrels)
##       MALE            FEMALE     
##  Min.   :0.2900   Min.   :0.320  
##  1st Qu.:0.4350   1st Qu.:0.410  
##  Median :0.5600   Median :0.495  
##  Mean   :0.5908   Mean   :0.518  
##  3rd Qu.:0.7025   3rd Qu.:0.605  
##  Max.   :1.1300   Max.   :0.850

Q1A

Plot the data using histogram(s) and boxplot(s).

Plot histogram

# Convert data from wide to long format
squirrels_long <- squirrels %>% 
  pivot_longer(
    c(MALE, FEMALE),   # Columns to convert
    names_to = "Sex",      # New column name for MALE/FEMALE
    values_to = "Weight"      # New column name for values
  )

# Plot histogram of squirrel weight
squirrels_long %>%
  ggplot(aes(x = Weight, fill = Sex)) +
  geom_histogram(binwidth = 0.05) +  # Set binwidth to see distri
  labs(x = "Weight", y = "Count") +
  scale_fill_manual(values = c("lightpink","lightblue"))

Fig 1. Histogram of female and male squirrel weights

Plot Box plot

# Plot box plot of squirrel weights
squirrels_long %>% 
  ggplot(aes(x = Sex, y = Weight)) +
  geom_boxplot(fill = c("lightpink","lightblue"))

Fig 2. Box plot of female and male squirrel weights

Q1B

Perform a two-sample t-test and a Mann-Whitney test.

Check if sample variance is equal

var.test(Weight ~ Sex, squirrels_long)
## 
##  F test to compare two variances
## 
## data:  Weight by Sex
## F = 0.44156, num df = 49, denom df = 49, p-value = 0.004977
## alternative hypothesis: true ratio of variances is not equal to 1
## 95 percent confidence interval:
##  0.2505735 0.7781079
## sample estimates:
## ratio of variances 
##          0.4415577

Since p-value = 0.00498 is less than 0.05, the variance between male and female squirrel weights is statistically significant hence are unequal. Thus, it must be accounted for when performing the t-test.

Two-sample t-test

t.test(Weight ~ Sex, squirrels_long, var.equal = FALSE)
## 
##  Welch Two Sample t-test
## 
## data:  Weight by Sex
## t = -2.1662, df = 85.212, p-value = 0.03309
## alternative hypothesis: true difference in means between group FEMALE and group MALE is not equal to 0
## 95 percent confidence interval:
##  -0.139618083 -0.005981917
## sample estimates:
## mean in group FEMALE   mean in group MALE 
##               0.5180               0.5908

Mann-Whitney test

wilcox.test(Weight ~ Sex, squirrels_long)
## 
##  Wilcoxon rank sum test with continuity correction
## 
## data:  Weight by Sex
## W = 1012, p-value = 0.1014
## alternative hypothesis: true location shift is not equal to 0

Q1C

What features of the plots suggest that the use of a parametric test may be unwise?

From the histogram, both female and male squirrel weights look right-skewed (Fig. 1), which do not meet the assumption of a parametric test that the data is normally distributed. Thus, a non-parametric test that does not have this assumption would be more suitable to compare between both groups. Additionally, the box plot shows the male squirrel weights have a larger variance than the female weights as seen from the larger interquartile range and longer whiskers (Fig. 2). This suggests unequal variances which do not fulfill the assumptions of a parametric test.

Q1D

From the results of the (parametric) t-test, report the t, dfs, and p-values, and interpret them.

t = -2.1662, df = 85.212, p-value = 0.03309

The difference in weight between female and male squirrels are statistically significant from the t-test with p-value = 0.0331 which is less than 0.05 and t1,85.21 = -2.17. However, as the weights of female and male are not normally distributed, they do not follow the assumptions of a t-test and further interpretation should not be conducted.

Q1E

From the results of the (non-parametric) Mann-Whitney test, report the W, and p-values, and interpret them.

W = 1012, p-value = 0.1014

From the Mann-Whitney test, the p-value = 0.101 is greater than 0.05 and W = 1012, which means the difference in weights between female and male squirrels are not statistically significant.

Q1F

What type of error do you get if you use a t-test?

Type 1 error. From the Mann-Whitney test, the difference in weights between both sexes are not statistically significant but the t-test shows that the difference is statistically significant, even though it was run under false assumptions that the weights are normally distributed. Thus, the t-test shows a statistically significant result even though it is not, meaning a type 1 error.

Question 2: Melons

Load and review data

melons <- read_excel("C:/Users/charl/OneDrive - Nanyang Technological University/ES Mods/ES3307 Experimental Design/Assignment 1/Assignment 1 data.xls", sheet = "Melons")

# Check data structure and summary
str(melons)
## tibble [22 × 2] (S3: tbl_df/tbl/data.frame)
##  $ YIELDM : num [1:22] 25.1 17.2 26.4 16.1 22.1 ...
##  $ VARIETY: num [1:22] 1 1 1 1 1 1 2 2 2 2 ...
summary(melons)
##      YIELDM         VARIETY     
##  Min.   :15.05   Min.   :1.000  
##  1st Qu.:22.26   1st Qu.:1.250  
##  Median :27.82   Median :2.000  
##  Mean   :27.66   Mean   :2.455  
##  3rd Qu.:32.90   3rd Qu.:3.750  
##  Max.   :43.32   Max.   :4.000

Q2A

What is/are the hypothesis/es?

The hypothesis is there is a statistically significant difference between the 4 different melon variety’s yield because of the variety’s genetics.

Q2B

Plot the data so that you can see the difference in yield among the varieties.

Plot a box plot of the yield by melon variety

# Convert VARIETY to factor variable
melons$VARIETY <- as.factor(melons$VARIETY)
str(melons)
## tibble [22 × 2] (S3: tbl_df/tbl/data.frame)
##  $ YIELDM : num [1:22] 25.1 17.2 26.4 16.1 22.1 ...
##  $ VARIETY: Factor w/ 4 levels "1","2","3","4": 1 1 1 1 1 1 2 2 2 2 ...
# Box plot of each variety's yield
melons %>% 
  ggplot(aes(x = VARIETY, y = YIELDM)) +
  labs(x = "Melon variety", y = "Yield") +
  geom_boxplot(fill = c("pink","lightblue","lightgreen","purple"))

Fig 3. Box plot of melon yield by variety

Q2C

Produce descriptive statistics that show you the means of the groups, and the confidence intervals.

Summarise by melon variety

melons_summary <- melons %>%
  group_by(VARIETY) %>%
  summarise(
    count = n(),
    mean_yield = mean(YIELDM, na.rm = TRUE),
    sd_yield = sd(YIELDM, na.rm = TRUE),
    
    # Calculate the 95% CI bounds
    se_yield   = sd_yield / sqrt(count), # standard error
    # Use t-distri to calc bounds
    ci_lower   = mean_yield - qt(1 - (0.05/2), count - 1) * se_yield,
    ci_upper   = mean_yield + qt(1 - (0.05/2), count - 1) * se_yield
  )

print(melons_summary)
## # A tibble: 4 × 7
##   VARIETY count mean_yield sd_yield se_yield ci_lower ci_upper
##   <fct>   <int>      <dbl>    <dbl>    <dbl>    <dbl>    <dbl>
## 1 1           6       20.5     4.69    1.92      15.6     25.4
## 2 2           6       37.4     3.95    1.61      33.3     41.5
## 3 3           4       20.5     4.76    2.38      12.9     28.0
## 4 4           6       29.9     2.23    0.910     27.6     32.2

Q2D

Produce an analysis of variance of YIELDM with respect to VARIETY.

melons_lm <- lm(YIELDM ~ VARIETY, melons)
anova(melons_lm)
## Analysis of Variance Table
## 
## Response: YIELDM
##           Df  Sum Sq Mean Sq F value    Pr(>F)    
## VARIETY    3 1115.28  371.76  23.798 1.735e-06 ***
## Residuals 18  281.19   15.62                      
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Q2E

Show and describe the model’s resulting diagnostic plots.

par(mfrow = c(2, 2))
plot(melons_lm)

Fig 4. Diagnostic plots of the model

The Residuals vs Fitted graph shows a straight line with points randomly arranged around 0. This suggests a linear relationship.

The Q-Q residuals graph shows the points generally follow the dotted line, suggesting a normal distribution of residuals.

The Scale-Location plot shows the data randomly scatter around the red line, suggesting homoscedasticity.

The Residuals vs Leverage graph shows there are no points lying outside of Cook’s distance (not drawn as out of graph bounds), which means there are no outliers in the data set.

Q2F

Report the F, dfs, and p-values, and interpret them.

F = 23.798, df (VARIETY) = 3, df (Residuals) = 18, p-value = 1.735e-06

Results of ANOVA shows that F3,18 = 23.8 and p-value = 1.74e-06, which p-value is lower than 0.05, suggesting there is a statistically significant difference between the yields of the 4 melon varieties. From Fig 3, the yield of variety 2 and 4 are higher than that of variety 1 and 3. However, as 2 plots of variety 3 were destroyed, the number of replicates across varieties were inconsistent and ranged from 4-6. Future studies can be conducted with the same number of replicates to allow for more accurate statistical results.

Question 3: Dioecious trees

Load and review data

trees <- read_excel("C:/Users/charl/OneDrive - Nanyang Technological University/ES Mods/ES3307 Experimental Design/Assignment 1/Assignment 1 data.xls", sheet = "Dioecious trees")

# Check data structure and summary
str(trees)
## tibble [50 × 3] (S3: tbl_df/tbl/data.frame)
##  $ SEX    : num [1:50] 2 1 2 1 2 1 2 2 1 1 ...
##  $ DBH    : num [1:50] 174 171 212 187 227 274 86 161 285 182 ...
##  $ FLOWERS: num [1:50] 234 110 512 133 1107 ...
summary(trees)
##       SEX           DBH           FLOWERS      
##  Min.   :1.0   Min.   : 68.0   Min.   :   1.0  
##  1st Qu.:1.0   1st Qu.:140.5   1st Qu.: 115.8  
##  Median :2.0   Median :188.0   Median : 253.0  
##  Mean   :1.6   Mean   :191.6   Mean   : 383.5  
##  3rd Qu.:2.0   3rd Qu.:242.0   3rd Qu.: 533.0  
##  Max.   :2.0   Max.   :316.0   Max.   :1482.0
# Convert SEX to factor variable
trees$SEX <- as.factor(trees$SEX)

str(trees)
## tibble [50 × 3] (S3: tbl_df/tbl/data.frame)
##  $ SEX    : Factor w/ 2 levels "1","2": 2 1 2 1 2 1 2 2 1 1 ...
##  $ DBH    : num [1:50] 174 171 212 187 227 274 86 161 285 182 ...
##  $ FLOWERS: num [1:50] 234 110 512 133 1107 ...

Q3A

Show graphically how the number of flowers differs between the sexes.

Plot histogram to see distribution between the sexes

trees %>% 
  ggplot(aes(x = FLOWERS, fill = SEX)) +
  geom_histogram(binwidth = 40) +
  labs(x = "Number of Flowers", y = "Count") +
  scale_fill_manual(
    name = "Sex", 
    values = c("lightblue", "lightpink"),
    labels = c("1" = "Male", "2" = "Female")
  )

Fig 5. Histogram of the number of flowers between the sexes

Q3B

Test the hypothesis that male and female trees produce different number of flowers.

From Fig 5, the number of flowers produced by female trees looks right skewed, while that of male trees looks bimodal. Both distributions do not look normal hence the normal distribution assumption to conduct a t-test does not hold and a Mann-Whitney test is more appropriate the compare whether the number of flowers produced between male and female trees are statistically different.

wilcox.test(FLOWERS ~ SEX, trees)
## Warning in wilcox.test.default(x = DATA[[1L]], y = DATA[[2L]], ...): cannot
## compute exact p-value with ties
## 
##  Wilcoxon rank sum test with continuity correction
## 
## data:  FLOWERS by SEX
## W = 298, p-value = 0.9763
## alternative hypothesis: true location shift is not equal to 0

Q3C

Report the appropriate test statistic, dfs, and p-values, and interpret them.

W = 298, p-value = 0.9763

From the Mann-Whitney test, W = 298 and p-value = 0.976. The difference between the number of flowers produced by male and female trees are not statistically different as p-value is greater than 0.05, hence the number of flowers produced by both sexes are around the same.