Question 1: Describe and summarise your assigned dataset.

# Loading relevant packages
library(tidyverse) # iris data set is found here
library (ggplot2)

#Load dataset
{data(iris)}

#Summarise dataset
head (iris) # view first 6 rows to get a brief understanding of the data set
##   Sepal.Length Sepal.Width Petal.Length Petal.Width Species
## 1          5.1         3.5          1.4         0.2  setosa
## 2          4.9         3.0          1.4         0.2  setosa
## 3          4.7         3.2          1.3         0.2  setosa
## 4          4.6         3.1          1.5         0.2  setosa
## 5          5.0         3.6          1.4         0.2  setosa
## 6          5.4         3.9          1.7         0.4  setosa
str (iris) # shows the internal strcture of the different objects
## 'data.frame':    150 obs. of  5 variables:
##  $ Sepal.Length: num  5.1 4.9 4.7 4.6 5 5.4 4.6 5 4.4 4.9 ...
##  $ Sepal.Width : num  3.5 3 3.2 3.1 3.6 3.9 3.4 3.4 2.9 3.1 ...
##  $ Petal.Length: num  1.4 1.4 1.3 1.5 1.4 1.7 1.4 1.5 1.4 1.5 ...
##  $ Petal.Width : num  0.2 0.2 0.2 0.2 0.2 0.4 0.3 0.2 0.2 0.1 ...
##  $ Species     : Factor w/ 3 levels "setosa","versicolor",..: 1 1 1 1 1 1 1 1 1 1 ...

From the basic summary table, we can see that the iris dataset consists of four physical characteristics of flowers, namely the sepal length, sepal width, petal length and petal width. They should all be measured and recorded in centimeters and are thus saved as numerical data types. In addition, this dataset has a species coloumn which is stored as a factor with 3 levels, where the flowers belong to 3 different species: setosa, versicolor and virginica. It is stored this way as species is a type of categorical variable. There are 50 flowers from each species, amounting to 150 records.

# Summary statistics table of all 4 variables by Species
iris_summary <- iris %>%
  #Convert the dataset from wide format to long format
  pivot_longer(
    cols = c(Sepal.Length, Sepal.Width, Petal.Length, Petal.Width),
    names_to = "Variable",     #name of measurement
    values_to = "Measurement"  #actual measurement value
  ) %>%
  
  #Group observations by Iris species and measurement
  group_by(Species, Variable) %>%
  
  #Calculate the mean and SD for each measurement within each species 
  summarise(
    Mean = mean(Measurement, na.rm = TRUE),
    SD = sd(Measurement, na.rm = TRUE)
  )

# Display the summary stats in a neat table
knitr::kable(iris_summary, 
             digits = 3,   #round to 3dp
             caption = "Summary Statistics of Petal and Sepal Dimensions by Iris Species")
Summary Statistics of Petal and Sepal Dimensions by Iris Species
Species Variable Mean SD
setosa Petal.Length 1.462 0.174
setosa Petal.Width 0.246 0.105
setosa Sepal.Length 5.006 0.352
setosa Sepal.Width 3.428 0.379
versicolor Petal.Length 4.260 0.470
versicolor Petal.Width 1.326 0.198
versicolor Sepal.Length 5.936 0.516
versicolor Sepal.Width 2.770 0.314
virginica Petal.Length 5.552 0.552
virginica Petal.Width 2.026 0.275
virginica Sepal.Length 6.588 0.636
virginica Sepal.Width 2.974 0.322

From the summary statistics, we can see that for each species the values for the four physical characteristics varies. The species, virginica, has the longest mean petal length, 5.55 cm (2 d.p.), largest mean petal width, 2.03 cm (2 d.p.), as well as the longest mean sepal length, 6.59 cm (2 d.p.). However, setosa, has the largest mean sepal width, 3.43 cm (2 d.p.).

# Histogram of the 4 variables 
library(patchwork) #used to combine multiple plots in one figure

#Plot 1: Sepal Length
p1 <- iris %>% 
  ggplot() +
  
  #Create histogram of sepal length and colour by species
  geom_histogram(
    aes(x = Sepal.Length, fill = Species)
  ) +
  
  #Separate the histogram into one row for each species
  facet_grid(rows = vars(Species)) +
  
  #Add axis labels and title 
  labs(
    x = "Sepal Length (cm)",
    y = "Count",
    title = "Sepal Length"
  ) +
  
  #Assign colours for the 3 species 
  scale_fill_manual(
    values = c("darkorange", "purple", "cyan4")
  ) +
  
  #Set graph theme and remove legend
  theme_minimal() +
  theme(legend.position = "none")

#Plot 2: Sepal Width
p2 <- iris %>% 
  ggplot() +
  
  #Create histogram for sepal width and colour by species
  geom_histogram(
    aes(x = Sepal.Width, fill = Species)
  ) +
  
  #Separate the histogram into one row for each species
  facet_grid(rows = vars(Species)) +
  
  #Add axis labels and title 
  labs(
    x = "Sepal Width (cm)",
    y = "Count",
    title = "Sepal Width"
  ) +
  
   #Assign colours for the 3 species 
  scale_fill_manual(
    values = c("darkorange", "purple", "cyan4")
  ) +
  
  #Set graph theme and remove legend 
  theme_minimal() +
  theme(legend.position = "none")

#Plot 3: Petal Length
p3 <- iris %>% 
  ggplot() +
  #Create histogram of petal length and colour by species
  geom_histogram(
    aes(x = Petal.Length, fill = Species)
  ) +

  #Separate the histogram into one row for each species
  facet_grid(rows = vars(Species)) +
  
  #Add axis labels and title 
  labs(
    x = "Petal Length (cm)",
    y = "Count",
    title = "Petal Length"
  ) +
  
   #Assign colours for the 3 species 
  scale_fill_manual(
    values = c("darkorange", "purple", "cyan4")
  ) +
  
  #Set graph theme and remove legend
  theme_minimal() +
  theme(legend.position = "none")

#Petal Width
p4 <- iris %>% 
  ggplot() +
  #Create histogram for petal width and colour by species
  geom_histogram(
    aes(x = Petal.Width, fill = Species)
  ) +
  
  #Separate the histogram into one row for each species
  facet_grid(rows = vars(Species)) +
  
  #Add axis and title 
  labs(
    x = "Petal Width (cm)",
    y = "Count",
    title = "Petal Width"
  ) +
  scale_fill_manual(
    values = c("darkorange", "purple", "cyan4")
  ) +
  
  #Set graph theme and remove legend
  theme_minimal() +
  theme(legend.position = "none")


#Combine all 4 plots
(p1 | p2) /
(p3 | p4)

The histograms show clear differences in flower dimensions among the three Iris species. Sepal measurements overlap considerably between species, while petal measurements show clearer separation. Setosa has the smallest petals, versicolor has intermediate-sized petals, and virginica has the largest petals.

Question 2: How do petal length and petal width vary across species?

#Side by side Box plot of petal length and petal width

# Convert petal measurements to long format
petal_long <- iris %>%
  pivot_longer(
    cols = c(Petal.Length, Petal.Width),
    names_to = "Petal_dimension",
    values_to = "Measurement"
  )

# Boxplots
ggplot(petal_long,
       aes(x = Species, y = Measurement, fill = Species)) +
  geom_boxplot() +
  facet_wrap(~Petal_dimension, ncol = 2) +
  labs(
    x = "Species",
    y = "Measurement (cm)",
    title = "Petal Dimensions Across Iris Species"
  ) +
  scale_fill_manual(
    values = c("darkorange", "purple", "cyan4")
  ) +
  theme_classic() +
  theme(legend.position = "none")

For petal length, the measurements increases sharply from setosa to virginica. Setosa has the lowest variation in petal length measurements, with 3 outliers, whereas versicolor shows slightly greater variability with 1 outlier. As for virginica, the petal length shows the greatest variability and no outliers. Overall, virginica has the longest petal length and setosa has the shortest.

Similarly,petal width increases from setosa to virginica. The spread in petal width for all 3 species are relatively smaller than that of petal length, and only setosa has 2 outliers. Overall, virginica has the widest petals, while setosa has the narrowest.

# Summary statistics of petal width and petal length by species
petal_summary <- iris %>%
  #group data by species
  group_by(Species) %>%
  #calculate mean and sd
  summarise(
    `Mean Petal Length` = mean(Petal.Length),
    `SD Petal Length` = sd(Petal.Length),
    `Mean Petal Width` = mean(Petal.Width),
    `SD Petal Width` = sd(Petal.Width)
  )

#Table
knitr::kable(
  petal_summary,
  digits = 3,  #round to 3dp
  caption = "Summary Statistics of Petal Dimensions by Iris Species"
)
Summary Statistics of Petal Dimensions by Iris Species
Species Mean Petal Length SD Petal Length Mean Petal Width SD Petal Width
setosa 1.462 0.174 0.246 0.105
versicolor 4.260 0.470 1.326 0.198
virginica 5.552 0.552 2.026 0.275

The mean values from the summary statistics complement our findings from the boxplot graphs.

Question 3: Graph the petal dimension across Iris species. Are sepal dimensions correlated with petal dimensions?

library (patchwork)
#Calculate overall correlation coefficients
r1 <- cor(iris$Sepal.Length, iris$Petal.Length)
r2 <- cor(iris$Sepal.Width, iris$Petal.Width)
r3 <- cor(iris$Sepal.Length, iris$Petal.Width)
r4 <- cor(iris$Sepal.Width, iris$Petal.Length)

#Assign species colours
species_colours <- c(
  "setosa" = "darkorange",
  "versicolor" = "purple",
  "virginica" = "cyan4"
)

#Plot 1: Sepal Length vs Petal Length
p1 <- ggplot(
  iris,
  aes(
    x = Sepal.Length, #x-axis
    y = Petal.Length, #y-axis
    colour = Species #colour points according to species
  )
) +
  #Add data points
  geom_point(size = 2) +
  
  #Add linear regression line 
  geom_smooth(
    aes(group = 1),   #treat all species as one group for the regression line
    method = "lm",    #fit linear regression model
    colour = "black"
  ) +
  
  #Add Person correlation coefficient (r) to the graph
  annotate(
    "text",
    x = Inf, y = -Inf,  #position text at bottom right of graph
    label = paste("r =", round(r1, 3)), #display r rounded to 3dp
    hjust = 1.2, vjust = -1,#horizontal and vertical position of text
    size = 3
  ) +
  
  #Set colours for the 3 species 
  scale_colour_manual(values = species_colours) +
  
  #Axis labels, graph title and legend title
  labs(
    x = "Sepal Length (cm)",  
    y = "Petal Length (cm)",  
    title = "Sepal Length vs Petal Length",
    colour = "Species"        
  ) +
  
  #Set theme 
  theme_classic()

#Plot 2: Sepal Width vs Petal Width
p2 <- ggplot(
  iris,
  aes(
    x = Sepal.Width,  #x-axis
    y = Petal.Width,  #y-axis
    colour = Species  #colour points according to species
  )
) +
  #Add data points
  geom_point(size = 2) +
  
  #Add linear regression line 
  geom_smooth(
    aes(group = 1), #treat all species as one group for the regression line
    method = "lm",  #fit linear regression model 
    colour = "black"
  ) +
  
  #Add r values on graph 
  annotate(
    "text",
    x = 4.3, y = 0.8, #position of text
    label = paste("r =", round(r2, 3)),  #r round to 3dp
    size = 3
  ) +
  
  #Set colours for the 3 species 
  scale_colour_manual(values = species_colours) +
  
  #Axis labels, graph title and legend title
  labs(
    x = "Sepal Width (cm)",
    y = "Petal Width (cm)",
    title = "Sepal Width vs Petal Width",
    colour = "Species"
  ) +
  theme_classic()


#Plot 3: Sepal Length vs Petal Width
p3 <- ggplot(
  iris,
  aes(
    x = Sepal.Length, #x-axis
    y = Petal.Width,  #y-axis
    colour = Species #colour according to species
  )
) +
  #Add data points
  geom_point(size = 2) +

  #Add linear regression line 
  geom_smooth(
    aes(group = 1), #treat all species as one group for regression line
    method = "lm",  #fit linear regression model 
    colour = "black"
  ) +
  
  #Add r value on graph 
  annotate(
    "text",
    x = Inf, y = -Inf,  #position at bottom right of graph
    label = paste("r =", round(r3, 3)), #r round to 3dp
    hjust = 1.2, vjust = -1, #horizontal and vertical position
    size = 3
  ) +
  
  #Set colours for the 3 species 
  scale_colour_manual(values = species_colours) +
  
  #Axis labels, graph title and legend title
  labs(
    x = "Sepal Length (cm)",
    y = "Petal Width (cm)",
    title = "Sepal Length vs Petal Width",
    colour = "Species"
  ) +
  theme_classic()


#Plot 4: Sepal Width vs Petal Length
p4 <- ggplot(
  iris,
  aes(
    x = Sepal.Width,  #xaxis
    y = Petal.Length, #yaxis
    colour = Species  #colour according to species
  )
) +
  #Add data points
  geom_point(size = 2) +
  
  #Add linear regression line
  geom_smooth(
    aes(group = 1),  #Treat all species as one group for regression line
    method = "lm",   #fit linear regression model 
    colour = "black"
  ) +
  
  #Add r value on graph 
  annotate(
    "text",
    x = 4.3, y = 3,  #position of text
    label = paste("r =", round(r4, 3)), #r round to 3 dp
    size = 3
  ) +
  
  #Set colours for the 3 species 
  scale_colour_manual(values = species_colours) +
  
  #Axis labels, graph title and legend title
  labs(
    x = "Sepal Width (cm)",
    y = "Petal Length (cm)",
    title = "Sepal Width vs Petal Length",
    colour = "Species"
  ) +
  theme_classic()

#Combine all 4 plots
(p1 | p2) /
(p3 | p4) +
  plot_layout(guides = "collect") & #Combine species legend into one single legend
  theme(legend.position = "bottom") #Position legend at bototm of figure

  1. Sepal length and petal length shows a strong positive correlation, with a r value of 0.87 (2 d.p.). This r value is the Pearson’s correlation coefficient, which is used to analyse the strength and direction of the linear relationship of 2 continuous variables. If the r value is closer to -1 (for negative linear relationships) or 1 (for positive linear relationships), the correlation is stronger.

  2. Sepal width and petal width shows a weak to moderate negative correlation, of an r value of -0.36 (2 d.p.).

  3. Sepal length and petal width shows a strong positive correlation, with a r value of 0.82 (2 d.p.).

  4. Sepal width and petal length show a moderate negative correlation, with a r value of -0.42 (2 d.p).

Overall, sepal length shows a strong positive correlation with both petal length and petal width, while sepal width shows weak-to-moderate negative correlations with petal dimensions. The clustering of observations by species suggests that these overall correlations may also be influenced by morphological differences among the three Iris species.

## Box plot for all 4 dimensions 
# Convert all flower dimensions to long format
iris_long <- iris %>%
  pivot_longer(
    cols = c(Sepal.Length, Sepal.Width,
             Petal.Length, Petal.Width),
    names_to = "Flower_dimension",
    values_to = "Measurement"
  )

# Boxplots for all 4 dimensions
ggplot(iris_long,
       aes(x = Species, y = Measurement, fill = Species)) +
  geom_boxplot() +
  facet_wrap(~Flower_dimension, ncol = 2) +
  labs(
    x = "Species",
    y = "Measurement (cm)",
    title = "Flower Dimensions Across Iris Species"
  ) +
  scale_fill_manual(
    values = c("darkorange", "purple", "cyan4")
  ) +
  theme_classic() +
  theme(legend.position = "none")

## Statistic Tests Do the petal and sepal dimension differ between species?

X variable –> Species (Single categorical variable) Y variable –> Petal width, Petal length, sepal width, sepal length (multiple continuous groups)

One-Way ANOVA

# ANOVA for Petal Length by Species
model_petal_len <- lm(Petal.Length ~ Species, data = iris) 
anova(model_petal_len)
## Analysis of Variance Table
## 
## Response: Petal.Length
##            Df Sum Sq Mean Sq F value    Pr(>F)    
## Species     2 437.10 218.551  1180.2 < 2.2e-16 ***
## Residuals 147  27.22   0.185                      
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
# ANOVA for Petal Width by Species
model_petal_wid <- lm(Petal.Width ~ Species, data = iris)
anova(model_petal_wid)
## Analysis of Variance Table
## 
## Response: Petal.Width
##            Df Sum Sq Mean Sq F value    Pr(>F)    
## Species     2 80.413  40.207  960.01 < 2.2e-16 ***
## Residuals 147  6.157   0.042                      
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
# ANOVA for Sepal Length by Species
model_sepal_len <- lm(Sepal.Length ~ Species, data = iris)
anova(model_sepal_len)
## Analysis of Variance Table
## 
## Response: Sepal.Length
##            Df Sum Sq Mean Sq F value    Pr(>F)    
## Species     2 63.212  31.606  119.26 < 2.2e-16 ***
## Residuals 147 38.956   0.265                      
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
# ANOVA for Sepal Width by Species
model_sepal_wid <- lm(Sepal.Width ~ Species, data = iris)
anova(model_sepal_wid)
## Analysis of Variance Table
## 
## Response: Sepal.Width
##            Df Sum Sq Mean Sq F value    Pr(>F)    
## Species     2 11.345  5.6725   49.16 < 2.2e-16 ***
## Residuals 147 16.962  0.1154                      
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Based on the ANOVA results, petal and sepal dimensions for all 4 groups were statistically significantly different between the 3 species (p <0.05).

However, based on the boxplots plotted earlier, it seems that the the 3 species do not have equal variances. Hence, we conduct Kruskal–Wallis test as a non-parametric approach to assess whether the flower dimensions differed among the three species.

Kruskal Wallis Test

# Kruskal Wallis Test 

# Petal Length by Species
kruskal.test(Petal.Length ~ Species, data = iris)
## 
##  Kruskal-Wallis rank sum test
## 
## data:  Petal.Length by Species
## Kruskal-Wallis chi-squared = 130.41, df = 2, p-value < 2.2e-16
# Petal Width by Species
kruskal.test(Petal.Width ~ Species, data = iris)
## 
##  Kruskal-Wallis rank sum test
## 
## data:  Petal.Width by Species
## Kruskal-Wallis chi-squared = 131.19, df = 2, p-value < 2.2e-16
# Sepal Length by Species
kruskal.test(Sepal.Length ~ Species, data = iris)
## 
##  Kruskal-Wallis rank sum test
## 
## data:  Sepal.Length by Species
## Kruskal-Wallis chi-squared = 96.937, df = 2, p-value < 2.2e-16
# Sepal Width by Species
kruskal.test(Sepal.Width ~ Species, data = iris)
## 
##  Kruskal-Wallis rank sum test
## 
## data:  Sepal.Width by Species
## Kruskal-Wallis chi-squared = 63.571, df = 2, p-value = 1.569e-14

The petal and sepal dimensions across the three species were statistically significantly different (p < 0.05). Hence, we conducted a post-hoc Tukey HSD test to determine which pairs of species differed significantly from each other for each floral dimension.

# Tukey HSD test

# Petal Length
tukey_petal_len <- TukeyHSD(aov(Petal.Length ~ Species, data = iris))
print(tukey_petal_len)
##   Tukey multiple comparisons of means
##     95% family-wise confidence level
## 
## Fit: aov(formula = Petal.Length ~ Species, data = iris)
## 
## $Species
##                       diff     lwr     upr p adj
## versicolor-setosa    2.798 2.59422 3.00178     0
## virginica-setosa     4.090 3.88622 4.29378     0
## virginica-versicolor 1.292 1.08822 1.49578     0
# Petal Width
tukey_petal_wid <- TukeyHSD(aov(Petal.Width ~ Species, data = iris))
print(tukey_petal_wid)
##   Tukey multiple comparisons of means
##     95% family-wise confidence level
## 
## Fit: aov(formula = Petal.Width ~ Species, data = iris)
## 
## $Species
##                      diff       lwr       upr p adj
## versicolor-setosa    1.08 0.9830903 1.1769097     0
## virginica-setosa     1.78 1.6830903 1.8769097     0
## virginica-versicolor 0.70 0.6030903 0.7969097     0
# Sepal Length
tukey_sepal_len <- TukeyHSD(aov(Sepal.Length ~ Species, data = iris))
print(tukey_sepal_len)
##   Tukey multiple comparisons of means
##     95% family-wise confidence level
## 
## Fit: aov(formula = Sepal.Length ~ Species, data = iris)
## 
## $Species
##                       diff       lwr       upr p adj
## versicolor-setosa    0.930 0.6862273 1.1737727     0
## virginica-setosa     1.582 1.3382273 1.8257727     0
## virginica-versicolor 0.652 0.4082273 0.8957727     0
# Sepal Width
tukey_sepal_wid <- TukeyHSD(aov(Sepal.Width ~ Species, data = iris))
print(tukey_sepal_wid)
##   Tukey multiple comparisons of means
##     95% family-wise confidence level
## 
## Fit: aov(formula = Sepal.Width ~ Species, data = iris)
## 
## $Species
##                        diff         lwr        upr     p adj
## versicolor-setosa    -0.658 -0.81885528 -0.4971447 0.0000000
## virginica-setosa     -0.454 -0.61485528 -0.2931447 0.0000000
## virginica-versicolor  0.204  0.04314472  0.3648553 0.0087802

The Tukey HSD test showed significant differences between all pairs of species for all four flower dimensions (p < 0.05).