ADwP - Mini Projekt

1 Przygotowania wstępne

Na początku mojego projektu załaduję wszystkie wykorzystywane podczas jego realizacji biblioteki. Przygotuję także dwie proste funkcje (kable oraz kable2) ułatwiające formatowanie i wizualizację danych.

library(kableExtra)
library(validate)
library(readr)
library(class)
library(rstatix)
library(ggpubr)
library(ggpmisc)
library(GGally)
library(tidyverse)

kable1 = function(data, digits=3) data %>% kable(digits = digits) %>% 
  kable_styling(bootstrap_options = 
                  c("striped", "hover", "condensed", "responsive"))

kable2 = function(data, digits=3, height="320px") data %>% 
  kable1(digits = digits) %>% 
  scroll_box(width = "100%", height = height)

2 Dane

2.1 Charakterystyka danych

W niniejszym projekcie będę się zajmować ramką danych zawierającą kilka wielkości fizycznych trzech gatunków pingwinów Chinstrap, Gentoo oraz Adelie.

Pingwiny badane były na trzech wyspach (Torgersen, Biscoe oraz Dream). Dokonano pomiarów wielkości fizycznych takich jak długość i szerokość dzioba, długość płetw, masa ciała oraz określono płeć badanych osobników. Drobnego wyjaśnienia wymaga określenie długości oraz szerokości dzioba. Sposób pomiaru tych wielkości został przedstawiony na poniższym rysunku.

Wszystkie wartości dotyczące długości są w milimetrach a masa jest w gramach.

2.2 Wczytanie danych

Dane pochodzą z dostarczonego do projektu pliku penguins.csv. Do ich wczytania użyję funkcji read_csv z pakietu readr. Podczas wczytania danych od razu zmienię także nazwy poszczególnych zmiennych.

pingwiny = read_csv(
  "penguins.csv", skip = 1, show_col_types = FALSE,
  col_names = c("gatunek", "wyspa", "długość.dzioba", "szerokość.dzioba", 
                "długość.płetw", "masa.ciała", "płeć"))
pingwiny %>% head(10) %>% kable2

gatunek	wyspa	długość.dzioba	szerokość.dzioba	długość.płetw	masa.ciała	płeć
Adelie	Torgersen	39.1	18.7	181	3750	MALE
Adelie	Torgersen	39.5	17.4	186	3800	FEMALE
Adelie	Torgersen	40.3	18.0	195	3250	FEMALE
Adelie	Torgersen	NA	NA	NA	NA	NA
Adelie	Torgersen	36.7	19.3	193	3450	FEMALE
Adelie	Torgersen	39.3	20.6	190	3650	MALE
Adelie	Torgersen	38.9	17.8	181	3625	FEMALE
Adelie	Torgersen	39.2	19.6	195	4675	MALE
Adelie	Torgersen	34.1	18.1	193	3475	NA
Adelie	Torgersen	42.0	20.2	190	4250	NA

2.3 Czyszczenie danych

Załadowane dane należy poddać procesowi oczyszczania. Poza standardowymi funkcjami z biblioteki dplyr żyję tu jednej własnej funkcji. Będzie to funkcja naKnn która wypełnia dowolnie wybraną zmienną zawierającą wartości NA wartościami wyznaczonymi metodą k najbliższych sąsiadów. Funkcja ta ma szereg zalet. Po pierwsze automatycznie konwertuje wszystkie wartości typu character w typ factor a następnie wszystkie zmienne typu factor konwertuje na typ liczbowy (wymagany przez funkcję knn). Tak przekonwertowane zmienne są następnie normalizowane. Funkcja ponadto dopuszcza aby w zmiennych obserwacyjnych występowały wartości NA, byle tylko nie występowały one w rekordach w których również i w zmiennej objaśnianej występują NA. Na koniec funkcja jest napisana w sposób przyjazny do użycia w potokach (pipe-friendly). Jako wynik zwracana jest ramka danych z uzupełnionymi brakującymi wartościami, która dodatkowo jest uzupełniana o dwa atrybuty. Atrybut knn.row zawiera numery wierszy w których dokonano wynaczynienia brakujących wartości. Natomiast atrybut knn.val zawiera wartości które zostały ustalone metodą knn.

Zdefiniuję także jedna drobną pomocniczą funkcję dsignif która zwraca maksymalną ilość cyfr znaczących z wszystkich wartości wektora x będącego jedynym argumentem tej prostej funkcji.

naKnn = function(df, var, ..., k = 5, l = 0){
  var = enquo(var)
  
  df.na = df %>% select(...) %>% na.omit() %>% attributes() %>% .$na.action
  if(df[df.na, quo_name(var)][[1]] %>% is.na %>% sum > 0) stop(
    paste("W wybranych zmiennych obserwacji wykryto warości NA.\n",
          "Zmienna objaśniana dla tych rekordów nie może być wyznaczona!\n"))

  cf.na = df[, quo_name(var)][[1]] %>% is.na() %>% which()
  if(length(cf.na)>0){
    dane <<- df %>%
      ungroup() %>%
      select(...) %>%
      mutate_if(is.character, factor) %>%
      mutate_if(is.factor, as.integer) %>%
      scale()
    
    cf = df[-c(df.na, cf.na), quo_name(var)][[1]] %>% factor()
    knn.value = class::knn(dane[-c(df.na, cf.na),], dane[cf.na,], cf, k=k, l=l)
    df[cf.na, quo_name(var)] = knn.value
    
    attributes(df)$knn.row = cf.na
    attributes(df)$knn.val = knn.value
  }
  df
}

dsignif = function(x) gsub(
  "(^0+|0+$)", "", sub("\\.", "", x|> na.omit() |> as.character())) |> 
  nchar() |> max()

Właściwe czyszczenie danych będzie wykonane w jednym krótkim potoku. Wykonam w nim następujące operacje:

• zmienne gatunek, wyspa oraz płeć (zmienne typu character) zostaną przekonwertowane na typ factor,

• dodatkowo etykiety zmiennej płeć zostaną zamienione na samiec oraz samica,

• dla wszystkich zmiennych liczbowych brakujące wartości zostaną uzupełnione wartościami średnimi wyliczonymi w obrębie danego gatunku i zaokrąglane do tylu cyfr znaczących ile mają wartości w danej zmiennej (dzięki funkcji dsignif),

• brakujące wartości płci będą wyznaczone metodą k-najbliższych sąsiadów co wykonam używając własnej funkcji naKnn.

pingwiny = pingwiny %>%
  mutate_if(is.character, factor) %>% 
  group_by(gatunek) %>%
  mutate_if(is.double, function(x){ 
    ifelse(is.na(x), mean(x, na.rm = T) %>% signif(dsignif(x)), x)}) %>% 
  mutate(płeć = płeć %>% factor(labels = c("samica", "samiec"))) %>%
  naKnn(płeć, gatunek, długość.dzioba:masa.ciała)

pingwiny %>% kable2

gatunek	wyspa	długość.dzioba	szerokość.dzioba	długość.płetw	masa.ciała	płeć
Adelie	Torgersen	39.1	18.7	181	3750	samiec
Adelie	Torgersen	39.5	17.4	186	3800	samica
Adelie	Torgersen	40.3	18.0	195	3250	samica
Adelie	Torgersen	38.8	18.3	190	3701	samiec
Adelie	Torgersen	36.7	19.3	193	3450	samica
Adelie	Torgersen	39.3	20.6	190	3650	samiec
Adelie	Torgersen	38.9	17.8	181	3625	samica
Adelie	Torgersen	39.2	19.6	195	4675	samiec
Adelie	Torgersen	34.1	18.1	193	3475	samica
Adelie	Torgersen	42.0	20.2	190	4250	samiec
Adelie	Torgersen	37.8	17.1	186	3300	samica
Adelie	Torgersen	37.8	17.3	180	3700	samica
Adelie	Torgersen	41.1	17.6	182	3200	samica
Adelie	Torgersen	38.6	21.2	191	3800	samiec
Adelie	Torgersen	34.6	21.1	198	4400	samiec
Adelie	Torgersen	36.6	17.8	185	3700	samica
Adelie	Torgersen	38.7	19.0	195	3450	samica
Adelie	Torgersen	42.5	20.7	197	4500	samiec
Adelie	Torgersen	34.4	18.4	184	3325	samica
Adelie	Torgersen	46.0	21.5	194	4200	samiec
Adelie	Biscoe	37.8	18.3	174	3400	samica
Adelie	Biscoe	37.7	18.7	180	3600	samiec
Adelie	Biscoe	35.9	19.2	189	3800	samica
Adelie	Biscoe	38.2	18.1	185	3950	samiec
Adelie	Biscoe	38.8	17.2	180	3800	samiec
Adelie	Biscoe	35.3	18.9	187	3800	samica
Adelie	Biscoe	40.6	18.6	183	3550	samiec
Adelie	Biscoe	40.5	17.9	187	3200	samica
Adelie	Biscoe	37.9	18.6	172	3150	samica
Adelie	Biscoe	40.5	18.9	180	3950	samiec
Adelie	Dream	39.5	16.7	178	3250	samica
Adelie	Dream	37.2	18.1	178	3900	samiec
Adelie	Dream	39.5	17.8	188	3300	samica
Adelie	Dream	40.9	18.9	184	3900	samiec
Adelie	Dream	36.4	17.0	195	3325	samica
Adelie	Dream	39.2	21.1	196	4150	samiec
Adelie	Dream	38.8	20.0	190	3950	samiec
Adelie	Dream	42.2	18.5	180	3550	samica
Adelie	Dream	37.6	19.3	181	3300	samica
Adelie	Dream	39.8	19.1	184	4650	samiec
Adelie	Dream	36.5	18.0	182	3150	samica
Adelie	Dream	40.8	18.4	195	3900	samiec
Adelie	Dream	36.0	18.5	186	3100	samica
Adelie	Dream	44.1	19.7	196	4400	samiec
Adelie	Dream	37.0	16.9	185	3000	samica
Adelie	Dream	39.6	18.8	190	4600	samiec
Adelie	Dream	41.1	19.0	182	3425	samiec
Adelie	Dream	37.5	18.9	179	2975	samica
Adelie	Dream	36.0	17.9	190	3450	samica
Adelie	Dream	42.3	21.2	191	4150	samiec
Adelie	Biscoe	39.6	17.7	186	3500	samica
Adelie	Biscoe	40.1	18.9	188	4300	samiec
Adelie	Biscoe	35.0	17.9	190	3450	samica
Adelie	Biscoe	42.0	19.5	200	4050	samiec
Adelie	Biscoe	34.5	18.1	187	2900	samica
Adelie	Biscoe	41.4	18.6	191	3700	samiec
Adelie	Biscoe	39.0	17.5	186	3550	samica
Adelie	Biscoe	40.6	18.8	193	3800	samiec
Adelie	Biscoe	36.5	16.6	181	2850	samica
Adelie	Biscoe	37.6	19.1	194	3750	samiec
Adelie	Biscoe	35.7	16.9	185	3150	samica
Adelie	Biscoe	41.3	21.1	195	4400	samiec
Adelie	Biscoe	37.6	17.0	185	3600	samica
Adelie	Biscoe	41.1	18.2	192	4050	samiec
Adelie	Biscoe	36.4	17.1	184	2850	samica
Adelie	Biscoe	41.6	18.0	192	3950	samiec
Adelie	Biscoe	35.5	16.2	195	3350	samica
Adelie	Biscoe	41.1	19.1	188	4100	samiec
Adelie	Torgersen	35.9	16.6	190	3050	samica
Adelie	Torgersen	41.8	19.4	198	4450	samiec
Adelie	Torgersen	33.5	19.0	190	3600	samica
Adelie	Torgersen	39.7	18.4	190	3900	samiec
Adelie	Torgersen	39.6	17.2	196	3550	samica
Adelie	Torgersen	45.8	18.9	197	4150	samiec
Adelie	Torgersen	35.5	17.5	190	3700	samica
Adelie	Torgersen	42.8	18.5	195	4250	samiec
Adelie	Torgersen	40.9	16.8	191	3700	samica
Adelie	Torgersen	37.2	19.4	184	3900	samiec
Adelie	Torgersen	36.2	16.1	187	3550	samica
Adelie	Torgersen	42.1	19.1	195	4000	samiec
Adelie	Torgersen	34.6	17.2	189	3200	samica
Adelie	Torgersen	42.9	17.6	196	4700	samiec
Adelie	Torgersen	36.7	18.8	187	3800	samica
Adelie	Torgersen	35.1	19.4	193	4200	samiec
Adelie	Dream	37.3	17.8	191	3350	samica
Adelie	Dream	41.3	20.3	194	3550	samiec
Adelie	Dream	36.3	19.5	190	3800	samiec
Adelie	Dream	36.9	18.6	189	3500	samica
Adelie	Dream	38.3	19.2	189	3950	samiec
Adelie	Dream	38.9	18.8	190	3600	samica
Adelie	Dream	35.7	18.0	202	3550	samica
Adelie	Dream	41.1	18.1	205	4300	samiec
Adelie	Dream	34.0	17.1	185	3400	samica
Adelie	Dream	39.6	18.1	186	4450	samiec
Adelie	Dream	36.2	17.3	187	3300	samica
Adelie	Dream	40.8	18.9	208	4300	samiec
Adelie	Dream	38.1	18.6	190	3700	samica
Adelie	Dream	40.3	18.5	196	4350	samiec
Adelie	Dream	33.1	16.1	178	2900	samica
Adelie	Dream	43.2	18.5	192	4100	samiec
Adelie	Biscoe	35.0	17.9	192	3725	samica
Adelie	Biscoe	41.0	20.0	203	4725	samiec
Adelie	Biscoe	37.7	16.0	183	3075	samica
Adelie	Biscoe	37.8	20.0	190	4250	samiec
Adelie	Biscoe	37.9	18.6	193	2925	samica
Adelie	Biscoe	39.7	18.9	184	3550	samiec
Adelie	Biscoe	38.6	17.2	199	3750	samica
Adelie	Biscoe	38.2	20.0	190	3900	samiec
Adelie	Biscoe	38.1	17.0	181	3175	samica
Adelie	Biscoe	43.2	19.0	197	4775	samiec
Adelie	Biscoe	38.1	16.5	198	3825	samica
Adelie	Biscoe	45.6	20.3	191	4600	samiec
Adelie	Biscoe	39.7	17.7	193	3200	samica
Adelie	Biscoe	42.2	19.5	197	4275	samiec
Adelie	Biscoe	39.6	20.7	191	3900	samica
Adelie	Biscoe	42.7	18.3	196	4075	samiec
Adelie	Torgersen	38.6	17.0	188	2900	samica
Adelie	Torgersen	37.3	20.5	199	3775	samiec
Adelie	Torgersen	35.7	17.0	189	3350	samica
Adelie	Torgersen	41.1	18.6	189	3325	samiec
Adelie	Torgersen	36.2	17.2	187	3150	samica
Adelie	Torgersen	37.7	19.8	198	3500	samiec
Adelie	Torgersen	40.2	17.0	176	3450	samica
Adelie	Torgersen	41.4	18.5	202	3875	samiec
Adelie	Torgersen	35.2	15.9	186	3050	samica
Adelie	Torgersen	40.6	19.0	199	4000	samiec
Adelie	Torgersen	38.8	17.6	191	3275	samica
Adelie	Torgersen	41.5	18.3	195	4300	samiec
Adelie	Torgersen	39.0	17.1	191	3050	samica
Adelie	Torgersen	44.1	18.0	210	4000	samiec
Adelie	Torgersen	38.5	17.9	190	3325	samica
Adelie	Torgersen	43.1	19.2	197	3500	samiec
Adelie	Dream	36.8	18.5	193	3500	samica
Adelie	Dream	37.5	18.5	199	4475	samiec
Adelie	Dream	38.1	17.6	187	3425	samica
Adelie	Dream	41.1	17.5	190	3900	samiec
Adelie	Dream	35.6	17.5	191	3175	samica
Adelie	Dream	40.2	20.1	200	3975	samiec
Adelie	Dream	37.0	16.5	185	3400	samica
Adelie	Dream	39.7	17.9	193	4250	samiec
Adelie	Dream	40.2	17.1	193	3400	samica
Adelie	Dream	40.6	17.2	187	3475	samiec
Adelie	Dream	32.1	15.5	188	3050	samica
Adelie	Dream	40.7	17.0	190	3725	samiec
Adelie	Dream	37.3	16.8	192	3000	samica
Adelie	Dream	39.0	18.7	185	3650	samiec
Adelie	Dream	39.2	18.6	190	4250	samiec
Adelie	Dream	36.6	18.4	184	3475	samica
Adelie	Dream	36.0	17.8	195	3450	samica
Adelie	Dream	37.8	18.1	193	3750	samiec
Adelie	Dream	36.0	17.1	187	3700	samica
Adelie	Dream	41.5	18.5	201	4000	samiec
Chinstrap	Dream	46.5	17.9	192	3500	samica
Chinstrap	Dream	50.0	19.5	196	3900	samiec
Chinstrap	Dream	51.3	19.2	193	3650	samiec
Chinstrap	Dream	45.4	18.7	188	3525	samica
Chinstrap	Dream	52.7	19.8	197	3725	samiec
Chinstrap	Dream	45.2	17.8	198	3950	samica
Chinstrap	Dream	46.1	18.2	178	3250	samica
Chinstrap	Dream	51.3	18.2	197	3750	samiec
Chinstrap	Dream	46.0	18.9	195	4150	samica
Chinstrap	Dream	51.3	19.9	198	3700	samiec
Chinstrap	Dream	46.6	17.8	193	3800	samica
Chinstrap	Dream	51.7	20.3	194	3775	samiec
Chinstrap	Dream	47.0	17.3	185	3700	samica
Chinstrap	Dream	52.0	18.1	201	4050	samiec
Chinstrap	Dream	45.9	17.1	190	3575	samica
Chinstrap	Dream	50.5	19.6	201	4050	samiec
Chinstrap	Dream	50.3	20.0	197	3300	samiec
Chinstrap	Dream	58.0	17.8	181	3700	samica
Chinstrap	Dream	46.4	18.6	190	3450	samica
Chinstrap	Dream	49.2	18.2	195	4400	samiec
Chinstrap	Dream	42.4	17.3	181	3600	samica
Chinstrap	Dream	48.5	17.5	191	3400	samiec
Chinstrap	Dream	43.2	16.6	187	2900	samica
Chinstrap	Dream	50.6	19.4	193	3800	samiec
Chinstrap	Dream	46.7	17.9	195	3300	samica
Chinstrap	Dream	52.0	19.0	197	4150	samiec
Chinstrap	Dream	50.5	18.4	200	3400	samica
Chinstrap	Dream	49.5	19.0	200	3800	samiec
Chinstrap	Dream	46.4	17.8	191	3700	samica
Chinstrap	Dream	52.8	20.0	205	4550	samiec
Chinstrap	Dream	40.9	16.6	187	3200	samica
Chinstrap	Dream	54.2	20.8	201	4300	samiec
Chinstrap	Dream	42.5	16.7	187	3350	samica
Chinstrap	Dream	51.0	18.8	203	4100	samiec
Chinstrap	Dream	49.7	18.6	195	3600	samiec
Chinstrap	Dream	47.5	16.8	199	3900	samica
Chinstrap	Dream	47.6	18.3	195	3850	samica
Chinstrap	Dream	52.0	20.7	210	4800	samiec
Chinstrap	Dream	46.9	16.6	192	2700	samica
Chinstrap	Dream	53.5	19.9	205	4500	samiec
Chinstrap	Dream	49.0	19.5	210	3950	samiec
Chinstrap	Dream	46.2	17.5	187	3650	samica
Chinstrap	Dream	50.9	19.1	196	3550	samiec
Chinstrap	Dream	45.5	17.0	196	3500	samica
Chinstrap	Dream	50.9	17.9	196	3675	samica
Chinstrap	Dream	50.8	18.5	201	4450	samiec
Chinstrap	Dream	50.1	17.9	190	3400	samica
Chinstrap	Dream	49.0	19.6	212	4300	samiec
Chinstrap	Dream	51.5	18.7	187	3250	samiec
Chinstrap	Dream	49.8	17.3	198	3675	samica
Chinstrap	Dream	48.1	16.4	199	3325	samica
Chinstrap	Dream	51.4	19.0	201	3950	samiec
Chinstrap	Dream	45.7	17.3	193	3600	samica
Chinstrap	Dream	50.7	19.7	203	4050	samiec
Chinstrap	Dream	42.5	17.3	187	3350	samica
Chinstrap	Dream	52.2	18.8	197	3450	samiec
Chinstrap	Dream	45.2	16.6	191	3250	samica
Chinstrap	Dream	49.3	19.9	203	4050	samiec
Chinstrap	Dream	50.2	18.8	202	3800	samiec
Chinstrap	Dream	45.6	19.4	194	3525	samica
Chinstrap	Dream	51.9	19.5	206	3950	samiec
Chinstrap	Dream	46.8	16.5	189	3650	samica
Chinstrap	Dream	45.7	17.0	195	3650	samica
Chinstrap	Dream	55.8	19.8	207	4000	samiec
Chinstrap	Dream	43.5	18.1	202	3400	samica
Chinstrap	Dream	49.6	18.2	193	3775	samiec
Chinstrap	Dream	50.8	19.0	210	4100	samiec
Chinstrap	Dream	50.2	18.7	198	3775	samica
Gentoo	Biscoe	46.1	13.2	211	4500	samica
Gentoo	Biscoe	50.0	16.3	230	5700	samiec
Gentoo	Biscoe	48.7	14.1	210	4450	samica
Gentoo	Biscoe	50.0	15.2	218	5700	samiec
Gentoo	Biscoe	47.6	14.5	215	5400	samiec
Gentoo	Biscoe	46.5	13.5	210	4550	samica
Gentoo	Biscoe	45.4	14.6	211	4800	samica
Gentoo	Biscoe	46.7	15.3	219	5200	samiec
Gentoo	Biscoe	43.3	13.4	209	4400	samica
Gentoo	Biscoe	46.8	15.4	215	5150	samiec
Gentoo	Biscoe	40.9	13.7	214	4650	samica
Gentoo	Biscoe	49.0	16.1	216	5550	samiec
Gentoo	Biscoe	45.5	13.7	214	4650	samica
Gentoo	Biscoe	48.4	14.6	213	5850	samiec
Gentoo	Biscoe	45.8	14.6	210	4200	samica
Gentoo	Biscoe	49.3	15.7	217	5850	samiec
Gentoo	Biscoe	42.0	13.5	210	4150	samica
Gentoo	Biscoe	49.2	15.2	221	6300	samiec
Gentoo	Biscoe	46.2	14.5	209	4800	samica
Gentoo	Biscoe	48.7	15.1	222	5350	samiec
Gentoo	Biscoe	50.2	14.3	218	5700	samiec
Gentoo	Biscoe	45.1	14.5	215	5000	samica
Gentoo	Biscoe	46.5	14.5	213	4400	samica
Gentoo	Biscoe	46.3	15.8	215	5050	samiec
Gentoo	Biscoe	42.9	13.1	215	5000	samica
Gentoo	Biscoe	46.1	15.1	215	5100	samiec
Gentoo	Biscoe	44.5	14.3	216	4100	samica
Gentoo	Biscoe	47.8	15.0	215	5650	samiec
Gentoo	Biscoe	48.2	14.3	210	4600	samica
Gentoo	Biscoe	50.0	15.3	220	5550	samiec
Gentoo	Biscoe	47.3	15.3	222	5250	samiec
Gentoo	Biscoe	42.8	14.2	209	4700	samica
Gentoo	Biscoe	45.1	14.5	207	5050	samica
Gentoo	Biscoe	59.6	17.0	230	6050	samiec
Gentoo	Biscoe	49.1	14.8	220	5150	samica
Gentoo	Biscoe	48.4	16.3	220	5400	samiec
Gentoo	Biscoe	42.6	13.7	213	4950	samica
Gentoo	Biscoe	44.4	17.3	219	5250	samiec
Gentoo	Biscoe	44.0	13.6	208	4350	samica
Gentoo	Biscoe	48.7	15.7	208	5350	samiec
Gentoo	Biscoe	42.7	13.7	208	3950	samica
Gentoo	Biscoe	49.6	16.0	225	5700	samiec
Gentoo	Biscoe	45.3	13.7	210	4300	samica
Gentoo	Biscoe	49.6	15.0	216	4750	samiec
Gentoo	Biscoe	50.5	15.9	222	5550	samiec
Gentoo	Biscoe	43.6	13.9	217	4900	samica
Gentoo	Biscoe	45.5	13.9	210	4200	samica
Gentoo	Biscoe	50.5	15.9	225	5400	samiec
Gentoo	Biscoe	44.9	13.3	213	5100	samica
Gentoo	Biscoe	45.2	15.8	215	5300	samiec
Gentoo	Biscoe	46.6	14.2	210	4850	samica
Gentoo	Biscoe	48.5	14.1	220	5300	samiec
Gentoo	Biscoe	45.1	14.4	210	4400	samica
Gentoo	Biscoe	50.1	15.0	225	5000	samiec
Gentoo	Biscoe	46.5	14.4	217	4900	samica
Gentoo	Biscoe	45.0	15.4	220	5050	samiec
Gentoo	Biscoe	43.8	13.9	208	4300	samica
Gentoo	Biscoe	45.5	15.0	220	5000	samiec
Gentoo	Biscoe	43.2	14.5	208	4450	samica
Gentoo	Biscoe	50.4	15.3	224	5550	samiec
Gentoo	Biscoe	45.3	13.8	208	4200	samica
Gentoo	Biscoe	46.2	14.9	221	5300	samiec
Gentoo	Biscoe	45.7	13.9	214	4400	samica
Gentoo	Biscoe	54.3	15.7	231	5650	samiec
Gentoo	Biscoe	45.8	14.2	219	4700	samica
Gentoo	Biscoe	49.8	16.8	230	5700	samiec
Gentoo	Biscoe	46.2	14.4	214	4650	samica
Gentoo	Biscoe	49.5	16.2	229	5800	samiec
Gentoo	Biscoe	43.5	14.2	220	4700	samica
Gentoo	Biscoe	50.7	15.0	223	5550	samiec
Gentoo	Biscoe	47.7	15.0	216	4750	samica
Gentoo	Biscoe	46.4	15.6	221	5000	samiec
Gentoo	Biscoe	48.2	15.6	221	5100	samiec
Gentoo	Biscoe	46.5	14.8	217	5200	samica
Gentoo	Biscoe	46.4	15.0	216	4700	samica
Gentoo	Biscoe	48.6	16.0	230	5800	samiec
Gentoo	Biscoe	47.5	14.2	209	4600	samica
Gentoo	Biscoe	51.1	16.3	220	6000	samiec
Gentoo	Biscoe	45.2	13.8	215	4750	samica
Gentoo	Biscoe	45.2	16.4	223	5950	samiec
Gentoo	Biscoe	49.1	14.5	212	4625	samica
Gentoo	Biscoe	52.5	15.6	221	5450	samiec
Gentoo	Biscoe	47.4	14.6	212	4725	samica
Gentoo	Biscoe	50.0	15.9	224	5350	samiec
Gentoo	Biscoe	44.9	13.8	212	4750	samica
Gentoo	Biscoe	50.8	17.3	228	5600	samiec
Gentoo	Biscoe	43.4	14.4	218	4600	samica
Gentoo	Biscoe	51.3	14.2	218	5300	samiec
Gentoo	Biscoe	47.5	14.0	212	4875	samica
Gentoo	Biscoe	52.1	17.0	230	5550	samiec
Gentoo	Biscoe	47.5	15.0	218	4950	samica
Gentoo	Biscoe	52.2	17.1	228	5400	samiec
Gentoo	Biscoe	45.5	14.5	212	4750	samica
Gentoo	Biscoe	49.5	16.1	224	5650	samiec
Gentoo	Biscoe	44.5	14.7	214	4850	samica
Gentoo	Biscoe	50.8	15.7	226	5200	samiec
Gentoo	Biscoe	49.4	15.8	216	4925	samiec
Gentoo	Biscoe	46.9	14.6	222	4875	samica
Gentoo	Biscoe	48.4	14.4	203	4625	samica
Gentoo	Biscoe	51.1	16.5	225	5250	samiec
Gentoo	Biscoe	48.5	15.0	219	4850	samica
Gentoo	Biscoe	55.9	17.0	228	5600	samiec
Gentoo	Biscoe	47.2	15.5	215	4975	samica
Gentoo	Biscoe	49.1	15.0	228	5500	samiec
Gentoo	Biscoe	47.3	13.8	216	4725	samica
Gentoo	Biscoe	46.8	16.1	215	5500	samiec
Gentoo	Biscoe	41.7	14.7	210	4700	samica
Gentoo	Biscoe	53.4	15.8	219	5500	samiec
Gentoo	Biscoe	43.3	14.0	208	4575	samica
Gentoo	Biscoe	48.1	15.1	209	5500	samiec
Gentoo	Biscoe	50.5	15.2	216	5000	samica
Gentoo	Biscoe	49.8	15.9	229	5950	samiec
Gentoo	Biscoe	43.5	15.2	213	4650	samica
Gentoo	Biscoe	51.5	16.3	230	5500	samiec
Gentoo	Biscoe	46.2	14.1	217	4375	samica
Gentoo	Biscoe	55.1	16.0	230	5850	samiec
Gentoo	Biscoe	44.5	15.7	217	4875	samiec
Gentoo	Biscoe	48.8	16.2	222	6000	samiec
Gentoo	Biscoe	47.2	13.7	214	4925	samica
Gentoo	Biscoe	47.5	15.0	217	5076	samiec
Gentoo	Biscoe	46.8	14.3	215	4850	samica
Gentoo	Biscoe	50.4	15.7	222	5750	samiec
Gentoo	Biscoe	45.2	14.8	212	5200	samica
Gentoo	Biscoe	49.9	16.1	213	5400	samiec

2.4 Walidacja danych

Po wyczyszczeniu danych warto wykonać także prostą walidację aby przekonać się, że w danych nie są ukryte jakieś ewidentnie niepoprawne wartości. Walidację będę wykonywał przy użyciu funkcji z pakietu validate.

Przyjąłem następujące warunki walidacyjne:

• gatunek musi zawierać jedna z trzech wartości Adelie, Chinstrap, Gentoo,

• wyspa musi zawierać jedna z trzech wartości Biscoe, Dream, Torgersen,

• długość.dzioba musi być typu numerycznego i mieścić się w zakresie [20, 100],

• szerokość.dzioba musi być typu numerycznego i mieścić się w zakresie [10, 50],

• długość.płet musi być typu numerycznego i mieścić się w zakresie [100, 300],

• masa.ciała musi być typu numerycznego i mieścić się w zakresie [1000, 10000],

• płeć musi zawierać jedna z dwóch wartości samica, samiec.

Ponadto w żadnej ze zmiennych nie powinny już występować wartości NA.

pingwiny %>% 
  confront(
    validator(
      gatunek.list = gatunek %in% c('Adelie', 'Chinstrap', 'Gentoo'),
      wyspa.list = wyspa %in% c('Biscoe', 'Dream', 'Torgersen'),
      długość.dzioba.num = is.numeric(długość.dzioba),
      długość.dzioba.zakres = długość.dzioba >= 20 & długość.dzioba <= 100,
      szerokość.dzioba.num = is.numeric(szerokość.dzioba),
      szerokość.dzioba.zakres = szerokość.dzioba >= 10 & szerokość.dzioba <= 50,
      długość.płetw.num = is.numeric(długość.płetw),
      długość.płetw.zakres = długość.płetw >= 100 & długość.płetw <= 300,
      masa.ciała.num = is.numeric(masa.ciała),
      masa.ciała.zakres = masa.ciała >= 1000 & masa.ciała <= 10000,
      płeć.list = płeć %in% c('samica', 'samiec')
    )
  ) %>% summary %>% select(-expression) %>% kable1

name	items	passes	fails	nNA	error	warning
gatunek.list	344	344	0	0	FALSE	FALSE
wyspa.list	344	344	0	0	FALSE	FALSE
długość.dzioba.num	1	1	0	0	FALSE	FALSE
długość.dzioba.zakres	344	344	0	0	FALSE	FALSE
szerokość.dzioba.num	1	1	0	0	FALSE	FALSE
szerokość.dzioba.zakres	344	344	0	0	FALSE	FALSE
długość.płetw.num	1	1	0	0	FALSE	FALSE
długość.płetw.zakres	344	344	0	0	FALSE	FALSE
masa.ciała.num	1	1	0	0	FALSE	FALSE
masa.ciała.zakres	344	344	0	0	FALSE	FALSE
płeć.list	344	344	0	0	FALSE	FALSE

Po upewnieniu się, że wszystkie zmienne zawierają dane poprawnego typu oraz z poprawnego zakresu można było przejść do ich analizy.

3 Analiza

3.1 Funkcja podsumowująca

Zanim przejdę do analizy danych, przygotuję kilka własnych funkcji ułatwiających i upraszczających wszystkie planowane operacje. Pierwszą będzie funkcja przygotowująca podsumowanie sum_stats. Funkcja ta poza podstawowymi statystykami takimi jak średnia, media, min, max, kurtoza, skośność itp od razu testuje podgrupy wartości pod kątem zgodności z rozkładem normalnym, zwracając statystykę W oraz wartość p testu Shapiro-Wilka.

sum_stats = function(data, value){
  ShapiroTest = function(d, alpha=0.05){
    if(length(d)>5000) d=sample(d, 5000)
    if(length(d)>=3 & length(d)<=5000){
      sw = shapiro.test(d)
      tibble(
        SW.stat = sw$statistic,
        SW.p = sw$p.value,
        SW.test = sw$p.value>alpha)
    } else tibble(
      SW.stat = NA,
      SW.p = NA,
      SW.test = NA)
  }
  
  value = enquo(value)
  data %>%
  summarise(
    n = n(),
    nout = length(boxplot.stats(!!value)$out),
    q1 = quantile(!!value,1/4,8),
    min = min(!!value),
    mean = mean(!!value),
    median = median(!!value),
    q3 = quantile(!!value,3/4,8),
    max = max(!!value),
    sd = sd(!!value),
    kurtosis = e1071::kurtosis(!!value),
    skewness = e1071::skewness(!!value),
    ShapiroTest(!!value),
    .groups = "drop"
  )
}

3.2 Wizualziacja

Kolejne dwie funkcje będą pomocne w przygotowaniu odpowiednich wykresów. Pierwszą z nich jest funkcja ggBoxplot która tworzy połączony wykres box-plot z wykresem wiolinowym. Dodatkowo funkcja wykonuje niejako wewnętrznie odpowiednie testy na istotność różnic. Testy mogą być wykonywane w dwóch wariantach. Jeżeli parametr parametric ma wartość TRUE wykonywany jest test ANOVA a dla każdej pary poszczególnych podgrup wartości jest wykonywany test t-Studenta z poprawką Benferroniego. Wyniki testów są nanoszone na wykres. Jeżeli jednak wartości nie spełniają założeń wymaganych dla testów parametrycznych wykonywane są testy nieparametryczne czyli test Kruskal-Wallis oraz pomiędzy parami wykonywany jest test Wilcoxona z poprawką Holma.

Druga funkcja (która bazuje na ggBoxplot) ggGroupBoxplot wykonuje podobne wykresy z tym, że generuje kilka wykresów dla wartości rozdzielonych na odpowiednie podgrupy. Oczywiście obie funkcje są potoko-przyjazne co zdecydowanie ułatwia ich późniejsze wykorzystanie.

ggBoxplot = function(data, 
                     x, y,
                     parametric = T,
                     add.test = T,
                     add.pwc = T,
                     step.increase = 0,
                     title = NULL,
                     caption = "M.F."){
  x = enquo(x)
  y = enquo(y)
  data = data %>% ungroup()
  
  p = data %>% ggplot(aes(!!x, !!y, fill=!!x))+  #wykres właściwy
    geom_violin(alpha=0.8)+
    geom_boxplot(outlier.shape = 23, outlier.size = 3, outlier.alpha = 1,
                 alpha = 0.6, width = .5, fill="white")+
  scale_colour_manual(values = c("darkorange","purple","cyan4")) +
  scale_fill_manual(values = c("darkorange","purple","cyan4"))
  
  if(data %>% distinct(!!x) %>% nrow()>1){
    formula = formula(paste(quo_name(y), "~", quo_name(x)))
    if(parametric) {
      pwc = data %>%  #przygotowanie wyników testu ANOVA 
        #oraz t-Studenta na wszystkich parach podgrup
        pairwise_t_test(formula, p.adjust.method = "bonferroni") %>%
        add_xy_position(x = quo_name(x)) %>%
        mutate(!!x := group1, lab = paste(p, " - ", p.adj.signif))
      res.test = data %>% anova_test(formula)
    } else {
      pwc = data %>%  #przygotowanie wyników testu Kruskal-Wallis 
        #oraz Wilcoxona na wszystkich parach podgrup
        pairwise_wilcox_test(formula) %>%
        add_xy_position(x = quo_name(x)) %>%
        mutate(!!x := group1, lab = paste(p, " - ", p.adj.signif))
      res.test = data %>% kruskal_test(formula)
    }
    #dodanie wyników testów do wykresu
    p = p + stat_pvalue_manual(pwc, step.increase=step.increase, label = "lab")
    if(add.test) p = p + labs(title = get_test_label(res.test, detailed = TRUE))
    if(add.pwc) p = p + labs(subtitle = get_pwc_label(pwc))
  } else {
    if(add.test) p = p + labs(title = "Brak testu")
    if(add.pwc) p = p + labs(subtitle = "pwc: ")
  }
  if(!is.null(title)) p = p %>% annotate_figure(
    top = text_grob(title, size = 16, x=0, just="left"))
  if(!is.null(caption)) p = p %>% annotate_figure(
    bottom = text_grob(caption, size = 10, just="right", x = 1))
  p
}

ggGroupBoxplot = function(data, 
                          group, x, y,
                          parametric = T,
                          add.test = F,
                          add.pwc = T,
                          step.increase = 0,
                          ncol = NULL,
                          nrow = NULL,
                          title = NULL,
                          caption = "M.F."){
  group = enquo(group)
  x = enquo(x)
  y = enquo(y)
  data = data %>% ungroup()
  
  plots = list()
  for(g in data %>% distinct(!!group) %>% pull(!!group)){
    plots[[g]] = data %>% filter(!!group == g) %>%
      ggBoxplot(!!x, !!y, parametric, add.test, add.pwc, 
                step.increase, title = NULL, caption = NULL) %>%
      facet(facet.by = quo_name(group))
  }
  
  p = ggarrange(plotlist=plots, ncol=ncol, nrow=nrow,
                common.legend = TRUE, legend="bottom") 
  if(!is.null(title)) p = p %>% annotate_figure(
    top = text_grob(title, size = 16, x=0, just="left"))
  if(!is.null(caption)) p = p %>% annotate_figure(
    bottom = text_grob(caption, size = 10, just="right", x = 1))
  p
}

3.3 Długość dzioba

Pierwszą wielkością którą będę analizował będzie długość dzioba w zależności od gatunku badanych pingwinów. Na początku przyjrzyjmy się podstawowym statystykom.

pingwiny %>% sum_stats(długość.dzioba) %>% 
  mutate_if(is.numeric, signif, 3) %>% kable1

gatunek	n	nout	q1	min	mean	median	q3	max	sd	kurtosis	skewness	SW.stat	SW.p	SW.test
Adelie	152	0	36.8	32.1	38.8	38.8	40.7	46.0	2.65	-0.21	0.159	0.994	0.742	TRUE
Chinstrap	68	0	46.3	40.9	48.8	49.6	51.1	58.0	3.34	-0.14	-0.087	0.975	0.194	TRUE
Gentoo	124	1	45.3	40.9	47.5	47.3	49.5	59.6	3.07	1.16	0.638	0.973	0.013	FALSE

Niestety, jak można zauważyć dane dla gatunku Gentoo nie są zgodne z rozkładem normalnym. Nie będzie więc można stosować testów parametrycznych.

pingwiny %>% ggBoxplot(
  gatunek, długość.dzioba, F, 
  title = "Rozkład dzługości dzioba dla poszczególnych gatunków pingwinów")

Jak widać różnice pomiędzy gatunkami są istotne statystycznie w każdym przypadku. Rzuca się jednak w oczy dość “pagórkowaty” rozkład widoczny bardzo wyraźnie w przypadku gatunku Chinstrap i nieco mniej wyraźnie w przypadku Gentoo. Można podejrzewać, że jest to spowodowane zebraniem w jedną grupę danych dla obu płci.

Sprawdźmy zatem co się stanie jeżeli wprowadzimy dodatkową zmienną grupującą którą będzie płeć.

pingwiny %>% group_by(gatunek, płeć) %>% 
  sum_stats(długość.dzioba) %>% 
  mutate_if(is.numeric, signif, 3) %>% kable1

gatunek	płeć	n	nout	q1	min	mean	median	q3	max	sd	kurtosis	skewness	SW.stat	SW.p	SW.test
Adelie	samica	77	0	35.9	32.1	37.2	37.0	38.7	42.2	2.01	-0.343	0.012	0.994	0.984	TRUE
Adelie	samiec	75	4	38.9	34.6	40.4	40.6	41.6	46.0	2.26	0.183	0.055	0.987	0.629	TRUE
Chinstrap	samica	34	3	45.4	40.9	46.6	46.3	47.4	58.0	3.11	3.360	1.270	0.883	0.002	FALSE
Chinstrap	samiec	34	1	50.0	48.5	51.1	51.0	52.0	55.8	1.56	0.736	0.769	0.955	0.177	TRUE
Gentoo	samica	61	0	44.0	40.9	45.6	45.5	46.9	50.5	2.02	-0.443	-0.052	0.991	0.927	TRUE
Gentoo	samiec	63	3	47.7	44.4	49.4	49.5	50.5	59.6	2.76	1.920	0.844	0.944	0.006	FALSE

Tym razem niezgodność z rozkładem normalnym jest w przypadku samców Gentoo oraz samic Chinstrap. Dodatkowo jest znacznie więcej wartości odstających. Porównajmy zatem, czy w obrębie jednego gatunku są istotne różnice w długości dzioba ze względu na płeć.

pingwiny %>% ggGroupBoxplot(
  gatunek, płeć, długość.dzioba, F, ncol=3, 
  title="Rozkład długości dzioba dla trzech gatunków w zależności od płci")

Jak widać na powyższym wykresie, rzeczywiście w każdym z trzech gatunków istnieje bardzo istotna różnica pomiędzy długością dziobów samców i samic. Samce mają z zwyczaj znacznie dłuższe dzioby niż samice choć w przypadku jednej samicy z gatunku Gentoo zaobserwowano jedną odstającą wartość ulokowaną daleko poza maksimum występującym u samców w tym gatunku.

Porównajmy zatem osobno samców i osobno samice dla poszczególnych gatunków.

pingwiny %>% ggGroupBoxplot(
  płeć, gatunek, długość.dzioba, F, T,
  title=paste("Rozkład długości dzioba dla poszczególnych płci",
              "w zależności od gatunku"))

W przypadku samców w każdym z gatunków występuje statystycznie różna wartość mediany długości dzioba. Dla samic różnice występują pomiędzy gatunkiem Adelie a pozostałymi gatunkami. Natomiast pomiędzy gatunkiem Chinstrap a Gentoo brak jest istotnych różnic.

Na koniec sprawdźmy jeszcze czy istnieją różnice w długości dzioba pomiędzy osobnikami tego samego gatunku, jednak zasiedlającymi różne wyspy. Najpierw jednak przekonajmy się który gatunek zasiedla którą wyspę.

pingwiny %>% group_by(gatunek, wyspa) %>% 
  summarise(n = n(), .groups="drop") %>% kable1

gatunek	wyspa	n
Adelie	Biscoe	44
Adelie	Dream	56
Adelie	Torgersen	52
Chinstrap	Dream	68
Gentoo	Biscoe	124

Okazuje się, że sensowne jest jedynie testowanie różnic w obrębie gatunku Adelie. Tylko bowiem ten gatunek występował na każdej z trzech wysp. Gatunek Chinstrap upodobał sobie wyspę Dream a gatunek Gentoo wyspę Biscoe. Mając jednak świadomość różnic pomiędzy samcami a samicami porównanie to wykonam osobno dla każdej z tych podgrup.

pingwiny %>% filter(gatunek == "Adelie") %>% 
  ggGroupBoxplot(
    płeć, wyspa, długość.dzioba, F, T, step.increase=.05, 
    title=paste("Rozkład długości dzioba gatunku Adelie dla samców oraz samic\n",
                "w zależności od zasiedlonej wyspy"))

Tym razem nie zaobserwowano żadnej istotnej różnicy. Tak więc długość dzioba gatunku Adelie nie zależy od wyspy którą zasiedlają dane osobniki.

3.4 Szerokość dzioba

Sprawdźmy zatem czy podobnie ma się rozkład szerokości dzioba.

pingwiny %>% sum_stats(szerokość.dzioba) %>% 
  mutate_if(is.numeric, signif, 3) %>% kable1

gatunek	n	nout	q1	min	mean	median	q3	max	sd	kurtosis	skewness	SW.stat	SW.p	SW.test
Adelie	152	1	17.5	15.5	18.3	18.4	19.0	21.5	1.210	-0.117	0.317	0.985	0.091	TRUE
Chinstrap	68	0	17.5	16.4	18.4	18.4	19.4	20.8	1.140	-0.960	0.007	0.973	0.142	TRUE
Gentoo	124	0	14.2	13.1	15.0	15.0	15.7	17.3	0.977	-0.628	0.317	0.977	0.031	FALSE

Jest bardzo zaskakujące, że podobnie jak w przypadku długości dzioba jego szerokość dla gatunku Gentoo również nie pochodzi z rozkładu normalnego.

pingwiny %>% ggBoxplot(
  gatunek, szerokość.dzioba, F, 
  title = "Rozkład szerokości dzioba dla poszczególnych gatunków pingwinów")

Gatunek ten również odznacza się dziobem o istotnie mniejszej szerokości. Warto jednak sprawdzić czy istnieją różnice ze względu na płeć.

pingwiny %>% group_by(gatunek, płeć) %>% sum_stats(szerokość.dzioba) %>% 
  mutate_if(is.numeric, signif, 3) %>% kable1

gatunek	płeć	n	nout	q1	min	mean	median	q3	max	sd	kurtosis	skewness	SW.stat	SW.p	SW.test
Adelie	samica	77	1	17.0	15.5	17.6	17.6	18.3	20.7	0.934	0.308	0.340	0.983	0.407	TRUE
Adelie	samiec	75	1	18.4	17.0	19.1	18.9	19.6	21.5	1.020	-0.276	0.463	0.965	0.036	FALSE
Chinstrap	samica	34	0	17.0	16.4	17.6	17.6	18.0	19.4	0.781	-0.841	0.293	0.959	0.230	TRUE
Chinstrap	samiec	34	0	18.8	17.5	19.3	19.3	19.8	20.8	0.761	-0.516	-0.104	0.983	0.863	TRUE
Gentoo	samica	61	0	13.8	13.1	14.2	14.3	14.6	15.5	0.530	-0.577	0.031	0.987	0.747	TRUE
Gentoo	samiec	63	0	15.2	14.1	15.7	15.7	16.1	17.3	0.735	-0.273	0.127	0.979	0.338	TRUE

Tu zgodność z rozkładem normalnym jest już znacznie lepsza. Co prawda na poziomie istotności \(\alpha=0.05\) należało by uznać, że szerokość dziobów samców gatunku Adelie nie pochodzi z rozkładu normalnego, to wartość p testu Shapiro-Wilka jest bardzo bliska naszemu poziomowi istotności. Zdecydowałem wiec aby w następnych analizach tej wielkości stosować mimo wszystko testy parametryczne.

pingwiny %>% ggGroupBoxplot(
  gatunek, płeć, szerokość.dzioba, T, ncol=3, 
  title="Rozkład szerokości dzioba dla trzech gatunków w zależności od płci")

Podobnie jak to było przypadku długości dzioba, również i w przypadku szerokości dzioba występują bardzo istotne różnice pomiędzy średnimi dla samców oraz samic.

Dokonajmy więc porównania dla samców i samic w zależności od gatunku.

pingwiny %>% ggGroupBoxplot(
  płeć, gatunek, szerokość.dzioba, T, T, 
  title=paste("Rozkład szerokości dzioba dla poszczególnych płci",
              "w zależności od gatunku"))

Mimo różnic występujących pomiędzy samcami i samicami w obrębie jednego gatunku, brak jest istonie statystycznych różnic średniej szerokości dzioba dla osobników tych samych płci pomiędzy gatunkiem Adelie a gatunkiem Chinstramp. Natomiast zarówno samce jak i samice gatunku Gentoo maja dzioby o zdecydowanie mniejszej szerokości niż osobniki tej samej płci pozostałych dwóch gatunków.

Podobnie też jak poprzednio sprawdźmy czy istnieją różnice w długości dzioba pomiędzy osobnikami gatunku Adelie zasiedlającymi różne wyspy.

pingwiny %>% filter(gatunek == "Adelie") %>% 
  ggGroupBoxplot(
    płeć, wyspa, szerokość.dzioba, T, T, step.increase=.05, 
    title=paste("Rozkład szerokości dzioba gatunku Adelie dla samców oraz samic\n",
                "w zależności od zasiedlonej wyspy"))

Tak jak w przypadku długości dzioba gatunku Adelie, szerokość dzioba nie jest zależna od tego którą wyspę zasiedla dany osobnik.

3.5 Długość płetw

Kolejnym parametrem do zbadania będzie długość płetw.

pingwiny %>% sum_stats(długość.płetw) %>% 
  mutate_if(is.numeric, signif, 3) %>% kable1

gatunek	n	nout	q1	min	mean	median	q3	max	sd	kurtosis	skewness	SW.stat	SW.p	SW.test
Adelie	152	2	186	172	190	190	195	210	6.52	0.26	0.086	0.993	0.677	TRUE
Chinstrap	68	0	191	178	196	196	201	212	7.13	-0.13	-0.009	0.989	0.811	TRUE
Gentoo	124	0	212	203	217	216	221	231	6.46	-0.62	0.388	0.963	0.002	FALSE

I znów, już po raz trzeci wartości dla gatunku Gentoo nie pochodzą z rozkładu normalnego!

pingwiny %>% ggBoxplot(
  gatunek, długość.płetw, F, 
  title = "Rozkład długości płetw dla poszczególnych gatunków pingwinów")

Chociaż gatunek Gentoo odznacza się najmniejszą szerokością dzioba, to nadrabia to długością płetw które są najdłuższe wśród trzech badanych gatunków. Z kolei różnice długości płetw pomiędzy gatunkami Adelie oraz Chinstrap nie są już tak duże, jednak nadal są one istotne statystycznie.

Mając na uwadze wcześniejsze obserwacje sprawdźmy także czy istnieją różnice pomiędzy długością płetw samców i samic.

pingwiny %>% group_by(gatunek, płeć) %>% sum_stats(długość.płetw) %>% 
  mutate_if(is.numeric, signif, 3) %>% kable1

gatunek	płeć	n	nout	q1	min	mean	median	q3	max	sd	kurtosis	skewness	SW.stat	SW.p	SW.test
Adelie	samica	77	3	185	172	188	187	191	202	5.64	0.206	-0.254	0.987	0.647	TRUE
Adelie	samiec	75	3	190	178	192	192	196	210	6.52	-0.006	0.067	0.984	0.446	TRUE
Chinstrap	samica	34	0	187	178	192	192	196	202	5.75	-0.475	-0.391	0.972	0.507	TRUE
Chinstrap	samiec	34	0	196	187	200	200	203	212	5.98	-0.645	0.172	0.975	0.620	TRUE
Gentoo	samica	61	0	210	203	213	213	216	222	3.85	-0.473	0.132	0.978	0.322	TRUE
Gentoo	samiec	63	0	217	208	221	221	225	231	5.64	-0.737	-0.045	0.962	0.052	TRUE

Tym razem mamy zgodność z rozkładem normalnym dla każdej z podgrupy danych. Bez oporów będę więc stosował testy parametryczne.

pingwiny %>% ggGroupBoxplot(
  gatunek, płeć, długość.płetw, T, ncol=3, 
  title="Rozkład długości płetw dla trzech gatunków w zależności od płci")

Kolejny raz możemy się przekonać o istotnych różnicach pomiędzy samcami a samicami.

pingwiny %>% ggGroupBoxplot(
  płeć, gatunek, długość.płetw, T, T, 
  title=paste("Rozkład długości płetw dla poszczególnych płci",
              "w zależności od gatunku"))

Ze względu na długość płetw można zaobserwować istotne różnice pomiędzy każdym z badanych gatunków niezależnie od płci osobników.

Sprawdźmy jeszcze czy może miejsce zasiedlenia osobników gatunku Adelie ma wpływ na długość płetw.

pingwiny %>% filter(gatunek == "Adelie") %>% 
  ggGroupBoxplot(
    płeć, wyspa, długość.płetw, T, T, step.increase=.05, 
    title=paste("Rozkład długości płetw gatunku Adelie dla samców oraz samic\n",
                "w zależności od zasiedlonej wyspy"))

I znów nie zaobserwowano żadnej istotnej różnicy.

3.6 Masa ciała

Ostatnim analizowanym parametrem fizycznym będzie masa ciała.

pingwiny %>% sum_stats(masa.ciała) %>% 
  mutate_if(is.numeric, signif, 3) %>% kable1

gatunek	n	nout	q1	min	mean	median	q3	max	sd	kurtosis	skewness	SW.stat	SW.p	SW.test
Adelie	152	0	3350	2850	3700	3700	4000	4780	457	-0.610	0.281	0.981	0.034	FALSE
Chinstrap	68	2	3490	2700	3730	3700	3950	4800	384	0.363	0.237	0.984	0.561	TRUE
Gentoo	124	0	4700	3950	5080	5020	5500	6300	502	-0.760	0.068	0.987	0.271	TRUE

Tym razem brak zgodności z rozkładem normalnym obserwujemy w przypadku gatunku Adelie choć również i tu wartość p testu Shapiro-Wilka jest tylko nieznacznie mniejsza od przyjętego przez nas poziomu istotności. Dlatego też zdecydowałem się na testy parametryczne.

pingwiny %>% ggBoxplot(
  gatunek, masa.ciała, T, 
  title = "Rozkład masy ciała dla poszczególnych gatunku pingwinów")

## Coefficient covariances computed by hccm()

Podobnie jak w przypadku długości płetw także i masa ciała osobników z gatunku Gentoo istotnie różni się od średniej masy osobników pozostałych dwóch gatunków, pomiędzy którymi nie ma w tym przypadku istotnych różnic.

Pora przejść do rozdziału również ze względu na płeć.

pingwiny %>% group_by(gatunek, płeć) %>% sum_stats(masa.ciała) %>% 
  mutate_if(is.numeric, signif, 3) %>% kable1

gatunek	płeć	n	nout	q1	min	mean	median	q3	max	sd	kurtosis	skewness	SW.stat	SW.p	SW.test
Adelie	samica	77	0	3180	2850	3370	3400	3550	3900	269	-0.897	-0.072	0.976	0.159	TRUE
Adelie	samiec	75	0	3800	3320	4040	4000	4290	4780	345	-0.711	0.159	0.983	0.416	TRUE
Chinstrap	samica	34	1	3360	2700	3530	3550	3690	4150	285	0.903	-0.585	0.963	0.306	TRUE
Chinstrap	samiec	34	1	3730	3250	3940	3950	4100	4800	362	-0.392	0.224	0.984	0.891	TRUE
Gentoo	samica	61	0	4450	3950	4670	4700	4880	5200	284	-0.449	-0.338	0.980	0.423	TRUE
Gentoo	samiec	63	0	5250	4750	5470	5500	5700	6300	321	-0.402	0.075	0.993	0.971	TRUE

I znów dla wszystkich podgrup mamy pełną zgodność z rozkładem normalnym.

pingwiny %>% ggGroupBoxplot(
  gatunek, płeć, masa.ciała, T, ncol=3, 
  title="Rozkład masu ciała dla trzech gatunków w zależności od płci")

Nie jest najmniejszym zaskoczeniem, że i w przypadku masy ciała dla samców i samic występują istotnie statystycznie różnice co ostatecznie przekonuje mnie o bardzo silnej determinacji płci każdej badanej cechy.

Rozdzielmy więc, po raz ostatni badaną cechę ze względu na płeć i zbadajmy różnice pomiędzy gatunkami.

pingwiny %>% ggGroupBoxplot(
  płeć, gatunek, masa.ciała, T, T, 
  title=paste("Rozkład masy ciała dla poszczególnych płci",
              "w zależności od gatunku"))

Jak widać w przypadku masy znów wybija się gatunek Gentoo przy braku różnic pomiędzy Adelie a Chinstrap.

Już wyłącznie formalnie sprawdźmy jeszcze czy masa ciała osobników gatunku Adelie będzie rożna ze względu na zasiedlaną wyspę.

pingwiny %>% filter(gatunek == "Adelie") %>% 
  ggGroupBoxplot(
    płeć, wyspa, masa.ciała, T, T, step.increase=.05, 
    title=paste("Rozkład szerokości dzioba gatunku Adelie na różnych wyspach",
                "z podziałem według płci"))

Tak jak można się było spodziewać nie widać tu żadnych istotnych różnic.

4 Korelacje

Analiza wartości poszczególnych zmiennych nasuwa wniosek, iż można się spodziewać dość silnej korelacji pomiędzy częścią z nich. Sprawdźmy zatem czy tak jest w istocie.

4.1 Wizualizacja

Do efektywnej wizualizacji wartości współczynników korelacji będę używał własnej funkcji ggCorr (bazującej na funkcji ggcorr z pakietu GGally) tworzącej wykres macierzy korelacji pomiędzy każdą parą zmiennych numerycznych znajdujących się w ramce danych będącej pierwszym argumentem tej funkcji oraz funkcje ggGroupCorr która rozpoznaje grupowanie zastosowane w ramce i jeżeli zostanie ono wykryte stworzy oddzielne wykresy dla każdej z podgrup danych. Funkcje te umożliwiają także swobodny wybór metody wyznaczenia współczynników korelacji. W tym zadaniu będę używał metody Pearsona.

ggCorr = function(data,
                  group = tibble(),
                  method = c("pairwise", "pearson"),
                  nbreaks = 10,
                  label_round = 2,
                  label_size = 4){
  group = group %>% mutate_if(is.factor, as.character)
  data %>% ungroup() %>%
    select_if(is.numeric) %>%
    ggcorr(method = method, nbreaks = nbreaks, label_round = label_round,
           hjust = 0.8, label = TRUE, label_size = label_size, color = "grey50") +
    labs(subtitle = paste(names(group), group, sep = ": ", collapse = ", "))
}

ggGroupCorr = function(data,
                   ncol = NULL,
                   nrow = NULL,
                   title = NULL,
                   caption = "M.F.",
                   common.legend = T,
                   legend = "bottom",
                   method = c("pairwise", "pearson"),
                   nbreaks = 10,
                   label_round = 2,
                   label_size = 4){
  p = data %>%
    group_map(~ggCorr(.x, .y, method, nbreaks, label_round, label_size)) %>%
    ggarrange(plotlist = ., common.legend = common.legend, legend=legend)
  if(!is.null(title)) p = p %>%
    annotate_figure(top = text_grob(title, size = 16, x=0, just="left"))
  if(!is.null(caption)) p = p %>% annotate_figure(
            bottom = text_grob(caption, size = 10, just="right", x = 1))
  p
}

4.2 Wszystkie obserwacje

Na początku zobaczmy czy istnieje korelacja pomiędzy zmiennymi bez dzielenia ich ze względy na jakiekolwiek podgrupy.

pingwiny %>% ggCorr(label_size=6)+
  labs(title = "Macierz korelacji zmiennych liczbowych", caption = "M.F")

Pierwsza, szybka analiza wyników wskazuje, że mamy bardzo silną, dodatnią korelację pomiędzy masą ciała oraz długością płetw. Również silna korelacja występuje pomiędzy długością i szerokością dzioba a długością płetw oraz masą ciała. Należy jednak pamiętać, o znaczących różnicach pomiędzy gatunkami oraz pomiędzy osobnikami różnych płci. Zobaczmy więc jak będą wyglądały współczynniki korelacji kiedy zgrupujemy dane ze względu na płeć.

4.3 Grupowanie wg płci

pingwiny %>%
  group_by(płeć) %>%
  ggGroupCorr(title = "Macierz korelacji dla poszczególnych płci")

Co może być zaskakujące, to że w wielu przypadkach wyniki wskazują na znacznie silniejszą korelację. Szczególnie w przypadku korelacji pomiędzy szerokością dzioba a długością płetw i masą ciała. Widocznie wzrosła też korelacja pomiędzy szerokością i długością dzioba.

4.4 Grupowanie wg gatunku

Zobaczmy w takim razie jakie otrzymamy wartości współczynników korelacji jeżeli będziemy grupować dane nie względem płci a względem gatunku. Tu przecież również występowały dość istotne różnice.

pingwiny %>%
  group_by(gatunek) %>%
  ggGroupCorr(title = "Macierz korelacji dla poszczególnych gatunków")

No i mamy kolejne spore zaskoczenie. Przy grupowaniu według płci oraz kiedy nie stosowałem żadnego grupowania, korelacja była zarówno dodatnia jak i ujemna. Tym czasem kiedy pogrupowałem dane ze względu na gatunek wszystkie współczynniki korelacji są dodatnie! Ponadto można zauważyć, że dla gatunku Gentoo pomiędzy każdą parą zmiennych istnieje dość silna korelacja.

4.5 Grupowanie wg gatunku oraz płci

Nie pozostaje nam więc nic innego jak tylko grupowanie zarówno względem płci jak i gatunku.

pingwiny %>%
  group_by(płeć, gatunek) %>%
  ggGroupCorr(title = "Macierz korelacji dla poszczególnych gatunków i płci")

Takie grupowanie danych spowodowało, że w wielu przypadkach współczynniki korelacji spadły do wartości średnich lub nawet bardzo niskich. W żadnym przypadku nie mamy już tak sielnej korelacji jak początkowe 0.87.

W zasadzie jedyna silna korelacja występuje dla samców Chinstrap pomiędzy długością płetw a masą ciała.

5 Regresja

Informacje z poprzedniego rozdziału prowadzą do wniosku, że tylko dla nielicznych par, odpowiednio grupowanych zmiennych będzie można znaleźć istotne funkcje regresji. Przekonajmy się zatem o słuszności takich oczekiwań. Biorąc jednak pod uwagę wszystkie wcześniej dostrzeżone zależności od razu będę grupował zmienne zarówno ze względu na gatunek jak i na płeć.

Na początku jednak przygotuję odpowiednią funkcję wizualizacyjną.

5.1 Wizualizacja

Moja funkcja ggSmooth przygotowuje wykres rozrzutu wraz z dopasowaną funkcją regresji liniowej. Funkcja ta również rozróżnia zastosowane na wejściowej ramce danych grupowanie. Ponadto prezentuje na wykresie równanie regresji wraz z wartościami testu F co ułatwia szybką ocenę jakości dopasowania.

ggSmooth = function(data, 
                    x, y, 
                    v.step = 1, 
                    y.lim = NULL,
                    eq.size = 3){
  x = enquo(x)
  y = enquo(y)
  
  group = data %>% group_vars()
  color = NULL
  if(length(group)>=1) color = sym(group[1])
  facet = NULL
  if(length(group)>=2) facet = sym(group[2])
  
  p = data %>%
    ggplot(aes(!!x, !!y, color= !!color))+
    stat_poly_eq(aes(
      label =  paste(
        after_stat(eq.label), "*\", \"*",
        after_stat(rr.label), "*\", \"*",
        after_stat(f.value.label), "*\", \"*",
        after_stat(p.value.label), "*\".\"",
        sep = "")),
      formula = y~x, geom = "label_npc", vstep = v.step, size = eq.size)+
    geom_point()+
    geom_smooth(formula = y~x, method = "lm") +
    theme(legend.position="bottom")+
    scale_colour_manual(values = c("darkorange","purple","cyan4")) +
    scale_fill_manual(values = c("darkorange","purple","cyan4"))
  
  if(!is.null(facet)) p = p + facet_wrap(vars(!!facet))
  if(!is.null(y.lim)) p = p + ylim(y.lim)
  p
}

5.2 Długość dzioba w funkcji jego szerokości

Długość dzioba oraz jego szerokość była dość słabo skorelowana. Sprawdźmy zatem czy pomiędzy tymi wielkościami można zaleźć jakąś zależność liniową.

pingwiny %>% group_by(gatunek, płeć) %>% 
  ggSmooth(szerokość.dzioba, długość.dzioba, .07, c(NA,70), 3.8) %>% 
  ggpar(title = "Zależność pomiędzy szerokością a długością dzioba", caption = "M.F.")

Można zauważyć, że tylko w przypadku samców z gatunków Chinstrap oraz Gentoo model regresji jest istotny. W przypadku samic dotyczy to wyłącznie gatunku Gentoo. Jednak współczynniki determinacji są bardzo niskie i nigdzie nie przekraczają 0.2. Trudno także wyciągnąć wnioski o tym, że mamy tu do czynienia z regresją nieliniową. Widać natomiast jak duże różnice występując pomiędzy gatunkami przy uwzględnieniu płci. Praktycznie żaden z obszarów rozrzutu nie nachodzi na innym obszar.

5.3 Długość płetw w funkcji długości dzioba

pingwiny %>% group_by(gatunek, płeć) %>% 
  ggSmooth(długość.dzioba, długość.płetw, .07, c(NA,250), 3.8) %>% 
  ggpar(title = "Zależność pomiędzy długością płetw a długością dzioba", caption = "M.F.")

Gdyby polegać tylko na współczynniku korelacji wyliczonym dla wszystkich danych można by oczekiwać dość dobrego dopasowania funkcji regresji pomiędzy tymi zmiennymi. Tymczasem po podzieleniu danych na podgrupy model okazał się istotny tylko w przypadku samców z gatunku Adelie oraz Gentoo, choć również i i tutaj współczynnik determinacji należy uznać za raczej bardzo niski. Również i w tym przypadku odpowiednie obszary rozrzutu praktycznie nie nachodzą na siebie.

5.4 Długość płetw w funkcji masy ciała

pingwiny %>% group_by(gatunek, płeć) %>% 
  ggSmooth(masa.ciała, długość.płetw, .07, c(NA, 250), 3.8) %>% 
  ggpar(title = "Zależność pomiędzy długością płetw a masą ciała", caption = "M.F.")

Jak można się było spodziewać masa ciała oraz długość płetw rzeczywiście są wielkościami gdzie występuje (poza samicami z gatunku Chinstrap) istotna zależność liniowa. Chociaż nawet tutaj współczynnik determinacji nigdzie nie przekroczył wartości 0.5.

Wnioski

Zrealizowany projekt pozwolił mi wyciągnąć kilka ciekawych wniosków oraz zastosować w praktyce zdobytą wiedzę. Szczególnie istotna była dla mnie, stworzona właśnie do tego projektu funkcja naKnn która z pewnością będzie mi pomocna w niejednym przyszłym projekcie.

Sama analiza danych pokazała mi zaś, jak bardzo istotne jest dokładne rozważenie wszystkich zależności kryjących się w badanych zmiennych. W przypadku sympatycznych pingwinów okazało się, że nie wystarczy brać pod uwagę jedynie samego gatunku. Różnice bowiem które występują pomiędzy samcami i samicami są bardzo istotne i występowały w przypadku każdego badanego parametru. Pochopne pominięcie wpływu takich znaczących zmiennych może prowadzić do błędny wniosków, co szczególnie widoczne było w przypadku analizowania współczynników korelacji.