Statystyka I z R
Warsztat 6. Parametry rozproszenia
+Tomasz Żółtak
+14 listopada 2016
+ +Na dzisiejszych zajęciach poznamy funkcje pozwalające obliczyć typowo wykorzystywane parametry parametry rozproszenia zmiennych statystycznych.
+Wczytanie danych
+Zacznijmy od wczytania danych, na których będziemy dalej pracować. Funkcja load() pozwala wczytać obiekty R zapisane w natywnym formacie R-a, czyli .RData (linijka wcześniej służy upewnieniu się, że bęziemy próbowali wczytać dane z odpowiedniego folderu). Funkcja load() zwraca nazwy wczytanych obiektów - w tym przypadku jest to 15 wektorów. Wektor o nazwie variablesDescription opisuje znaczenie pozostałych wektorów, które zawierają dane - zmienne z dodatkowej próby badawczej uczniów szkół pogimnazjalnych, zrealizowanej w ramach badania PISA 2009.
> try(setwd("warsztat 2016.11.14"), silent = TRUE)
+> nazwyObiektow = load("dane_2016.11.14.RData")
+> nazwyObiektow
+ [1] "schoolType" "sex" "age"
+ [4] "scoreTMR" "scoreKKS" "scoreKNS"
+ [7] "scorePISAMath" "scorePISARead" "scorePISAScie"
+[10] "noPersHous" "wealth" "income"
+[13] "parEdu" "hisei" "variablesDescription"
+> variablesDescription
+ schoolType
+ "School type (track)"
+ sex
+ "Sex of student"
+ age
+ "How old student was in 2009"
+ scoreTMR
+ "Score on Raven's Progressive Matrices test [0-60]"
+ scoreKKS
+ "Score on Cultural Competences Scale [60-240]"
+ scoreKNS
+ "Score on Snyder's Adult Hope Scale [8-64]"
+ scorePISAMath
+"Score on PISA mathematics test (pseudoEAP estimator)"
+ scorePISARead
+ "Score on PISA reading test (pseudoEAP estimator)"
+ scorePISAScie
+ "Score on PISA science test (pseudoEAP estimator)"
+ noPersHous
+ "Number of persons in a household"
+ wealth
+ "Highest parental occupational status"
+ income
+ "Highest parental education in years"
+ parEdu
+ "Wealth"
+ hisei
+ "Houshold monthly income" Parametry rozproszenia
+Rozstęp i odchylenie ćwiartkowe
+W R nie mamy funkcji pozwalającej bezpośrednio obliczyć rozstęp, niemniej zrobienie tego samemu jest oczywiście banalnie proste. Z kolei funkcja IQR(), służąca obliczeniu odchylenia ćwiartkowego robi to nieco inaczej, niż w typowych (stosowanych w Polsce) definicjach (nie dzieli różnicy wartości 3. i 1. kwartyla przez dwa). Niemniej obliczenie samemu odchylenia ćwiartkowego również nie nastręcza trudności (jeśli umiemy obliczyć kwartyle).
> # rozstęp
+> max(scoreTMR, na.rm = TRUE) - min(scoreTMR, na.rm = TRUE)
+[1] 49
+> # odchylenie ćwiartkowe
+> IQR(scoreTMR, na.rm = TRUE) / 2
+[1] 4.5
+> # równoważnie
+> q = quantile(scoreTMR, na.rm = TRUE)
+> q
+ 0% 25% 50% 75% 100%
+ 11 45 50 54 60
+> (q[4] - q[2]) / 2
+75%
+4.5 +
Wariancja i odchylenie standardowe
+Do obliczenia wariancji i odchylenia standardowego domyślnie wykorzystue się w R funkcje odpowiednio var() i sd(). W praktyce jest jednak z nimi pewien problem, bowiem zwracają one wartości, które są nieobciążonymi przewidywaniami (odpowiednio wariancji i odchylenia standardowego), jeśli zadany im wektor traktować jako prostą próbę losową z populacji, dla której chcemy oszacować wartość danego parametru.
Od strony technicznej oznacza to, że w przypadku wariancji suma kwadratów różnic od średniej jest dzielona nie przez liczbę elmentów wektora, lecz przez liczbę elementów pomniejszoną o jeden. Jeśli wektor jest długi, robi to niewielką różnicę, jednak jeśli jest krótki, będzie ona bardzo wyraźna.
+> x = 0:10 # niewielka liczba obserwacji
+> y = rep(x, 10) # spora liczba obserwacji
+> z = rep(x, 100) # duża liczba obserwacji
+> srednie = c(x = mean(x), y = mean(y), z = mean(z))
+> # wariancja obliczona jako oszacowanie na podstawie prostej próby losowej
+> warOszac = c(x = var(x), y = var(y), z = var(z))
+> # wariancja obliczona na całej populacji (jeden ze sposobów)
+> warPopul = c(x = mean(x^2) - mean(x)^2, y = mean(y^2) - mean(y)^2,
++ z = mean(z^2) - mean(z)^2)
+> # porównajmy
+> round(rbind(srednie, warOszac, warPopul), 1)
+ x y z
+srednie 5 5.0 5
+warOszac 11 10.1 10
+warPopul 10 10.0 10
+> # to się w ogóle nie uda!
+> var(1)
+[1] NAJeśli chcemy uzyskać normalną wariancję/odchylenie standardowe musimy więc obliczyć ją sobie sami, lub skorygować wynik działania funkcji var() lub sd():
-
+
- Wariancja policzona na piechotę ze wzoru definicyjnego: +
> varTMR = mean((scoreTMR - mean(scoreTMR, na.rm = TRUE))^2, na.rm = TRUE)
+> sdTMR = varTMR^0.5
+> varTMR
+[1] 54.03407
+> sdTMR
+[1] 7.350787
+> # dla porównania
+> var(scoreTMR, na.rm = TRUE)
+[1] 54.04975
+> sd(scoreTMR, na.rm = TRUE)
+[1] 7.351854-
+
- Wariancja policzona na piechotę z alternatywnego wzoru, jako różnica średniej kwadratów wartości zmiennej i kwadratu średniej tej zmiennej: +
> var2TMR = mean(scoreTMR^2, na.rm = TRUE) - mean(scoreTMR, na.rm = TRUE)^2
+> sd2TMR = var2TMR^0.5
+> var2TMR
+[1] 54.03407
+> sd2TMR
+[1] 7.350787-
+
- Skorygowanie (przeliczenie) wyników działania funkcji
var()lubsd():
+
> poprawka = (length(scoreTMR) - 1) / length(scoreTMR)
+> poprawka
+[1] 0.9997366
+> var(scoreTMR^2, na.rm = TRUE) * poprawka
+[1] 423791.9
+> sd(scoreTMR, na.rm = TRUE) * poprawka^0.5
+[1] 7.350885Współczynnik zmienności
+Współczynniki zmienności są często bardziej użyteczną miarą zróżnicowania w przypadku zmiennych mierzonych na skalach ilorazowych, które przyjmują tylko wartości nieujemne, a więc mają zakres wartości ograniczony z jednej strony, ale otwarty z drugiej. Pozwalają bowiem uwzględnić, że w takim przypadku w ramach grupy o wysokiej średniej wartości zmiennej jest znacznie więcej miejsca na różnicowanie się wartości, niż w ramach grupy, w której średnia wartości tej zmiennej jest niska, a więc bliska granicy zakresu możliwych do przyjęcia wartości.
+Najszerzej wykorzystywany typ wskaźnika zmieności obliczamy dzieląc wartość odchylenia standardowego przez wartość średniej.
+Aby prześledzić opisane wyżej własności współczynika zmienności w stosunku do odchylenia standardowego, rozpatrzmy dwie zmienne: x i y, które mogą przyjmować tylko wartości nieujemne i są określone w 10-cio elementowej zbiorowości. Przyjmiemy przy tym arbitralnie, że średnia zmiennej x wynosi 1, a średnia zmiennej y wynosi 10. Następnie dobierzemy wartości tych zmiennych dla poszczególnych jednostek zbiorowości w ten sposób, aby zmaksymalizować wartość odchylenia standardowego (w ramach przyjętych założeń, że zmienne mają opisane wyżej średnie i nie mogą przyjmować wartości ujemnych).
> # dla 10-cio elementowej zbiorowości odch. std. zm. x i y maksymalizuje
+> # takie przypisanie wartości tych zmiennych:
+> x = c(rep(0, 9), 10)
+> y = c(rep(0, 9), 100)
+> # obejrzyjmy rozkłady (pionowa niebieska linia wskazuje średnią)
+> par(mfcol = c(1, 2))
+> hist(x, xlim = c(0, 100), breaks = 0:100, col = 2)
+> abline(v = mean(x), lwd = 2, col = 4)
+> hist(y, xlim = c(0, 100), breaks = 0:100, col = 2)
+> abline(v = mean(y), lwd = 2, col = 4)
> # obliczmy parametry rozkładów tych zmiennych
+> srednie = c(x = mean(x), y = mean(y))
+> odchStd = c(x = mean(x^2) - mean(x)^2, y = mean(y^2) - mean(y)^2)^0.5
+> round(rbind(srednie, odchStd), 1)
+ x y
+srednie 1 10
+odchStd 3 30
+> # jak widać, odchylenia standardowe dramatycznie się różnią
+> wspZm = odchStd / srednie
+> round(rbind(srednie, odchStd, wspZm), 1)
+ x y
+srednie 1 10
+odchStd 3 30
+wspZm 3 3
+> # ale współczynniki zmienności są już takie same+
Zadanie
+Jaka zmienna/zmienne spośród wczytanych na początku zajęć spełnia warunki umożliwiające policzenie współczynnika zmienności? Odpowiedz na podstawie opisu zmiennych i wartości parametrów wartości (kod poniżej).
+> variablesDescription
+ schoolType
+ "School type (track)"
+ sex
+ "Sex of student"
+ age
+ "How old student was in 2009"
+ scoreTMR
+ "Score on Raven's Progressive Matrices test [0-60]"
+ scoreKKS
+ "Score on Cultural Competences Scale [60-240]"
+ scoreKNS
+ "Score on Snyder's Adult Hope Scale [8-64]"
+ scorePISAMath
+"Score on PISA mathematics test (pseudoEAP estimator)"
+ scorePISARead
+ "Score on PISA reading test (pseudoEAP estimator)"
+ scorePISAScie
+ "Score on PISA science test (pseudoEAP estimator)"
+ noPersHous
+ "Number of persons in a household"
+ wealth
+ "Highest parental occupational status"
+ income
+ "Highest parental education in years"
+ parEdu
+ "Wealth"
+ hisei
+ "Houshold monthly income"
+> summary(cbind(age, scoreTMR, scoreKKS, scoreKNS, scorePISAMath, scorePISARead,
++ scorePISAScie, noPersHous, wealth, income, parEdu, hisei))
+ age scoreTMR scoreKKS scoreKNS
+ Min. :15.79 Min. :11.00 Min. : 60.0 Min. : 8.00
+ 1st Qu.:17.17 1st Qu.:45.00 1st Qu.:162.0 1st Qu.:43.00
+ Median :17.45 Median :50.00 Median :178.0 Median :48.00
+ Mean :17.47 Mean :48.54 Mean :178.3 Mean :47.74
+ 3rd Qu.:17.72 3rd Qu.:54.00 3rd Qu.:195.0 3rd Qu.:53.00
+ Max. :20.21 Max. :60.00 Max. :240.0 Max. :64.00
+ NA's :322 NA's :350 NA's :325 NA's :337
+ scorePISAMath scorePISARead scorePISAScie noPersHous
+ Min. :199.2 Min. :231.7 Min. :226.8 Min. : 2.000
+ 1st Qu.:457.2 1st Qu.:467.3 1st Qu.:473.9 1st Qu.: 4.000
+ Median :528.2 Median :534.8 Median :538.6 Median : 4.000
+ Mean :523.5 Mean :524.5 Mean :531.0 Mean : 4.536
+ 3rd Qu.:588.2 3rd Qu.:589.0 3rd Qu.:596.2 3rd Qu.: 5.000
+ Max. :793.7 Max. :749.1 Max. :756.3 Max. :14.000
+ NA's :86
+ wealth income parEdu hisei
+ Min. :-2.8256 Min. : 1.136 Min. : 3.00 Min. :16.00
+ 1st Qu.:-0.7216 1st Qu.: 5.305 1st Qu.:11.00 1st Qu.:33.00
+ Median :-0.2500 Median : 7.494 Median :11.00 Median :43.00
+ Mean :-0.1678 Mean : 9.598 Mean :12.34 Mean :45.86
+ 3rd Qu.: 0.2512 3rd Qu.:10.818 3rd Qu.:12.00 3rd Qu.:55.00
+ Max. : 2.7703 Max. :68.456 Max. :16.00 Max. :90.00
+ NA's :9 NA's :9 NA's :43 NA's :111 Oblicz wartość współczynnika zmienności dla tej zmiennej (zmienych).
+> # to jest miejsce na TWój kod+
Odchylenie przeciętne od mediany
++
Zadanie
+Oblicz odchylenie przeciętne od mediany zmiennej scoreTMR, pamiętając że:
-
+
- Odchylenie przeciętne od mediany to średnia modułóW różnic pomiędzy wartością zmiennej dla danej jednostki obserwacji a wartością mediany. +
- Do obliczenia modułu liczby można wykorzystać funkcję
abs().
+
> # to jest miejsce na Twój kod+
Na następne zajęcia
+Praca domowa
+Zostanie nadesłana mailem.
+Do przeczytania na następne zajęcia
+Zostanie nadesłana mailem.
+