Metody probabilistyczne i statytyka - laboratorium 6 (realizacja)

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

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

1 Realizacja 7 maja 2020r.

1.1 Grupa 1 i 2

1.1.1 Podczytanie biblioteki i danych

rm(list = ls())

library(tidyverse)

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()
)

1.1.2 Przekształcenie zmienncyh napisowych na czynnikowe

daneMieszkania %>% summary

      cena            pokoi       powierzchnia    dzielnica         typ budynku       
 Min.   : 83280   Min.   :1.00   Min.   :17.00   Length:200         Length:200        
 1st Qu.:143304   1st Qu.:2.00   1st Qu.:31.15   Class :character   Class :character  
 Median :174935   Median :3.00   Median :43.70   Mode  :character   Mode  :character  
 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                                        

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

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                                   

1.1.3 Budowa modeli regresji i ANOVA (testowanie równości średnich w podgrupach)

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

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

Coefficients:
         (Intercept)       dzielnicaKrzyki  dzielnicaSrodmiescie  
              189494                -21321                -18351  

Interpretacja wyniku!

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

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

Coefficients:
   dzielnicaBiskupin       dzielnicaKrzyki  dzielnicaSrodmiescie  
              189494                168173                171143  

Interpretacja wyniku!

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!

anova(model_1)

Analysis of Variance Table

Response: cena
           Df     Sum Sq    Mean Sq F value    Pr(>F)    
dzielnica   3 6.2086e+12 2.0695e+12  1160.5 < 2.2e-16 ***
Residuals 197 3.5130e+11 1.7833e+09                      
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Interpretacja wyniku!

model$residuals  %>% head

        1         2         3         4         5         6 
 58125.20 -61372.64 -52810.03  30920.71 -37869.97 -33252.10 

1.1.3.1 Tesowanie założeń regresji

library(normtest)

jb.norm.test(model$residuals)

    Jarque-Bera test for normality

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

Interpretacja wyniku!

library(lmtest)

dwtest(daneMieszkania$cena ~ daneMieszkania$dzielnica)

    Durbin-Watson test

data:  daneMieszkania$cena ~ daneMieszkania$dzielnica
DW = 2.1565, p-value = 0.8655
alternative hypothesis: true autocorrelation is greater than 0

Interpretacja wyniku!

gqtest(daneMieszkania$cena ~ daneMieszkania$dzielnica)

    Goldfeld-Quandt test

data:  daneMieszkania$cena ~ daneMieszkania$dzielnica
GQ = 1.0691, df1 = 97, df2 = 97, p-value = 0.3713
alternative hypothesis: variance increases from segment 1 to 2

Interpretacja wyniku!

1.1.3.2 Wykres pudełko-wąsy

boxplot(daneMieszkania$cena ~ daneMieszkania$dzielnica)

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

daneMieszkania %>%
  ggplot(aes(dzielnica, cena)) +
  geom_boxplot() +
  coord_flip()

detach(package:tidyverse)
detach(package:ggplot2)
detach(package:tibble)
detach(package:tidyr)
detach(package:readr)
detach(package:purrr)
detach(package:dplyr)
detach(package:stringr)
detach(package:forcats)

detach(package:normtest)
detach(package:lmtest)

1.2 Grupa 5

1.2.1 Generowanie ciągóW liczb pseudolosowych z zdanego rozkładu

rm(list = ls())

library(MASS)

set.seed(20200507)
x.norm <- rnorm(100, mean = 1, sd =2)

set.seed(20200507)
x.lnorm <- rlnorm(100, meanlog = .1, sdlog = 2)

set.seed(20200507)
x.gamma <- rgamma(100, 2, 3)

#MLE

1.2.2 Metoda największej wiarogodności (bibliotek MASS) - wyznaczanie parametrów rozkładu

fitdistr(x.norm, "normal")

     mean         sd    
  1.1618982   2.0215240 
 (0.2021524) (0.1429433)

x.norm.fitt <- fitdistr(x.norm, "normal")
x.norm.fitt$estimate

    mean       sd 
1.161898 2.021524 

x.lnorm.fitt <- fitdistr(x.lnorm, "log-normal")
x.lnorm.fitt$estimate

  meanlog     sdlog 
0.2618982 2.0215240 

(x.lnorm.fitt <- fitdistr(x.lnorm, "log-normal"))

    meanlog      sdlog  
  0.2618982   2.0215240 
 (0.2021524) (0.1429433)

(x.lnormGamma.fitt <- fitdistr(x.lnorm, "gamma", list(shape = 3, rate = 2)))

wyprodukowano warto㤼㹣ci NaNwyprodukowano warto㤼㹣ci NaNwyprodukowano warto㤼㹣ci NaNwyprodukowano warto㤼㹣ci NaNwyprodukowano warto㤼㹣ci NaNwyprodukowano warto㤼㹣ci NaNwyprodukowano warto㤼㹣ci NaNwyprodukowano warto㤼㹣ci NaNwyprodukowano warto㤼㹣ci NaNwyprodukowano warto㤼㹣ci NaNwyprodukowano warto㤼㹣ci NaNwyprodukowano warto㤼㹣ci NaNwyprodukowano warto㤼㹣ci NaNwyprodukowano warto㤼㹣ci NaNwyprodukowano warto㤼㹣ci NaNwyprodukowano warto㤼㹣ci NaNwyprodukowano warto㤼㹣ci NaN

     shape         rate   
  0.36475915   0.04496287 
 (0.04138815) (0.00901971)

1.2.3 Podczytanie biblioteki i danych

library(readr)
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()
)

library(tidyverse)

-- Attaching packages --------------------------------------- tidyverse 1.3.0 --
<U+221A> ggplot2 3.3.0     <U+221A> dplyr   0.8.5
<U+221A> tibble  3.0.1     <U+221A> stringr 1.4.0
<U+221A> tidyr   1.0.2     <U+221A> forcats 0.5.0
<U+221A> purrr   0.3.4     
-- Conflicts ------------------------------------------ tidyverse_conflicts() --
x dplyr::filter() masks stats::filter()
x dplyr::lag()    masks stats::lag()
x dplyr::select() masks MASS::select()

daneMieszkania %>% summary

      cena            pokoi       powierzchnia    dzielnica         typ budynku       
 Min.   : 83280   Min.   :1.00   Min.   :17.00   Length:200         Length:200        
 1st Qu.:143304   1st Qu.:2.00   1st Qu.:31.15   Class :character   Class :character  
 Median :174935   Median :3.00   Median :43.70   Mode  :character   Mode  :character  
 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                                        

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

daneMieszkania %>% summary

1.2.4 Testowanie równości średnich (t.test)

t.test((daneMieszkania[daneMieszkania$dzielnica == "Biskupin",])$cena, 
       (daneMieszkania[daneMieszkania$dzielnica == "Krzyki",])$cena)

    Welch Two Sample t-test

data:  (daneMieszkania[daneMieszkania$dzielnica == "Biskupin", ])$cena and (daneMieszkania[daneMieszkania$dzielnica == "Krzyki", ])$cena
t = 2.9793, df = 140.82, p-value = 0.003404
alternative hypothesis: true difference in means is not equal to 0
95 percent confidence interval:
  7173.093 35468.945
sample estimates:
mean of x mean of y 
   189494    168173 

Interpretacja wyniku!

daneMieszkania[daneMieszkania$dzielnica == "Krzyki",]

1.2.5 Regresja i testowanie równości średnich (ANOVA)

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

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

Coefficients:
         (Intercept)       dzielnicaKrzyki  dzielnicaSrodmiescie  
              189494                -21321                -18351  

Interpretacja wyniku!

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!

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

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

Coefficients:
   dzielnicaBiskupin       dzielnicaKrzyki  dzielnicaSrodmiescie  
              189494                168173                171143  

Interpretacja wyniku!

anova(model2)

Analysis of Variance Table

Response: cena
           Df     Sum Sq    Mean Sq F value    Pr(>F)    
dzielnica   3 6.2086e+12 2.0695e+12  1160.5 < 2.2e-16 ***
Residuals 197 3.5130e+11 1.7833e+09                      
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Interpretacja wyniku!

1.2.6 Coś na boku - użycie attach do ramki danych

# Użycie attach
attach(daneMieszkania)

Nast攼㹡puj戼㸹ce obiekty zosta戼㸳y zakryte z daneMieszkania (pos = 13):

    cena, dzielnica, pokoi, powierzchnia, typ budynku

boxplot(cena~dzielnica)

detach(daneMieszkania)

Było 28 ostrzeżenie (użyj 'warnings()' aby je zobaczyć)

1.2.7 Regresja i testowanie równości średnich (ANOVA) c.d.

colnames(daneMieszkania)[5] <- "typ"

model3 <- lm(cena ~ typ, data = daneMieszkania)
model3

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

Coefficients:
  (Intercept)  typniski blok   typwiezowiec  
       178318          10473         -14955  

Interpretacja wyniku!

anova(model3)

Analysis of Variance Table

Response: cena
           Df     Sum Sq    Mean Sq F value   Pr(>F)   
typ         2 2.2770e+10 1.1385e+10  6.4725 0.001895 **
Residuals 197 3.4653e+11 1.7590e+09                    
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Interpretacja wyniku!

model4 <- lm(cena ~ typ - 1, data = daneMieszkania)
model4

Call:
lm(formula = cena ~ typ - 1, data = daneMieszkania)

Coefficients:
 typkamienica  typniski blok   typwiezowiec  
       178318         188791         163363  

Interpretacja wyniku!

anova(model4)

Analysis of Variance Table

Response: cena
           Df     Sum Sq    Mean Sq F value    Pr(>F)    
typ         3 6.2133e+12 2.0711e+12  1177.4 < 2.2e-16 ***
Residuals 197 3.4653e+11 1.7590e+09                      
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Interpretacja wyniku!

1.2.8 Testowanie założeń regresji

model4$residuals %>% length()

[1] 200

ks.test(model4$residuals, "pnorm")

    One-sample Kolmogorov-Smirnov test

data:  model4$residuals
D = 0.505, p-value < 2.2e-16
alternative hypothesis: two-sided

Interpretacja wyniku!

nortest::ad.test(model4$residuals)

    Anderson-Darling normality test

data:  model4$residuals
A = 0.63027, p-value = 0.09918

Interpretacja wyniku!

tseries::jarque.bera.test(model4$residuals)

    Jarque Bera Test

data:  model4$residuals
X-squared = 4.0672, df = 2, p-value = 0.1309

Interpretacja wyniku!

detach(package:tidyverse)
detach(package:ggplot2)
detach(package:tibble)
detach(package:tidyr)
detach(package:readr)
detach(package:purrr)
detach(package:dplyr)
detach(package:stringr)
detach(package:forcats)

detach(package:MASS)

2 Realizacja 14 maja 2020r.

2.1 Grupa 1 i 2

2.1.1 Kilka uwag o pozostałych testach nieparametrycznych

rm(list = ls())

2.1.1.1 Funkcja sample

set.seed(908)
probka <- sample(c("x", "y", "z"), 100, replace = TRUE, prob = c(.3, .5, .2))

summary(probka)

   Length     Class      Mode 
      100 character character 

probka

  [1] "z" "z" "x" "y" "x" "y" "y" "x" "y" "y" "y" "x" "y" "y" "x" "z" "z" "y" "y" "y" "y" "y" "y"
 [24] "x" "z" "y" "x" "x" "y" "x" "y" "x" "y" "x" "y" "y" "x" "y" "x" "y" "y" "y" "z" "y" "y" "x"
 [47] "y" "y" "y" "x" "x" "y" "y" "x" "y" "y" "y" "x" "x" "y" "y" "z" "x" "z" "x" "y" "y" "z" "y"
 [70] "y" "x" "x" "y" "y" "x" "y" "y" "y" "y" "y" "z" "x" "x" "z" "z" "x" "y" "x" "x" "z" "y" "x"
 [93] "z" "y" "y" "y" "z" "x" "z" "y"

2.1.1.2 Funkcja data.frame

probka2 <- sample(c("a", "b"), 100, replace = TRUE, prob = c(.8, .2))
df <- data.frame(probka = probka, probka2 = probka2)

2.1.1.3 Funkcja table - tablica kontyngencji

tablica <- table(df$probka, df$probka2)
tablica

   
     a  b
  x 21 10
  y 41 12
  z 10  6

2.1.1.4 Test niezależności chi-kwadrat

chisq.test(tablica)

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

    Pearson's Chi-squared test

data:  tablica
X-squared = 1.7499, df = 2, p-value = 0.4169

Interpretacja wyniku!

set.seed(908)
x.norm<-rnorm(1000,5,2)
set.seed(908)
x.gamma<-rgamma(1000,2,3)

2.1.1.5 Testowanie niezależności przy użyciu współczynika korelacji (Pearsona i Spearmana)

Niewłaściwe stosowanie (niespełnione założenie o normalności)

cor.test(x.norm,x.gamma)

    Pearson's product-moment correlation

data:  x.norm and x.gamma
t = -1.9241, df = 998, p-value = 0.05462
alternative hypothesis: true correlation is not equal to 0
95 percent confidence interval:
 -0.122326536  0.001203135
sample estimates:
        cor 
-0.06079448 

Interpretacja wyniku!

cor.test(x.norm,x.gamma,method="spearman")

    Spearman's rank correlation rho

data:  x.norm and x.gamma
S = 177858726, p-value = 0.03374
alternative hypothesis: true rho is not equal to 0
sample estimates:
        rho 
-0.06715342 

Interpretacja wyniku!

cor.test(x.norm[1:100],x.norm[901:1000])

    Pearson's product-moment correlation

data:  x.norm[1:100] and x.norm[901:1000]
t = 1.5964, df = 98, p-value = 0.1136
alternative hypothesis: true correlation is not equal to 0
95 percent confidence interval:
 -0.0384156  0.3448386
sample estimates:
      cor 
0.1592038 

Interpretacja wyniku!

2.1.2 Kilka uwag o pozostałych testach parametrycznych

2.1.2.1 Testowanie równości wariancji

set.seed(908)
x.norm2 <- rnorm(1000,3,2)

Pamiętaj o założeniach!

var.test(x.norm, x.norm2)

    F test to compare two variances

data:  x.norm and x.norm2
F = 1, num df = 999, denom df = 999, p-value = 1
alternative hypothesis: true ratio of variances is not equal to 1
95 percent confidence interval:
 0.8832987 1.1321198
sample estimates:
ratio of variances 
                 1 

Interpretacja wyniku!

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

    F test to compare two variances

data:  x.norm[1:100] and x.norm[901:1000]
F = 1.1009, num df = 99, denom df = 99, p-value = 0.6334
alternative hypothesis: true ratio of variances is not equal to 1
95 percent confidence interval:
 0.7407362 1.6362053
sample estimates:
ratio of variances 
          1.100907 

Interpretacja wyniku!

mood.test(x.norm, x.norm2)

    Mood two-sample test of scale

data:  x.norm and x.norm2
Z = -0.15617, p-value = 0.8759
alternative hypothesis: two.sided

Interpretacja wyniku!

ansari.test(x.norm, x.norm2)

    Ansari-Bradley test

data:  x.norm and x.norm2
AB = 501504, p-value = 0.8764
alternative hypothesis: true ratio of scales is not equal to 1

Interpretacja wyniku!

2.1.2.2 Testowanie frakcji

prop.test(399, 1000, .5)

    1-sample proportions test with continuity correction

data:  399 out of 1000, null probability 0.5
X-squared = 40.401, df = 1, p-value = 2.068e-10
alternative hypothesis: true p is not equal to 0.5
95 percent confidence interval:
 0.3685995 0.4301862
sample estimates:
    p 
0.399 

Interpretacja wyniku!

2.1.2.2.1 Funkcja sample

set.seed(908)
probka3 <- sample(c(0,1), 100, replace=TRUE, prob = c(.4, .6))

sum(probka3)

[1] 63

prop.test(sum(probka3), length(probka3), .6)

    1-sample proportions test with continuity correction

data:  sum(probka3) out of length(probka3), null probability 0.6
X-squared = 0.26042, df = 1, p-value = 0.6098
alternative hypothesis: true p is not equal to 0.6
95 percent confidence interval:
 0.5271463 0.7227373
sample estimates:
   p 
0.63 

Interpretacja wyniku!

prop.test(sum(probka3), length(probka3), .7)

    1-sample proportions test with continuity correction

data:  sum(probka3) out of length(probka3), null probability 0.7
X-squared = 2.0119, df = 1, p-value = 0.1561
alternative hypothesis: true p is not equal to 0.7
95 percent confidence interval:
 0.5271463 0.7227373
sample estimates:
   p 
0.63 

Interpretacja wyniku!

2.1.3 Miary zależności oparte na statystyce chi-kwadrat

DescTools::CramerV(tablica)

[1] 0.1322852

library(DescTools)

ContCoef(tablica)

[1] 0.1311428

Phi(tablica)

[1] 0.1322852

TschuprowT(tablica)

[1] 0.1112382

2.1.4 Regresja, ANOVA i testowanie ich założeń

library(tidyverse)

-- 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()

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 <- daneMieszkania %>%
  mutate_if(is.character, list(factor))

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                                   

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

2.1.4.1 Testowanie autokorelacji

lmtest::dwtest(model)

    Durbin-Watson test

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

# brak podstaw do odrzucenia hioptezy o braku autokorelacji

Interpretacja wyniku!

model$residuals %>% head(10)

         1          2          3          4          5          6          7          8 
 58125.195 -61372.644 -52810.033  30920.707 -37869.975 -33252.104  -2761.984 -69247.315 
         9         10 
 22043.257 -51866.395 

2.1.4.2 Testowanie normalności

nortest::cvm.test(model$residuals)

    Cramer-von Mises normality test

data:  model$residuals
W = 0.102, p-value = 0.1048

Interpretacja wyniku!

nortest::sf.test(model$residuals)

    Shapiro-Francia normality test

data:  model$residuals
W = 0.98656, p-value = 0.05273

Interpretacja wyniku!

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

    Jarque-Bera test for normality

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

Interpretacja wyniku!

detach(package:tidyverse)
detach(package:ggplot2)
detach(package:tibble)
detach(package:tidyr)
detach(package:readr)
detach(package:purrr)
detach(package:dplyr)
detach(package:stringr)
detach(package:forcats)

detach(package:DescTools)

2.2 Grupa 6

2.2.1 Pozostałe testy parametryczne i nieparametryczne

rm(list = ls())

set.seed(1735)
x.norm <- rnorm(1000, mean = 5, sd = 1)

set.seed(1735)
x.gamma <- rgamma(1000, 3, 5)

2.2.2 Testowanie wariancji

var.test(x.norm[1:100], x.norm[301:400])

    F test to compare two variances

data:  x.norm[1:100] and x.norm[301:400]
F = 1.0601, num df = 99, denom df = 99, p-value = 0.7721
alternative hypothesis: true ratio of variances is not equal to 1
95 percent confidence interval:
 0.7132824 1.5755629
sample estimates:
ratio of variances 
          1.060104 

Interpretacja wyniku!

mood.test(x.norm[1:100], x.norm[901:1000])

    Mood two-sample test of scale

data:  x.norm[1:100] and x.norm[901:1000]
Z = 0.004164, p-value = 0.9967
alternative hypothesis: two.sided

Interpretacja wyniku!

ansari.test(x.norm[501:600], x.norm[701:800])

    Ansari-Bradley test

data:  x.norm[501:600] and x.norm[701:800]
AB = 5116, p-value = 0.747
alternative hypothesis: true ratio of scales is not equal to 1

Interpretacja wyniku!

2.2.3 Testowanie frakcji

prop.test(295, 750, .45)

    1-sample proportions test with continuity correction

data:  295 out of 750, null probability 0.45
X-squared = 9.503, df = 1, p-value = 0.002051
alternative hypothesis: true p is not equal to 0.45
95 percent confidence interval:
 0.3583488 0.4294256
sample estimates:
        p 
0.3933333 

Interpretacja wyniku!

prop.test(295, 750, .4)

    1-sample proportions test with continuity correction

data:  295 out of 750, null probability 0.4
X-squared = 0.1125, df = 1, p-value = 0.7373
alternative hypothesis: true p is not equal to 0.4
95 percent confidence interval:
 0.3583488 0.4294256
sample estimates:
        p 
0.3933333 

Interpretacja wyniku!

2.2.3.1 Funkcja sample

proba <- sample(0:1, 1000, replace=TRUE, prob = c(.2, .8))
sum(proba)

[1] 804

prop.test(sum(proba), length(proba), .78)

    1-sample proportions test with continuity correction

data:  sum(proba) out of length(proba), null probability 0.78
X-squared = 3.2182, df = 1, p-value = 0.07282
alternative hypothesis: true p is not equal to 0.78
95 percent confidence interval:
 0.7777307 0.8278955
sample estimates:
    p 
0.804 

Interpretacja wyniku!

z1 <- sample(c("x","y","z"), 100, replace = TRUE, prob = c(.3, .4, .3))
z2 <- sample(c("a","b"), 100, replace = TRUE, prob = c(.7, .3))

2.2.3.2 Tablica kontyngencji i funkcja table

tablica <- table(z1, z2)
tablica 

   z2
z1   a  b
  x 22  9
  y 32 13
  z 17  7

2.2.4 Test niezależności chi-kwadrat

chisq.test(tablica)

    Pearson's Chi-squared test

data:  tablica
X-squared = 0.00060927, df = 2, p-value = 0.9997

Interpretacja wyniku!

2.2.4.1 Funkcja data.frame

df <- data.frame(z1=z1, z2=z2)

2.2.5 Testowanie współczynnika korelacji Pearsona i rang Spearmana

Tak stosowany test chi-kwadrat obowiązuje tylko dla wielowymiarowego rozkładu normalnego

cor.test(x.norm, x.gamma)

    Pearson's product-moment correlation

data:  x.norm and x.gamma
t = 0.70572, df = 998, p-value = 0.4805
alternative hypothesis: true correlation is not equal to 0
95 percent confidence interval:
 -0.03971457  0.08420999
sample estimates:
       cor 
0.02233349 

# tak stosowany test pearsona obowiązuje tylko dla wielowymiarowego

Interpretacja wyniku!

cor.test(x.norm, x.gamma, method = "spearman")

    Spearman's rank correlation rho

data:  x.norm and x.gamma
S = 165662512, p-value = 0.8491
alternative hypothesis: true rho is not equal to 0
sample estimates:
        rho 
0.006023934 

Interpretacja wyniku!

2.2.6 Miary zależności oparte na statystyce chi-kwadrat

DescTools::CramerV(tablica)

[1] 0.002468333

CramerV(tablica)

[1] 0.002468333

DescTools::Phi(tablica)

[1] 0.002468333

DescTools::ContCoef(tablica)

[1] 0.002468325

DescTools::TschuprowT(tablica)

[1] 0.002075612

2.2.7 Testowanie równości średnich (t.test i ANOVA)

library(readr)
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()
)

library(tidyverse)

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

daneMieszkania.Biskupin <- daneMieszkania %>% 
  filter(dzielnica=="Biskupin")

daneMieszkania.Biskupin %>% summary

      cena            pokoi        powierzchnia         dzielnica      typ budynku
 Min.   :120290   Min.   :1.000   Min.   :17.10   Biskupin   :65   kamienica :26  
 1st Qu.:156655   1st Qu.:2.000   1st Qu.:35.20   Krzyki     : 0   niski blok:17  
 Median :189291   Median :3.000   Median :45.10   Srodmiescie: 0   wiezowiec :22  
 Mean   :189494   Mean   :2.585   Mean   :47.05                                   
 3rd Qu.:214462   3rd Qu.:3.000   3rd Qu.:61.20                                   
 Max.   :295762   Max.   :4.000   Max.   :87.70                                   

daneMieszkania.Biskupin %>% dim

[1] 65  5

t.test(daneMieszkania.Biskupin$cena)

    One Sample t-test

data:  daneMieszkania.Biskupin$cena
t = 37.585, df = 64, p-value < 2.2e-16
alternative hypothesis: true mean is not equal to 0
95 percent confidence interval:
 179422 199566
sample estimates:
mean of x 
   189494 

Interpretacja wyniku!

daneMieszkania.krzyki <- daneMieszkania %>% 
  filter(dzielnica=="Krzyki")

t.test(daneMieszkania.krzyki$cena, daneMieszkania.Biskupin$cena)

    Welch Two Sample t-test

data:  daneMieszkania.krzyki$cena and daneMieszkania.Biskupin$cena
t = -2.9793, df = 140.82, p-value = 0.003404
alternative hypothesis: true difference in means is not equal to 0
95 percent confidence interval:
 -35468.945  -7173.093
sample estimates:
mean of x mean of y 
   168173    189494 

Interpretacja wyniku!

daneMieszkania.srodmiescie <- daneMieszkania %>% 
  filter(dzielnica=="Srodmiescie")

t.test(daneMieszkania.Biskupin$cena, daneMieszkania.srodmiescie$cena)

    Welch Two Sample t-test

data:  daneMieszkania.Biskupin$cena and daneMieszkania.srodmiescie$cena
t = 2.5079, df = 117.13, p-value = 0.01351
alternative hypothesis: true difference in means is not equal to 0
95 percent confidence interval:
  3859.764 32841.317
sample estimates:
mean of x mean of y 
 189494.0  171143.5 

Interpretacja wyniku!

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

model

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

Coefficients:
         (Intercept)       dzielnicaKrzyki  dzielnicaSrodmiescie  
              189494                -21321                -18351  

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!

detach(package:tidyverse)
detach(package:ggplot2)
detach(package:tibble)
detach(package:tidyr)
detach(package:readr)
detach(package:purrr)
detach(package:dplyr)
detach(package:stringr)
detach(package:forcats)

# detach(package:normtest)
# detach(package:lmtest)