Zgodnie z zapowiedziami dzisiaj będzie kilka technicznych komentarzy dotyczących wtorkowego wpisu. Pojawi się temat ważenia obserwacji, analizy wariancji, pojawi się bardzo przydatna funkcja by() z pakietu R i oczywiście opis krok po kroku jak zrobić wykres prezentowany w poprzednim odcinku.
Zakładam, że dane z badania PISA są już wczytane i dostępne w ramce danych o nazwie dane. Będzie mnie interesowała średnia wyników gimnazjalistów z testów językowych i testów matematycznych w podziale na płeć i na liczbę książek w domu.
W kolumnach PV1MATH i PV1READ znajdują się wartości, które będę traktował jako wyniki z, odpowiednio, testów matematycznych i językowych. W kolumnie ST04Q01 znajduje się zakodowana płeć osobnika (wartości 1 – kobieta, 2 – mężczyzna) a w kolumnie ST22Q01 zakodowana jest liczba książek w domu.
Zanim zaczniemy analizę dwa komentarze do tych danych.
Zmienna ST22Q01 została zakodowana jako zmienna jakościowa na 6 poziomach (wartość 1 odpowiada '0-10 books’, …, wartość 6 odpowiada 'More than 5000 books’). Nie dla wszystkich uczniów informacja o liczbie książek jest dostępna. Uczeń mógł nie wypełnić tego pola, zaznaczyć opcję 'odmawiam odpowiedzi’
lub z innego powodu nie podać liczby książek. Zwyczajowo jeżeli zmienną można zakodować na co najwyżej 6. poziomach to brakujące obserwacje są kodowane wartościami 7, 8 i 9 (w zależności od powodu niedostępności, nie dotyczy, niepoprawna odpowiedź, brak odpowiedzi). Tak jest też w tym przypadku. W kolumnie ST22Q01 występują wartości większe niż 6, które będziemy musieli w dalszym etapie odpowiednio uwzględnić (usunąć).
Drugi komentarz dotyczy ważenia. Z przyczyn technicznych podczas badania PISA nie wylosowano z tym samym prawdopodobieństwem z populacji wszystkich 15 latków (około 440 tys osób) próby 5 tys osób. Zdecydowanie łatwiej było wylosować pewną grupę szkół z której następnie losowano uczniów. Takie losowanie ma jednak tę wadę, że struktura wylosowanych szkół, lub klas, lub innych zmiennych, może nie odpowiadać strukturze całej populacji. Może nam się wylosować więcej szkół z dużych miast, albo szkół prywatnych itp. Aby móc próbować odnosić te wyniki do całej populacji potrzebne jest ważenie przypadków uwzględniający takie nierównomierne losowanie. Opis metodologii liczenia wag przedstawiony jest w raporcie technicznym PISA. Dla nas ważne będzie, że w kolumnie W_FSTUWT przedstawione są wagi, które w uproszczeniu można interpretować jako liczbę wszystkich 15-latków, którą reprezentuje ta konkretna obserwacja. W przypadku uczniów z Polski wagi przyjmują wartości od 17 do 293, czyli niektórzy uczniowie będą ponad 15 razy ”ważniejsi” niż inni. Ale 50% obserwacji ma wagi od 84 do 100 więc rozkład wag jest w miarę skupiony wokół mediany 91.
Ok, jesteśmy gotowi by policzyć średnie. Poniżej wykorzystujemy funkcję by() aby policzyć ważone średnie (trzeci argument to wywołanie funkcji weighted.mean()) dla cechy PV1MATH z wagami W_FSTUWT w grupach wyznaczonych przez dwie zmienne grupujące: płeć ST04Q01 i liczbę książek w domu ST22Q01. Wynikiem jest macierz o wymiarach 2×9 (2 płcie 9 grup zmiennej ST22Q01) z ważonymi średnimi z matematyki.
Następnie robimy to samo dla wyników z testów językowych i dla sumy wag.
1 2 3 4 5 6 7 8 9 | resMath <- by(dane[, c("PV1MATH", "W_FSTUWT")], list(dane[, "ST04Q01"], dane[, "ST22Q01"]), function(x) weighted.mean(x[, 1], x[, 2])) resRead <- by(dane[, c("PV1READ", "W_FSTUWT")], list(dane[, "ST04Q01"], dane[, "ST22Q01"]), function(x) weighted.mean(x[, 1], x[, 2])) resCount <- by(dane[, "W_FSTUWT"], list(dane[, "ST04Q01"], dane[, "ST22Q01"]), sum) |
Mamy policzone średnie, czas je narysować.
Zaczniemy od starannego przygotowania tła. Funkcją plot() narysujemy szare kropki odpowiadające wynikom poszczególnych uczniów. Kropki będą w kolorze „#33333333” czyli półprzezroczysty jasno-szary. Zakres osi ustawiamy na 400-600 co pozwala na narysowanie 3/4 uczniów a jednocześnie ułatwi porównywanie średnich. Następnie funkcją abline() rysujemy szare poziome i pionowe linie dla pomocniczych siatek. Znacznie to ułatwia odczytywanie wartości z wykresu. Dorysowujemy też przekątną, ci co ponad nią mają lepsze wyniki z testów językowych, ci co poniżej niej mają lepsze wyniki z matematyki.
10 11 12 13 14 15 | plot(dane$PV1MATH, dane$PV1READ, col = "#33333333", cex = 0.2, pch = 19, xlim = c(400, 600), ylim = c(400, 600), las = 1, bty = "n", ylab = "Results in Reading", xlab = "Results in Mathematics") abline(0, 1, lty = 3, col = "grey") abline(h = seq(400, 600, 25), lty = 3, col = "grey") abline(v = seq(400, 600, 25), lty = 3, col = "grey") |
Teraz czas dorysować do wykresu uprzednio wyliczone średnie. Średnie w płciach są opisane przez kolejne wiersze tabel resMath i resRead. Aby nie kopiować kodu dodajemy pętlę for iterującą po wierszach tych tabel. Pierwszy wiersz rysowany będzie trójkątami (pch=17), drugi okręgami (pch=19). Rysujemy tylko kolumny od 1 do 6, a więc nie zaznaczamy średnich dla osób które nie odpowiedziały ile mają książek w domu (kolumny 7-9). Punkty rysujemy z różną wielkością, w zależności od liczebności reprezentowanej populacji, po to nam macierz rozmiar.
16 17 18 | rozmiar <- sqrt(resCount[, 1:6])/100 for (i in 1:2) points(resMath[i, 1:6], resRead[i, 1:6], pch = c(17, 19)[i], cex = rozmiar[i,], type = "o") |
Pozostała kosmetyka. Dodajemy do wykresu legendę i napisy, oznaczające które punkty odpowiadają którym odpowiedziom.
19 20 21 22 23 24 | etykiety <- c("0-10 books", "11-25 books", "26-100 books", "101-200 books", "201-500 books", "More than 500 books") text(resMath[i, 1:6], resRead[i, 1:6], etykiety, adj = c(0, 1)) legend("topleft", c("Male", "Female"), bty = "n", lwd = 2, pch = c(17, 19), col = "red4", cex = 1.8) |
I oto gotowy wykres.
Na koniec użyjemy jeszcze funkcji summary.lm(), a dokładniej właściwości r.squared z wyniku tej funkcji, by pokazać jaki procent zmienności jest wyjaśnione przez zmienne płeć i liczba książek w domu. Zauważmy, że i tutaj uwzględniamy wagi poszczególnych obserwacji.
25 26 27 28 29 30 31 | summary(lm(PV1READ ~ factor(ST04Q01) * factor(ST22Q01), data = dane, weights = dane$W_FSTUWT))$r.squared ## [1] 0.2476 summary(lm(PV1MATH ~ factor(ST04Q01) * factor(ST22Q01), data = dane, weights = dane$W_FSTUWT))$r.squared ## [1] 0.1646 |
I jak opis techniczny? Nie za szczegółowo / za mało szczegółowo?
Dzięki za opis. Jak dla mnie szczegółów w sam raz 🙂
Zastanawiam się nad ostatnim przykładem i ważoną regresją. W szczególności, czy argument 'weights’ do funkcj 'lm’ rzeczywiście robi to, czego oczekujemy. Poniżej mały eksperyment, który porównuje dwie metody szacowania regresji z wagami:
1. via 'lm’ i argument 'weights’ („weightedLm”),
2. używając funkcji z pakietu „survey” („weightedSurv”)
3. nieuwzględniając wag („unweighted”)
Z eksperymentu wynika, że o ile oszacowania parametrów są poprawne w przypadkach (1) i (2), to obie metody przeszacowują błędy standardowe. Metoda (1) obliczna poprawne R-kwadrat.
Zastanawiam się, czemu te ważone 'lm’ nie liczy poprawnych błędów standardowych i dlaczego obie metody dają różne wyniki. Może Ty masz jakiś pomysł/wytłumaczenie?
—e-k-s-p-e-r-y-m-e-n-t—
# generujemy dane z 'y’, 'x’ i wagami 'w’
set.seed(1234)
d <- data.frame(x=1:5)
d$y <- with(d, 2*x + 5 + rnorm(5))
w <- sample(1:5, 30, replace=TRUE)
d$w <- table(w)
# dane2 zmultiplikowane wagą 'w'
d2 <- d[w,c("x", "y")]
# Dla pewnosci: frekwencje 'x' w 'd2' powinny byc rowne wagom w 'd'
d[,c("x", "w")]
table(d2$x)
# modele
models <- list()
# "prawdziwy" dla danych rozagregowanych
models$true <- lm(y~x, data=d2)
# z gruntu błędny model dla danych zagregowanych bez wagi
models$unweighted <- lm(y ~ x, data=d)
# dla danych zagregowanych via 'lm(.., weights=)'
models$weightedLm <- lm(y~x, data=d, weights=w)
# dla danych zagregowanych via pakiet "survey"
des <- svydesign(~1, data=d, weights=as.numeric(d$w))
models$weightedSurv <- svyglm(y ~ x, design=des)
s <- lapply(models, summary)
# porownanie
r <- sapply(s, function(m) coef(m)[,1:2])
rownames(r) <- c("int", "x", "se.int", "se.x")
r
dotchart(t(r))
# 'lm' z wagami daje poprawne R^2
sapply(s, "[[", "r.squared")
Ciekawe pytanie.
Afaik, w tych dwóch funkcjach parametr weights oznacza coś innego, w przypadku lm wagi dotyczą odwrotności wariancji y, która jest pochodną tego, że y_i jest traktowana jako średnia z w_i obserwacji (czyli w_i = 4 odpowiada temu, że ten wiersz to średnia z czterech obserwacji). W funkcji svydesign wagi są proporcjonalne do prawdopodobieństw wylosowania danej obserwacji (sampling weights). Czyli tutaj wariancja 'y’ jest taka sama, bo każdy 'y’ to jedna wylosowana obserwacja.
Dla oceny efektów to bez znaczenia, dlatego w Twoim przykładzie wszystkie efekty przy x są równe sobie a ja mogłem wykorzystać funkcję lm do obliczeń (jest szybsza niż svyglm).
Dla oceny sd to ma znaczenie, funkcja lm traktuje obserwacje z dużą wagą jako 'pewniejsze’ a funkcja svydesign jako przedstawicieli 'częstszych’. W tym przypadku gdyby chcieć wykorzystywać informację o błędach standardowych ocen to pewnie lepiej skorzystać z svyglm (ale ta nie liczy R2).
Ok, dzięki. To chyba rzeczywiście ma sens (wagi dla wariancji y|x vs wagi zwiazane z sampliing probs).
Swoją drogą czy nie jest tak, że przeważnie to własnie wielkości efektów są interesujące (i przedziały ufności wokół nich) a nie wartości R^2? W końcu nie przewidujemy pogody tylko chcemy poznać mechanizm…
Jeszcze co do eksperymentu, zastanawia mnie cały czas, czemu wyniki otrzymane za pomocą 'svyglm’ nie są identyczne z wynikami 'lm’ dla niezagregowanych danych. Oba zbiory danych zawierają identyczną informację, a wagi zostały użyte jako „replication weights”… Może ponapastuję Thomasa Lumleya… 😉
Co do wielkości efektów to pełna zgoda, co do przedziałów ufności… też, choć wolę te oparte o bootstrap.
Zgadzam się też że R^2 raczej trzeba odradzać bo jest nagminnie nadużywany, ale jest dosyć wygodny w automatycznym przeszukiwaniu małych modeli.
Dlatego na wykresach pokazywane są średnie, ale potrzebowałem R^2 by móc porównać ,,ważność” efektu liczby książek (6 poziomów o podobnej liczebności) i efektu wykształcenia (5 poziomów z czego 2 bardzo nieliczne).
Swoją drogą ponieważ wybierając zmienne do pokazania opierałem się na R^2, więc nie opisałem efektu ”dwójki rodziców w domu”, który może i jest ciekawy ale ma ,,małe R^2”..
Jeżeli porównywać dzieci, które mieszkały z dwójką rodziców vs. dzieci, które mieszkały co najwyżej z jednym rodzicem (wliczając przybranych) to ta druga grupa ma niższą średnią wyników i z testów matematycznych i językowych. Ale zdecydowana większość dzieci mieszka jednak z dwójką rodziców. Gdyby wybierać jeden z efektów 'płeć’ lub 'dwójka rodziców w domu’ to za ,,większy” R^2 wskaże płeć, a same współczynniki wskażą dwójkę rodziców. Pewnie to kwestia wyboru, ale ja za większy ,w skali kraju’ uznałbym efekt płci. Być może to pozostałość z moich doświadczeń z rolnictwem, gdzie genetyczny efekt rzadkiego genu jest mniej ciekawy bo nie ma ,,przemysłowej” wartości.
W każdym razie trafia do mnie porównanie tego do przewidywania pogody, więc do PISA i efektów różnych zmiennych jeszcze wrócę.
Co do różnic pomiędzy lm i svyglm, to nie wiem. Jeżeli to nieduże różnice to może wynikają z tego, że lm to zwykła jednokrokowa algebra, a svyglm to iteracyjna procedura maksymalizująca numerycznie funkcję wiarogodności.
Przyjrzę się temu jak wrócę z urlopu 😉
Cześć, pozwólcie, że podsunę temacik.
Pojawiło się takie coś:
http://www.mapawydatkow.pl/ (można tam pobrać mapę w pliku pdf lub jpg)
Z jednej strony super fajna inicjatywa, ale z drugiej strony jest to materiał do krytycznej analizy sposobu wizualizacji danych liczbowych.
Co mi się nie podoba, to że prezentacja danych liczbowych w formie mapy jest bez sensu. To powinien być np. wykres kołowy, gdzie do każdej cząstki tortu można zrobić zoom i zobaczyć podział na pod-cząstki. W ten sposób łatwiej porównać poszczególne liczby między sobą. Co prawda wielkość kółeczek mówi coś wielkościach liczbowych. Ale nie wiem, czy są one proporcjonalne do średnicy kółka, czy do pola powierzchni. Taka wizualizacja jest moim zdaniem bardzo nieczytelna.
Dzięki za link, pisał mi o nim też Błażej O, wysłałem pytanie do autorów o dostęp do ich danych, będzie można pokusić się o alternatywną wizualizacje.
Bo i temat bardzo ciekawy.
Z tymi cząsteczkami tortu to nie wiem, mnie na tym wykresie najbardziej przeszkadza:
– czasem wydatki są per agencja wydająca (Narodowy Bank Polski), czasem per aktywność (obsługa długu / misja w Afganistanie). Wolałbym wszystko per aktywność, byłoby wiadomo na co pieniądze są wydawane a nie kto wydaje. Ciekawe też ile kosztują takie aktywności jak budowa stadionów na Euro lub budowa autostrad.
– porozrzucane koła nie pozwalają zobaczyć różnicy w skali, np. wydatki na emerytury i renty są 20x większe niż na gospodarkę. Na kołach tej skali 20x nie widać.
– brak zestawienia przychodów i wydatków. A aż się prosi o porównanie przychód/wydatek: składki na NFZ – wydatki na zdrowie, składki na KRUS – wydatki na krus, składki na ZUS – ZUS, Akcyza – infrastruktura.
Bardzo bardzo bardzo chciałbym też zobaczyć takie dane dla kilku lat, można zobaczyć na co wydatki % rosną a na co nie.
W każdym razie fajnie, że takie dane są gromadzone. A nawet jeżeli ta wizualizacja jest nieczytelna to oby każda kolejna była coraz lepsza.
Co do wizualizacji wydatków, to myślę, że nie jest aż taka zła zakładając, że przede wszystkim chcemy pokazać strukturę i względne wielkości poszczególnych części. No i oczywiście zakładając, że powierzchnie kółek wszędzie są proporcjonalne do danych, które pokazują, czego nie sprawdzałem.
Swego czasu Centrum Cyfrowe Polska ogłosiło konkurs na wizualizację prototypowego budżetu zadaniowego dla Polski. Jest kilka ciekawych wizualizacji, ale żadna IMHO nie jest bez wad http://www.otwartybudzet.pl/
Dla porównania NYT zaprezentował interaktywną wizualizację wydatków budżetowych USA. Można sobie klikać i „zoomować”, czego brakowało, jak rozumiem, anuszce http://www.nytimes.com/interactive/2010/02/01/us/budget.html
Porozrzucane kółka nie pozwalają na porównanie wielkości czegokolwiek (różnice 3x są ledwie zauważalne jeżeli kółka są daleko od siebie).
Pod tym względem mapa NYT jest o niebo lepsza, pola prostokątów są lepsze niż pola kół, nie ma zbędnych ozdobników, a zmieściła się jeszcze jedna ciekawa zmienna tzn. % zmiana w stosunku do poprzedniego roku.
Ale największą zaletą mapy NYT jest to, że jest za nią historia, duże wydatki na obronę narodową (powinni napisać obronę plus agresję), rosnące wydatki na edukację.
Ma wrażenie że problem z wizualizacją polskiego budżetu / wydatków bierze się stąd że jest strasznie bałaganiarski. Skoro nie wiadomo na co wydajemy to jak to pokazać graficznie?
Btw, zgodnie z tą wizualizacją wydatki na instytucje d.s. nauki (NCBiR, NCN, PAN, PAA, CNK, MN) w 2011 były o 3% wyższe niż koszt budowy stadionu narodowego (1 973mln vs 1 915mln).
Super, tak trzymać.
Bardzo podobna mi się pomysł przedstawiania „silnika” analizy.
I też zdecydowanie popieram inicjatywę rewizualizacji wapy wydatków.