Wykład 12 - Statystyka matematyczna (kierunek matematyka specjalność matematyka finansowa)

Praca w środowisku R

¹ Zakład Bioinformatyki, Instytut Informatyki, Uniwersytet w Białymstoku

^✉ Correspondence: Jarosław Kotowicz <j.kotowicz@uwb.edu.pl>

1 Przygotowanie do pracy środowiska pracy

1.1 Instalacja oprogramowania

1. R strona domowa www.r-project.org

• wybieramy download
• wybieramy mirror np. serwer niemiecki
• ściągamy plik instalacyjny

2. Rstudio (GUI) dla R RStudio

• należy pobrać system opensource

3. Przeglądarka plików pdf natywna dla RStudio (pod Windows) sumatra (opensource)

• pobieramy ze strony Sumatra

4. Instalujemy w kolejności R, sumatra, RStudio

1.2 Przygotowanie do pracy

1. W katalagu dokumenty tworzymy folder (nazwa zgodnie z zasadą Inf_xxxxx, gdzie xxxxx to numer albumu).
2. Uruchamiamy RSudio
3. Zachowujemy porządek. Tworzymy nowy projekt

• Files -> New Project… -> Existing Directory -> Browse… i wybieramy folder, który utworzyliśmy oraz zaznaczamy Open in new session.
  Uwaga: Zawsze na koniec pracy/zajęć należy zamknąć projekt. Files -> Close Project
  Uwaga: Zawsze na początku pracy/zajęć należy otworzyć projekt. Files -> Recent Projects i wybrać swój (będzie miał on taką nazwę, jak nasz katalog).

4. Dostosowujem RSudio do następujących wymagań:

• Tools -> Project Options ->
  1. Code Editing -> Text encoding wybór UTF-8
  2. Sweave -> PDF Generaion wybór knitr
     Uwaga: Jeżeli chcemy korzystać z TeX (LaTeX), to musimy mieć zaistalowny system/program np. MikTeX.
• Tools -> Global Options ->
  1. Appearance wybór wyglądu RSudion np. ambiance (ciemny motyw).

---

5. W domu instalujemy konieczne biblioteki R nazywane package

• dolne prawe okienko - okno Packages -> Install wpisujemy tidyverse

6. Otwieramy pliki

• R Script: Files -> New File -> R Script
• R Notebook: Files -> New File -> R Notebook
• R Markdown: Files -> New File -> R Markdown i wybieramy rodzaj dokumentu np. Document lub Presentation.
  

---

7. Będziemy zawsze pracować albo w pliku skryptowym R albo w R Notebook.
8. Informacje dla pliku skryptowego

• wywołanie aktualnej linii w konsoli CTRL+ENTER
• odwołanie się do funkcji z danej biblioteki package::function np. stats::filter
• podczytanie biblioteki library(nazwa bibilioteki)
• znak komentarza hashtag tzn. #

---

9. Import danych

• korzystamy z danych dostępnych na stronie Jareda P. Landera
• plik csv (textowe) import: File -> Import Dataset -> From Text (readr) (w polu File/URL kopiujemy adres https://jaredlander.com/data/ i dopisujemy nazwę pliku np. acs_ny.csv) -> Update
• wykorzystujemy opcje Import Options aby wybrać np. znak rozdzielający dane średnik to semicolon itd.
• możemy zmienić też typy kolumn wybierając dla kolumny z menu rozwijalnego wartość Uwaga: W przypadku typu factor (zmienna czynnikowa) należy wskazać wszystkie poziomy (wartości).
• kopiujemy z okienka Code Preview do pliku skryptowego polecenie poczytując wcześniej bibliotekę readr, ja tutaj podczytuję tidyverse.
• braki danych NA (skrót od not available), nieliczby NaN (not a number)
• typ numeric w R identyczny jak double w C++ i co więcej używane też jest równoważnie numeric i double

library(tidyverse)

Registered S3 method overwritten by 'dplyr':
  method           from
  print.rowwise_df     
Registered S3 methods overwritten by 'dbplyr':
  method         from
  print.tbl_lazy     
  print.tbl_sql      
-- Attaching packages --------------------------------------- tidyverse 1.3.0 --
<U+221A> ggplot2 3.3.0     <U+221A> purrr   0.3.4
<U+221A> tibble  3.0.1     <U+221A> dplyr   0.8.5
<U+221A> tidyr   1.0.2     <U+221A> stringr 1.4.0
<U+221A> readr   1.3.1     <U+221A> forcats 0.5.0
-- Conflicts ------------------------------------------ tidyverse_conflicts() --
x dplyr::filter() masks stats::filter()
x dplyr::lag()    masks stats::lag()

Przykład (skopiowany kod i następnie uporządkowany)

acsNY <- read_csv("https://jaredlander.com/data/acs_ny.csv", 
                  col_types = cols(ElectricBill = col_integer(), 
                                   FamilyIncome = col_integer(), 
                                   FamilyType = col_factor(levels = c("Married", "Female Head", "Male Head")), 
                                   HouseCosts = col_integer(), 
                                   Insurance = col_integer(), 
                                   NumBedrooms = col_integer(), 
                                   NumChildren = col_integer(), 
                                   NumPeople = col_integer(), 
                                   NumRooms = col_integer(), 
                                   NumVehicles = col_integer(), 
                                   NumWorkers = col_integer()))

---

10. Przetwarzanie potokowe tzw. pipe (biblioteka magrittr) Uwaga: Użycie biblioteki dplyr powoduje możliwość korzystania z przetwarzania potokowego. Używając biblioteki tidyverse podczytuję również dplyr. co wysyłam %>% gdzie wysyłam Zasada przetwarzania potkowego f(a) jest to a %>% f

f(a,b) jest to a %>% f(b)

f(a,b) jest to b %>% f(a,.) f(a,b) Uwaga: Kropka wskazuje miejsca wstawienia zmiennej.

g(f(a)) jest to a %>% f %>% g

---

11. Reguła wcięć poprawny zapis zgodny z regułą wcięć to

acsNY %>% 
  summary()

    Acres            FamilyIncome           FamilyType     NumBedrooms     NumChildren     
 Length:22745       Min.   :     50   Married    :18326   Min.   :0.000   Min.   : 0.0000  
 Class :character   1st Qu.:  52540   Female Head: 3266   1st Qu.:3.000   1st Qu.: 0.0000  
 Mode  :character   Median :  87000   Male Head  : 1153   Median :3.000   Median : 0.0000  
                    Mean   : 110281                       Mean   :3.385   Mean   : 0.9012  
                    3rd Qu.: 133800                       3rd Qu.:4.000   3rd Qu.: 2.0000  
                    Max.   :1605000                       Max.   :8.000   Max.   :12.0000  
   NumPeople        NumRooms        NumUnits          NumVehicles      NumWorkers   
 Min.   : 2.00   Min.   : 1.000   Length:22745       Min.   :0.000   Min.   :0.000  
 1st Qu.: 2.00   1st Qu.: 6.000   Class :character   1st Qu.:2.000   1st Qu.:1.000  
 Median : 3.00   Median : 7.000   Mode  :character   Median :2.000   Median :2.000  
 Mean   : 3.39   Mean   : 7.175                      Mean   :2.113   Mean   :1.745  
 3rd Qu.: 4.00   3rd Qu.: 8.000                      3rd Qu.:3.000   3rd Qu.:2.000  
 Max.   :18.00   Max.   :21.000                      Max.   :6.000   Max.   :3.000  
   OwnRent           YearBuilt           HouseCosts    ElectricBill  FoodStamp        
 Length:22745       Length:22745       Min.   :   4   Min.   :  1   Length:22745      
 Class :character   Class :character   1st Qu.: 650   1st Qu.:100   Class :character  
 Mode  :character   Mode  :character   Median :1200   Median :150   Mode  :character  
                                       Mean   :1480   Mean   :175                     
                                       3rd Qu.:2000   3rd Qu.:220                     
                                       Max.   :7090   Max.   :580                     
 HeatingFuel          Insurance        Language        
 Length:22745       Min.   :   0.0   Length:22745      
 Class :character   1st Qu.: 400.0   Class :character  
 Mode  :character   Median : 720.0   Mode  :character  
                    Mean   : 960.9                     
                    3rd Qu.:1200.0                     
                    Max.   :6600.0                     

---

2 Funkcje wywołujące rozkłady w R

Budowa funkcji prefix NazwaRozkładu

Rodzaje prefiksów

1. d - funkcja gęstości w przypadku rozkładu ciągłego lub funkcja prawdopodobieństwa w przypadku rozkłądu dyskretnego
2. p - dystrybuanta rozkładu
3. q - kwantyle rozkładu
4. r - liczba (liczby) pseudolosowe z rozkładu

Nazwy rozkładów

1. norm - noramlny,

2. gamma

3. beta

4. chisq

5. itd

6. Wyznaczamy wartość gęstości rozkładu normalnego standardowego dla argumentu 0.

dnorm(0)

[1] 0.3989423

2. Wyznaczamy wartość gęstości rozkładu N(2,5) dla argumentu 0.

dnorm(0, mean = 2, sd = 5)

[1] 0.07365403

3. Wyznaczamy wartość dystrybuanty rozkładu normalnego standardowego dla argumentu od 0 do 1 z krokiem .2.

pnorm(seq(0, 1, by = .2))

[1] 0.5000000 0.5792597 0.6554217 0.7257469 0.7881446 0.8413447

2.1 Liczby pseudolosowe z zadanego rozkładu prawdopodobieństwa

set.seed(20200512)
x.norm32 <- rnorm(1000, mean = 3, sd = 2)
set.seed(20200512)
x.norm01 <- rnorm(1000, mean = 0, sd = 1)

2.2 Dystrybuanta

pnorm(0)

[1] 0.5

2.3 Gęstość

dnorm(0)

[1] 0.3989423

2.4 Kwantyle

qnorm(c(.25, .5, .75))

[1] -0.6744898  0.0000000  0.6744898

---

3 Statystyka opisowa

3.1 Szereg rozdzielczy przedziałowy (prosty, skumulowany)

3.1.1 Zmienna jednowymiarowa

Registered S3 methods overwritten by 'htmltools':
  method               from         
  print.html           tools:rstudio
  print.shiny.tag      tools:rstudio
  print.shiny.tag.list tools:rstudio
Registered S3 method overwritten by 'GGally':
  method from   
  +.gg   ggplot2
Registered S3 method overwritten by 'gdata':
  method         from     
  reorder.factor DescTools
Registered S3 method overwritten by 'htmlwidgets':
  method           from         
  print.htmlwidget tools:rstudio

table(cut(x.norm01, seq(min(x.norm01), max(x.norm01)+.2, by = 1), include.lowest = TRUE))

 [-3.17,-2.17]  (-2.17,-1.17] (-1.17,-0.169] (-0.169,0.831]   (0.831,1.83]    (1.83,2.83] 
            15            109            293            381            175             24 
   (2.83,3.83] 
             3 

sjmisc::frq(sample(c("x", "y", "z"), size = 1000, replace = TRUE, prob = c(.3, .5, .2)))

x <character>
# total N=1000  valid N=1000  mean=1.88  sd=0.68

Value |   N | Raw % | Valid % | Cum. %
--------------------------------------
x     | 297 | 29.70 |   29.70 |  29.70
y     | 525 | 52.50 |   52.50 |  82.20
z     | 178 | 17.80 |   17.80 | 100.00
<NA>  |   0 |  0.00 |    <NA> |   <NA>

3.1.2 Zmienna wielowymiarowa (ramka danych data.frame)

df01 <- data.frame(zm1 = sample(c("a", "b"), 100, replace = TRUE),
                   zm2 = sample(c("c", "d"), 100, replace = TRUE))

sjmisc::frq(df01)

zm1 <character>
# total N=100  valid N=100  mean=1.38  sd=0.49

Value |  N | Raw % | Valid % | Cum. %
-------------------------------------
a     | 62 |    62 |      62 |     62
b     | 38 |    38 |      38 |    100
<NA>  |  0 |     0 |    <NA> |   <NA>


zm2 <character>
# total N=100  valid N=100  mean=1.55  sd=0.50

Value |  N | Raw % | Valid % | Cum. %
-------------------------------------
c     | 45 | 45.00 |   45.00 |     45
d     | 55 | 55.00 |   55.00 |    100
<NA>  |  0 |  0.00 |    <NA> |   <NA>

3.2 Tablica kontyngencji

(tabela <- table(df01$zm1, df01$zm2))

   
     c  d
  a 24 38
  b 21 17

3.3 Średnia

mean(x.norm01)

[1] -0.009051554

mean(x.norm32)

[1] 2.981897

3.4 Wariancja

var(x.norm01)

[1] 0.9785055

3.5 Odchylenie standradowe

sd(x.norm32)

[1] 1.978389

3.6 Tzw. fivenum (minimum, kwartyle, maksimum)

fivenum(x.norm01)

[1] -3.1686257 -0.6939557  0.0476431  0.6494921  3.6518932

3.7 Kwantyle

quantile(x.norm01, seq(0,1, by = .1))

        0%        10%        20%        30%        40%        50%        60%        70% 
-3.1686257 -1.2769306 -0.8797193 -0.5415714 -0.2266632  0.0476431  0.2510401  0.4988340 
       80%        90%       100% 
 0.8378308  1.2253506  3.6518932 

3.8 Zakres zmienności

range(x.norm32)

[1] -3.337251 10.303786

3.9 Mediana

median(x.norm32)

[1] 3.095286

3.10 Podsumowanie i summary vs. fivenum

summary(x.norm32)

   Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
 -3.337   1.612   3.095   2.982   4.293  10.304 

fivenum(x.norm32)

[1] -3.337251  1.612089  3.095286  4.298984 10.303786

3.11 Podstawowe miary statystyki opisowej

psych::describe(x.norm32, quant = c(.25, .75), IQR = TRUE)

Hmisc::describe(x.norm32)

x.norm32 
       n  missing distinct     Info     Mean      Gmd      .05      .10      .25      .50 
    1000        0     1000        1    2.982    2.229  -0.3510   0.4461   1.6121   3.0953 
     .75      .90      .95 
  4.2929   5.4507   6.0952 

lowest : -3.337251 -3.304396 -2.868932 -2.750761 -1.974559
highest:  8.327035  8.470047  8.957606  9.994734 10.303786

3.12 Rozstęp międzyćwiartkowy (czyli podwojone odchylenie ćwiartkowe)

IQR(x.norm32)

[1] 2.68075

3.13 Skośność (wywołanie na dwa sposoby i trzech bibliotek)

e1071::skewness(x.norm01)

[1] -0.03115604

library(e1071)
skewness(x.norm32)

[1] -0.03115604

detach(package:e1071)

agricolae::skewness(x.norm01)

[1] -0.03124973

psych::skew(x.norm01)

[1] -0.03115604

3.14 Moda

RVAideMemoire::mod(x.norm32)

[1] 3.371455

3.15 Kurtoza

e1071::kurtosis(x.norm01)

[1] 0.1129802

agricolae::kurtosis(x.norm01)

[1] 0.1258378

psych::kurtosi(x.norm01)

[1] 0.1129802

3.16 Momenty

3.16.1 Z domyślnymi parametrami

e1071::moment(x.norm32)

[1] 2.981897

3.16.2 Odpowiedniego rzędu

e1071::moment(x.norm32, order = 2)

[1] 12.80182

3.16.3 Centralne

e1071::moment(x.norm32, order = 2, center = TRUE)

[1] 3.910108

3.16.4 Moment centralny rzędu 2 vs. wariancja (var)

e1071::moment(x.norm32, order = 2, center = TRUE)

[1] 3.910108

var(x.norm32)

[1] 3.914022

var(x.norm32)/length(x.norm32)*(length(x.norm32)-1)

[1] 3.910108

3.17 Korelacja i kowariancja (macierze korelacji i kowariancji)

df <- data.frame(norm01 = x.norm01, norm32 = x.norm32, ga = rgamma(1000, 2, 1)) %>% as_tibble()

3.17.1 Kowariancja

cov(x.norm01, df$ga)

[1] -0.02882743

cov(df)

            norm01      norm32          ga
norm01  0.97850552  1.95701104 -0.02882743
norm32  1.95701104  3.91402208 -0.05765486
ga     -0.02882743 -0.05765486  2.19773905

3.17.2 Korelacja

cor(df)

            norm01      norm32          ga
norm01  1.00000000  1.00000000 -0.01965786
norm32  1.00000000  1.00000000 -0.01965786
ga     -0.01965786 -0.01965786  1.00000000

Hmisc::rcorr(as.matrix(df), type = "pearson")

       norm01 norm32    ga
norm01   1.00   1.00 -0.02
norm32   1.00   1.00 -0.02
ga      -0.02  -0.02  1.00

n= 1000 


P
       norm01 norm32 ga    
norm01        0.0000 0.5347
norm32 0.0000        0.5347
ga     0.5347 0.5347       

Hmisc::rcorr(as.matrix(df), type = "spearman")

       norm01 norm32    ga
norm01   1.00   1.00 -0.01
norm32   1.00   1.00 -0.01
ga      -0.01  -0.01  1.00

n= 1000 


P
       norm01 norm32 ga    
norm01        0.0000 0.8288
norm32 0.0000        0.8288
ga     0.8288 0.8288       

---

4 Statystyki graficzne z ggplot2

4.1 Dystrybuanta empiryczna

df %>%
ggplot(aes(x = norm01)) +
  stat_ecdf()

4.2 Gęstość

df %>%
  ggplot(aes(x.norm01)) +
  geom_density(kernel = "gaussian")

4.3 Histogram

df %>%
  ggplot(aes(x.norm01)) +
  geom_histogram()

df %>%
  ggplot(aes(x.norm01)) +
  geom_freqpoly()

4.4 Wykres kwantyli

df %>%
  ggplot() +
  geom_qq(aes(sample= x.norm01))

GGally::ggpairs(df)

---

5 Estymacja punktowa metodą największej wiarogodnosci (biblioteka MASS)

MASS::fitdistr(x.norm01, "normal")

       mean            sd     
  -0.009051554    0.988699659 
 ( 0.031265428) ( 0.022107996)

MASS::fitdistr(df$ga, "gamma")

     shape         rate   
  1.81814044   0.93017351 
 (0.07504241) (0.04415700)

---

6 Estymacja przedziałowa

6.1 Różne funkcje wyznaczające przedział ufności dla średniej

Rmisc::CI(x.norm01)

       upper         mean        lower 
 0.052332592 -0.009051554 -0.070435701 

DescTools::MeanCI(x.norm01)

        mean       lwr.ci       upr.ci 
-0.009051554 -0.070435701  0.052332592 

t.test(x.norm01)

    One Sample t-test

data:  x.norm01
t = -0.28936, df = 999, p-value = 0.7724
alternative hypothesis: true mean is not equal to 0
95 percent confidence interval:
 -0.07043570  0.05233259
sample estimates:
   mean of x 
-0.009051554 

gmodels::ci(x.norm01)

No class or unkown class.  Using default calcuation.

    Estimate     CI lower     CI upper   Std. Error 
-0.009051554 -0.070435701  0.052332592  0.031281073 

6.2 Przedział ufności dla frakcji

set.seed(20200512)
x.bin <-rbinom(1000, size = 1, prob = .3)
gmodels::ci.binom(x.bin)

     Estimate  CI lower  CI upper Std. Error
[1,]    0.294 0.2659037 0.3233147 0.01440708

prop.test(sum(x.bin), length(x.bin), .3)

    1-sample proportions test with continuity correction

data:  sum(x.bin) out of length(x.bin), null probability 0.3
X-squared = 0.14405, df = 1, p-value = 0.7043
alternative hypothesis: true p is not equal to 0.3
95 percent confidence interval:
 0.2661099 0.3234946
sample estimates:
    p 
0.294 

6.3 Budowa własnej funkcji wyznaczającej przedział ufności dla wariancji

CI_var_N <- function(dane, alpha = .05) {
  konieclewy <- (length(dane) -1) * var(dane) / qchisq(1- alpha/2, df = length(dane) - 1)
  koniecprawy <- (length(dane) -1) * var(dane) / qchisq(alpha/2, df = length(dane) - 1)
  cat(paste("Przedział ufności dla wariancji wynosi:\n [", konieclewy, ",", koniecprawy, "]\n"))
}

CI_var_N(x.norm01)

Przedział ufności dla wariancji wynosi:
 [ 0.898060293361267 , 1.07032294630049 ]

---

7 Testowanie hipotez istotoności

7.1 Testy parametryczne

7.1.1 Test dla średniej

t.test(x.norm01)

    One Sample t-test

data:  x.norm01
t = -0.28936, df = 999, p-value = 0.7724
alternative hypothesis: true mean is not equal to 0
95 percent confidence interval:
 -0.07043570  0.05233259
sample estimates:
   mean of x 
-0.009051554 

Interpretacja wyniku!

t.test(x.norm32)

    One Sample t-test

data:  x.norm32
t = 47.663, df = 999, p-value < 2.2e-16
alternative hypothesis: true mean is not equal to 0
95 percent confidence interval:
 2.859129 3.104665
sample estimates:
mean of x 
 2.981897 

Interpretacja wyniku!

t.test(x.norm32, mu = 3)

    One Sample t-test

data:  x.norm32
t = -0.28936, df = 999, p-value = 0.7724
alternative hypothesis: true mean is not equal to 3
95 percent confidence interval:
 2.859129 3.104665
sample estimates:
mean of x 
 2.981897 

Interpretacja wyniku!

7.1.1.1 Różne alternatywy

t.test(x.norm32, mu = 3, alternative = "less")

    One Sample t-test

data:  x.norm32
t = -0.28936, df = 999, p-value = 0.3862
alternative hypothesis: true mean is less than 3
95 percent confidence interval:
     -Inf 3.084898
sample estimates:
mean of x 
 2.981897 

Interpretacja wyniku!

t.test(x.norm32, mu = 3, alternative = "greater")

    One Sample t-test

data:  x.norm32
t = -0.28936, df = 999, p-value = 0.6138
alternative hypothesis: true mean is greater than 3
95 percent confidence interval:
 2.878896      Inf
sample estimates:
mean of x 
 2.981897 

Interpretacja wyniku!

7.1.2 Testy dwóch średnich

t.test(x.norm32, x.norm01)

    Welch Two Sample t-test

data:  x.norm32 and x.norm01
t = 42.76, df = 1469.1, p-value < 2.2e-16
alternative hypothesis: true difference in means is not equal to 0
95 percent confidence interval:
 2.853743 3.128154
sample estimates:
   mean of x    mean of y 
 2.981896891 -0.009051554 

Interpretacja wyniku!

t.test(x.norm32[1:500], x.norm32[501:1000])

    Welch Two Sample t-test

data:  x.norm32[1:500] and x.norm32[501:1000]
t = 0.46174, df = 997.23, p-value = 0.6444
alternative hypothesis: true difference in means is not equal to 0
95 percent confidence interval:
 -0.1878366  0.3034311
sample estimates:
mean of x mean of y 
 3.010796  2.952998 

Interpretacja wyniku!

7.1.3 Testy dla dwóch wariancji

var.test(x.norm01, x.norm32)

    F test to compare two variances

data:  x.norm01 and x.norm32
F = 0.25, num df = 999, denom df = 999, p-value < 2.2e-16
alternative hypothesis: true ratio of variances is not equal to 1
95 percent confidence interval:
 0.2208247 0.2830300
sample estimates:
ratio of variances 
              0.25 

Interpretacja wyniku!

var.test(x.norm01[1:100], x.norm01[901:100])

    F test to compare two variances

data:  x.norm01[1:100] and x.norm01[901:100]
F = 0.97939, num df = 99, denom df = 801, p-value = 0.9215
alternative hypothesis: true ratio of variances is not equal to 1
95 percent confidence interval:
 0.7406962 1.3406822
sample estimates:
ratio of variances 
          0.979391 

Interpretacja wyniku!

var.test(x.norm01, x.norm32)

    F test to compare two variances

data:  x.norm01 and x.norm32
F = 0.25, num df = 999, denom df = 999, p-value < 2.2e-16
alternative hypothesis: true ratio of variances is not equal to 1
95 percent confidence interval:
 0.2208247 0.2830300
sample estimates:
ratio of variances 
              0.25 

Interpretacja wyniku!

mood.test(x.norm01, x.norm32)

    Mood two-sample test of scale

data:  x.norm01 and x.norm32
Z = -5.3831, p-value = 7.323e-08
alternative hypothesis: two.sided

Interpretacja wyniku!

ansari.test(x.norm01, x.norm32)

    Ansari-Bradley test

data:  x.norm01 and x.norm32
AB = 541388, p-value = 2.408e-10
alternative hypothesis: true ratio of scales is not equal to 1

Interpretacja wyniku!

7.1.4 Test równości frakcji

prop.test(470, 1000, .5)

    1-sample proportions test with continuity correction

data:  470 out of 1000, null probability 0.5
X-squared = 3.481, df = 1, p-value = 0.06208
alternative hypothesis: true p is not equal to 0.5
95 percent confidence interval:
 0.4387437 0.5014896
sample estimates:
   p 
0.47 

Interpretacja wyniku!

7.2 Testy niemaprametryczne

7.2.1 Testy zgodności z rozkładem normalnym z biblioteki nortest

library(nortest)

cvm.test(x.norm01)

    Cramer-von Mises normality test

data:  x.norm01
W = 0.10747, p-value = 0.08869

Interpretacja wyniku!

ad.test(x.norm01)

    Anderson-Darling normality test

data:  x.norm01
A = 0.52378, p-value = 0.1819

Interpretacja wyniku!

shapiro.test(x.norm01)

    Shapiro-Wilk normality test

data:  x.norm01
W = 0.99829, p-value = 0.428

Interpretacja wyniku!

lillie.test(x.norm01)

    Lilliefors (Kolmogorov-Smirnov) normality test

data:  x.norm01
D = 0.028205, p-value = 0.05865

Interpretacja wyniku!

pearson.test(x.norm01)

    Pearson chi-square normality test

data:  x.norm01
P = 31.04, p-value = 0.3635

Interpretacja wyniku!

detach(package:nortest)

7.2.2 Inne testy zgosności z rozkładem normalnym

fBasics::dagoTest(x.norm01)

Title:
 D'Agostino Normality Test

Test Results:
  STATISTIC:
    Chi2 | Omnibus: 0.8855
    Z3  | Skewness: -0.4057
    Z4  | Kurtosis: 0.8491
  P VALUE:
    Omnibus  Test: 0.6423 
    Skewness Test: 0.685 
    Kurtosis Test: 0.3959 

Description:
 Thu May 14 22:08:00 2020 by user: user

7.2.3 Zodności z rozkładem jednostajnym (test chi-kwdrat)

set.seed(20200512)
x.unif <- runif(1000, 1, 3)

chisq.test(x.unif)

Aproksymacja chi-kwadrat mo戼㹦e by攼㸶 niepoprawna

    Chi-squared test for given probabilities

data:  x.unif
X-squared = 166.57, df = 999, p-value = 1

Interpretacja wyniku!

chisq.test(x.norm32 - min(x.norm32))

    Chi-squared test for given probabilities

data:  x.norm32 - min(x.norm32)
X-squared = 618.77, df = 999, p-value = 1

Interpretacja wyniku!

7.2.4 Zgodności z rozkładem wzorcowym (Kołmogorowa-Smirnowa)

ks.test(x.unif, "punif")

    One-sample Kolmogorov-Smirnov test

data:  x.unif
D = 1, p-value < 2.2e-16
alternative hypothesis: two-sided

Interpretacja wyniku!

ks.test(x.unif, "punif", 1, 3)

    One-sample Kolmogorov-Smirnov test

data:  x.unif
D = 0.022248, p-value = 0.7054
alternative hypothesis: two-sided

Interpretacja wyniku!

ks.test(x.norm32, "pnorm")

    One-sample Kolmogorov-Smirnov test

data:  x.norm32
D = 0.70315, p-value < 2.2e-16
alternative hypothesis: two-sided

Interpretacja wyniku!

ks.test(x.norm32, "pnorm", 3, 2)

    One-sample Kolmogorov-Smirnov test

data:  x.norm32
D = 0.024589, p-value = 0.581
alternative hypothesis: two-sided

Interpretacja wyniku!

7.2.5 Zgodności dwóch rozkładów (Kołmogorowa-Smirnowa)

ks.test(x.norm01, x.norm32)

    Two-sample Kolmogorov-Smirnov test

data:  x.norm01 and x.norm32
D = 0.714, p-value < 2.2e-16
alternative hypothesis: two-sided

Interpretacja wyniku!

ks.test(x.norm01, x.norm32, alternative = "greater")

    Two-sample Kolmogorov-Smirnov test

data:  x.norm01 and x.norm32
D^+ = 0.714, p-value < 2.2e-16
alternative hypothesis: the CDF of x lies above that of y

Interpretacja wyniku!

ks.test(x.norm01, x.norm32, alternative = "less")

    Two-sample Kolmogorov-Smirnov test

data:  x.norm01 and x.norm32
D^- = 0.002, p-value = 0.996
alternative hypothesis: the CDF of x lies below that of y

Interpretacja wyniku!

7.2.6 Testy niezależności

cor.test(x.norm01, x.norm32)

    Pearson's product-moment correlation

data:  x.norm01 and x.norm32
t = Inf, df = 998, p-value < 2.2e-16
alternative hypothesis: true correlation is not equal to 0
95 percent confidence interval:
 1 1
sample estimates:
cor 
  1 

Interpretacja wyniku!

cor.test(x.norm01, df$ga)

    Pearson's product-moment correlation

data:  x.norm01 and df$ga
t = -0.62113, df = 998, p-value = 0.5347
alternative hypothesis: true correlation is not equal to 0
95 percent confidence interval:
 -0.08155156  0.04238688
sample estimates:
        cor 
-0.01965786 

Interpretacja wyniku!

cor.test(x.norm01, df$ga, method = "spearman")

    Spearman's rank correlation rho

data:  x.norm01 and df$ga
S = 167807514, p-value = 0.8288
alternative hypothesis: true rho is not equal to 0
sample estimates:
         rho 
-0.006846091 

Interpretacja wyniku!

chisq.test(tabela)

    Pearson's Chi-squared test with Yates' continuity correction

data:  tabela
X-squared = 1.9825, df = 1, p-value = 0.1591

Interpretacja wyniku!

7.2.7 Miary zależności oparte na statystyce chi kwadrat

7.2.7.1 V Cramera

lsr::cramersV(tabela)

[1] 0.1408003

DescTools::CramerV(tabela)

[1] 0.1615063

7.2.7.2 Czuprowa

DescTools::TschuprowT(tabela)

[1] 0.1615063

7.2.7.3 Phi

DescTools::Phi(tabela)

[1] 0.1615063

7.2.7.4 Kontyngencji C Pearsona

DescTools::ContCoef(tabela)

[1] 0.1594402

7.2.7.5 Yule’a

DescTools::YuleQ(tabela)

[1] -0.3233831

DescTools::YuleY(tabela)

[1] -0.1661555

---

8 Analiza regresji i ANOVA

daneMieszkania <- read_delim("http://www.biecek.pl/R/dane/daneMieszkania.csv", 
                             ";", escape_double = FALSE, trim_ws = TRUE)

Parsed with column specification:
cols(
  cena = col_double(),
  pokoi = col_double(),
  powierzchnia = col_double(),
  dzielnica = col_character(),
  `typ budynku` = col_character()
)

daneMieszkania %>% DT::datatable()

8.1 Zmiana zmiennych dzielnica i typ budynku na zmienne czynnikowe

daneMieszkania <- daneMieszkania %>%
  mutate_if(is.character, list(factor))

8.2 Podsumowanie obserwacji

daneMieszkania %>% summary

      cena            pokoi       powierzchnia         dzielnica      typ budynku
 Min.   : 83280   Min.   :1.00   Min.   :17.00   Biskupin   :65   kamienica :61  
 1st Qu.:143304   1st Qu.:2.00   1st Qu.:31.15   Krzyki     :79   niski blok:63  
 Median :174935   Median :3.00   Median :43.70   Srodmiescie:56   wiezowiec :76  
 Mean   :175934   Mean   :2.55   Mean   :46.20                                   
 3rd Qu.:208741   3rd Qu.:3.00   3rd Qu.:61.40                                   
 Max.   :295762   Max.   :4.00   Max.   :87.70                                   

8.3 Budowa modelu liniowego (regresja)

model <- lm(cena~dzielnica, data = daneMieszkania)

model %>% summary

Call:
lm(formula = cena ~ dzielnica, data = daneMieszkania)

Residuals:
   Min     1Q Median     3Q    Max 
-84893 -31892   -880  29611 106268 

Coefficients:
                     Estimate Std. Error t value Pr(>|t|)    
(Intercept)            189494       5238  36.178  < 2e-16 ***
dzielnicaKrzyki        -21321       7072  -3.015  0.00291 ** 
dzielnicaSrodmiescie   -18351       7699  -2.383  0.01810 *  
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 42230 on 197 degrees of freedom
Multiple R-squared:  0.04873,   Adjusted R-squared:  0.03907 
F-statistic: 5.046 on 2 and 197 DF,  p-value: 0.007294

8.4 Sprawdzanie równości średnich w podgrupach

anova(model)

Analysis of Variance Table

Response: cena
           Df     Sum Sq    Mean Sq F value   Pr(>F)   
dzielnica   2 1.7995e+10 8997691613  5.0456 0.007294 **
Residuals 197 3.5130e+11 1783263361                    
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Interpretacja wyniku!

8.5 Wykres

daneMieszkania %>%
  ggplot(aes(powierzchnia, cena)) +
  geom_point(aes(color = dzielnica)) +
  geom_smooth()

8.6 Sprawdzanie założeń regresji

8.6.1 Zależności liniowa

lmtest::harvtest(model)

    Harvey-Collier test

data:  model
HC = 1.3246, df = 196, p-value = 0.1868

Interpretacja wyniku!

lmtest::raintest(model)

    Rainbow test

data:  model
Rain = 0.98189, df1 = 100, df2 = 97, p-value = 0.5364

Interpretacja wyniku!

8.6.2 Normalność reszt

normtest::jb.norm.test(model$residuals)

    Jarque-Bera test for normality

data:  model$residuals
JB = 5.2583, p-value = 0.0555

Interpretacja wyniku!

8.6.3 Brak autokorelacji

lmtest::dwtest(model)

    Durbin-Watson test

data:  model
DW = 2.1565, p-value = 0.8655
alternative hypothesis: true autocorrelation is greater than 0

Interpretacja wyniku!

8.6.4 Homoskedatyczność (równość wariancji)

lmtest::gqtest(model)

    Goldfeld-Quandt test

data:  model
GQ = 1.0691, df1 = 97, df2 = 97, p-value = 0.3713
alternative hypothesis: variance increases from segment 1 to 2

Interpretacja wyniku!

lmtest::bptest(model)

    studentized Breusch-Pagan test

data:  model
BP = 2.2201, df = 2, p-value = 0.3295

Interpretacja wyniku!

lmtest::hmctest(model)

    Harrison-McCabe test

data:  model
HMC = 0.48313, p-value = 0.362

Interpretacja wyniku!

bartlett.test(cena~dzielnica, data = daneMieszkania)

    Bartlett test of homogeneity of variances

data:  cena by dzielnica
Bartlett's K-squared = 1.3075, df = 2, p-value = 0.5201

Interpretacja wyniku!

fligner.test(cena~dzielnica, data = daneMieszkania)

    Fligner-Killeen test of homogeneity of variances

data:  cena by dzielnica
Fligner-Killeen:med chi-squared = 1.7358, df = 2, p-value = 0.4198

Interpretacja wyniku!

car::leveneTest(residuals(model) ~ dzielnica, data = daneMieszkania)

Levene's Test for Homogeneity of Variance (center = median)
       Df F value Pr(>F)
group   2  0.6934 0.5011
      197               

lawstat::levene.test(residuals(model), daneMieszkania$dzielnica)

    Modified robust Brown-Forsythe Levene-type test based on the absolute deviations from
    the median

data:  residuals(model)
Test Statistic = 0.6934, p-value = 0.5011

lawstat::levene.test(residuals(model), daneMieszkania$dzielnica, location = "mean")

    Classical Levene's test based on the absolute deviations from the mean ( none not
    applied because the location is not set to median )

data:  residuals(model)
Test Statistic = 0.72784, p-value = 0.4842

mood.test(cena~dzielnica, data = daneMieszkania)

Błąd w poleceniu 'mood.test.formula(cena ~ dzielnica, data = daneMieszkania)':
  grupujący czynnik musi mieć dokładnie 2 poziomy

Interpretacja wyniku!

8.7 Wykresy diagnostyczne

plot(model)

8.8 Wykres pudełko-wąsy

daneMieszkania %>%
  ggplot(aes(x = dzielnica, y = cena)) +
  geom_boxplot()

8.9 Test post-hoc

TukeyHSD(aov(model))

  Tukey multiple comparisons of means
    95% family-wise confidence level

Fit: aov(formula = model)

$dzielnica
                           diff       lwr        upr     p adj
Krzyki-Biskupin      -21321.019 -38021.10 -4620.9333 0.0081457
Srodmiescie-Biskupin -18350.541 -36532.88  -168.2053 0.0473579
Srodmiescie-Krzyki     2970.478 -14450.28 20391.2340 0.9145465

agricolae::scheffe.test(aov(model), "dzielnica", group=TRUE, console=TRUE)

Study: aov(model) ~ "dzielnica"

Scheffe Test for cena 

Mean Square Error  : 1783263361 

dzielnica,  means

Alpha: 0.05 ; DF Error: 197 
Critical Value of F: 3.041753 

Groups according to probability of means differences and alpha level( 0.05 )

Means with the same letter are not significantly different.

agricolae::LSD.test(aov(model), "dzielnica", p.adj="none", group=FALSE, console = TRUE)

Study: aov(model) ~ "dzielnica"

LSD t Test for cena 

Mean Square Error:  1783263361 

dzielnica,  means and individual ( 95 %) CI

Alpha: 0.05 ; DF Error: 197
Critical Value of t: 1.972079 

Comparison between treatments means

agricolae::LSD.test(aov(model), "dzielnica", p.adj="bonferroni", group=FALSE, console = TRUE) 

Study: aov(model) ~ "dzielnica"

LSD t Test for cena 
P value adjustment method: bonferroni 

Mean Square Error:  1783263361 

dzielnica,  means and individual ( 95 %) CI

Alpha: 0.05 ; DF Error: 197
Critical Value of t: 2.414597 

Comparison between treatments means

agricolae::duncan.test(aov(model), "dzielnica", group=TRUE, console=TRUE) 

Study: aov(model) ~ "dzielnica"

Duncan's new multiple range test
for cena 

Mean Square Error:  1783263361 

dzielnica,  means

Groups according to probability of means differences and alpha level( 0.05 )

Means with the same letter are not significantly different.