library(ggplot2)
library(dplyr)
library(tidyr)
library(readxl)
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
Plot the data using histogram(s) and boxplot(s).
# 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 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
Perform a two-sample t-test and a Mann-Whitney test.
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.
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
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
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.
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.
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.
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.
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
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.
Plot the data so that you can see the difference in yield among the varieties.
# 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
Produce descriptive statistics that show you the means of the groups, and the confidence intervals.
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
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
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.
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.
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 ...
Show graphically how the number of flowers differs 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
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
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.