Przygotowanie danych

Dane pobrano z bazy danych NOAA. Zawierają one dane meteorologiczne dla wybranej stacji. W tym przypadku została wybrana stacja w Zakopanem. Anova będzie stosowana dla grup danych temperatury w zależności od roku (2017, 2018, 2019, 2020), w którym zostały wykonane pomiary.

library(tidyverse)
library(worldmet)
library(lubridate)
library(openair)
library(e1071)
library(lawstat)
library(ggplot2)
library(DT)

#dane1 <- read.csv(file = "C:/Users/domka/Desktop/AGH/II semestr/PNOZ/projekt1_sz/Projekt_1_PNOZ/isd-history.csv",header = T,sep = ",",dec = ".")
#dane2 <- read.csv(file = "C:/Users/domka/Desktop/AGH/II semestr/PNOZ/projekt1_sz/Projekt_1_PNOZ/isd-inventory.csv",header = T,sep = ",",dec = ".")

#dane_PL <- dane1 %>% 
#  filter(.,CTRY == "PL")

#dane2 %>% 
#  filter(.,USAF == "125100")

#dane_NOAA <- importNOAA(code = "125100-99999",year = c(2017, 2018, 2019, 2020),hourly = T, quiet = FALSE)
dane_NOAA <- readRDS(file="noaa.rds")

dane_NOAA <- cutData(dane_NOAA,type="season")
dane_NOAA <- cutData(dane_NOAA,type="month")
dane_NOAA <- cutData(dane_NOAA,type="year")

dane_NOAA %>% select(air_temp, season, date, month, year) %>%
  mutate(day = day(.$date)) %>%
  group_by(day,month,year,season) %>% 
  summarise(mean_temp = round(mean(air_temp, na.rm=TRUE),2)) %>% 
  filter(season == "winter (DJF)")-> dane_NOAA_daily

dane_NOAA_daily_by_year <- dane_NOAA_daily[c(3,5)]

dane_NOAA_daily_by_year_2017 <- dane_NOAA_daily_by_year %>% filter(year == "2017")
dane_NOAA_daily_by_year_2018 <- dane_NOAA_daily_by_year %>% filter(year == "2018")
dane_NOAA_daily_by_year_2019 <- dane_NOAA_daily_by_year %>% filter(year == "2019")
dane_NOAA_daily_by_year_2020 <- dane_NOAA_daily_by_year %>% filter(year == "2020")

Ze wszystkich danych wybrano 30-elementową próbkę średniodniowej temperatury w sezonie zimowym.

#temperature_2017 <- sample(dane_NOAA_daily_by_year_2017$mean_temp, 30)
#temperature_2018 <- sample(dane_NOAA_daily_by_year_2018$mean_temp, 30)
#temperature_2019 <- sample(dane_NOAA_daily_by_year_2019$mean_temp, 30)
#temperature_2020 <- sample(dane_NOAA_daily_by_year_2020$mean_temp, 30)

#temperature <- append(temperature_2017,temperature_2018,length(temperature_2017))
#temperature <- append(temperature,temperature_2019,length(temperature))
#temperature <- append(temperature,temperature_2020,length(temperature))
temperature <- readRDS(file="temp.rds")

rok = c(rep("2017",30), rep("2018",30), rep("2019",30), rep("2020",30))
df <- data.frame(rok,temperature)

Wstępny przegląd danych za pomocą histogramu

hist(df$temperature,col="blue")

Podstawowe statystyki

stats <- data.frame(round(mean(df$temperature),3), round(median(df$temperature),3), round(var(df$temperature),3), round(sd(df$temperature),3), round(kurtosis(df$temperature),3), round(skewness(df$temperature),3))
colnames(stats) <- c("Srednia","Mediana","Wariancja","Odchylenie std","Kurtoza","Skosnosc")
datatable(stats)

Sprawdzanie założeń

Test Shapiro-Wilka

P-value jest większe niż 0.05, zatem przyjmuje się hipotezę zerową, czyli wykres jest podobny do rozkładu normalnego.

shapiro.test(df$temperature)
## 
##  Shapiro-Wilk normality test
## 
## data:  df$temperature
## W = 0.98891, p-value = 0.4417

Testowanie homogeniczności wariancji

W obu testach p-value jest większe od 0.05, zatem przyjmuje się hipotezę zerową, czyli jest zachowana homogeniczność wariancji.

levene.test(df$temperature, df$rok)
## 
##  Modified robust Brown-Forsythe Levene-type test based on the absolute
##  deviations from the median
## 
## data:  df$temperature
## Test Statistic = 0.28617, p-value = 0.8353
bartlett.test(df$temperature, df$rok)
## 
##  Bartlett test of homogeneity of variances
## 
## data:  df$temperature and df$rok
## Bartlett's K-squared = 2.1031, df = 3, p-value = 0.5513

Anova

Tworzenie Anovy

W Anovie p-value jest mniejsze od 0.05, czyli przyjmuje się hipotezę alternatywną, zatem przynajmniej jedna z grup różni się od pozostałych.

res.aov = aov(temperature~rok,data = df)
summary(res.aov)
##              Df Sum Sq Mean Sq F value  Pr(>F)   
## rok           3  319.2  106.40   5.446 0.00154 **
## Residuals   116 2266.5   19.54                   
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Test Tukey’a

Patrząc na p adj, dwie pary różnią się między sobą - 2018/2017 i 2020/2018, ponieważ p adj < 0.05

TukeyHSD(res.aov)
##   Tukey multiple comparisons of means
##     95% family-wise confidence level
## 
## Fit: aov(formula = temperature ~ rok, data = df)
## 
## $rok
##                 diff        lwr         upr     p adj
## 2018-2017 -3.0233333 -5.9983541 -0.04831257 0.0448735
## 2019-2017 -0.7563333 -3.7313541  2.21868743 0.9108539
## 2020-2017  1.4943333 -1.4806874  4.46935410 0.5589046
## 2019-2018  2.2670000 -0.7080208  5.24202077 0.1991003
## 2020-2018  4.5176667  1.5426459  7.49268743 0.0007450
## 2020-2019  2.2506667 -0.7243541  5.22568743 0.2045811

Boxplot

Z wykresu widać, że wartość średnia grup lat 2017 i 2020 nie nachodzi na pudełko roku 2018. Dodatkowo wartość średnia grupy z roku 2018 też nie nachodzi na pudełko 2020, co zgadza się z największą różnicą pomiędzy tymi grupami.

ggplot(df,aes(x=rok,y=temperature,fill=rok))+
  geom_boxplot()

Jitter

W 2018 roku widać, że temperatura była średnio niższa niż w innych latach i że jedna wartość przekroczyła temperaturę poniżej -20 stopni.

ggplot(df,aes(x=rok,y=temperature,fill=rok,color=rok))+
  geom_jitter()