install.packages("palmerpenguins")
library(palmerpenguins)
## Warning: package 'palmerpenguins' was built under R version 4.3.3
library(tidyverse)
## Warning: package 'tidyverse' was built under R version 4.3.3
## Warning: package 'ggplot2' was built under R version 4.3.3
## Warning: package 'tidyr' was built under R version 4.3.3
## Warning: package 'readr' was built under R version 4.3.2
## Warning: package 'purrr' was built under R version 4.3.2
## Warning: package 'dplyr' was built under R version 4.3.2
## Warning: package 'stringr' was built under R version 4.3.2
## Warning: package 'lubridate' was built under R version 4.3.2
## ── Attaching core tidyverse packages ──────────────────────── tidyverse 2.0.0 ──
## ✔ dplyr     1.1.4     ✔ readr     2.1.5
## ✔ forcats   1.0.0     ✔ stringr   1.5.1
## ✔ ggplot2   3.5.1     ✔ tibble    3.2.1
## ✔ lubridate 1.9.3     ✔ tidyr     1.3.1
## ✔ purrr     1.0.2     
## ── Conflicts ────────────────────────────────────────── tidyverse_conflicts() ──
## ✖ dplyr::filter() masks stats::filter()
## ✖ dplyr::lag()    masks stats::lag()
## ℹ Use the conflicted package (<http://conflicted.r-lib.org/>) to force all conflicts to become errors
library(ggplot2)
library(plotly)
## Warning: package 'plotly' was built under R version 4.3.3
## 
## Attaching package: 'plotly'
## 
## The following object is masked from 'package:ggplot2':
## 
##     last_plot
## 
## The following object is masked from 'package:stats':
## 
##     filter
## 
## The following object is masked from 'package:graphics':
## 
##     layout
library(GGally)
## Warning: package 'GGally' was built under R version 4.3.3
## Registered S3 method overwritten by 'GGally':
##   method from   
##   +.gg   ggplot2
data(penguins)

# Remove rows with missing values for a complete dataset
penguins_complete <- na.omit(penguins)

# Basic summary statistics
summary(penguins)
##       species          island    bill_length_mm  bill_depth_mm  
##  Adelie   :152   Biscoe   :168   Min.   :32.10   Min.   :13.10  
##  Chinstrap: 68   Dream    :124   1st Qu.:39.23   1st Qu.:15.60  
##  Gentoo   :124   Torgersen: 52   Median :44.45   Median :17.30  
##                                  Mean   :43.92   Mean   :17.15  
##                                  3rd Qu.:48.50   3rd Qu.:18.70  
##                                  Max.   :59.60   Max.   :21.50  
##                                  NA's   :2       NA's   :2      
##  flipper_length_mm  body_mass_g       sex           year     
##  Min.   :172.0     Min.   :2700   female:165   Min.   :2007  
##  1st Qu.:190.0     1st Qu.:3550   male  :168   1st Qu.:2007  
##  Median :197.0     Median :4050   NA's  : 11   Median :2008  
##  Mean   :200.9     Mean   :4202                Mean   :2008  
##  3rd Qu.:213.0     3rd Qu.:4750                3rd Qu.:2009  
##  Max.   :231.0     Max.   :6300                Max.   :2009  
##  NA's   :2         NA's   :2

Exploratory Data Analysis (EDA)

Basic Scatterplot

We start with a simple scatterplot to visualize the relationship between flipper length and body mass:

ggplot(data = penguins) +
  geom_point(mapping = aes(x = flipper_length_mm, y = body_mass_g)) +
  labs(title = "Penguins: Body Mass vs. Flipper Length") 
## Warning: Removed 2 rows containing missing values or values outside the scale range
## (`geom_point()`).

Detailed Scatterplot

Now, we enhance the scatterplot to incorporate species and island information:

ggplot(penguins_complete, aes(x = flipper_length_mm, y = body_mass_g, color = species, shape = island)) +
  geom_point(size = 3) +
  labs(title = "Body Mass vs. Flipper Length by Species and Island",
       x = "Flipper Length (mm)", y = "Body Mass (g)") +
  theme_minimal() +
  scale_color_viridis_d()

Bill Dimension Distributions

Next, we examine the distributions of bill length and depth, first with boxplots and then with violin plots enhanced with points:

# Boxplot
ggplot(penguins_complete, aes(x = species, y = bill_length_mm, fill = species)) +
  geom_boxplot() +
  labs(title = "Bill Length Distribution by Species",
       x = "Species", y = "Bill Length (mm)") +
  theme_classic() +
  scale_fill_brewer(palette = "Set2")

# Violin plot with overlaid points
ggplot(penguins_complete, aes(x = species, y = bill_depth_mm, fill = species)) +
  geom_violin(alpha = 0.6) +
  geom_jitter(width = 0.2, alpha = 0.8) +
  labs(title = "Bill Depth Distribution by Species",
       x = "Species", y = "Bill Depth (mm)") +
  theme_light() +
  scale_fill_manual(values = c("#E69F00", "#56B4E9", "#009E73"))

Faceting by Island

To further understand how bill dimensions vary across species and islands, we create faceted scatterplots:

ggplot(penguins_complete, aes(x = bill_length_mm, y = bill_depth_mm, color = species)) +
  geom_point() +
  facet_wrap(~ island) +
  labs(title = "Bill Dimensions by Species and Island",
       x = "Bill Length (mm)", y = "Bill Depth (mm)") +
  theme_bw()

Advanced Analysis & Statistical Tests

Interactive Plots with Trend Lines, Means, and Error Bars

We can enhance our plots by making them interactive using the plotly library and adding trend lines, mean points, and error bars:

# Scatterplot with Trend Line (geom_smooth())
p1 <- ggplot(penguins_complete, aes(x = flipper_length_mm, y = body_mass_g, color = species)) +
  geom_point(size = 3) +
  geom_smooth(method = "lm", se = TRUE) + # Linear regression with confidence interval
  labs(title = "Body Mass vs. Flipper Length with Trend Lines",
       x = "Flipper Length (mm)", y = "Body Mass (g)") +
  theme_minimal()

# Boxplot with Individual Points and Means
p2 <- ggplot(penguins_complete, aes(x = species, y = bill_length_mm, fill = species)) +
  geom_boxplot() +
  stat_summary(fun = mean, geom = "point", shape = 23, size = 4, fill = "white") + # Mean points
  labs(title = "Bill Length Distribution by Species with Means",
       x = "Species", y = "Bill Length (mm)") +
  theme_classic()

# Violin Plot with Individual Points and Error Bars
p3 <- ggplot(penguins_complete, aes(x = species, y = bill_depth_mm, fill = species)) +
  geom_violin() +
  geom_point(position = position_jitter(width = 0.2)) + # Add jitter to avoid overlap
  stat_summary(fun.data = mean_se, geom = "errorbar", width = 0.2) + # Error bars (mean +/- SE)
  labs(title = "Bill Depth Distribution by Species with Error Bars",
       x = "Species", y = "Bill Depth (mm)") +
  theme_light()

# Convert static ggplot2 plots to interactive plotly plots
ggplotly(p1)
## `geom_smooth()` using formula = 'y ~ x'
ggplotly(p2)
ggplotly(p3)

Scatterplot Matrix and Pairwise Relationships

Next, we use a scatterplot matrix (ggpairs) to visualize pairwise relationships among multiple variables and also create separate scatterplots for each species:

# Scatterplot Matrix with Correlation Coefficients
ggpairs(penguins_complete[, c("bill_length_mm", "bill_depth_mm", "flipper_length_mm", "body_mass_g")],
       upper = list(continuous = wrap("cor", method = "pearson")),
       lower = list(continuous = "points"))

# Pairwise Scatterplots with Trend Lines by Species
ggplot(penguins_complete, aes(x = flipper_length_mm, y = body_mass_g, color = species)) +
  geom_point() +
  geom_smooth(method = "lm") +
  facet_wrap(~species) +
  labs(title = "Pairwise Relationships by Species",
       x = "Flipper Length (mm)", y = "Body Mass (g)") +
  theme_bw()
## `geom_smooth()` using formula = 'y ~ x'

Statistical Tests: T-test and ANOVA

To assess whether observed differences are statistically significant, we perform a t-test for bill length between two species and an ANOVA for bill depth across all species:

# T-test for Bill Length Difference Between Two Species
t.test(bill_length_mm ~ species, data = penguins_complete, subset = species %in% c("Adelie", "Gentoo"))
## 
##  Welch Two Sample t-test
## 
## data:  bill_length_mm by species
## t = -24.286, df = 233.51, p-value < 2.2e-16
## alternative hypothesis: true difference in means between group Adelie and group Gentoo is not equal to 0
## 95 percent confidence interval:
##  -9.453448 -8.034741
## sample estimates:
## mean in group Adelie mean in group Gentoo 
##             38.82397             47.56807
# ANOVA for Bill Depth Differences Across All Species
bill_depth_anova <- aov(bill_depth_mm ~ species, data = penguins_complete)
summary(bill_depth_anova)
##              Df Sum Sq Mean Sq F value Pr(>F)    
## species       2  870.8   435.4   344.8 <2e-16 ***
## Residuals   330  416.7     1.3                   
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
# Post-hoc Tukey Test if ANOVA is significant
TukeyHSD(bill_depth_anova)
##   Tukey multiple comparisons of means
##     95% family-wise confidence level
## 
## Fit: aov(formula = bill_depth_mm ~ species, data = penguins_complete)
## 
## $species
##                         diff       lwr        upr     p adj
## Chinstrap-Adelie  0.07332796 -0.315078  0.4617339 0.8968734
## Gentoo-Adelie    -3.35062162 -3.677347 -3.0238962 0.0000000
## Gentoo-Chinstrap -3.42394958 -3.826113 -3.0217860 0.0000000

Conclusion

This comprehensive analysis provides insights into the relationships between penguin measurements, the variations across species and islands, and the statistical significance of these differences. The interactive visualizations enable deeper exploration, while the statistical tests quantify the evidence supporting the observed patterns.